00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 #include <limits.h>
00030 #include <float.h>
00031 #include "CubicLS.hh"
00032
00033 #define VERBOSE
00034
00035 CubicLineSearch::CubicLineSearch(ObjectiveFunction* f, int it)
00036 : LineSearch(f)
00037 { iterMax = it;
00038 }
00039
00040 CubicLineSearch::CubicLineSearch(ObjectiveFunction* f, int it, Vector<double>* interval)
00041 : LineSearch(f,interval)
00042 { iterMax = it;
00043 }
00044
00045 CubicLineSearch::~CubicLineSearch() { if (fp != NULL) fp = NULL;}
00046
00047
00048 Model<double> CubicLineSearch::search(Model<double>& current_solution, Vector<double>& p, double slope, double lambda)
00049 {
00050 int tst = 0;
00051 double alpha2=0, alpha_tmp=0, alpha_prev=0;
00052 double alpha_prev2=0, alpha=0;
00053 double f1=0, f2=0, fprev=0;
00054 double a=0, b=0;
00055 double c=0, cm11=0, cm12=0, cm21=0, cm22=0;
00056 double disc=0;
00057 double new_m=0, old_m=0;
00058 Model<double> new_solution(current_solution);
00059
00060
00061 old_m = fp->performance(current_solution);
00062
00063 iterNum = 0; iterNum++;
00064 alpha = 1.;
00065
00066 new_solution = current_solution.update(1,alpha,p);
00067 new_m = fp->performance(new_solution);
00068
00069 iterNum++;
00070
00071
00072 while (new_m < old_m + (1. - lambda)*alpha*slope && iterNum< iterMax)
00073 {
00074 alpha *= 3;
00075 new_solution = current_solution.update(1, alpha, p);
00076 new_m = fp->performance(new_solution);
00077 iterNum++;
00078 }
00079 if (iterNum == iterMax)
00080 cerr << "Alpha over flowed! \n";
00081
00082
00083 alpha_prev = alpha;
00084 while (new_m > old_m + lambda*alpha*slope && iterNum < iterMax)
00085 {
00086 alpha2 = alpha * alpha;
00087 f1 = new_m - old_m - slope * alpha;
00088
00089 if (tst == 0)
00090 {
00091 alpha_tmp = -slope * alpha2 / (f1 * 2.);
00092
00093
00094 tst = 1;
00095 }
00096 else
00097 {
00098 alpha_prev2 = alpha_prev * alpha_prev;
00099 f2 = fprev - old_m - alpha_prev * slope;
00100
00101 c = 1. / (alpha - alpha_prev);
00102 cm11 = 1. / alpha2;
00103 cm12 = -1. / alpha_prev2;
00104 cm21 = -alpha_prev / alpha2;
00105 cm22 = alpha / alpha_prev2;
00106
00107 a = c * (cm11 * f1 + cm12 * f2);
00108 b = c * (cm21 * f1 + cm22 * f2);
00109 disc = b * b - 3. * a * slope;
00110
00111 if ((Abs(a) > FLT_MIN) && (disc > FLT_MIN))
00112 alpha_tmp = (-b + sqrt(disc)) / (3. * a);
00113 else
00114 alpha_tmp = slope * alpha2 / (2. * f1);
00115
00116 if (alpha_tmp >= .5 * alpha)
00117 alpha_tmp = .5 * alpha;
00118 }
00119 alpha_prev = alpha;
00120 fprev = new_m;
00121
00122 if (alpha_tmp < .1 * alpha)
00123 alpha *= .1;
00124 else
00125 alpha = alpha_tmp;
00126
00127 new_solution = current_solution.update(1, alpha, p);
00128 new_m = fp->performance(new_solution);
00129 iterNum++;
00130 }
00131 #ifdef COOOL_VERBOSE
00132 if (iterNum == iterMax)
00133 cerr << "Alpha under flowed! \n";
00134 #endif
00135 value = new_m;
00136 appendSearchNumber();
00137 return (new_solution);
00138 }
00139
00140 Model<long> CubicLineSearch::search(Model<long>& current_solution, Vector<double>& p, double slope, double lambda)
00141 {
00142 Model<double> temp(current_solution);
00143 temp = search(temp, p, slope, lambda);
00144 Model<long> optm(temp);
00145 return optm;
00146 }