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>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>
 <TR> <TD>recweight</TD> <TD> </TD> <TD>The receiver function class has the possibility to calculate absolute velocity information from the receiver function.
 This and the following weight determine in how far the traditional receiver function and the absolute velocity transformation determine the misfit
 for the receiver function. </TD></TR>
 <TR> <TD>absvelweight</TD> <TD> </TD> <TD> See above </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$\Omega_0\f$ for gaussian filter </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 anytime 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 <fstream>
#include <algorithm>
#include <numeric>
#include <sstream>
#include <string>
#include "GA.h"
#include "Iso1DMTObjective.h"
#include "MTStation.h"
#include "C1dInvGaConf.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 <boost/bind.hpp>
#include <boost/shared_ptr.hpp>
#include <boost/filesystem.hpp>
#include <boost/date_time/posix_time/posix_time.hpp>
#include "Util.h"
#include "MTRoughness.h"
#include "SurfaceWaveData.h"
#include "FileAbort.h"
#include "SeisTools.h"
enum tgatype
  { pareto,anneal};

typedef boost::shared_ptr<GeneralObjective> pCGeneralObjective;
using namespace std;

void SetupAnnealingGA(boost::shared_ptr<GeneralGA> &GA,
    const C1dInvGaConf 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 SetupRecObjective(const string recinfofile,
    Multi1DRecObjective &RecObjective)
  {
    ifstream recinfo(recinfofile.c_str());
    double slowness, omega, sigma, cc, recerror, recweight, absweight;
    int shift;
    string recinputdata, absinputdata;
    char wantabsvel; // do we want absolute velocities as well
    while (recinfo.good())
      {
        wantabsvel = 'r';
        recinfo >> slowness >> omega >> sigma >> cc >> recerror >> shift
            >> wantabsvel >> recinputdata;
        if (wantabsvel == 'a') // if we want absolute velocities we have to read in additional information
          {
            recinfo >> recweight >> absweight >> absinputdata;
          }
        if (recinfo.good())
          {
            boost::shared_ptr<const SeismicDataComp>
                RecData(new SeismicDataComp(recinputdata, SeismicDataComp::sac));
            Normalize(const_cast<SeismicDataComp*>(RecData.get())->GetData());
            if (wantabsvel == 'a')
              {
                SurfaceWaveData AbsVelData;
                AbsVelData.ReadAscii(absinputdata);
                RecObjective.AddAbsVelFunction(RecData, AbsVelData, shift,
                    omega, sigma, cc, slowness, RecCalc::iterdecon, recerror,
                    absweight, recweight);
              }
            else
              {
                RecObjective.AddRecFunction(RecData, shift, omega, sigma, cc,
                    slowness, RecCalc::iterdecon, recerror);
              }
          }
      }
  }

int main(int argc, char *argv[])
  {
    boost::posix_time::ptime starttime =
        boost::posix_time::microsec_clock::local_time();
    string version = "$Id: 1dinvga.cpp 1696 2008-05-20 12:08:38Z 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;

    C1dInvGaConf Configuration;
    MTStation MTData;
    //SeismicDataComp RecData;
    SurfaceWaveData RFAbsVel;
    SurfaceWaveData DispData;
    UniformRNG Random;
    Configuration.GetData("1dinvga.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.resbase.size())
        || (Configuration.thickbase.size() != Configuration.svelbase.size()))
      {
        cerr << "Configuration of genes is not consistent ! ";
        return 100;
      }
    const int nparam = 3;
    const int genes = (Configuration.thickbase.size()
        + Configuration.resbase.size() + Configuration.svelbase.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.resbase.begin(), Configuration.resbase.end(),
        basevalues.begin());
    copy(Configuration.resstep.begin(), Configuration.resstep.end(),
        stepsizes.begin());
    copy(Configuration.ressizes.begin(), Configuration.ressizes.end(),
        genesizes.begin());

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

    copy(Configuration.svelbase.begin(), Configuration.svelbase.end(),
        basevalues.begin()+2*paramlength);
    copy(Configuration.svelstep.begin(), Configuration.svelstep.end(),
        stepsizes.begin()+2*paramlength);
    copy(Configuration.svelsizes.begin(), Configuration.svelsizes.end(),
        genesizes.begin()+2*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);

    //read in data
    try
      {
        boost::shared_ptr<Multi1DRecObjective> RecObjective(new Multi1DRecObjective());
        //we only read the data that has positive weights in the inversion
        if (Configuration.weights.at(0)> 0)
          {
            MTData.GetData(Configuration.mtinputdata);
          }
        if (Configuration.weights.at(1)> 0)
          {
            SetupRecObjective(Configuration.recinfofile,*RecObjective.get());
          }
        if (Configuration.weights.at(2)> 0)
          {
            DispData.ReadFile(Configuration.dispdata);
          }
        //setup MT Objective function
        boost::shared_ptr<Iso1DMTObjective> MTObjective(new Iso1DMTObjective(MTData));
        SetupMTFitParameters(Configuration,*MTObjective);
        MTObjective->SetFitExponent(Configuration.mtfitexponent);

        //setup RF objective function
        RecCalc::trfmethod rfmethod = RecCalc::specdiv;
        if (Configuration.recmethod == "iterdecon")
          {
            rfmethod = RecCalc::iterdecon;
          }

        RecObjective->SetPoisson(Configuration.poisson);
        RecObjective->SetTimeWindow(Configuration.starttime,Configuration.endtime);

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

        //setup Regularization
        boost::shared_ptr<GeneralObjective> RecRegularization;
        if (!Configuration.usevrefmodel)
          {
            RecRegularization = boost::shared_ptr<GeneralObjective>(new CombinedRoughness);
          }
        else
          {
            if (!boost::filesystem::exists(Configuration.vrefmodel))
              {
                cerr << "Cannot find reference model for regularization ! Exiting !" << endl;
                return 100;
              }
            ResPkModel VRefModel;
            VRefModel.GetData(Configuration.vrefmodel);
            RecRegularization = boost::shared_ptr<GeneralObjective>(new SeismicModelDiff(VRefModel));
          }
        boost::shared_ptr<GeneralObjective> MTRegularization (new MTRoughness());
        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;
        if (Configuration.gatype == "pareto")
          {
            GAtype = pareto;
            boost::shared_ptr<ParetoGA> Pareto(new ParetoGA(&Propagator,&Population,&Transcribe,ObjVector));
            GA = Pareto;
          }
        else if (Configuration.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(Configuration.weights);
        //Setting up indices for forward modelling
        GeneralGA::tparamindv ParamIndices;
        const int nmtparam = 2 * Configuration.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 * Configuration.thickbase.size();
        LocalIndices.clear();
        // generate Indice thickbase .. thickbase+nseisparam for rf and surface wave data
        generate_n(back_inserter(LocalIndices),nseisparam,IntSequence(Configuration.thickbase.size()));
        ParamIndices.push_back(LocalIndices);
        ParamIndices.push_back(LocalIndices);
        const int nmodelparam = 3 * Configuration.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->SetThreads(Configuration.threads);
        GA->SetElitist(Configuration.elitist);

        if (GAtype == anneal)
          {
            SetupAnnealingGA(GA,Configuration);
            Population.InitPop();
          }
        //Population.PrintPopulation(cout);
        for (int i = 0; i < maxgen; ++i)
          {
            GA->DoIteration(i, i == (maxgen-1));
            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);
            if (WantAbort()) //if we want to exit prematurely

              {
                RemoveAbort();
                i = maxgen; // set counter to last generation, so we do everything that follows
              }
          }

        int bestcount = 0;
        UniquePop FinalUnique;
        tfitvec bestfit(ObjVector.size());
        tpopmember member(Population.GetGenesize());

        vector<int> BestIndices(GA->GetBestModelIndices());
        const int noutmodels = GA->GetNBestmodels();
        const int nobjective = ObjVector.size();
        GeneralGA::tparamvector LocalParameters(nobjective);
        for (int i = 0; i < nobjective; ++i)
          {
            LocalParameters.at(i).resize(ParamIndices.at(i).size(),false);
          }
        ofstream finalmodels((outputbase+".final").c_str());
        for (int h= 0; h < noutmodels; ++h)
          {
            member = row(Population.GetPopulation(),BestIndices.at(h));
            cout << "Best: " << h << " Index : " << BestIndices.at(h) << " ";

            transcribed = Transcribe.GetValues(member);
            GA->SetupParams(transcribed,LocalParameters);
            for (int f = 0; f < nobjective; ++f)
              {
                if (Configuration.weights.at(f)> 0)
                  {
                    bestfit(f) = ObjVector.at(f)->CalcPerformance(LocalParameters.at(f));
                  }
                else
                  {
                    bestfit(f) = 0;
                  }
                cout << bestfit(f) << " ";
              }

            if (FinalUnique.Insert(bestfit,transcribed))
              {
                copy(bestfit.begin(),bestfit.end(),ostream_iterator<double>(finalmodels," "));
                cout << " " << h;
                finalmodels << " " << bestcount << endl;
                ostringstream filename;
                filename << "best_" << bestcount;
                cout << " Saved as : " << filename.str();
                if (Configuration.weights.at(0)> 0)
                  {
                    MTObjective->WriteModel(filename.str()+"_mt.mod");
                    MTObjective->WritePlot(filename.str()+"_mt.plot");
                    MTObjective->WriteData(filename.str());
                  }
                if (Configuration.weights.at(1)> 0)
                  {
                    RecObjective->WriteData(filename.str()+".rec");
                    RecObjective->WriteModel(filename.str()+"_rec.mod");
                    RecObjective->WritePlot(filename.str()+"_rec");
                  }
                if (Configuration.weights.at(2)> 0)
                  {
                    SurfObjective->WriteData(filename.str()+".disp");
                    SurfObjective->WriteModel(filename.str()+"_disp.mod");
                    SurfObjective->WritePlot(filename.str()+"_disp.plot");
                  }
                filename.str("");
                bestcount++;
              }
            cout << endl;
          }
        FinalUnique.PrintAll(bestfile);
        GA->PrintUniquePop(uniquepop);
        boost::posix_time::ptime endtime = boost::posix_time::microsec_clock::local_time();
        std::cout << "Total Runtime: " << endtime - starttime << std::endl;
      }
    catch(CFatalException &e)
      {
        cerr << e.what() << endl;
        return -1;
      }
  }

---