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 (fabs(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());
}

---