anisogajoint.cpp

Go to the documentation of this file.

#include <iostream>
#include <fstream>
#include <algorithm>
#include <numeric>
#include <sstream>
#include <string>
#include "SimpleSelect.h"
#include "BinaryTournamentSelect.h"
#include "UniformRNG.h"
#include "gentypes.h"
#include "BinaryPopulation.h"
#include "StandardPropagation.h"
#include "GrayTranscribe.h"
#include "AnnealingGA.h"
#include "ParetoGA.h"
#include "Aniso1DMTObjective.h"
#include "MultiAnisoSurfaceWaveObjective.h"
#include "ParkSurfaceWaveData.h"
#include "PTensor1DMTObjective.h"
#include "MTStation.h"
#include "Adaptors.h"
#include "MTFitSetup.h"
#include "PTensorMTStation.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 "MTAnisoRoughness.h"
#include "SWAnisoRoughness.h"
#include "AnisoJointConf.h"
#include "MTInvConf.h"
#include "GAConf.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 GAConf 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);
  }

void SetupSurfaceWaveObjective(const string surfinfofile,
    MultiAnisoSurfaceWaveObjective &SurfObjective)
  {
    ifstream surfinfo(surfinfofile.c_str());
    double backazimuth; //first column in swaveinfo
    double avelratio; //third column in swaveinfo
    string inputdata; //second column in swaveinfo
    while (surfinfo.good())
      {

        surfinfo >> backazimuth >> inputdata >> avelratio;
        if (surfinfo.good())
          {
            ParkSurfaceWaveData Data;
            Data.ReadAscii(inputdata);
            SurfObjective.AddMeasurement(Data, backazimuth, avelratio);
          }
      }
  }

//! Program to invert MT and SW 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 1572 2007-12-17 18:20:43Z mmoorkamp $";
    cout << endl << endl;
    cout << "Program " << version << endl;
    cout
        << "Genetic algorithm to jointly invert MT and SW data for 1D anisotropy"
        << endl;
    cout << "look at anisogajoint.conf for configuration and setup " << endl;
    cout << "The root name of log-files can be overwritten by using " << endl;
    cout << "1dinvga newrootname " << endl;

    MTInvConf MTConf;
    AnisoJointConf JointConf;
    GAConf GaConf;

    std::ifstream conffile("anisogajoint.conf");

    MTConf.GetData(conffile);
    conffile.clear();
    conffile.seekg(0, std::ios::beg);
    JointConf.GetData(conffile);
    conffile.clear();
    conffile.seekg(0, std::ios::beg);
    GaConf.GetData(conffile);

    UniformRNG Random;
    try
      {
        ttranscribed transcribed;
        tgatype GAtype = pareto;

        string outputbase;

        if (argc > 1)
          outputbase = argv[1];
        else
          outputbase = JointConf.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 (JointConf.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 ((JointConf.thickbase.size() != JointConf.resbase.size())
            || (JointConf.thickbase.size() != JointConf.aresbase.size())
            || (JointConf.thickbase.size()
                != JointConf.strikebase.size())
            || (JointConf.thickbase.size()
                != JointConf.dstrikebase.size())
            || (JointConf.thickbase.size() != JointConf.avelbase.size()))
          {
            cerr << "Configuration of genes is not consistent ! ";
            return 100;
          }
        const int nparam = 7;
        const int genes = nparam * JointConf.thickbase.size();
        const int popsize = GaConf.popsize;
        const int maxgen = GaConf.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(JointConf.resbase.begin(), JointConf.resbase.end(),
            basevalues.begin());
        copy(JointConf.resstep.begin(), JointConf.resstep.end(),
            stepsizes.begin());
        copy(JointConf.ressizes.begin(), JointConf.ressizes.end(),
            genesizes.begin());

        copy(JointConf.thickbase.begin(), JointConf.thickbase.end(),
            basevalues.begin() + paramlength);
        copy(JointConf.thickstep.begin(), JointConf.thickstep.end(),
            stepsizes.begin() + paramlength);
        copy(JointConf.thicksizes.begin(), JointConf.thicksizes.end(),
            genesizes.begin() + paramlength);

        copy(JointConf.velbase.begin(), JointConf.velbase.end(),
            basevalues.begin() + 2 * paramlength);
        copy(JointConf.velstep.begin(), JointConf.velstep.end(),
            stepsizes.begin() + 2 * paramlength);
        copy(JointConf.velsizes.begin(), JointConf.velsizes.end(),
            genesizes.begin() + 2 * paramlength);

        copy(JointConf.strikebase.begin(), JointConf.strikebase.end(),
            basevalues.begin() + 3 * paramlength);
        copy(JointConf.strikestep.begin(), JointConf.strikestep.end(),
            stepsizes.begin() + 3 * paramlength);
        copy(JointConf.strikesizes.begin(),
            JointConf.strikesizes.end(), genesizes.begin() + 3
                * paramlength);

        copy(JointConf.dstrikebase.begin(),
            JointConf.dstrikebase.end(), basevalues.begin() + 4
                * paramlength);
        copy(JointConf.dstrikestep.begin(),
            JointConf.dstrikestep.end(), stepsizes.begin() + 4
                * paramlength);
        copy(JointConf.dstrikesizes.begin(),
            JointConf.dstrikesizes.end(), genesizes.begin() + 4
                * paramlength);

        copy(JointConf.aresbase.begin(), JointConf.aresbase.end(),
            basevalues.begin() + 5 * paramlength);
        copy(JointConf.aresstep.begin(), JointConf.aresstep.end(),
            stepsizes.begin() + 5 * paramlength);
        copy(JointConf.aressizes.begin(), JointConf.aressizes.end(),
            genesizes.begin() + 5 * paramlength);

        copy(JointConf.avelbase.begin(), JointConf.avelbase.end(),
            basevalues.begin() + 6 * paramlength);
        copy(JointConf.avelstep.begin(), JointConf.avelstep.end(),
            stepsizes.begin() + 6 * paramlength);
        copy(JointConf.avelsizes.begin(), JointConf.avelsizes.end(),
            genesizes.begin() + 6 * 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 - std::pow(1.0
            - GaConf.mutationprob, 1. / totalgenesize); // we specify the overall probability for mutation in config file
        cout << "Individual Mutation set to: " << indmutationprob << endl;
        Propagator.SetParams(indmutationprob, GaConf.crossoverprob);
        GrayTranscribe Transcribe(basevalues, stepsizes, genesizes);

        //setup MT Objective function
        boost::shared_ptr<PlottableObjective> MTObjective;

        if (MTConf.mtfit == "ptensor")
          {
            PTensorMTStation PTData;
            PTData.GetData(JointConf.ptensordata);
            boost::shared_ptr<PTensor1DMTObjective> PTObjective(
                new PTensor1DMTObjective(PTData));
            PTObjective->SetErrorLevel(MTConf.phaseerror);
            MTObjective = PTObjective;
          }
        else
          {
            MTStation MTData;
            MTData.GetData(MTConf.mtinputdata);
            boost::shared_ptr<Aniso1DMTObjective> AnisoObjective(
                new Aniso1DMTObjective(MTData));
            SetupMTFitParameters(MTConf, *AnisoObjective);
            MTObjective = AnisoObjective;
          }
        MTObjective->SetFitExponent(MTConf.mtfitexponent);

        //Setup Surface wave objective
        boost::shared_ptr<MultiAnisoSurfaceWaveObjective>
            MultiSurfaceObjective(new MultiAnisoSurfaceWaveObjective);
        SetupSurfaceWaveObjective("swaveinfo", *MultiSurfaceObjective.get());

        //setup MT regularization
        boost::shared_ptr<MTAnisoRoughness> MTRegularization(
            new MTAnisoRoughness());
        MTRegularization->SetCondDiffWeight(JointConf.conddiffweight);
        MTRegularization->SetAnisotropyWeight(JointConf.anisotropyweight);
        MTRegularization->SetStrikeDiffWeight(JointConf.strikediffweight);

        //setup SW regularization
        boost::shared_ptr<SWAnisoRoughness> SWRegularization(
            new SWAnisoRoughness());
        SWRegularization->SetVelDiffWeight(JointConf.veldiffweight);
        SWRegularization->SetAnisovelWeight(JointConf.anisovelweight);
        SWRegularization->SetStrikeDiffWeight(JointConf.strikediffweight);
        SWRegularization->SetDeltaStrikeDiffWeight(
            JointConf.deltastrikediffweight);

        vector<pCGeneralObjective> ObjVector;
        ObjVector.push_back(MTObjective);
        ObjVector.push_back(MultiSurfaceObjective);
        ObjVector.push_back(SWRegularization);
        ObjVector.push_back(MTRegularization);

        boost::shared_ptr<GeneralGA> GA;
        if (GaConf.gatype == "pareto")
          {
            GAtype = pareto;
            boost::shared_ptr<ParetoGA> Pareto(new ParetoGA(&Propagator,
                &Population, &Transcribe, ObjVector));
            GA = Pareto;
          }
        else if (GaConf.gatype == "anneal")
          {
            GAtype = anneal;
            boost::shared_ptr<AnnealingGA> Annealing(new AnnealingGA(
                &Propagator, &Population, &Transcribe, ObjVector));
            GA = Annealing;
          }
        else
          {
            cerr << "Invalid GA type in configuration file." << endl;
            cerr << " Has to be either 'pareto' or 'anneal' !" << endl;
            return 100;
          }
        GA->SetWeights(JointConf.weights);

        //Setting up indices for forward modelling
        GeneralGA::tparamindv ParamIndices;
        const int nmtparam = 4 * JointConf.resbase.size();
        std::vector<int> MTIndices, SurfaceWaveIndices;
        generate_n(back_inserter(MTIndices), 2 * paramlength, IntSequence(0)); // generate indices for rho and t
        generate_n(back_inserter(MTIndices), paramlength, IntSequence(5
            * paramlength));
        generate_n(back_inserter(MTIndices), paramlength, IntSequence(3
            * paramlength));

        ParamIndices.push_back(MTIndices);//add indices for MT

        generate_n(back_inserter(SurfaceWaveIndices), 6 * paramlength,
            IntSequence(paramlength)); // generate indices for Surface Wave objective
        ParamIndices.push_back(SurfaceWaveIndices); //work on all parts of the model vector
        ParamIndices.push_back(SurfaceWaveIndices);// regularization for SW
        ParamIndices.push_back(MTIndices);//regularization takes same parameters as MT, because only resistivity is regularized

        GA->SetParameterIndices(ParamIndices);
        GA->SetElitist(GaConf.elitist);
        //we have to do a few things that are specific to the type of genetic algorithm
        if (GAtype == anneal)
          {
            SetupAnnealingGA(GA, GaConf);
            Population.InitPop();
          }
        //now comes the core loop
        for (int i = 0; i < maxgen; ++i) //iterate over the specified number of generations

          {
            GA->DoIteration(i, i == (maxgen - 1)); // pass iteration number to GA and tell whether it's the final iteration
            if (JointConf.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);
          }

        //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 (JointConf.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();
                if (JointConf.weights.at(0) > 0)
                  {
                    MTObjective->WriteModel(filename.str() + "_mt.mod"); //write model, plotting file and synthetic data
                    MTObjective->WritePlot(filename.str() + "_mt.plot");
                    MTObjective->WriteData(filename.str());
                  }
                //we also have to write out surface wave data
                if (JointConf.weights.at(1) > 0)
                  {
                    MultiSurfaceObjective->WriteModel(filename.str()
                        + "_sw.mod"); //write model, plotting file and synthetic data
                    MultiSurfaceObjective->WritePlot(filename.str()
                        + "_sw.plot");
                    MultiSurfaceObjective->WriteData(filename.str()
                        + "_sw.cvel"); //write synthetic data as ascii file
                  }
                filename.str("");
                bestcount++;
                cout << endl;
              }

          }

        if (JointConf.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;
      }
  }