CGeneralLinInversion.cpp

Go to the documentation of this file.

#include "CGeneralLinInversion.h"
#include <iostream>
#include <limits>
using namespace std;

CGeneralLinInversion::CGeneralLinInversion(const unsigned int ndatapoints, const unsigned int nmodel)
{
        Iterationnumber = 0;
        CurrentRMS = 1e90;
        nparams = nmodel;
        ndata = ndatapoints;
        InversionResult.resize(ndata,false); //The vector containing the best fitting synthetic data after the current iteration
        CurrentData.resize(ndata,false);
        CurrentMisfit.resize(ndata,false);
        Gradient.resize(nparams,false);
        CurrentModel.resize(nparams,false);
        DataCovar.resize(ndata,ndata,false);
        ModelCovar.resize(nparams,nparams,false);
        PartialDerivs.resize(ndata,nparams,false);
        epsilon = sqrt(numeric_limits<double>::epsilon());
}

CGeneralLinInversion::~CGeneralLinInversion()
{
}

void CGeneralLinInversion::CalcGradient()
{
        
        const double delta = 0.0001;
        const double poschange = 1. + delta;
        const double negchange = 1. - delta;
        double denominator;
        CurrentRMS = 0;
        double posrms, negrms;
        
    CurrentRMS = MisfitCalculator->CalcPerformance(CurrentModel);
    CurrentMisfit = MisfitCalculator->GetMisfit();
    CurrentData = MisfitCalculator->GetSynthData();
    //cout << "Model: " << CurrentModel << " Misfit:" << CurrentMisfit << endl;
    posmodel = CurrentModel;
    negmodel = CurrentModel;
        
    for (unsigned int i = 0; i < nparams; ++i)
    {
        if (i > 0) //revert the change of the last iteration
        {
                posmodel(i-1) = CurrentModel(i-1);
                negmodel(i-1) = CurrentModel(i-1);
        }
        posmodel(i) = CurrentModel(i) * poschange;
        negmodel(i) = CurrentModel(i) * negchange;
        posrms = MisfitCalculator->CalcPerformance(posmodel);
        negrms = MisfitCalculator->CalcPerformance(negmodel);
        
        //cout << "Pos Model: " << posmodel << " Data:" << posdata << endl;
        //cout << "Neg Model: " << negmodel << " Data:" << negdata << endl;
        denominator = 1./(2.* delta * CurrentModel(i));
        Gradient(i) = (posrms - negrms)*denominator;
    }
}

void CGeneralLinInversion::CalcPartialDerivs()
{
        unsigned int i,j;
        
        const double delta = 0.0001;
        const double poschange = 1. + delta;
        const double negchange = 1. - delta;
        double denominator;
        CurrentRMS = 0;
        
        
    CurrentRMS = MisfitCalculator->CalcPerformance(CurrentModel);
    CurrentMisfit = MisfitCalculator->GetMisfit();
    CurrentData = MisfitCalculator->GetSynthData();
    //cout << "Model: " << CurrentModel << " Misfit:" << CurrentMisfit << endl;
    posmodel = CurrentModel;
    negmodel = CurrentModel;
        
    for (i = 0; i < nparams; ++i)
    {
        if (i > 0)
        {
                posmodel(i-1) = CurrentModel(i-1);
                negmodel(i-1) = CurrentModel(i-1);
        }
        posmodel(i) = CurrentModel(i) * poschange;
        negmodel(i) = CurrentModel(i) * negchange;
        MisfitCalculator->CalcPerformance(posmodel);
        posdata = MisfitCalculator->GetSynthData();
        MisfitCalculator->CalcPerformance(negmodel);
        negdata = MisfitCalculator->GetSynthData();
        //cout << "Pos Model: " << posmodel << " Data:" << posdata << endl;
        //cout << "Neg Model: " << negmodel << " Data:" << negdata << endl;
        denominator = 1./(2.* delta * CurrentModel(i));
        for (j = 0; j < ndata; ++j)
        {
                                PartialDerivs(j,i) = (posdata(j) - negdata(j))*denominator;
        }
    }
    //cout << "CurrentModel: " << CurrentModel << endl;
        //std::cout << "Partial Deriv: " <<  PartialDerivs <<  std::endl;
}