1dinvga.cpp

Go to the documentation of this file.

/*!
 * \file
 * The program 1dinvga is the core joint inversion program for MT, receiver functions
 * and surface wave data. It reads in its settings from the configuration file <CODE> 1dinvga.conf </CODE>, described
 * below, and produces several log files and the final data, model(s) and plot file(s) for each
 * inverted dataset.
 *
 <H2> General Configuration in 1dinvga.conf </H2>
 The options in <CODE> 1dinvga.conf </CODE> are as follows:
 <H3> General program options </H3>
 <TABLE>
 <TR> <TD> <B>Name</B> </TD> <TD><B> Default</B></TD> <TD><B>Description</B></TD>
 <TR> <TD>verbose</TD> <TD> false </TD> <TD> Print detailed information about all models used in the invbersion, only useful for debugging</TD></TR>
 <TR> <TD>outputbase</TD> <TD> </TD> <TD> The start of the name of all logfiles. This can also be specified on the command lin.e</TD> </TR>
 <TR> <TD>threads</TD><TD> 1 </TD> <TD>If the program has been compiled with openmp support, this specifies the number of simultaneous threads. </TD></TR>
 </TABLE>

 <H3> Genetic algorithm options </H3>
 <TABLE>
 <TR> <TD> <B>Name</B> </TD> <TD><B> Default</B></TD> <TD><B>Description</B></TD>
 <TR> <TD>gatype</TD> <TD> </TD> <TD> Can be either <CODE>anneal</CODE> or <CODE> pareto</CODE>. If set to <CODE>anneal</CODE> we use an annealing type GA that minimizes
 a weighted sum of the objective functions with the weights specified below. For this type of GA we also have to specify inittemp and coolingratio.
 If set to <CODE> pareto</CODE>, we use NSGA-II as described by Deb et al. (2002). In this case inittemp and coolingratio are ignored. </TD></TR>
 <TR> <TD>inittemp</TD><TD> </TD> <TD>(For GA type <CODE>anneal</CODE>)The initial temperature in the cooling schedule. High values mean that there will be no big differences in
 probability between well fitting models and bad fitting models. Low values make those differences much more severe. If set to negative values the temperature
 will be determined from the average misfit values in the first iteration.   </TD> </TR>
 <TR> <TD>annealinggeneration</TD><TD> </TD><TD>(For GA type <CODE>anneal</CODE>) The generation number when the initial tempature is divided by coolingratio to enforce stronger
 discrimination between good and bad fitting models.</TD></TR>
 <TR> <TD>coolingratio</TD><TD> </TD><TD>(For GA type <CODE>anneal</CODE>) The value by which inittemp is divided at generation <CODE>annealinggeneration</CODE>
 Higher values mean stronger discrimination.</TD>/TR>
 <TR> <TD>popsize</TD><TD> </TD><TD>The population size for the genetic algorithm. Typical values 50...1000 </TD></TR>
 <TR> <TD>generations</TD><TD> </TD><TD>The number of generations for the GA to run.Typical values 50...500.  </TD>/TR>
 <TR> <TD>mutationprob</TD><TD> </TD><TD>The probability for a mutation within one member, typically 0.1...0.2. The probability for mutation of
 individual bits is calculated from this number depending on the length of the encoding specified below. </TD></TR>
 <TR> <TD>crossoverprob</TD><TD> </TD><TD>The probability for two model to perform crossover. Typically 0.4...0.7. </TD></TR>
 <TR> <TD>elitist</TD><TD> </TD><TD>Do we want to keep the best model(s)  from the last generation (true), or fully depend
 on crossover and mutation with the risk of loosing good models (false).  </TD></TR>
 </TABLE>
 <H3> Data information</H3>

 <TABLE>
 <TR> <TD> <B>Name</B> </TD> <TD><B> Default</B></TD> <TD><B>Description</B></TD>
 <TR> <TD>recinfofile</TD><TD> </TD><TD>The file containing information about the receiver function data. This file is decribed below. </TD></TR>
 <TR> <TD>mtinputdata</TD><TD> </TD><TD>The name of the file containing the MT impedances </TD></TR>
 <TR> <TD>dispdata</TD><TD> </TD><TD>The name of the file containing the Rayleigh wave dispersion data</TD></TR>
 </TABLE>

 <H3> Error floor information </H3>
 <TABLE>
 <TR> <TD> <B>Name</B> </TD> <TD><B> Default</B></TD> <TD><B>Description</B></TD>
 <TR> <TD>tensorerror</TD> <TD> 0.02</TD> <TD>Error floor for MT impedances, if selected below </TD></TR>
 <TR> <TD>reserror  </TD> <TD>0.02 </TD> <TD>Error floor for MT apparent resistivities, if selected below </TD></TR>
 <TR> <TD>phaseerror </TD> <TD>0.02 </TD><TD>Error floor for MT phase, if selected below </TD></TR>
 <TR> <TD>surferror </TD><TD>0.02 </TD><TD>Error floor for Rayleigh wave dispersion data </TD></TR>
 </TABLE>

 <H3> Global receiver function parameters </H3>
 <TABLE>
 <TR> <TD> <B>Name</B> </TD> <TD><B> Default</B></TD> <TD><B>Description</B></TD>
 <TR> <TD>poisson</TD><TD> </TD><TD> Poissons ratio, we calculate Vp from Vs using this ratio</TD></TR>
 <TR> <TD>normrec</TD><TD> true </TD><TD> Are we inverting receiver functions that have been normalized to a maximum amplitude of 1 ?</TD></TR>
 <TR> <TD>starttime</TD><TD> </TD><TD>Start of the fit window in seconds. Any data in the receiver function before this time is ignored. </TD></TR>
 <TR> <TD>endtime</TD><TD> </TD><TD>End of the fit window in seconds. Any data in the receiver function after this time is ignored. </TD></TR>
 <TR> <TD>recmethod</TD><TD> </TD><TD>Method used to calculate the measured data. Can be <CODE>specdiv</CODE> for spectral division,
 or <CODE>iterdecon</CODE> for iterative deconvolution. </TD></TR>
 </TABLE>

 <H3> MT data parameters </H3>
 <TABLE>
 <TR> <TD> <B>Name</B> </TD> <TD><B> Default</B></TD> <TD><B>Description</B></TD>
 <TR> <TD>mode</TD><TD> </TD><TD>Which component of the tensor to fit. Can be <CODE>xy</CODE> or <CODE>yx</CODE>. </TD></TR>
 <TR> <TD>mtfit</TD><TD> </TD><TD>Which aspect of the MT data to fit. Possible values are shown in the table below.  </TD></TR>
 </TABLE>

 <H3> Possible mtfit values </H3>
 <TABLE>
 <TR> <TD> <B>Name</B> </TD>  <TD><B>Description</B></TD> </TR>
 <TR> <TD> appres </TD>  <TD>Fit only apparent resistivity of selected mode</TD> </TR>
 <TR> <TD> phase </TD>  <TD>Fit only phase of selected mode</TD> </TR>
 <TR> <TD> resphase </TD>  <TD>Fit apparent resistivity  and phase of selected mode</TD> </TR>
 <TR> <TD> berd </TD>  <TD>Fit Berdichevskiy invariant. This is independent of selected mode</TD> </TR>
 <TR> <TD> det</TD>  <TD>Fit determinant of impedance tensor. This is independent of selected mode.</TD> </TR>
 </TABLE>

 <H3> Regularization and fit information </H3>
 <TABLE>
 <TR> <TD> <B>Name</B> </TD> <TD><B> Default</B></TD> <TD><B>Description</B></TD>
 <TR> <TD>usevrefmodel </TD><TD> false </TD><TD>Do we want to regularize the seismic parameters by the difference to a reference model (true/false) ? </TD></TR>
 <TR> <TD>vrefmodel</TD><TD> </TD><TD>If true above, the name of the reference model file. </TD></TR>
 <TR> <TD>mtfitexponent </TD><TD> 2</TD><TD>The exponent for the difference  between measured and calculated data. Higher exponents penalize
 any departure from the measured data. This ensures a tighter fit, but also increases the influence of noise.</TD></TR>
 <TR> <TD>recfitexponent </TD><TD>2 </TD><TD>Same for receiver function data. </TD></TR>
 <TR> <TD>surffitexponent </TD><TD> 2</TD><TD>Same for dispersion data. </TD></TR>
 </TABLE>


 <H3> Parameter encoding </H3>
 The number of layers is determined by the number of entries, it has to be equal for all
 encoding parameters.
 <TABLE>
 <TR> <TD> <B>Name</B> </TD> <TD><B> Default</B></TD> <TD><B>Description</B></TD>
 <TR> <TD>thickbase</TD><TD> </TD><TD>Minimum thickness for each layer in km </TD></TR>
 <TR> <TD>thickstep</TD><TD> </TD><TD>Stepsize for each layer in km </TD></TR>
 <TR> <TD>thicksizes</TD><TD> </TD><TD>Number of bits used to encode the thickness of each layer </TD></TR>
 <TR> <TD>resbase</TD><TD> </TD><TD>Base 10 logarithm of the minimum resistivity, i.e. 0 corresponds to 1 Ohm.m </TD></TR>
 <TR> <TD>resstep</TD><TD> </TD><TD>Stepsize for base 10 logarithm of resistivity  </TD></TR>
 <TR> <TD>ressizes</TD><TD> </TD><TD>Number of bits used to encode resistivity </TD></TR>
 <TR> <TD>svelbase</TD><TD> </TD><TD>Minimum S-Wave velocity for each layer in km/s </TD></TR>
 <TR> <TD>svelstep</TD><TD> </TD><TD>Stepsize for S_Wave velocity for each layer </TD></TR>
 <TR> <TD>svelsizes</TD><TD> </TD><TD>Number of bits used to encode Vs </TD></TR>
 </TABLE>

 <H3> Weighting of datasets </H3>
 <TABLE>
 <TR> <TD> <B>Name</B> </TD> <TD><B> Default</B></TD> <TD><B>Description</B></TD>
 <TR> <TD>weights</TD> <TD> </TD>  <TD>These are the weights for the different objective functions. The current version needs exactly 5 weights.
 They correspond to MT misfit, Receiver function misfit, Surface Wave Misfit, Seismic Regularization and MT regularization, respectively.
 If a weight is set to zero the calculation of the misfit is completely disabled and the inversion thus runs faster. For pareto type inversion
 the value of weight is irrelevant as long as it is different from zero. For annealing it determines the influence of the corresponding dataset. </TD></TR>
 </TABLE>

 <H2> Recinfo file </H2>
 The recinfo file contains information about the different receiver functions we are trying to fit. Each
 line in the file corresponds to one receiver function and requires the following entries
 <TABLE>
 <TR> <TD>Ray parameter in s/km </TD>
 <TD>\f$ \sigma\f$ for gaussian filter </TD> <TD> Waterlevel </TD> <TD>Relative error </TD> <TD>Shift in s </TD> <TD> Type of data </TD>
 <TD>Filename </TD></TR>
 </TABLE>
 Type of data is a single character. If it is 'a' we are also fitting absolute velocity information as calculated by the
 program rfvel and described in Svenningsen, GJI 170, 2007. In this case the following
 additional fields are needed
 <TABLE>
 <TR> <TD>Weight for RF </TD> <TD>Weight for absolute velocities </TD> <TD>Filename for absolute velocities </TD> </TR>
 </TABLE>
 If data type is any other character we only fit receiver function information and the absolute velocity information cannot be there.
 <H2> Other features</H2>
 - You can gracefully stop the inversion any time by creating a file called "abort" in the working directory.
 The inversion will stop after the current iteration and will write out the usual information about best models etc.
 This is useful when the chosen number of generations was too high and the algorithm has already reached an acceptable solution.
 The file is automatically deleted.
 */

#include <iostream>
#include <iomanip>
#include <fstream>
#include <algorithm>
#include <numeric>
#include <sstream>
#include <string>
#include <boost/bind.hpp>
#include <boost/shared_ptr.hpp>
#include <boost/filesystem.hpp>
#include <boost/date_time/posix_time/posix_time.hpp>
#include "GA.h"
#include "Iso1DMTObjective.h"
#include "MTStation.h"
#include "AbsVelRecObjective.h"
#include "Multi1DRecObjective.h"
#include "SeismicDataComp.h"
#include "CombinedRoughness.h"
#include "SeismicModelDiff.h"
#include "ResPkModel.h"
#include "SurfaceWaveObjective.h"
#include "SurfaceWaveData.h"
#include "Adaptors.h"
#include "MTFitSetup.h"
#include "Util.h"
#include "MTRoughness.h"
#include "SurfaceWaveData.h"
#include "FileAbort.h"
#include "SeisTools.h"
#include "IsoJointConf.h"
#include "MTInvConf.h"
#include "SurfInvConf.h"
#include "RecInvConf.h"
#include "GAConf.h"

using namespace std;
using namespace gplib;

enum tgatype
  {
  pareto, anneal
  };

typedef boost::shared_ptr<GeneralObjective> pCGeneralObjective;

//! Set the annealing GA temperature parameter, we should only call this function if we are SURE that it is an AnnealinGA
void SetupAnnealingGA(boost::shared_ptr<GeneralGA> &GA, const double inittemp,
    const int annealinggeneration, const double coolingratio)
  {
    // we get the temperature from the configuration file
    double InitTemperature = inittemp;
    // if it is negative we determine it by performing one iteration
    if (inittemp < 0)
      {
        double avgcombfit = 0;
        GA->DoIteration(0, false);
        //we calculate the average fit across all objective functions
        avgcombfit = accumulate(GA->GetAvgFit().begin(), GA->GetAvgFit().end(),
            0.0);
        // for negative temperatures we calculate the actual temperature from the average fit
        InitTemperature = avgcombfit * abs(inittemp);
      }
    AnnealingGA *AnnealGA = dynamic_cast<AnnealingGA*> (GA.get());
    //we set the annealing parameters for the GA
    AnnealGA->SetParams(InitTemperature, annealinggeneration, coolingratio);
  }

//! Use the information stored in recinforfile to setup the multi receiver function objective function
void SetupRecObjective(const string recinfofile,
    Multi1DRecObjective &RecObjective, const bool normrec,
    RecCalc::trfmethod rfmethod)
  {
    ifstream recinfo(recinfofile.c_str());
    double slowness, sigma, cc, recerror, recweight, absweight;
    int shift;
    string recinputdata, absinputdata;
    char wantabsvel; // do we want absolute velocities as well
    //as long as we can read from recinfofile
    while (recinfo.good())
      {
        //assume we do not want absolute velocity information
        wantabsvel = 'r';
        //read the mandatory values from the file
        recinfo >> slowness >> sigma >> cc >> recerror >> shift >> wantabsvel
            >> recinputdata;
        // if we want absolute velocities we have to read in additional information
        if (wantabsvel == 'a')
          {
            recinfo >> recweight >> absweight >> absinputdata;
          }
        //if reading all values succeeded
        if (recinfo.good())
          {
            //create a shared pointer to a new receiver function data object
            boost::shared_ptr<const SeismicDataComp> RecData(
                new SeismicDataComp(recinputdata, SeismicDataComp::sac));
            //if we need absolute velocity information
            if (wantabsvel == 'a')
              {
                //read in data from the input file
                SurfaceWaveData AbsVelData;
                AbsVelData.ReadAscii(absinputdata);
                // add an objective function that fits receiver functions and absolute velocities
                RecObjective.AddAbsVelFunction(RecData, AbsVelData, shift,
                    sigma, cc, slowness, rfmethod, recerror, normrec,
                    absweight, recweight);
              }
            //if we don't want absolute velocity information
            else
              {
                ResPkModel::WaveType InWave = ResPkModel::PWave;
                if (wantabsvel == 's')
                  {
                    InWave = ResPkModel::SWave;
                  }
                //add an objective function that only considers standard receiver functions
                RecObjective.AddRecFunction(RecData, shift, sigma, cc,
                    slowness, rfmethod, recerror, normrec, InWave);
              }
          }
      }
  }

int main(int argc, char *argv[])
  {
    //we report the total run time at the end so we record when we started
    boost::posix_time::ptime starttime =
        boost::posix_time::microsec_clock::local_time();
    //output some information about the program
    string version = "$Id: 1dinvga.cpp 1888 2015-04-09 13:40:29Z mmoorkamp $";
    cout << endl << endl;
    cout << "Program " << version << endl;
    cout << "Genetic algorithm to jointly invert seismic and MT data " << endl;
    cout << "look at 1dinvga.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;
    RecInvConf RecConf;
    SurfInvConf SurfConf;
    IsoJointConf JointConf;
    GAConf GaConf;

    MTStation MTData;
    SurfaceWaveData RFAbsVel;
    SurfaceWaveData DispData;
    UniformRNG Random;

    std::ifstream conffile("1dinvga.conf");

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

    ttranscribed transcribed;
    tgatype GAtype = pareto;

    string outputbase;
    // we can specify the start of the names of the log files on the command line or in the configuration file
    if (argc > 1)
      outputbase = argv[1];
    else
      outputbase = JointConf.outputbase;
    //initialize the streams for the log files
    //the ones that are always written out get the filename
    ofstream fitstat((outputbase + ".ftst").c_str());
    ofstream rankfile((outputbase + ".rank").c_str());
    ofstream frontfile((outputbase + ".front").c_str());
    //the optional ones are left without a filename
    ofstream misfitfile;
    ofstream valuefile;
    ofstream probfile;
    ofstream distancefile;
    // in most cases we do not need all of these values
    //so we only write them if the user insists
    if (JointConf.verbose)
      {
        misfitfile.open((outputbase + ".msft").c_str());
        valuefile.open((outputbase + ".vals").c_str());
        probfile.open((outputbase + ".probs").c_str());
        distancefile.open((outputbase + ".dist").c_str());
      }
    //we do some very basic consistency checking
    if ((JointConf.thickbase.size() != JointConf.resbase.size())
        || (JointConf.thickbase.size() != JointConf.svelbase.size()))
      {
        cerr << "Configuration of genes is not consistent ! ";
        return 100;
      }
    // we have three different types of inversion parameters
    const int nparam = 3;
    //determine the sizes of a number of GA structures from the configuration file
    const int genes = (JointConf.thickbase.size()
        + JointConf.resbase.size() + JointConf.svelbase.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 the vectors from the congiguration object to the GA
    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());
    //the GA vectors are just concatenations of the different parameters
    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.svelbase.begin(), JointConf.svelbase.end(),
        basevalues.begin() + 2 * paramlength);
    copy(JointConf.svelstep.begin(), JointConf.svelstep.end(),
        stepsizes.begin() + 2 * paramlength);
    copy(JointConf.svelsizes.begin(), JointConf.svelsizes.end(),
        genesizes.begin() + 2 * paramlength);

    const int totalgenesize = accumulate(genesizes.begin(), genesizes.end(), 0);
    //create the objects for the genetic algorithm
    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);

    //read in data
    try
      {
        boost::shared_ptr<Multi1DRecObjective> RecObjective(
            new Multi1DRecObjective());
        //we only read the data that has positive weights in the inversion
        if (JointConf.weights.at(0) > 0)
          {
            MTData.GetData(MTConf.mtinputdata);
          }
        if (JointConf.weights.at(1) > 0)
          {
            RecCalc::trfmethod rfmethod = RecCalc::specdiv;
            if (RecConf.recmethod == "iterdecon")
              {
                rfmethod = RecCalc::iterdecon;
              }
            //use the information in the recinfofile to initialize the RF objective function
            SetupRecObjective(RecConf.recinfofile, *RecObjective.get(),
                RecConf.normrec, rfmethod);
          }
        if (JointConf.weights.at(2) > 0)
          {
            DispData.ReadFile(SurfConf.dispdata);
          }
        //setup MT Objective function
        boost::shared_ptr<Iso1DMTObjective> MTObjective(new Iso1DMTObjective(
            MTData));
        SetupMTFitParameters(MTConf, *MTObjective);
        MTObjective->SetFitExponent(MTConf.mtfitexponent);

        //setup RF objective function
        RecObjective->SetPoisson(JointConf.poisson);
        RecObjective->SetTimeWindow(RecConf.starttime,
            RecConf.endtime);

        //setup surface wave data
        boost::shared_ptr<SurfaceWaveObjective> SurfObjective(
            new SurfaceWaveObjective(DispData));
        SurfObjective->SetFitExponent(SurfConf.surffitexponent);
        SurfObjective->SetPoisson(JointConf.poisson);
        SurfObjective->SetErrorLevel(SurfConf.surferror);

        //setup Regularization
        boost::shared_ptr<GeneralObjective> RecRegularization;
        //if we don't want a reference model regularization
        // we regularize by smoothness over resistivity and velocity
        if (!JointConf.usevrefmodel)
          {
            RecRegularization = boost::shared_ptr<GeneralObjective>(
                new CombinedRoughness);
          }
        //otherwise we regularize with a reference model
        else
          {
            //we check if the file exists
            if (!boost::filesystem::exists(JointConf.vrefmodel))
              {
                cerr
                    << "Cannot find reference model for regularization ! Exiting !"
                    << endl;
                return 100;
              }
            //get the model from the file
            ResPkModel VRefModel;
            VRefModel.ReadModel(JointConf.vrefmodel);
            // and initialize the regularization object
            RecRegularization = boost::shared_ptr<GeneralObjective>(
                new SeismicModelDiff(VRefModel));
          }
        //MT is always regularized by smoothness
        //this can create some redundancy and should be changed in the future
        boost::shared_ptr<GeneralObjective> MTRegularization(new MTRoughness());
        //create the vector of objective functions for the GA
        vector<pCGeneralObjective> ObjVector;
        ObjVector.push_back(MTObjective);
        ObjVector.push_back(RecObjective);
        ObjVector.push_back(SurfObjective);
        ObjVector.push_back(RecRegularization);
        ObjVector.push_back(MTRegularization);

        boost::shared_ptr<GeneralGA> GA;
        const size_t nthreads =1;
        if (GaConf.gatype == "pareto")
          {
            GAtype = pareto;
            boost::shared_ptr<ParetoGA> Pareto(new ParetoGA(&Propagator,
                &Population, &Transcribe, ObjVector, nthreads));
            GA = Pareto;
          }
        else if (GaConf.gatype == "anneal")
          {
            GAtype = anneal;
            boost::shared_ptr<AnnealingGA> Annealing(new AnnealingGA(
                &Propagator, &Population, &Transcribe, ObjVector,
                nthreads));
            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 = 2 * JointConf.resbase.size();
        std::vector<int> LocalIndices;
        //generate Indices 0 .. nmtparam for MT Objective function
        generate_n(back_inserter(LocalIndices), nmtparam, IntSequence(0));
        ParamIndices.push_back(LocalIndices);
        const int nseisparam = 2 * JointConf.thickbase.size();
        LocalIndices.clear();
        // generate Indice thickbase .. thickbase+nseisparam for rf and surface wave data
        generate_n(back_inserter(LocalIndices), nseisparam, IntSequence(
            JointConf.thickbase.size()));
        ParamIndices.push_back(LocalIndices);
        ParamIndices.push_back(LocalIndices);
        const int nmodelparam = 3 * JointConf.thickbase.size();
        LocalIndices.clear();
        //generate Indices 0 .. nmodelparam for Regularization Objective function
        generate_n(back_inserter(LocalIndices), nmodelparam, IntSequence(0));
        ParamIndices.push_back(LocalIndices);
        ParamIndices.push_back(LocalIndices); // we push back twice, because both Regularization functionals work with all parameters
        GA->SetParameterIndices(ParamIndices);
        GA->SetElitist(GaConf.elitist);

        if (GAtype == anneal)
          {
            SetupAnnealingGA(GA,GaConf.inittemp,GaConf.annealinggeneration,GaConf.coolingratio );
            Population.InitPop();
          }
        //Population.PrintPopulation(cout);
        for (int i = 0; i < maxgen; ++i)
          {
            GA->DoIteration(i, i == (maxgen - 1));
            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);
            if (WantAbort()) //if we want to exit prematurely
              {
                RemoveAbort();
                // set counter to two before the last generation
                // this way we do the last iteration for the GA which might have some special treatment
                // and do everything that follows below to store the models and plots
                i = maxgen - 2;
              }
          }

        int bestcount = 0;
        // at the end we only want to output the unique models
        // so we initialize a new UniquePop object that keeps track
        UniquePop FinalUnique;
        tfitvec bestfit(ObjVector.size());
        tpopmember member(Population.GetGenesize());
        //Get the indices of the models in the pareto front
        vector<int> BestIndices(GA->GetBestModelIndices());
        const int noutmodels = GA->GetNBestmodels();
        const int nobjective = ObjVector.size();
        //We have to translate between the parameters in the GA and the parameters
        // for the objective functions
        GeneralGA::tparamvector LocalParameters(nobjective);
        //we allocate memory for each parameter set and objective function
        for (int i = 0; i < nobjective; ++i)
          {
            LocalParameters.at(i).resize(ParamIndices.at(i).size(), false);
          }
        // we write the parameters for the unique final models into a file for overview
        ofstream finalmodels((outputbase + ".final").c_str());
        //for all good models including duplicates
        for (int h = 0; h < noutmodels; ++h)
          {
            //get the genetic string for the current best model
            member = row(Population.GetPopulation(), BestIndices.at(h));
            //translate into numerical values
            transcribed = Transcribe.GetValues(member);
            //copy the right portion of the vector for each objective function
            GA->SetupParams(transcribed, LocalParameters);
            //if we haven't calculated the model before
            if (!FinalUnique.Find(transcribed, bestfit))
              {
                //write out some information
                cout << "Best: " << h << " Index : " << BestIndices.at(h)
                    << " ";
                //recalculate the data, we didn't store this during the inversion
                //we only calculate if we have a weight greater 0
                for (int f = 0; f < nobjective; ++f)
                  {
                    if (JointConf.weights.at(f) > 0)
                      {
                        bestfit(f) = ObjVector.at(f)->CalcPerformance(
                            LocalParameters.at(f));
                      }
                    else
                      {
                        bestfit(f) = 0;
                      }
                    cout << setw(6) << setfill(' ') << setprecision(4)
                        << bestfit(f) << " ";
                  }
                // write out the fit information and the index of the model
                //so we can find the right files when looking at tradeoff etc.
                copy(bestfit.begin(), bestfit.end(), ostream_iterator<double> (
                    finalmodels, " "));
                cout << " " << h;
                finalmodels << " " << bestcount << endl;
                ostringstream filename;
                filename << "best_" << bestcount;
                //save all models, data and plots if the weight for the dataset
                // is different from zero, we have to treat data and regularization differently
                // so we cannot condense it with the for/if combination above that goes over all objectives
                cout << " Saved as : " << filename.str();
                if (JointConf.weights.at(0) > 0)
                  {
                    MTObjective->WriteModel(filename.str() + "_mt.mod");
                    MTObjective->WritePlot(filename.str() + "_mt.plot");
                    MTObjective->WriteData(filename.str());
                  }
                if (JointConf.weights.at(1) > 0)
                  {
                    RecObjective->WriteData(filename.str() + ".rec");
                    RecObjective->WriteModel(filename.str() + "_rec.mod");
                    RecObjective->WritePlot(filename.str() + "_rec");
                  }
                if (JointConf.weights.at(2) > 0)
                  {
                    SurfObjective->WriteData(filename.str() + ".disp");
                    SurfObjective->WriteModel(filename.str() + "_disp.mod");
                    SurfObjective->WritePlot(filename.str() + "_disp.plot");
                  }
                filename.str("");
                //store the model and fit, so we do not recalculate
                FinalUnique.Insert(bestfit, transcribed);
                bestcount++;
                cout << endl;
              }
          }
        //if we want verbose information
        if (JointConf.verbose)
          {
            //print all the good models in raw form to a file
            ofstream bestfile((outputbase + ".best").c_str());
            FinalUnique.PrintAll(bestfile);
            // and print the unique population members in the whole run
            ofstream uniquepop((outputbase + ".unpop").c_str());
            GA->PrintUniquePop(uniquepop);
          }
        //write out some run time information
        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;
      }
  }