ReClaM: CMAOptimizer.cpp Source File

00001 //===========================================================================
00040 //===========================================================================
00041 
00042 
00043 #include <ReClaM/CMAOptimizer.h>
00044 
00045 
00046 CMAOptimizer::CMAOptimizer(int verbosity)
00047 {
00048     parents = NULL;
00049     offspring = NULL;
00050 
00051     this->verbosity = verbosity;
00052 }
00053 
00054 CMAOptimizer::~CMAOptimizer()
00055 {
00056     if (parents != NULL) delete parents;
00057     if (offspring != NULL) delete offspring;
00058 }
00059 
00060 
00061 void CMAOptimizer::init(Model& model)
00062 {
00063     init(model, 0.01);
00064 }
00065 
00066 void CMAOptimizer::init(Model& model, double sigma, eMode mode, bool best, int lambda, int mu)
00067 {
00068     Array<double> s(model.getParameterDimension());
00069     s = sigma;
00070     init(model, s, mode, best, lambda, mu);
00071 }
00072 
00073 void CMAOptimizer::init(Model& model, const Array<double>& sigma, eMode mode, bool best, int lambda, int mu)
00074 {
00075     if (parents != NULL) delete parents;
00076     if (offspring != NULL) delete offspring;
00077 
00078     cmaMode = mode;
00079     returnBestIndividual = best;
00080 
00081     int i;
00082     int dim = model.getParameterDimension();
00083 
00084     if ((cmaMode == modeRankMuUpdate) || (cmaMode == modeRankOneUpdate))
00085     {
00086         if (lambda <= 0) lambda = cma.suggestLambda(dim);
00087         if (mu <= 0) mu = cma.suggestMu(lambda);
00088         if (verbosity > 0) printf("[CMA] lambda=%d mu=%d\n", lambda, mu);
00089 
00090         ChromosomeT<double> chrom_w(dim);
00091         ChromosomeT<double> chrom_0(dim);
00092         std::vector<double> stdv(dim);
00093         for (i = 0; i < dim; i++)
00094         {
00095             chrom_w[i] = model.getParameter(i);
00096             chrom_0[i] = 0.0;
00097             stdv[i] = sigma(i);
00098         }
00099 
00100         parents = new Population(mu, chrom_w, chrom_0);
00101         parents->setMinimize();
00102         offspring = new Population(lambda, chrom_w, chrom_0);
00103         offspring->setMinimize();
00104 
00105         if (cmaMode == modeRankMuUpdate)
00106             cma.init(dim, stdv, 1.0, *parents, CMA::superlinear, CMA::rankmu);
00107         else if (cmaMode == modeRankOneUpdate)
00108             cma.init(dim, stdv, 1.0, *parents, CMA::superlinear, CMA::rankone);
00109     }
00110     else if (cmaMode == modeOnePlusOne)
00111     {
00112         ChromosomeCMA chrom(dim);
00113         for (i = 0; i < dim; i++)
00114         {
00115             chrom[i] = model.getParameter(i);
00116         }
00117 
00118         parents = new PopulationCT<ChromosomeCMA>(1, chrom);
00119         parents->setMinimize();
00120         offspring = new PopulationCT<ChromosomeCMA>(1, chrom);
00121         offspring->setMinimize();
00122 
00123         ecma.init((*(PopulationCT<ChromosomeCMA>*)parents)[0], sigma, 1);
00124     }
00125 
00126     bestFitness = 1e100;
00127     bestParameters.resize(dim, false);
00128 
00129     uncertaintyHandling = false;
00130     bFirstIteration = true;
00131 }
00132 
00133 void CMAOptimizer::initUncertainty(Model& model, double sigma, unsigned int maxEvals, double alpha, double theta, eMode mode, int lambda, int mu)
00134 {
00135     init(model, sigma, mode, false, lambda, mu);
00136 
00137     uncertaintyHandling = true;
00138 
00139     this->maxEvals = maxEvals;
00140     this->alpha = alpha;
00141     this->theta = theta;
00142 
00143     objective.resetCount();
00144 }
00145 
00146 void CMAOptimizer::initUncertainty(Model& model, const Array<double>& sigma, unsigned int maxEvals, double alpha, double theta, eMode mode, int lambda, int mu)
00147 {
00148     init(model, sigma, mode, false, lambda, mu);
00149 
00150     uncertaintyHandling = true;
00151 
00152     this->maxEvals = maxEvals;
00153     this->alpha = alpha;
00154     this->theta = theta;
00155 
00156     objective.resetCount();
00157 }
00158 
00159 double CMAOptimizer::optimize(Model& model, ErrorFunction& errorfunction, const Array<double>& input, const Array<double>& target)
00160 {
00161     objective.Set(model, errorfunction, input, target);
00162 
00163     int dim = model.getParameterDimension();
00164     int i, o;
00165     Individual* pI;
00166     ChromosomeCMA* pC;
00167 
00168     if (bFirstIteration)
00169     {
00170         // ensure that the parents have a valid fitness value
00171         int i, ic = parents->size();
00172         for (i = 0; i < ic; i++)
00173         {
00174 //          Ind2Model((*parents)[i], model);
00175 //          (*parents)[i].setFitness(errorfunction.error(model, input, target));
00176             pI = &((*offspring)[0]);
00177             pC = (ChromosomeCMA*)(&((*pI)[0]));
00178             (*parents)[i].setFitness(objective.fitness(*pC));
00179         }
00180         bFirstIteration = false;
00181         objective.resetCount();
00182     }
00183 
00184     if ((cmaMode == modeRankMuUpdate) || (cmaMode == modeRankOneUpdate))
00185     {
00186         int lambda = offspring->size();
00187 
00188         // create lambda feasible offspring
00189         for (o = 0; o < lambda; o++)
00190         {
00191             pI =& (offspring->operator [](o));
00192             do
00193             {
00194                 cma.create(*pI);
00195                 Ind2Model(*pI, model);
00196             }
00197             while (! model.isFeasible());
00198 //          pI->setFitness(errorfunction.error(model, input, target));
00199             pC = (ChromosomeCMA*)(&((*pI)[0]));
00200             pI->setFitness(objective.fitness(*pC));
00201         }
00202 
00203         if (uncertaintyHandling)
00204         {
00205             double uncertainty = UncertaintyQuantification(*offspring, objective, theta);
00206             unsigned int n = objective.getN();
00207             if (verbosity > 1) { printf("[n=%d]", n); fflush(stdout); }
00208             unsigned int new_n = n;
00209             if (uncertainty > 0.0)
00210             {
00211                 new_n = (unsigned int)(alpha * n);
00212                 if (new_n == n) new_n++;
00213                 if (new_n > maxEvals) new_n = maxEvals;
00214             }
00215             if (uncertainty < 0.0)
00216             {
00217                 new_n = (unsigned int)(n / alpha);
00218                 if (new_n == n) new_n--;
00219                 if (new_n == 0) new_n++;
00220             }
00221             objective.setN(new_n);
00222         }
00223 
00224         // selection
00225         parents->selectMuLambda(*offspring, 0);
00226 
00227         // strategy adaptation
00228         cma.updateStrategyParameters(*parents);
00229     }
00230     else if (cmaMode == modeOnePlusOne)
00231     {
00232         // create one feasible offspring
00233         do
00234         {
00235             ecma.Mutate((*(PopulationCT<ChromosomeCMA>*)parents)[0], *(PopulationCT<ChromosomeCMA>*)offspring);
00236             pI = &((*(PopulationCT<ChromosomeCMA>*)offspring)[0]);
00237             Ind2Model(*pI, model);
00238         }
00239         while (! model.isFeasible());
00240         pC = (ChromosomeCMA*)(&((*pI)[0]));
00241 //      pI->setFitness(errorfunction.error(model, input, target));
00242         pI->setFitness(objective.fitness(*pC));
00243 
00244         // selection and strategy adaptation
00245         ecma.SelectAndUpdateStrategyParameters((*(PopulationCT<ChromosomeCMA>*)parents)[0], *(PopulationCT<ChromosomeCMA>*)offspring);
00246     }
00247     else
00248     {
00249         throw SHARKEXCEPTION("[CMA::optimize] invalid mode");
00250     }
00251 
00252     // If this is the best solution ever seen
00253     // remember the parameters and the fitness
00254     pI = &((*parents)[0]);
00255     pC = (ChromosomeCMA*)(&((*pI)[0]));
00256     double f = pI->fitnessValue();
00257 
00258     if (returnBestIndividual)
00259     {
00260         if (verbosity > 0) printf("[CMA] current fitness: %g best fitness: %g\n", f, bestFitness);
00261         if (f < bestFitness)
00262         {
00263             bestFitness = f;
00264             for (i = 0; i < dim; i++) bestParameters(i) = (*pC)[i];
00265         }
00266 
00267         // Copy the globally best solution to the model
00268         // and return the corresponding error value.
00269         for (i = 0; i < dim; i++) model.setParameter(i, bestParameters(i));
00270         return bestFitness;
00271     }
00272     else
00273     {
00274         if (verbosity > 0) printf("[CMA] current fitness: %g\n", f);
00275 
00276         // Copy the best current solution to the model
00277         // and return the corresponding error value.
00278         for (i = 0; i < dim; i++) model.setParameter(i, (*pC)[i]);
00279         return f;
00280     }
00281 }
00282 
00283 int CMAOptimizer::getLambda()
00284 {
00285     return offspring->size();
00286 }
00287 
00288 void CMAOptimizer::Ind2Model(Individual& ind, Model& model)
00289 {
00290     int i, ic = model.getParameterDimension();
00291     ChromosomeCMA* pC = (ChromosomeCMA*)(&ind[0]);
00292 
00293     SIZE_CHECK(ic == (int)pC->size());
00294 
00295     for (i = 0; i < ic; i++) model.setParameter(i, (*pC)[i]);
00296 }
00297 
00298 
00300 
00301 
00302 void CMAOptimizer::ModelFitness::Set(Model& model, ErrorFunction& errorfunction, const Array<double>& input, const Array<double>& target)
00303 {
00304     m = &model;
00305     e = &errorfunction;
00306     i = &input;
00307     t = &target;
00308 }
00309 
00310 double CMAOptimizer::ModelFitness::fitness(const std::vector<double>& v)
00311 {
00312     unsigned int p, pc = v.size();
00313     for (p=0; p<pc; p++) m->setParameter(p, v[p]);
00314     double ret = 0.0;
00315     for (p=0; p<evals; p++) ret += e->error(*m, *i, *t);
00316     return ret / (double)evals;
00317 }