mtanisoenum.cpp

Go to the documentation of this file.

#include <iostream>
#include <fstream>
#include <algorithm>
#include <numeric>
#include <sstream>
#include <string>
#include "UniformRNG.h"
#include "gentypes.h"
#include "PTensor1DMTObjective.h"
#include "Aniso1DMTObjective.h"
#include "PTensorMTStation.h"
#include "MTStation.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 <complex>
#include "Util.h"
#include "MTAnisoRoughness.h"
#include "MTAnisoGAConf.h"
#include "convert.h"
enum tgatype {pareto,anneal};
const int nobjective = 2;
typedef boost::shared_ptr<GeneralObjective> pCGeneralObjective;
using namespace std;



int main(int argc, char *argv[])
{
        string version =  "$Id: 1dinvga.cpp 1272 2007-07-03 23:17:36Z max $";
        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;
        
        MTAnisoGAConf Configuration;
        
        try 
        {
                Configuration.GetData("mtanisoga.conf");
        

        string outputbase;
        
        if (argc > 1)
                outputbase = argv[1];
        else
                outputbase = Configuration.outputbase;
        
        
                        
        if ( (Configuration.thickbase.size() != Configuration.resbase.size()) || 
                (Configuration.thickbase.size() != Configuration.anisobase.size()) ||
                (Configuration.thickbase.size() != Configuration.strikebase.size()))
                {
                        cerr << "Configuration of genes is not consistent ! ";
                        return 100;
                }
        
        string modelfile;
        cout << "Model file name: ";
        cin >> modelfile;
        
        C1DAnisoMTSynthData StartModel;
        StartModel.ReadModel(modelfile);
        ttranscribed origmodel = StartModel.GetModelVector();
        const unsigned int nmodelparameters = origmodel.size();
        const unsigned int nsections = 4;
        const unsigned int nparmspersection = nmodelparameters/4;
        ttranscribed  basevalues(nmodelparameters), steps(nmodelparameters),sizes(nmodelparameters);
        
        copy(Configuration.resbase.begin(),Configuration.resbase.end(),basevalues.begin());
        copy(Configuration.thickbase.begin(),Configuration.thickbase.end(),basevalues.begin()+nparmspersection);
        copy(Configuration.anisobase.begin(),Configuration.anisobase.end(),basevalues.begin()+2*nparmspersection);
        copy(Configuration.strikebase.begin(),Configuration.strikebase.end(),basevalues.begin()+3*nparmspersection);
        
        copy(Configuration.resstep.begin(),Configuration.resstep.end(),steps.begin());
        copy(Configuration.thickstep.begin(),Configuration.thickstep.end(),steps.begin()+nparmspersection);
        copy(Configuration.anisostep.begin(),Configuration.anisostep.end(),steps.begin()+2*nparmspersection);
        copy(Configuration.strikestep.begin(),Configuration.strikestep.end(),steps.begin()+3*nparmspersection);
                
        copy(Configuration.ressizes.begin(),Configuration.ressizes.end(),sizes.begin());
        copy(Configuration.thicksizes.begin(),Configuration.thicksizes.end(),sizes.begin()+nparmspersection);
        copy(Configuration.anisosizes.begin(),Configuration.anisosizes.end(),sizes.begin()+2*nparmspersection);
        copy(Configuration.strikesizes.begin(),Configuration.strikesizes.end(),sizes.begin()+3*nparmspersection);
        //setup MT Objective function
        boost::shared_ptr<PlottableObjective> MTObjective;
                
        if(Configuration.mtfit == "ptensor")
        {
                PTensorMTStation PTData;
                PTData.GetData(Configuration.ptensordata);
                boost::shared_ptr<PTensor1DMTObjective> PTObjective(new PTensor1DMTObjective(PTData));
                MTObjective = PTObjective;
        }
        else
        {
                MTStation MTData;
                MTData.GetData(Configuration.mtinputdata);
                boost::shared_ptr<Aniso1DMTObjective> AnisoObjective(new Aniso1DMTObjective(MTData));
                SetupMTFitParameters(Configuration,*AnisoObjective);
                MTObjective = AnisoObjective;
        }
        
        ofstream misfitfile;
        ofstream gmtfile("plotenum");
        double misfit = 0;
        double delta = 0.1;
        gmtfile << "gmtset PAPER_MEDIA A0" << endl;
        gmtfile << "psbasemap -JX20 -R0/20/0/20 -G254 -K > plotenum.ps" << endl;
        for (unsigned int i= 0; i < nmodelparameters; ++i)
        {
                ttranscribed currmodel = origmodel;
                string misfitfilename = "misfit"+stringify(i);
                misfitfile.open(misfitfilename.c_str());
                currmodel(i) = basevalues(i);
                const unsigned int currsize = pow(2,sizes(i));
                double maxmisfit = -1;
                double minmisfit = 1e50;
                for (unsigned int j = 0; j < currsize; ++j)
                {       
                        misfit = MTObjective->CalcPerformance(currmodel);
                        misfitfile << currmodel(i) << " " << misfit << endl;
                        currmodel(i) += steps(i);
                        maxmisfit = max(maxmisfit,misfit);
                        minmisfit = min(minmisfit,misfit);
                }
                double minx = basevalues(i)-delta;
                double maxx = basevalues(i)+currsize*steps(i)+delta;
                string region = " -R"+stringify(minx)+"/"+stringify(maxx)+"/"+stringify(minmisfit*0.9)+"/"+stringify(maxmisfit*1.1);
                int xticks = ceil((maxx - minx)/10.);
                int yticks = ceil(maxmisfit/10);
                string ticks = " -B"+stringify(xticks)+"/"+stringify(yticks)+" ";
                string shiftstring;
                double xshift = 10;
                double yshift = 10;
                if ((i % nparmspersection) == 0 && i!=0)
                        shiftstring = " -X-"+stringify(xshift*(nparmspersection-1))+" -Y"+stringify(yshift);
                else
                        if (i!=0)
                                shiftstring = " -X"+stringify(xshift);
                gmtfile << "psxy " + misfitfilename + " " + shiftstring+ " -JX7/7l " + region + ticks + " -O -K >> plotenum.ps" << endl;
                misfitfile.close();
        }
        }
        catch(CFatalException &e)
        {
                        cerr << e.what() << endl; 
                        return -1;
                }
}