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 <gsl/gsl_math.h>
#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"
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);
        DataCovar *= 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) = gsl_pow_2(ErrorFunctions.at(e)(&MTData.GetMTData().at(i))); 
                        }
        }
        
        for (size_t i = nmtdata*nerrorfuncs; i < nmtdata*nerrorfuncs+nrecdata; ++i)
        {
                        DataCovar(i,i) = gsl_pow_2(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"); 
}

---