recanalys.cpp

Go to the documentation of this file.

#include "C1DRecObjective.h"
#include "GradFunc.h"
#include "convert.h"
#include "SeismicDataComp.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 "NumUtil.h"

using namespace std;

int main()
  {
    SeismicDataComp RecData;

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

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

    C1DRecObjective RecObjective(RecData, Configuration.shift,
        Configuration.omega, Configuration.sigma, Configuration.cc,
        Configuration.slowness);

    string modelnamebase;
    cout << "Model name base: ";
    cin >> modelnamebase;
    ResPkModel RecModel;
    RecModel.GetData(modelnamebase + "_rec.mod");

    const size_t nmodelparms = RecModel.SVelocity.size() * 2;
    gplib::rvec RecModelVector(nmodelparms);
    for (int i = 0; i < nmodelparms / 2; ++i)
      {
        RecModelVector(i) = RecModel.Thickness.at(i);
        RecModelVector(i + nmodelparms / 2) = RecModel.SVelocity.at(i);
      }

    gplib::rvec varfactor(nmodelparms);
    RecObjective.SetPoisson(Configuration.poisson);
    RecObjective.SetTimeWindow(Configuration.starttime, Configuration.endtime);
    RecObjective.SetErrorLevel(Configuration.recerror);
    RecObjective.SetFitExponent(Configuration.recfitexponent);
    double RecOldRMS = RecObjective.CalcPerformance(RecModelVector);
    cout << "RMS: " << RecOldRMS << endl;
    gplib::rvec OldData(RecObjective.GetSynthData());
    for (size_t i = 0; i < nmodelparms; ++i)
      {
        varfactor(i) = 1.0;
      }

    ModelAnalysis RecAnalyser(RecObjective);
    RecAnalyser.SetModel(RecModelVector);
    RecAnalyser.SetVariationFactor(varfactor);
    ofstream evals("mtevals.plot");
    gplib::rvec svals = RecAnalyser.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);
    fill_n(modelerrors.begin(), nmodelparms, 0.0);
    const double errorratio = 0.1;
    gplib::rvec parameterthresholds = RecModelVector;
    for (size_t i = 0; i < nmodelparms; ++i)
      parameterthresholds(i) *= varfactor(i);

    parameterthresholds *= errorratio;
    size_t eigenvindex = 0;

    const size_t ndata = RecData.GetData().size();
    gplib::rmat DataCovar(ndata, ndata);
    fill_n(DataCovar.data().begin(), DataCovar.size1() * DataCovar.size2(), 0.0);

    for (int i = 0; i < ndata; ++i)
      {
        DataCovar(i, i) = gplib::pow2(RecData.GetData().at(i)
            * Configuration.recerror);
      }

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

    while (sum(modelerrors) < sum(parameterthresholds))
      {
        cout << "Ratio: " << svalratios(eigenvindex) << endl;
        RecAnalyser.SetRelativeSVDThreshold(svalratios(eigenvindex));
        cout << "Accum Model Errors: " << sum(modelerrors) << endl;
        cout << "Accum Threshold Errors: " << sum(parameterthresholds) << endl;

        ModelCovar = RecAnalyser.GetModelCovariance(DataCovar);
        double currerr = 0;
        for (size_t i = 0; i < nmodelparms; ++i)
          {
            currerr = sqrt(ModelCovar(i, i));
            recerrors << i << " : " << currerr << endl;
            modelerrors(i) = currerr;
          }
        ++eigenvindex;
      }
    std::vector<double> thickerr, svelerr;
    for (int i = 0; i < nmodelparms / 2; ++i)
      {
        thickerr.push_back(modelerrors(i));
        svelerr.push_back(modelerrors(i + nmodelparms / 2));
      }
    RecModel.SetThickErrors(thickerr);
    RecModel.SetSVelErrors(svelerr);
    cout << "Final Accum Model Errors: " << sum(modelerrors) << endl;
    cout << "Accum Threshold Errors: " << sum(parameterthresholds) << endl;
    WriteMatrixAsNetCDF("reccovar.nc", ModelCovar);
    WriteMatrixAsNetCDF("reccorr.nc",
        RecAnalyser.GetModelCorrelation(DataCovar));
    WriteMatrixAsNetCDF("recsens.nc", RecAnalyser.GetSensitivityMatrix());
    WriteMatrixAsNetCDF("recres.nc", RecAnalyser.GetModelResolution());
    WriteMatrixAsNetCDF("recnull.nc", RecAnalyser.GetModelNullspace());
    WriteMatrixAsNetCDF("recsinv.nc", RecAnalyser.GetSensitivityInverse());
    RecModel.PlotVelWithErrors("recerr.plot");
  }