mtanalys.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"
#include "NumUtil.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);
    SetupMTFitParameters(Configuration, MTObjective);

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

    gplib::rvec MTModelVector = MTModel.GetModelVector();
    const size_t nmodelparms = MTModelVector.size();
    gplib::rvec varfactor(nmodelparms);
    double MTOldRMS = MTObjective.CalcPerformance(MTModelVector);
    cout << "RMS: " << MTOldRMS << endl;
    gplib::rvec OldData(MTObjective.GetSynthData());
    for (size_t i = 0; i < nmodelparms; ++i)
      {
        varfactor(i) = 1.0;
      }
    ModelAnalysis MTAnalyser(MTObjective);
    MTAnalyser.SetModel(MTModelVector);
    MTAnalyser.SetVariationFactor(varfactor);
    ofstream evals("mtevals.plot");
    gplib::rvec svals = MTAnalyser.GetSingularValues();
    gplib::rvec svalratios = svals;
    svalratios /= svals(0); // determine ratios to larges eigenvalue
    cout << "MT Singular Values : " << svals << endl;
    //cout << "Relative SVD Threshold: ";
    //double svdthreshold = 0.0;
    //cin >> svdthreshold;

    gplib::rvec modelerrors(nmodelparms);
    modelerrors *= 0.0;
    const double errorratio = 0.1;
    gplib::rvec parameterthresholds = MTModelVector;
    for (size_t i = 0; i < nmodelparms / 2; ++i)
      parameterthresholds(i) = pow(10, parameterthresholds(i));
    parameterthresholds *= errorratio;
    size_t eigenvindex = 0;

    C1DMTObjective::datafuncvector_t ErrorFunctions =
        MTObjective.GetErrorFunctions();
    const size_t nerrorfuncs = ErrorFunctions.size();
    const size_t ndata = MTData.GetMTData().size();
    gplib::rmat DataCovar(nerrorfuncs * ndata, nerrorfuncs * ndata);
    fill_n(DataCovar.data().begin(), DataCovar.size1() * DataCovar.size2(), 0.0);
    ;
    for (size_t e = 0; e < nerrorfuncs; ++e)
      {
        for (int i = 0; i < ndata; ++i)
          {
            DataCovar(i + e * ndata, i + e * ndata) = gplib::pow2(
                ErrorFunctions.at(e)(&MTData.GetMTData().at(i)));
          }
      }
    ofstream mterrors("mterrors");
    gplib::rmat ModelCovar = MTAnalyser.GetModelCovariance(DataCovar);
    std::vector<double> reserr(nmodelparms / 2, 0.0), thickerr(nmodelparms / 2,
        0.0);
    const unsigned int nevals = svalratios.size();
    while (sum(modelerrors) < sum(parameterthresholds) && eigenvindex < nevals)
      {
        cout << "Accum Model Errors: " << sum(modelerrors) << endl;
        cout << "Accum Threshold Errors: " << sum(parameterthresholds) << endl;
        MTAnalyser.SetRelativeSVDThreshold(svalratios(eigenvindex));
        ModelCovar = MTAnalyser.GetModelCovariance(DataCovar);
        double currerr = 0;
        for (size_t i = 0; i < MTModel.GetThicknesses().size(); ++i)
          {
            currerr = 2.3025 * MTModel.GetResistivities().at(i) * sqrt(
                ModelCovar(i, i));
            mterrors << "Res " << i << " : " << currerr << endl;
            reserr.at(i) = currerr;
            modelerrors(i) = currerr;
          }
        for (size_t i = MTModel.GetThicknesses().size(); i < 2
            * MTModel.GetThicknesses().size(); ++i)
          {
            currerr = sqrt(ModelCovar(i, i));
            mterrors << "Th " << i - MTModel.GetThicknesses().size() << " : "
                << currerr << endl;
            modelerrors(i) = currerr;
            thickerr.at(i - MTModel.GetThicknesses().size()) = currerr;
          }
        ++eigenvindex;
      }
    cout << "Final Accum Model Errors: " << sum(modelerrors) << endl;
    cout << "Accum Threshold Errors: " << sum(parameterthresholds) << endl;
    WriteMatrixAsNetCDF("mtcovar.nc", ModelCovar);
    WriteMatrixAsNetCDF("mtcorr.nc", MTAnalyser.GetModelCorrelation(DataCovar));
    WriteMatrixAsNetCDF("mtsens.nc", MTAnalyser.GetSensitivityMatrix());
    WriteMatrixAsNetCDF("mtres.nc", MTAnalyser.GetModelResolution());
    WriteMatrixAsNetCDF("mtnull.nc", MTAnalyser.GetModelNullspace());
    WriteMatrixAsNetCDF("mtsinv.nc", MTAnalyser.GetSensitivityInverse());
    WriteMatrixAsNetCDF("senscheck.nc", prod(
        MTAnalyser.GetSensitivityInverse(), MTAnalyser.GetSensitivityMatrix()));

    MTModel.SetResistivityErrors(reserr);
    MTModel.SetThicknessErrors(thickerr);
    MTModel.writeplot("mterr.plot");
  }