perturbmodels.cpp

Go to the documentation of this file.

#include "C1DMTObjective.h"
#include "MTRecObjective.h"
#include "SeismicDataComp.h"
#include "GradFunc.h"
#include "UniformRNG.h"
#include "convert.h"
#include "C1DMTSynthData.h"
#include "ResPkModel.h"
#include "netcdfcpp.h"
#include <iostream>
#include <fstream>
#include <boost/bind.hpp>
#include "VecMat.h"
#include <algorithm>
#include <string>
#include "C1dInvGaConf.h"
#include "MTFitSetup.h"
#include "NetCDFTools.h"
using namespace std;


int main()
{
        
        CMTStation MTData;
        SeismicDataComp RecData;
        string datanamebase;
        cout << "Data name base: ";
        cin >> datanamebase;
        RecData.GetData(datanamebase+".rec");
        MTData.GetData(datanamebase+".mtt");
        
        
        C1dInvGaConf Configuration;
        Configuration.GetData("1dinvga.conf");
        
        C1DMTObjective MTObjective(MTData);
        SetupMTFitParameters(Configuration,MTObjective);
        
        C1DRecObjective RecObjective(RecData,Configuration.shift,Configuration.omega,
                                                                Configuration.sigma,Configuration.cc,Configuration.slowness);
        RecObjective.SetPoisson(Configuration.poisson);
        RecObjective.SetTimeWindow(Configuration.starttime,Configuration.endtime);
        RecObjective.SetErrorLevel(Configuration.recerror);
        RecObjective.SetFitExponent(Configuration.recfitexponent);
        string modelnamebase;
        cout << "Model name base: ";
        cin >> modelnamebase;
        ResPkModel RecModel;
        C1DMTSynthData MTModel;
        RecModel.GetData(modelnamebase+"_rec.mod");
        MTModel.readmodel(modelnamebase+"_mt.mod");
        
        gplib::rvec MTModelVector = MTModel.GetModelVector();
        
        
        
        
        double MTOldRMS = MTObjective.CalcPerformance(MTModelVector);
        gplib::rvec OldData(MTObjective.GetSynthData());
        ModelAnalysis MTAnalyser(MTObjective);
        MTAnalyser.SetModel(MTModelVector);
        
        
        cout << "MT Singular Values : " << MTAnalyser.GetSingularValues() << endl;
        cout << "Relative SVD Threshold: ";
        double svdthreshold = 0.0;
        cin >> svdthreshold;
        MTAnalyser.SetRelativeSVDThreshold(svdthreshold);
        
        
        const size_t nrecmodelparms = RecModel.SVelocity.size()*2;
        const size_t nlayers = RecModel.SVelocity.size();
        gplib::rvec RecModelVector(nrecmodelparms);
        for (int i = 0; i < nlayers; ++i)
        {
                RecModelVector(i) = RecModel.Thickness.at(i);
                RecModelVector(i+nlayers) = RecModel.SVelocity.at(i);
        }
        double RecOldRMS = RecObjective.CalcPerformance(RecModelVector);
        RecObjective.WriteData("old");
        cout << "Saved old model " << endl << flush;
        ModelAnalysis RecAnalyser(RecObjective);
        RecAnalyser.SetModel(RecModelVector);
        RecAnalyser.SetRelativeSVDThreshold(svdthreshold);
        
        cout << "Rec Singular Values : " << RecAnalyser.GetSingularValues() << endl;
        //cout << "Rec Analysis R: " << RecAnalyser.GetModelResolution() << endl;
        //cout << "Rec Analysis Null: " << RecAnalyser.GetModelNullspace() << endl;
        
        
        gplib::rvec MTProjection(MTModelVector), MTPerturb(MTModelVector), MTNewModel(MTModelVector);
        gplib::rvec RecProjection(RecModelVector), RecPerturb(RecModelVector), RecNewModel(RecModelVector);
        UniformRNG Random;
        double MTNewRMS = MTOldRMS;
        
        const int maxit = 20;
        int currit = 0;
        ublas::vector<double> Perturb(MTModelVector.size()/2);
        MTPerturb *=0.0;
        RecPerturb *=0.0;
        for (int i = 0; i < Perturb.size(); ++i)
                        Perturb(i) = Random.GetNumber(-1.,1.);
        
        for (size_t i = 0; i < nlayers; ++i)
        {
                MTPerturb(i+nlayers) = Perturb(i);
                RecPerturb(i) = -Perturb(i);
        }

        MTNewModel = MTModelVector;
        const double rmsincrease = 0.05;
        while (currit < maxit &&  fabs(MTNewRMS - MTOldRMS)/MTOldRMS <  rmsincrease )
        {
                for (size_t i = nlayers; i < MTPerturb.size(); ++i)
                        MTPerturb(i) = Random.GetNumber(-1.,1.);
                MTProjection = prod(MTAnalyser.GetModelNullspace(),MTPerturb);
                MTNewModel = MTModelVector + MTProjection;
                MTNewRMS = MTObjective.CalcPerformance(MTNewModel);
                MTObjective.WritePlot("new_mt"+stringify(currit));
                cout << "MT New RMS: " << MTNewRMS << endl;
                ++currit;
        }
        currit = 0;
        double RecNewRMS = RecOldRMS;
        RecNewModel = RecModelVector;
        while (currit < maxit && fabs(RecNewRMS - RecOldRMS)/RecOldRMS <  rmsincrease)
        {
                for (size_t i = 0; i < RecPerturb.size(); ++i)
                        RecPerturb(i) = Random.GetNumber(-1.,1.);
                RecProjection = prod(RecAnalyser.GetModelNullspace(),RecPerturb);
                RecNewModel = RecModelVector + RecProjection;
                //RecNewModel(RecNewModel.size()/2-1) = 1.0;
                RecNewRMS = RecObjective.CalcPerformance(RecNewModel);
                RecObjective.WritePlot("new_rec"+stringify(currit));
                cout << "Rec New RMS: " << RecNewRMS << endl;
                ++currit;
        }
        
        
        RecObjective.WriteData("new.rec");
        MTObjective.WriteData("new.mtt");
        
        MTObjective.WritePlot("new_mt.plot");
        RecObjective.WriteModel("new_rec.mod");
        MTObjective.WriteModel("new_mt.mod");
        cout << "MT OldRMS: " << MTOldRMS <<  " MT NewRMS: " << MTNewRMS << endl;
        cout << "Rec OldRMS: " << RecOldRMS << " Rec NewRMS: " << RecNewRMS << endl;
        cout << "MT Old Model: " << MTModelVector << endl;
        cout << "MT New Model: " << MTNewModel << endl;
        cout << "Rec Old Model: " << RecModelVector << endl;
        cout << "Rec New Model: " << RecNewModel << endl;
        
        WriteMatrixAsNetCDF("mtres.nc",MTAnalyser.GetModelResolution());
        WriteMatrixAsNetCDF("mtnull.nc",MTAnalyser.GetModelNullspace());
        WriteMatrixAsNetCDF("recres.nc",RecAnalyser.GetModelResolution());
        WriteMatrixAsNetCDF("recnull.nc",RecAnalyser.GetModelNullspace());
        
}

---