GeneralGA.cpp

Go to the documentation of this file.

#include "GeneralGA.h"
#include <cmath>
#include <iostream>
#include <numeric>
#include <sstream>
#include <map>
#include "statutils.h"
#include "VecMat.h"
#include "FatalException.h"
#include <boost/bind.hpp>
#include <algorithm>
#include <cassert>
#include <unistd.h>
#ifdef _OPENMP
#include <omp.h>
#endif
using namespace std;

namespace gplib
  {
    //!Copy the objective function within the shared pointer
    struct CopyFromPointer
      {
      boost::shared_ptr<GeneralObjective> operator()(boost::shared_ptr<
          GeneralObjective> param)
        {
          return boost::shared_ptr<GeneralObjective>(param->clone());
        }
      };

    //! Generate a new copy of the Objective function vector
    struct GenObjective
      {
      GeneralGA::tObjectiveVector operator()(
          const GeneralGA::tObjectiveVector &IndObjective)
        {
          GeneralGA::tObjectiveVector result(IndObjective.size()); //allocate space
          transform(IndObjective.begin(), IndObjective.end(), result.begin(),
              CopyFromPointer());//copy each objective function into result
          return result;
        }
      };

    //! Generate a unique ID for filenames etc. so we can call generated scripts in parallel
    std::string MakeParallelID(const int j, const int i,
        const int iterationnumber, const int Programnum)
      {
        ostringstream filename;
        filename << "j" << j << "i" << i << "g" << iterationnumber << "p"
            << Programnum;
        return filename.str();
      }

    GeneralGA::GeneralGA(GeneralPropagation* const LocalPropagation,
        GeneralPopulation* const LocalPopulation,
        GeneralTranscribe* const LocalTranscribe,
        const tObjectiveVector &IndObjective, const int nthreads) :
      CombMisFit(LocalPopulation->GetPopsize()), AvgFit(IndObjective.size()),
          MaxFit(IndObjective.size()), MinFit(IndObjective.size()), Weights(
              IndObjective.size(), 1), nobjective(IndObjective.size()),
          Transcribed(LocalPopulation->GetPopsize(),
              LocalTranscribe->GetNparams()), MisFit(IndObjective.size(),
              LocalPopulation->GetPopsize()), Objective(IndObjective),
          Propagation(LocalPropagation), Population(LocalPopulation),
          Transcribe(LocalTranscribe)
      {
        Threads = nthreads;
        Elitist = true;
        Programnum = getpid();
      }

    GeneralGA::~GeneralGA()
      {
      }

    //Prepare the appropriate parameters for each objective function
    void GeneralGA::SetupParams(const ttranscribed &member,
        tparamvector &params)
      {
        for (unsigned int i = 0; i < nobjective; ++i)
          for (unsigned int j = 0; j < ParameterIndices.at(i).size(); ++j)
            {
              params.at(i)(j) = member(ParameterIndices.at(i).at(j));
            }
      }

    void GeneralGA::SetWeights(const std::vector<double> &LocalWeights)
      {
        if (LocalWeights.size() != Weights.size()) //Weights has been pre allocated based on the number of objective functions
          throw FatalException(
              "Configured number of weights does not match number of objective functions !");
        Weights = LocalWeights;
      }

    void GeneralGA::SetParameterIndices(const tparamindv &Indices)
      {
        ParameterIndices = Indices;
      }

    void GeneralGA::PrintFitStat(std::ostream &output)
      {
        for (unsigned int i = 0; i < nobjective; ++i)
          {
            output << MaxFit.at(i) << " " << MinFit.at(i) << " "
                << AvgFit.at(i) << "   ";
          }
        output << CombMaxFit << " " << CombMinFit << " " << CombAvgFit << "   ";
        output << endl << flush;
      }

    void GeneralGA::PrintMisfit(std::ostream &output)
      {
        const unsigned int isize = MisFit.size1();
        const unsigned int jsize = MisFit.size2();

        for (unsigned int i = 0; i < isize; ++i)
          {
            for (unsigned int j = 0; j < jsize; ++j)
              output << MisFit(i, j) << " ";
            output << endl;
          }
        copy(CombMisFit.begin(), CombMisFit.end(), ostream_iterator<double> (
            output, " "));
        output << endl << endl;
      }

    void GeneralGA::PrintTranscribed(std::ostream &output)
      {
        const unsigned int isize = Transcribed.size1();
        const unsigned int jsize = Transcribed.size2();
        unsigned int i, j;
        for (i = 0; i < isize; ++i)
          {
            for (j = 0; j < jsize; ++j)
              output << Transcribed(i, j) << " ";
            output << endl;
          }
        output << endl << flush;
      }

    void GeneralGA::PrintBestMisfit(std::ostream &output)
      {
        std::vector<int> Indices(GetBestModelIndices());
        const unsigned int isize = Indices.size();
        const unsigned int jsize = MisFit.size1();

        for (unsigned int i = 0; i < isize; ++i)
          {
            const int currindex = Indices.at(i);
            for (unsigned int j = 0; j < jsize; ++j)
              output << MisFit(j, currindex) << " ";
            output << endl;
          }
        output << endl << endl;
      }

    void GeneralGA::CalcMisfit(const int iterationnumber)
      {
        int calculatecount = 0;
        int newcount = 0;
        // popsize cannot be unsigned because loop variables for openmp have to be signed
        const int popsize = Population->GetPopsize();
        ublas::vector<bool> AlreadyCalculated(popsize);
        tfitvec fitvec(nobjective);

#pragma omp parallel for default(shared)
        for (int i = 0; i < popsize; ++i)
          {
            tparamvector LocalParameters(nobjective);
            for (unsigned int k = 0; k < nobjective; ++k)
              {
                LocalParameters.at(k).resize(ParameterIndices.at(k).size(),
                    false);
              }
            row(Transcribed, i) = Transcribe->GetValues(row(
                Population->GetPopulation(), i));
            AlreadyCalculated(i) = UniquePopHist.Find(row(Transcribed, i),
                fitvec);
            if (!AlreadyCalculated(i))
              {
                SetupParams(row(Transcribed, i), LocalParameters);
                tObjectiveVector LocalObjective(GenObjective()(Objective));
                for (unsigned int j = 0; j < nobjective; ++j)
                  {
                    LocalObjective.at(j)->SetParallelID(MakeParallelID(j, i,
                        iterationnumber, Programnum));
                    if (Weights.at(j) != 0)
                      {
#pragma omp critical
                          {
                            LocalObjective.at(j)->PreParallel(
                                LocalParameters.at(j));
                          }
                        LocalObjective.at(j)->SafeParallel(
                            LocalParameters.at(j));
#pragma omp critical
                          {
                            MisFit(j, i) = Weights.at(j)
                                * LocalObjective.at(j)->PostParallel(
                                    LocalParameters.at(j));
                            CombMisFit.at(i) += MisFit(j, i);
                          }

                      }
                    else
                      {
                        MisFit(j, i) = 0;
                      }
                  }
                tfitvec fitvector(column(MisFit, i));
#pragma omp critical
                  {
                    UniquePopHist.Insert(fitvector, ublas::matrix_row<
                        gplib::rmat>(Transcribed, i));
                  }
                newcount++;
              }
            else // if already calculated
              {
                column(MisFit, i) = fitvec;
                CombMisFit.at(i) = ublas::sum(fitvec);
                calculatecount++;
              }
          }
        cout << "New models: " << newcount << " Re-used models: "
            << calculatecount << endl;
      }

    void GeneralGA::DoIteration(const int iterationnumber, const bool last)
      {
        cout << endl << endl << "Iteration: " << iterationnumber + 1 << endl;
        cout << "Calculating Misfit ...";
        CalcMisfit(iterationnumber);
        cout << " done" << endl << flush;
        //cout << "Before Elitism " << endl;
        //PrintMisfit(cout);
        cout << endl;
        if ((iterationnumber >= 1) && Elitist)
          {
            Elitism(iterationnumber);
          }
        CalcMisfit(iterationnumber); // we have to update the misfit after elitism
        //because we cache the misfit values, the cost is low for this
        CalcProbabilities(iterationnumber, MisFit, *Population);
        for (unsigned int i = 0; i < nobjective; ++i) // do some statistics on the misfit

          {
            ublas::matrix_row<gplib::rmat> mr(MisFit, i);
            AvgFit.at(i) = Mean(mr.begin(), mr.end());
            MaxFit.at(i) = *max_element(mr.begin(), mr.end());
            MinFit.at(i) = *min_element(mr.begin(), mr.end());
          }
        CombAvgFit = Mean(CombMisFit.begin(), CombMisFit.end());
        CombMaxFit = *max_element(CombMisFit.begin(), CombMisFit.end());
        CombMinFit = *min_element(CombMisFit.begin(), CombMisFit.end());
        if (!last) // if we're not in the last iteration, we store the last population for elitism

          {
            Population->StoreOldPopulation();
            OldMisFit = MisFit;
            Propagation->NextGeneration();
          }
      }
  }