GPLIB++
perturbmt.cpp
Go to the documentation of this file.
1 #include "GradFunc.h"
2 #include "UniformRNG.h"
3 #include "convert.h"
4 #include "C1DMTSynthData.h"
5 #include "C1DMTObjective.h"
6 #include "C1dInvGaConf.h"
7 #include "NetCDFTools.h"
8 #include "netcdfcpp.h"
9 #include "MTFitSetup.h"
10 #include <iostream>
11 #include <fstream>
12 #include <boost/bind.hpp>
13 #include "VecMat.h"
14 #include <algorithm>
15 #include <string>
16 
17 using namespace std;
18 
19 int main()
20  {
21  C1dInvGaConf Configuration;
22  Configuration.GetData("1dinvga.conf");
23 
24  CMTStation MTData;
25  string datanamebase;
26  cout << "Data name base: ";
27  cin >> datanamebase;
28  MTData.GetData(datanamebase + ".mtt");
29 
30  C1DMTObjective MTObjective(MTData);
31  SetupMTFitParameters(Configuration, MTObjective);
32  string modelnamebase;
33  cout << "Model name base: ";
34  cin >> modelnamebase;
35  C1DMTSynthData MTModel;
36  MTModel.readmodel(modelnamebase + "_mt.mod");
37 
38  ublas::vector<double> MTModelVector = MTModel.GetModelVector();
39 
40  double MTOldRMS = MTObjective.CalcPerformance(MTModelVector);
41  ublas::vector<double> OldData(MTObjective.GetSynthData());
42  ModelAnalysis MTAnalyser(MTObjective);
43  MTAnalyser.SetModel(MTModelVector);
44 
45  cout << "MT Singular Values : " << MTAnalyser.GetSingularValues() << endl;
46  //cout << "MT Analysis R: " << MTAnalyser.GetModelResolution() << endl;
47  //cout << "MT Analysis Null: " << MTAnalyser.GetModelNullspace() << endl;
48  cout << "Relative SVD Threshold: ";
49  double svdthreshold = 0.0;
50  cin >> svdthreshold;
51  MTAnalyser.SetRelativeSVDThreshold(svdthreshold);
52  ublas::vector<double> MTProjection(MTModelVector),
53  MTPerturb(MTModelVector), MTNewModel(MTModelVector);
54 
55  UniformRNG Random;
56  double MTNewRMS = MTOldRMS;
57  cout << "Original RMS : " << MTOldRMS << endl;
58  const int maxit = 20;
59  int currit = 0;
60  MTPerturb *= 0.0;
61  for (int i = 0; i < MTPerturb.size(); ++i)
62  MTPerturb(i) = Random.GetNumber(-1., 1.);
63  MTNewModel = MTModelVector;
64  while (std::abs(MTNewRMS - MTOldRMS) / MTOldRMS < 0.05 && currit < maxit)
65  {
66  MTProjection = prod(MTAnalyser.GetModelNullspace(), MTPerturb);
67  MTNewModel += MTProjection;
68  MTNewRMS = MTObjective.CalcPerformance(MTNewModel);
69  cout << "MT New RMS: " << MTNewRMS << endl;
70  ++currit;
71  }
72 
73  MTObjective.WriteData("new.mtt");
74 
75  MTObjective.WritePlot("new_mt.plot");
76 
77  MTObjective.WriteModel("new_mt.mod");
78  cout << "MT OldRMS: " << MTOldRMS << " MT NewRMS: " << MTNewRMS << endl;
79 
80  cout << "MT Old Model: " << MTModelVector << endl;
81  cout << "MT New Model: " << MTNewModel << endl;
82 
83  WriteMatrixAsNetCDF("mtres.nc", MTAnalyser.GetModelResolution());
84  WriteMatrixAsNetCDF("mtnull.nc", MTAnalyser.GetModelNullspace());
85  }
void SetupMTFitParameters(const tConfObject &Configuration, C1DMTObjective &Objective)
Definition: MTFitSetup.h:8
int main()
Definition: perturbmt.cpp:19
boost::shared_ptr< PlottableObjective > MTObjective
Definition: cadianiso.cpp:16
CMTStation MTData
Definition: cadijoint.cpp:15
The class ModelAnalysis is used to calculate resolution matrix, nullspace and other parameters for mo...
Definition: GradFunc.h:20
CLevanisoConf Configuration
Definition: cadianiso.cpp:17