ParetoGA.cpp

Go to the documentation of this file.
00001 #include "ParetoGA.h"
00002 #include <iostream>
00003 #include <algorithm>
00004 #include <numeric>
00005 #include <gsl/gsl_math.h>
00006 #include <boost/algorithm/minmax_element.hpp>
00007 #include <boost/numeric/ublas/io.hpp>
00008 
00009 using namespace std;
00010 
00011 //! Determines whether one vector of misfit values is partially less than the other
00012 /*! This function object implements the concept of dominance used in NSGA-II and other 
00013  * pareto based optimization methods. The vector fit1 dominates fit2 (or is partially less)
00014  * \f$ fit1 \prec_p fit2 \Leftrightarrow \forall i fit1_i \leq fit2_i \:\wedge\: \exists i: fit1_i < fit2_i \f$
00015  */
00016 class dominates {
00017         public:
00018                 bool operator() (const gplib::rvec &fit1, const gplib::rvec &fit2) const
00019                 {
00020                         typedef vector<double>::const_iterator tit;
00021                         bool smaller = false;
00022                         const size_t length = fit2.size();
00023                         for (size_t i = 0; i < length; ++i)
00024                         {
00025                                 if (fit1(i) > fit2(i))
00026                                         return false;
00027                                 if (fit1(i) < fit2(i))
00028                                         smaller = true;
00029                         }
00030                         return smaller;
00031                 } 
00032 };
00033 
00034 void ParetoGA::CalcCrowdingDistance(gplib::rmat &LocalMisFit,GeneralPopulation &LocalPopulation)
00035 {
00036         const double NearInfinity = 1e50;
00037         const double tolerance = 1e-10;
00038         const unsigned int popsize = LocalPopulation.GetPopsize();
00039                 
00040         tcrowddistv CrowdingDistances(popsize);
00041         for (unsigned int i = 0; i < popsize; ++i)
00042                 CrowdingDistances(i) = 0;
00043         for (unsigned int i = 0; i < nobjective; ++i)
00044         {
00045                  if (gsl_fcmp(GetWeights().at(i),0,tolerance) != 0)
00046                  { 
00047                          pair<tMisfitIterator,tMisfitIterator> MinMax = boost::minmax_element(ublas::row(LocalMisFit,i).begin(),ublas::row(LocalMisFit,i).end());
00048                          double normalize = 1;
00049                          if (gsl_fcmp(*(MinMax.second) - *(MinMax.first),0,tolerance) != 0) // if minimum not equal maximum
00050                                 normalize = *(MinMax.second) - *(MinMax.first);
00051                          for (unsigned int j = 0; j < Ranks.size(); ++j)
00052                          {
00053                                 tIndexMap IndexMap;     
00054                                 for (unsigned int k = 0; k < Ranks.at(j).size(); ++k)
00055                                 {
00056                                         IndexMap.insert(make_pair(LocalMisFit(i,Ranks.at(j).at(k)),Ranks.at(j).at(k)));
00057                                 }
00058                                 tIndexMap::iterator end = --IndexMap.end();
00059                                 tIndexMap::iterator pos= IndexMap.begin();
00060                                 CrowdingDistances(pos->second) = NearInfinity; // set the crowding distance very large
00061                                 if ( pos != end)
00062                                 {
00063                                         CrowdingDistances(end->second) = NearInfinity; // for first and last element 
00064                                         tIndexMap::iterator previous = IndexMap.begin(); //the previous element is the first element
00065                                         tIndexMap::iterator next = pos; // we start at the second, not the first element
00066                                         advance(pos,1);
00067                                         advance(next,2);
00068                                         for (; pos != end; pos++) 
00069                                         {
00070                                                 CrowdingDistances(pos->second) += (next->first - previous->first)/normalize;
00071                                                 ++previous;
00072                                                 ++next;
00073                                         }
00074                                 }
00075                          }
00076                  }
00077         }
00078         LocalPopulation.SetCrowdingDistances(CrowdingDistances);
00079 }
00080 
00081 
00082 void ParetoGA::CalcProbabilities(const int iterationnumber,gplib::rmat &LocalMisFit,GeneralPopulation &LocalPopulation)
00083 {
00084         const unsigned int size = LocalPopulation.GetPopsize();
00085         gplib::rvec currmisfit(nobjective),compmisfit(nobjective); //local storage 
00086         vector<int> Indices; 
00087         unsigned int rankedelements = 0;
00088         
00089         vector<int> CurrRanks; //the indices of the members in the current rank
00090         unsigned int i,j,k = 0;
00091         vector<bool> dominated(size,false); //for each member whether it is dominated
00092         unsigned int maxindex = size; //the maximum index in the beginning is the population size
00093         
00094         if (!Ranks.empty()) //empty old storage
00095                 Ranks.clear();
00096         Indices.reserve(size); //intialize population index
00097         for (i = 0; i < size; ++i)
00098                 Indices.push_back(i);
00099         while (rankedelements < size) // while we haven't ranked everything
00100         {
00101                 CurrRanks.clear(); //empty the current rank
00102                 for (i = 0; i < maxindex; ++i) // for all population members that haven't been ranked
00103                 {
00104                         for (j = 0; j < nobjective; ++j) //get the misfit values
00105                                 currmisfit(j) = LocalMisFit(j,Indices.at(i));
00106                         j = 0;
00107                         dominated.at(i) = false; // assume it is not dominated
00108                         while ( !dominated.at(i) && ( j < maxindex )) //while we haven't found it dominated or reached the maximum index
00109                         {
00110                                 for (k = 0; k < nobjective; ++k)
00111                                         compmisfit(k) = LocalMisFit(k,Indices.at(j)); //get the misfit values for comparison
00112                                 if ( dominates()(compmisfit,currmisfit)) //if current member is dominated
00113                                 {
00114                                         dominated.at(i) = true; // set domination property in vector
00115                                 }
00116                                 ++j;
00117                         }
00118                 } //end for all population members that haven't been ranked
00119                 
00120                 for (i = 0; i < maxindex; ++i) //remove non dominated members from comparison
00121                 {
00122                         if (!dominated.at(i)) //if the current member is not dominated
00123                         { 
00124                                 CurrRanks.push_back(Indices.at(i)); // add it to the non-dominated front
00125                                 Indices.at(i) = Indices.at(maxindex-1); //copy the last active index in the current position
00126                                 dominated.at(i) = dominated.at(maxindex-1); // same with domination property
00127                                 ++rankedelements;//we have now ranked one element
00128                                 --maxindex; // and have one less to go, the previous last index is now copied
00129                                 --i; // we have to examine the element in this position again, because we copied
00130                         }
00131                 }
00132                 Ranks.push_back(CurrRanks); //add current rank to overall ranks
00133         } //now we have ranked everything
00134         const unsigned int nranks = Ranks.size();
00135         cout << "NRanks : " << nranks << endl; 
00136         double sum = 0.0; //now we calculate the probabilities
00137         double prob = 0.0;
00138         tprobabilityv Probabilities(LocalPopulation.GetPopsize()); //allocate an array
00139         for (i = 0; i < nranks; ++i) //for all ranks
00140         {
00141                 for (j = 0; j < Ranks.at(i).size(); ++j) //and members in the rank
00142                 {
00143                         prob = nranks - i;  //high rank is low probability
00144                         Probabilities(Ranks.at(i).at(j)) = prob; //Assign to the right member
00145                         sum += prob; //sum for normalization
00146                 }
00147         }
00148         Probabilities /= sum; //make values probabilities
00149         LocalPopulation.SetProbabilities(Probabilities); //assign to population object
00150         
00151         CalcCrowdingDistance(LocalMisFit,LocalPopulation);
00152 }
00153 
00154 void ParetoGA::Elitism(const int iterationnumber)
00155 {
00156         const unsigned int popsize = Population->GetPopsize();
00157         const unsigned int ngenes = Population->GetGenesize();
00158         GeneralPopulation IntermediatePop(Population->GetPopulation(),Population->GetOldPopulation());
00159         gplib::rmat IntermediateMisFit(nobjective, 2*popsize);
00160         ublas::matrix_range<gplib::rmat > MisFitFirstHalf(IntermediateMisFit, ublas::range(0,nobjective), ublas::range(0,popsize));
00161         ublas::matrix_range<gplib::rmat > MisFitSecondHalf(IntermediateMisFit, ublas::range(0,nobjective), ublas::range(popsize,2*popsize));
00162         MisFitFirstHalf = MisFit;
00163         MisFitSecondHalf = OldMisFit;
00164         CalcProbabilities(iterationnumber,IntermediateMisFit,IntermediatePop);
00165         
00166         unsigned int included = 0;
00167         unsigned int currentrank = 0;
00168         tpopulation Newpopulation(popsize,ngenes);
00169         while (Ranks.at(currentrank).size() <= (popsize - included)) // as long as we can fit the current rank into the population
00170         {
00171                 for (unsigned int j = 0; j < Ranks.at(currentrank).size(); ++j)
00172                 {
00173                         ublas::row(Newpopulation,included) = ublas::row(IntermediatePop.GetPopulation(),Ranks.at(currentrank).at(j));
00174                         ++included;
00175                 }
00176                 ++currentrank; 
00177         }
00178         //Now we have to fill up the rest sorted by distance
00179         tIndexMap IndexMap; //sorts the current rank by crowding distance and stores the associated index 
00180         for (unsigned int j = 0; j < Ranks.at(currentrank).size(); ++j)
00181                 IndexMap.insert(make_pair(IntermediatePop.GetCrowdingDistances()(Ranks.at(currentrank).at(j)),Ranks.at(currentrank).at(j))); 
00182         tIndexMap::iterator pos= IndexMap.begin();
00183         advance(pos, Ranks.at(currentrank).size()-(popsize - included)); // we have to take the last members with the biggest distance
00184         while (pos != IndexMap.end())
00185         {
00186                 ublas::row(Newpopulation,included) = ublas::row(IntermediatePop.GetPopulation(),pos->second);
00187                 ++included;
00188                 ++pos;
00189         }
00190         Population->SetPopulation(Newpopulation);
00191 }
00192 
00193 void ParetoGA::PrintRanks(std::ostream &output)
00194 {
00195         const unsigned int nobjective = MisFit.size1();
00196         const unsigned int nranks = Ranks.size(); 
00197         for (unsigned int i = 0; i < nranks; ++i)
00198         {
00199                 output << "Rank: " << i << endl;
00200                 for (unsigned int j = 0; j < Ranks.at(i).size(); ++j)
00201                 {
00202                         for (unsigned int k = 0; k < nobjective; ++k)   
00203                                 output << MisFit(k,Ranks.at(i).at(j)) << " ";
00204                         output << endl;
00205                 }
00206                 output << endl;
00207         }
00208         output << endl;
00209 }
00210 
00211 void ParetoGA::PrintFront(std::ostream &output)
00212 {
00213         const unsigned int nobjective = MisFit.size1();
00214         const unsigned int frontsize = Ranks.front().size(); 
00215         for (unsigned int i = 0; i < frontsize; ++i)
00216         {
00217                 for (unsigned int j = 0; j < nobjective; ++j)   
00218                         output << MisFit(j,Ranks.front().at(i)) << " ";
00219                 output << endl;
00220         }               
00221         output << endl << endl;
00222 }
00223 
00224 ParetoGA::~ParetoGA()
00225 {}
00226 

Generated on Fri Jul 4 15:30:21 2008 for GPLIB++ by  doxygen 1.5.5