mtcovar.cpp

Go to the documentation of this file.

#include "C1DMTObjective.h"
#include "GradFunc.h"
#include "convert.h"
#include "C1DMTSynthData.h"
#include "netcdfcpp.h"
#include <iostream>
#include <fstream>
#include <boost/bind.hpp>
#include "VecMat.h"
#include "Adaptors.h"
#include <algorithm>
#include <string>
#include "NetCDFTools.h"
#include "MTFitSetup.h"
#include "C1dInvGaConf.h"
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);
    SetupFitParameters(Configuration, MTObjective);

    string modelnamebase;
    cout << "Model name base: ";
    cin >> modelnamebase;
    C1DMTSynthData MTModel;
    MTModel.readmodel(modelnamebase + "_mt.mod");

    gplib::rvec MTModelVector(2 * MTModel.GetThicknesses().size());
    for (int i = 0; i < MTModel.GetThicknesses().size(); ++i)
      {
        MTModelVector(i) = log10(MTModel.GetResistivities().at(i));
        MTModelVector(i + MTModel.GetThicknesses().size())
            = MTModel.GetThicknesses().at(i);
      }

    double MTOldRMS = MTObjective.CalcPerformance(MTModelVector);
    cout << "RMS: " << MTOldRMS << endl;
    gplib::rvec OldData(MTObjective.GetSynthData());

    ModelAnalysis MTAnalyser(MTObjective);
    MTAnalyser.SetModel(MTModelVector);
    ofstream evals("mtevals.plot");
    for (int i = 0; i < MTAnalyser.GetSingularValues().size(); ++i)
      evals << i << " " << MTAnalyser.GetSingularValues()(i) << endl;
    cout << "MT Singular Values : " << MTAnalyser.GetSingularValues() << endl;
    cout << "Relative SVD Threshold: ";
    double svdthreshold = 0.0;
    cin >> svdthreshold;

    MTAnalyser.SetRelativeSVDThreshold(svdthreshold);

    gplib::rmat DataCovar(2 * MTData.GetMTData().size(), 2
        * MTData.GetMTData().size());
    fill_n(DataCovar.data().begin(), DataCovar.size1() * DataCovar.size2(), 0.0);
    const size_t length = MTData.GetMTData().size();
    for (int i = 0; i < length; ++i)
      {
        DataCovar(i, i) = MTData.GetMTData().at(i).GetdBerd()
            * MTData.GetMTData().at(i).GetdBerd();
        DataCovar(i + length, i + length) = MTData.GetMTData().at(i).GetdBerd()
            * MTData.GetMTData().at(i).GetdBerd();
      }
    gplib::rmat ModelCovar = MTAnalyser.GetModelCovariance(DataCovar);
    WriteMatrixAsNetCDF("mtcovar.nc", ModelCovar);
    WriteMatrixAsNetCDF("mtsens.nc", MTAnalyser.GetSensitivityMatrix());
    ofstream mterrors("mterrors");
    for (int i = 0; i < MTModel.GetThicknesses().size(); ++i)
      {
        //mterrors << sqrt(2.3025 * pow(10,MTModel.GetResistivities().at(i)) * ModelCovar(i,i)) << endl;
        mterrors << sqrt(ModelCovar(i, i)) << endl;
      }
    for (int i = MTModel.GetThicknesses().size(); i < 2
        * MTModel.GetThicknesses().size(); ++i)
      mterrors << sqrt(ModelCovar(i, i)) << endl;

  }