ParetoGA.cpp

Go to the documentation of this file.

#include "ParetoGA.h"
#include <iostream>
#include <algorithm>
#include <numeric>
#include "NumUtil.h"
#include <boost/algorithm/minmax_element.hpp>
#include <boost/numeric/ublas/io.hpp>

using namespace std;
namespace gplib
  {
    //! Determines whether one vector of misfit values is partially less than the other
    /*! This function object implements the concept of dominance used in NSGA-II and other
     * pareto based optimization methods. The vector fit1 dominates fit2 (or is partially less)
     * \f$ fit1 \prec_p fit2 \Leftrightarrow \forall i fit1_i \leq fit2_i \:\wedge\: \exists i: fit1_i < fit2_i \f$
     */
    class dominates
      {
    public:
      bool operator()(const gplib::rvec &fit1, const gplib::rvec &fit2) const
        {
          typedef vector<double>::const_iterator tit;
          bool smaller = false;
          const size_t length = fit2.size();
          for (size_t i = 0; i < length; ++i)
            {
              if (fit1(i) > fit2(i))
                return false;
              if (fit1(i) < fit2(i))
                smaller = true;
            }
          return smaller;
        }
      };

    void ParetoGA::CalcCrowdingDistance(gplib::rmat &LocalMisFit,
        GeneralPopulation &LocalPopulation)
      {
        const double NearInfinity = 1e50;
        const double tolerance = 1e-10;
        const unsigned int popsize = LocalPopulation.GetPopsize();

        tcrowddistv CrowdingDistances(popsize);
        for (unsigned int i = 0; i < popsize; ++i)
          CrowdingDistances(i) = 0;
        for (unsigned int i = 0; i < nobjective; ++i)
          {
            if (fcmp(GetWeights().at(i), 0, tolerance) != 0)
              {
                pair<tMisfitIterator, tMisfitIterator> MinMax =
                    boost::minmax_element(ublas::row(LocalMisFit, i).begin(),
                        ublas::row(LocalMisFit, i).end());
                double normalize = 1;
                if (fcmp(*(MinMax.second) - *(MinMax.first), 0, tolerance)
                    != 0) // if minimum not equal maximum
                  normalize = *(MinMax.second) - *(MinMax.first);
                for (unsigned int j = 0; j < Ranks.size(); ++j)
                  {
                    tIndexMap IndexMap;
                    for (unsigned int k = 0; k < Ranks.at(j).size(); ++k)
                      {
                        IndexMap.insert(make_pair(LocalMisFit(i,
                            Ranks.at(j).at(k)), Ranks.at(j).at(k)));
                      }
                    tIndexMap::iterator end = --IndexMap.end();
                    tIndexMap::iterator pos = IndexMap.begin();
                    CrowdingDistances(pos->second) = NearInfinity; // set the crowding distance very large
                    if (pos != end)
                      {
                        CrowdingDistances(end->second) = NearInfinity; // for first and last element
                        tIndexMap::iterator previous = IndexMap.begin(); //the previous element is the first element
                        tIndexMap::iterator next = pos; // we start at the second, not the first element
                        advance(pos, 1);
                        advance(next, 2);
                        for (; pos != end; pos++)
                          {
                            CrowdingDistances(pos->second) += (next->first
                                - previous->first) / normalize;
                            ++previous;
                            ++next;
                          }
                      }
                  }
              }
          }
        LocalPopulation.SetCrowdingDistances(CrowdingDistances);
      }

    void ParetoGA::CalcProbabilities(const int iterationnumber,
        gplib::rmat &LocalMisFit, GeneralPopulation &LocalPopulation)
      {
        const unsigned int size = LocalPopulation.GetPopsize();
        gplib::rvec currmisfit(nobjective), compmisfit(nobjective); //local storage
        vector<int> Indices;
        unsigned int rankedelements = 0;

        vector<int> CurrRanks; //the indices of the members in the current rank
        unsigned int i, j, k = 0;
        vector<bool> dominated(size, false); //for each member whether it is dominated
        unsigned int maxindex = size; //the maximum index in the beginning is the population size

        if (!Ranks.empty()) //empty old storage
          Ranks.clear();
        Indices.reserve(size); //intialize population index
        for (i = 0; i < size; ++i)
          Indices.push_back(i);
        while (rankedelements < size) // while we haven't ranked everything
          {
            CurrRanks.clear(); //empty the current rank
            for (i = 0; i < maxindex; ++i) // for all population members that haven't been ranked
              {
                for (j = 0; j < nobjective; ++j) //get the misfit values
                  currmisfit(j) = LocalMisFit(j, Indices.at(i));
                j = 0;
                dominated.at(i) = false; // assume it is not dominated
                while (!dominated.at(i) && (j < maxindex)) //while we haven't found it dominated or reached the maximum index
                  {
                    for (k = 0; k < nobjective; ++k)
                      compmisfit(k) = LocalMisFit(k, Indices.at(j)); //get the misfit values for comparison
                    if (dominates()(compmisfit, currmisfit)) //if current member is dominated
                      {
                        dominated.at(i) = true; // set domination property in vector
                      }
                    ++j;
                  }
              } //end for all population members that haven't been ranked

            for (i = 0; i < maxindex; ++i) //remove non dominated members from comparison
              {
                if (!dominated.at(i)) //if the current member is not dominated
                  {
                    CurrRanks.push_back(Indices.at(i)); // add it to the non-dominated front
                    Indices.at(i) = Indices.at(maxindex - 1); //copy the last active index in the current position
                    dominated.at(i) = dominated.at(maxindex - 1); // same with domination property
                    ++rankedelements;//we have now ranked one element
                    --maxindex; // and have one less to go, the previous last index is now copied
                    --i; // we have to examine the element in this position again, because we copied
                  }
              }
            Ranks.push_back(CurrRanks); //add current rank to overall ranks
          } //now we have ranked everything
        const unsigned int nranks = Ranks.size();
        cout << "NRanks : " << nranks << endl;
        double sum = 0.0; //now we calculate the probabilities
        double prob = 0.0;
        tprobabilityv Probabilities(LocalPopulation.GetPopsize()); //allocate an array
        for (i = 0; i < nranks; ++i) //for all ranks
          {
            for (j = 0; j < Ranks.at(i).size(); ++j) //and members in the rank
              {
                prob = nranks - i; //high rank is low probability
                Probabilities(Ranks.at(i).at(j)) = prob; //Assign to the right member
                sum += prob; //sum for normalization
              }
          }
        Probabilities /= sum; //make values probabilities
        LocalPopulation.SetProbabilities(Probabilities); //assign to population object

        CalcCrowdingDistance(LocalMisFit, LocalPopulation);
      }

    void ParetoGA::Elitism(const int iterationnumber)
      {
        const unsigned int popsize = Population->GetPopsize();
        const unsigned int ngenes = Population->GetGenesize();
        GeneralPopulation IntermediatePop(Population->GetPopulation(),
            Population->GetOldPopulation());
        gplib::rmat IntermediateMisFit(nobjective, 2 * popsize);
        ublas::matrix_range<gplib::rmat> MisFitFirstHalf(IntermediateMisFit,
            ublas::range(0, nobjective), ublas::range(0, popsize));
        ublas::matrix_range<gplib::rmat> MisFitSecondHalf(IntermediateMisFit,
            ublas::range(0, nobjective), ublas::range(popsize, 2 * popsize));
        MisFitFirstHalf = MisFit;
        MisFitSecondHalf = OldMisFit;
        CalcProbabilities(iterationnumber, IntermediateMisFit, IntermediatePop);

        unsigned int included = 0;
        unsigned int currentrank = 0;
        tpopulation Newpopulation(popsize, ngenes);
        while (Ranks.at(currentrank).size() <= (popsize - included)) // as long as we can fit the current rank into the population
          {
            for (unsigned int j = 0; j < Ranks.at(currentrank).size(); ++j)
              {
                ublas::row(Newpopulation, included) = ublas::row(
                    IntermediatePop.GetPopulation(),
                    Ranks.at(currentrank).at(j));
                ++included;
              }
            ++currentrank;
          }
        //Now we have to fill up the rest sorted by distance
        tIndexMap IndexMap; //sorts the current rank by crowding distance and stores the associated index
        for (unsigned int j = 0; j < Ranks.at(currentrank).size(); ++j)
          IndexMap.insert(make_pair(IntermediatePop.GetCrowdingDistances()(
              Ranks.at(currentrank).at(j)), Ranks.at(currentrank).at(j)));
        tIndexMap::iterator pos = IndexMap.begin();
        advance(pos, Ranks.at(currentrank).size() - (popsize - included)); // we have to take the last members with the biggest distance
        while (pos != IndexMap.end())
          {
            ublas::row(Newpopulation, included) = ublas::row(
                IntermediatePop.GetPopulation(), pos->second);
            ++included;
            ++pos;
          }
        Population->SetPopulation(Newpopulation);
      }

    void ParetoGA::PrintRanks(std::ostream &output)
      {
        const unsigned int nobjective = MisFit.size1();
        const unsigned int nranks = Ranks.size();
        for (unsigned int i = 0; i < nranks; ++i)
          {
            output << "Rank: " << i << endl;
            for (unsigned int j = 0; j < Ranks.at(i).size(); ++j)
              {
                for (unsigned int k = 0; k < nobjective; ++k)
                  output << MisFit(k, Ranks.at(i).at(j)) << " ";
                output << endl;
              }
            output << endl;
          }
        output << endl;
      }

    void ParetoGA::PrintFront(std::ostream &output)
      {
        const unsigned int nobjective = MisFit.size1();
        const unsigned int frontsize = Ranks.front().size();
        for (unsigned int i = 0; i < frontsize; ++i)
          {
            for (unsigned int j = 0; j < nobjective; ++j)
              output << MisFit(j, Ranks.front().at(i)) << " ";
            output << endl;
          }
        output << endl << endl;
      }

    ParetoGA::~ParetoGA()
      {
      }
  }