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 && std::abs(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 && std::abs(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());

  }