surfinvga.cpp

Go to the documentation of this file.

#include <iostream>
#include <fstream>
#include <algorithm>
#include <numeric>
#include <sstream>
#include <string>
#include "GA.h"
#include "C1dInvGaConf.h"
#include "Adaptors.h"
#include <boost/bind.hpp>
#include <boost/shared_ptr.hpp>
#include <boost/filesystem.hpp>
#include <boost/date_time/posix_time/posix_time.hpp>
#include <complex>
#include "Util.h"
#include "SurfaceWaveObjective.h"
#include "SurfaceWaveData.h"
#include "SurfInvGaConf.h"
#include "SeismicModelDiff.h"
#include "ResPkModel.h"


using namespace std;
using namespace gplib;

enum tgatype
  {
  pareto, anneal
  };
typedef boost::shared_ptr<GeneralObjective> pCGeneralObjective;


void SetupAnnealingGA(boost::shared_ptr<GeneralGA> &GA,
    const SurfInvGaConf Configuration)
  {
    double InitTemperature = Configuration.inittemp;
    if (Configuration.inittemp < 0)
      {
        double avgcombfit = 0;
        GA->DoIteration(0, false);
        /*for (int i = 0; i < GA->GetAvgFit().size(); ++i)
         avgcombfit += GA->GetAvgFit().at(i);*/
        avgcombfit = accumulate(GA->GetAvgFit().begin(), GA->GetAvgFit().end(),
            0.0);
        InitTemperature = avgcombfit * abs(Configuration.inittemp);
      }
    AnnealingGA *AnnealGA = dynamic_cast<AnnealingGA*> (GA.get());
    AnnealGA->SetParams(InitTemperature, Configuration.annealinggeneration,
        Configuration.coolingratio);
  }

//! Program to invert MT data for 1D anisotropic structure with a genetic algorithm
int main(int argc, char *argv[])
  {
    boost::posix_time::ptime starttime =
        boost::posix_time::microsec_clock::local_time();
    string version = "$Id: mtanisoga.cpp 1446 2007-11-09 11:05:47Z mmoorkamp $";
    cout << endl << endl;
    cout << "Program " << version << endl;
    cout << "Genetic algorithm to jointly invert MT data for 1D anisotropy"
        << endl;
    cout << "look at mtanisoga.conf for configuration and setup " << endl;
    cout << "The root name of log-files can be overwritten by using " << endl;
    cout << "1dinvga newrootname " << endl;

    SurfInvGaConf Configuration;

    UniformRNG Random;
    try
      {
        Configuration.GetData("surfinvga.conf");

        ttranscribed transcribed;
        tgatype GAtype = pareto;

        string outputbase;

        if (argc > 1)
          outputbase = argv[1];
        else
          outputbase = Configuration.outputbase;

        ofstream misfitfile;
        ofstream valuefile;
        ofstream probfile;
        ofstream distancefile;
        ofstream fitstat((outputbase + ".ftst").c_str());
        ofstream uniquepop((outputbase + ".unpop").c_str());
        ofstream rankfile((outputbase + ".rank").c_str());
        ofstream frontfile((outputbase + ".front").c_str());
        ofstream bestfile((outputbase + ".best").c_str());

        if (Configuration.verbose) // in most cases we do not need all of these values
          {
            misfitfile.open((outputbase + ".msft").c_str());
            valuefile.open((outputbase + ".vals").c_str());
            probfile.open((outputbase + ".probs").c_str());
            distancefile.open((outputbase + ".dist").c_str());
          }

        if (Configuration.thickbase.size() != Configuration.svelbase.size())
          {
            cerr << "Configuration of genes is not consistent ! ";
            return 100;
          }
        const int nparam = 2;
        const int genes = 2 * Configuration.thickbase.size();
        const int popsize = Configuration.popsize;
        const int maxgen = Configuration.generations;
        const int paramlength = genes / nparam;
        ttranscribed basevalues(genes);
        ttranscribed stepsizes(genes);
        tsizev genesizes(genes);

        misfitfile << popsize << " " << maxgen << " " << endl;
        valuefile << genes << " " << popsize << " " << maxgen << endl;

        copy(Configuration.thickbase.begin(), Configuration.thickbase.end(),
            basevalues.begin());
        copy(Configuration.thickstep.begin(), Configuration.thickstep.end(),
            stepsizes.begin());
        copy(Configuration.thicksizes.begin(), Configuration.thicksizes.end(),
            genesizes.begin());

        copy(Configuration.svelbase.begin(), Configuration.svelbase.end(),
            basevalues.begin() + paramlength);
        copy(Configuration.svelstep.begin(), Configuration.svelstep.end(),
            stepsizes.begin() + paramlength);
        copy(Configuration.svelsizes.begin(), Configuration.svelsizes.end(),
            genesizes.begin() + paramlength);

        const int totalgenesize = accumulate(genesizes.begin(),
            genesizes.end(), 0);
        BinaryPopulation Population(popsize, totalgenesize, Random, true);
        BinaryTournamentSelect Select(Random, boost::bind(
            &GeneralPopulation::GetProbabilities, boost::ref(Population)),
            boost::bind(&GeneralPopulation::GetCrowdingDistances, boost::ref(
                Population)));
        StandardPropagation Propagator(&Select, &Population, &Random);
        double indmutationprob = 1.0 - pow(1.0 - Configuration.mutationprob, 1.
            / totalgenesize); // we specify the overall probability for mutation in config file
        cout << "Individual Mutation set to: " << indmutationprob << endl;
        Propagator.SetParams(indmutationprob, Configuration.crossoverprob);
        GrayTranscribe Transcribe(basevalues, stepsizes, genesizes);

        SurfaceWaveData Data;
        Data.ReadFile(Configuration.inputdata);
        boost::shared_ptr<SurfaceWaveObjective> SeisObjective(
            new SurfaceWaveObjective(Data));
        SeisObjective->SetFitExponent(Configuration.fitexponent);
        SeisObjective->SetPoisson(Configuration.poisson);
        SeisObjective->SetErrorLevel(Configuration.errorlevel);

        //setup regularization
        ResPkModel VRefModel;
        VRefModel.ReadModel(Configuration.vrefmodel);
        boost::shared_ptr<SeismicModelDiff> Regularization(
            new SeismicModelDiff(VRefModel));

        vector<pCGeneralObjective> ObjVector;
        ObjVector.push_back(SeisObjective);
        ObjVector.push_back(Regularization);

        boost::shared_ptr<GeneralGA> GA;
        if (Configuration.gatype == "pareto")
          {
            GAtype = pareto;
            boost::shared_ptr<ParetoGA> Pareto(new ParetoGA(&Propagator,
                &Population, &Transcribe, ObjVector, Configuration.threads));
            GA = Pareto;
          }
        else if (Configuration.gatype == "anneal")
          {
            GAtype = anneal;
            boost::shared_ptr<AnnealingGA> Annealing(new AnnealingGA(
                &Propagator, &Population, &Transcribe, ObjVector,
                Configuration.threads));
            GA = Annealing;
          }
        else
          {
            cerr << "Invalid GA type in configuration file." << endl;
            cerr << " Has to be either 'pareto' or 'anneal' !" << endl;
            return 100;
          }
        GA->SetWeights(Configuration.weights);
        //Setting up indices for forward modelling
        GeneralGA::tparamindv ParamIndices;
        std::vector<int> LocalIndices;
        generate_n(back_inserter(LocalIndices), genes, IntSequence(0)); //generate Indices 0 .. nmtparam for MT Objective function
        ParamIndices.push_back(LocalIndices);//both data misfit and regularization
        ParamIndices.push_back(LocalIndices); //work on all parts of the model vector

        GA->SetParameterIndices(ParamIndices);
        GA->SetElitist(Configuration.elitist);
        //we have to do a few things that are specific to the type of genetic algorithm
        if (GAtype == anneal)
          {
            SetupAnnealingGA(GA, Configuration);
            Population.InitPop();
          }
        //now comes the core loop
        for (int i = 0; i < maxgen; ++i) //iterate over the specified number of generattions
          {
            GA->DoIteration(i, i == (maxgen - 1)); // pass iteration number to GA and tell whether it's the final iteration
            if (Configuration.verbose) // in most cases we do not need all of these values
              {
                GA->PrintMisfit(misfitfile);
                GA->PrintTranscribed(valuefile);
                Population.PrintProbabilities(probfile);
                Population.PrintDistances(distancefile);
              }
            GA->PrintBestMisfit(frontfile);

            fitstat << i << " ";
            GA->PrintFitStat(fitstat);
            //system("rm -r j*");
          }
        boost::posix_time::ptime partime =
            boost::posix_time::microsec_clock::local_time();
        std::cout << "Parallel Runtime: " << partime - starttime << std::endl;
        //the following is to organize the output of the best models
        int bestcount = 0;
        UniquePop FinalUnique; //we only want to write each model once, even if it is several times in the population
        tfitvec bestfit(ObjVector.size());
        tpopmember member(Population.GetGenesize());

        vector<int> BestIndices(GA->GetBestModelIndices()); // store the indices of the best models
        const unsigned int noutmodels = GA->GetNBestmodels(); // and how many of them we have
        assert(BestIndices.size() == noutmodels);
        const unsigned int nobjective = ObjVector.size();
        GeneralGA::tparamvector LocalParameters(nobjective);
        for (unsigned int i = 0; i < nobjective; ++i)
          LocalParameters.at(i).resize(ParamIndices.at(i).size(), false);
        ofstream finalmodels((outputbase + ".final").c_str());
        for (unsigned int h = 0; h < noutmodels; ++h) //for each best model
          {
            member = row(Population.GetPopulation(), BestIndices.at(h)); //get the genetic string


            transcribed = Transcribe.GetValues(member); //transcribe the genetic string to model parameters
            GA->SetupParams(transcribed, LocalParameters); //copy the right parameters for each objective function

            for (unsigned int f = 0; f < nobjective; ++f) //calculate the performance for each function
              {
                if (Configuration.weights.at(f) > 0)
                  bestfit(f) = ObjVector.at(f)->CalcPerformance(
                      LocalParameters.at(f));
                else
                  bestfit(f) = 0;

              }

            if (FinalUnique.Insert(bestfit, transcribed)) //if we can insert it, it is a new model that we didn't output yet
              {
                copy(bestfit.begin(), bestfit.end(), ostream_iterator<double> (
                    finalmodels, " ")); //print the misfits in the .final file
                finalmodels << "  " << bestcount << endl; // write the index in the last column for easier identification
                ostringstream filename;
                filename << "best_" << bestcount;
                cout << "Best: " << h << " Index : " << BestIndices.at(h)
                    << " "; //write some info to screen
                for (unsigned int f = 0; f < nobjective; ++f)
                  cout << bestfit(f) << " ";
                cout << " Saved as : " << filename.str();

                SeisObjective->WriteModel(filename.str()); //write model, plotting file and synthetic data
                SeisObjective->WritePlot(filename.str());
                SeisObjective->WriteData(filename.str());

                filename.str("");
                bestcount++;
                cout << endl;
              }

          }

        if (Configuration.verbose) // in most cases we do not need all of these values
          {
            GA->PrintUniquePop(uniquepop);
            FinalUnique.PrintAll(bestfile);
          }
        boost::posix_time::ptime endtime =
            boost::posix_time::microsec_clock::local_time();
        std::cout << "Total Runtime: " << endtime - starttime << std::endl;
      } catch (FatalException &e)
      {
        cerr << e.what() << endl;
        return -1;
      }
  }