perturbmt.cpp

Go to the documentation of this file.

#include "GradFunc.h"
#include "UniformRNG.h"
#include "convert.h"
#include "C1DMTSynthData.h"
#include "C1DMTObjective.h"
#include "C1dInvGaConf.h"
#include "NetCDFTools.h"
#include "netcdfcpp.h"
#include "MTFitSetup.h"
#include <iostream>
#include <fstream>
#include <boost/bind.hpp>
#include "VecMat.h"
#include <algorithm>
#include <string>

using namespace std;

int main()
  {
    C1dInvGaConf Configuration;
    Configuration.GetData("1dinvga.conf");

    CMTStation MTData;
    string datanamebase;
    cout << "Data name base: ";
    cin >> datanamebase;
    MTData.GetData(datanamebase + ".mtt");

    C1DMTObjective MTObjective(MTData);
    SetupMTFitParameters(Configuration, MTObjective);
    string modelnamebase;
    cout << "Model name base: ";
    cin >> modelnamebase;
    C1DMTSynthData MTModel;
    MTModel.readmodel(modelnamebase + "_mt.mod");

    ublas::vector<double> MTModelVector = MTModel.GetModelVector();

    double MTOldRMS = MTObjective.CalcPerformance(MTModelVector);
    ublas::vector<double> OldData(MTObjective.GetSynthData());
    ModelAnalysis MTAnalyser(MTObjective);
    MTAnalyser.SetModel(MTModelVector);

    cout << "MT Singular Values : " << MTAnalyser.GetSingularValues() << endl;
    //cout << "MT Analysis R: " << MTAnalyser.GetModelResolution() << endl;
    //cout << "MT Analysis Null: " << MTAnalyser.GetModelNullspace() << endl;
    cout << "Relative SVD Threshold: ";
    double svdthreshold = 0.0;
    cin >> svdthreshold;
    MTAnalyser.SetRelativeSVDThreshold(svdthreshold);
    ublas::vector<double> MTProjection(MTModelVector),
        MTPerturb(MTModelVector), MTNewModel(MTModelVector);

    UniformRNG Random;
    double MTNewRMS = MTOldRMS;
    cout << "Original RMS : " << MTOldRMS << endl;
    const int maxit = 20;
    int currit = 0;
    MTPerturb *= 0.0;
    for (int i = 0; i < MTPerturb.size(); ++i)
      MTPerturb(i) = Random.GetNumber(-1., 1.);
    MTNewModel = MTModelVector;
    while (std::abs(MTNewRMS - MTOldRMS) / MTOldRMS < 0.05 && currit < maxit)
      {
        MTProjection = prod(MTAnalyser.GetModelNullspace(), MTPerturb);
        MTNewModel += MTProjection;
        MTNewRMS = MTObjective.CalcPerformance(MTNewModel);
        cout << "MT New RMS: " << MTNewRMS << endl;
        ++currit;
      }

    MTObjective.WriteData("new.mtt");

    MTObjective.WritePlot("new_mt.plot");

    MTObjective.WriteModel("new_mt.mod");
    cout << "MT OldRMS: " << MTOldRMS << " MT NewRMS: " << MTNewRMS << endl;

    cout << "MT Old Model: " << MTModelVector << endl;
    cout << "MT New Model: " << MTNewModel << endl;

    WriteMatrixAsNetCDF("mtres.nc", MTAnalyser.GetModelResolution());
    WriteMatrixAsNetCDF("mtnull.nc", MTAnalyser.GetModelNullspace());
  }