ParetoGA.cpp

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

Generated on Tue May 4 16:52:14 2010 for GPLIB++ by  doxygen 1.5.8