jointanalys.cpp

Go to the documentation of this file.

#include "MTRecObjective.h"
#include "GradFunc.h"
#include "convert.h"
#include "SeismicDataComp.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 "ResPkModel.h"
#include "C1dInvGaConf.h"
#include "MatNorm.h"
#include "MTFitSetup.h"
#include "NumUtil.h"

using namespace std;

void DetermineEvalThreshold(ModelAnalysis &JointAnalyser,
    C1DMTSynthData &MTModel, gplib::rmat &DataCovar, gplib::rvec &modelerrors,
    gplib::rvec &parameterthresholds)
  {
    const size_t nmodelparms = parameterthresholds.size();
    const size_t nlayers = nmodelparms / 4;
    size_t eigenvindex = 0;
    gplib::rvec svals = JointAnalyser.GetSingularValues();
    gplib::rvec svalratios = svals;
    svalratios /= svals(0); // determine ratios to larges eigenvalue
    cout << "Joint Singular Values : " << svals << endl;
    gplib::rmat ModelCovar;
    cout << "Model errors: " << modelerrors << endl;
    while (sum(modelerrors) < sum(parameterthresholds))
      {
        cout << "Ratio: " << svalratios(eigenvindex) << endl;
        JointAnalyser.SetRelativeSVDThreshold(svalratios(eigenvindex));
        cout << "Accum Model Errors: " << sum(modelerrors) << endl;
        cout << "Accum Threshold Errors: " << sum(parameterthresholds) << endl;

        ModelCovar = JointAnalyser.GetModelCovariance(DataCovar);
        double currerr = 0;
        for (size_t i = 0; i < nlayers; ++i)
          modelerrors(i) = 2.3025 * MTModel.GetResistivities().at(i) * sqrt(
              ModelCovar(i, i));
        for (size_t i = nlayers; i < nmodelparms; ++i)
          {
            currerr = sqrt(ModelCovar(i, i));
            modelerrors(i) = currerr;
          }
        ++eigenvindex;
      }
  }

int main()
  {
    SeismicDataComp RecData;

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

    C1dInvGaConf Configuration;
    Configuration.GetData("1dinvga.conf");

    MTRecObjective JointObjective(MTData, RecData, Configuration.shift,
        Configuration.omega, Configuration.sigma, Configuration.cc,
        Configuration.slowness);
    double rfweight = 0;
    cout << "RF weight: ";
    cin >> rfweight;
    SetupMTFitParameters(Configuration, JointObjective.GetMTObjective());
    JointObjective.SetRecWeight(rfweight);
    JointObjective.GetRecObjective().SetPoisson(Configuration.poisson);
    JointObjective.GetRecObjective().SetTimeWindow(Configuration.starttime,
        Configuration.endtime);
    JointObjective.GetRecObjective().SetErrorLevel(Configuration.recerror);
    JointObjective.GetRecObjective().SetFitExponent(
        Configuration.recfitexponent);
    string modelnamebase;
    cout << "Model name base: ";
    cin >> modelnamebase;
    ResPkModel RecModel;
    RecModel.GetData(modelnamebase + "_rec.mod");
    C1DMTSynthData MTModel;
    MTModel.readmodel(modelnamebase + "_mt.mod");

    const size_t nmodelparms = RecModel.SVelocity.size() * 3;
    const size_t nlayers = RecModel.SVelocity.size();
    gplib::rvec JointModelVector(nmodelparms);
    for (size_t i = 0; i < nlayers; ++i)
      {
        JointModelVector(i) = log10(MTModel.GetResistivities().at(i));
        JointModelVector(i + nlayers) = RecModel.Thickness.at(i);
        JointModelVector(i + 2 * nlayers) = RecModel.SVelocity.at(i);
      }

    gplib::rvec varfactor(nmodelparms);
    double errorratio = 0.1;
    cout << "Relative Error for Model parameters; ";
    cin >> errorratio;
    double JointOldRMS = JointObjective.CalcPerformance(JointModelVector);

    cout << "RMS: " << JointOldRMS << endl;

    for (size_t i = 0; i < nmodelparms; ++i)
      {
        varfactor(i) = 1.0;
      }

    ModelAnalysis JointAnalyser(JointObjective);
    JointAnalyser.SetModel(JointModelVector);
    JointAnalyser.SetVariationFactor(varfactor);
    gplib::rvec svals = JointAnalyser.GetSingularValues();
    ofstream evals("jointevals.plot");
    for (size_t i = 0; i < svals.size(); ++i)
      evals << svals(i) << endl;
    //cout << "Relative SVD Threshold: ";
    //double svdthreshold = 0.0;
    //cin >> svdthreshold;

    gplib::rvec modelerrors(nmodelparms);
    for (size_t i = 0; i < nmodelparms; ++i)
      modelerrors(i) = 0.0;

    gplib::rvec parameterthresholds = JointModelVector;
    for (size_t i = 0; i < nlayers; ++i)
      parameterthresholds(i) = pow(10, parameterthresholds(i)) * varfactor(i);
    for (size_t i = nlayers; i < nmodelparms; ++i)
      parameterthresholds(i) *= varfactor(i);

    parameterthresholds *= errorratio;

    const size_t nrecdata = RecData.GetData().size();
    const size_t nmtdata = MTData.GetMTData().size();

    C1DMTObjective::datafuncvector_t ErrorFunctions =
        JointObjective.GetMTObjective().GetErrorFunctions();
    const size_t nerrorfuncs = ErrorFunctions.size();
    gplib::rmat DataCovar(nrecdata + nerrorfuncs * nmtdata, nrecdata
        + nerrorfuncs * nmtdata);
    fill_n(DataCovar.data().begin(), DataCovar.size1() * DataCovar.size2(), 0.0);
    for (size_t e = 0; e < nerrorfuncs; ++e)
      {
        for (size_t i = 0; i < nmtdata; ++i)
          {
            DataCovar(i + e * nmtdata, i + e * nmtdata) = gplib::pow2(
                ErrorFunctions.at(e)(&MTData.GetMTData().at(i)));
          }
      }

    for (size_t i = nmtdata * nerrorfuncs; i < nmtdata * nerrorfuncs + nrecdata; ++i)
      {
        DataCovar(i, i) = gplib::pow2(RecData.GetData().at(i - nmtdata
            * nerrorfuncs) * Configuration.recerror);
      }

    ofstream recerrors("jointerrors");
    gplib::rmat ModelCovar = JointAnalyser.GetModelCovariance(DataCovar);

    double svdthreshold;
    cout << "SVD Threshold (negative for automatic): ";
    cin >> svdthreshold;
    if (svdthreshold < 0.0)
      {
        DetermineEvalThreshold(JointAnalyser, MTModel, DataCovar, modelerrors,
            parameterthresholds);
      }
    else
      {
        JointAnalyser.SetRelativeSVDThreshold(svdthreshold);
      }

    ModelCovar = JointAnalyser.GetModelCovariance(DataCovar);
    double currerr = 0;
    for (size_t i = 0; i < nlayers; ++i)
      modelerrors(i) = 2.3025 * MTModel.GetResistivities().at(i) * sqrt(
          ModelCovar(i, i));
    for (size_t i = nlayers; i < nmodelparms; ++i)
      {
        currerr = sqrt(ModelCovar(i, i));
        modelerrors(i) = currerr;
      }

    std::vector<double> thickerr, svelerr, reserr;

    for (int i = 0; i < nlayers; ++i)
      {
        reserr.push_back(modelerrors(i));
        thickerr.push_back(modelerrors(i + nlayers));
        svelerr.push_back(modelerrors(i + nlayers * 2));
      }

    cout << "Final Accum Model Errors: " << sum(modelerrors) << endl;
    cout << "Accum Threshold Errors: " << sum(parameterthresholds) << endl;
    WriteMatrixAsNetCDF("jointcovar.nc", ModelCovar);
    WriteMatrixAsNetCDF("jointcorr.nc", JointAnalyser.GetModelCorrelation(
        DataCovar));
    WriteMatrixAsNetCDF("jointsens.nc", JointAnalyser.GetSensitivityMatrix());
    WriteMatrixAsNetCDF("jointres.nc", JointAnalyser.GetModelResolution());
    WriteMatrixAsNetCDF("jointnull.nc", JointAnalyser.GetModelNullspace());
    WriteMatrixAsNetCDF("jointsinv.nc", JointAnalyser.GetSensitivityInverse());

    RecModel.SetThickErrors(thickerr);
    RecModel.SetSVelErrors(svelerr);
    RecModel.PlotVelWithErrors("recerr.plot");
    MTModel.SetResistivityErrors(reserr);
    MTModel.SetThicknessErrors(thickerr);
    MTModel.writeplot("mterr.plot");
  }