CVariableMetric.cpp

Go to the documentation of this file.

#include "CVariableMetric.h"
#include <gsl/gsl_math.h>
#include <iostream>


using namespace std;

CVariableMetric::CVariableMetric(const unsigned int ndatapoints, const unsigned int nmodel)  : CGeneralLinInversion(ndatapoints,nmodel)
{
        gamma.resize(nparams,false);
        oldgamma.resize(nparams,false);
        oldmodel.resize(nparams,false);
        phi.resize(nparams,false);
        PartialDerivs.resize(ndata,nparams,false);
        DataCovarInv.resize(ndata,ndata);
        ModelCovarInv.resize(nparams,nparams);
        InversionResult.resize(ndata,false);
        CMu.resize(nparams,false);
        deltagamma.resize(nparams,false);
        deltamodel.resize(nparams,false);
        deltamodelF.resize(nparams,false);
        CMGT.resize(nparams,ndata);
        CDinvG.resize(ndata,nparams);
        CDMisfit.resize(ndata,true);
        b.resize(ndata,true);
        F.resize(nparams,nparams,false);
        Hinv.resize(nparams,nparams,false);
        deltaF.resize(nparams,nparams,false);
        currentstepsize = 5.;
        Wolfe2Factor = 0.1;
        FirstIteration = true;
        currentiteration = 0;
}

CVariableMetric::~CVariableMetric()
{
}

void CVariableMetric::Prepare()
{
        currentiteration = 0;
        ublas::identity_matrix<double> Iparams(nparams);
        CurrentModel = StartModel;
        //InvertMatrix(DataCovar,DataCovarInv);
        //InvertMatrix(ModelCovar,ModelCovarInv);
        DataCovarInv = DataCovar;
        ModelCovarInv = ModelCovar;
        CalcPartialDerivs();
        CalcGradient();
        F.assign(Iparams);
        axpy_prod(ModelCovar,trans(PartialDerivs),CMGT,true); // CMGT = Cm * G^T
        cout << "ModelCovar: " << ModelCovar << endl;
        cout << "PartialDerivs: " << PartialDerivs << endl;
        axpy_prod(DataCovarInv,CurrentMisfit,CDMisfit,true);  //CDMisfit = Cd * misfit
        axpy_prod(CMGT,CDMisfit,gamma,true); //gamma = CMGT * CDMisfit
        cout << "Gamma: " << gamma << endl;
        //gamma = Gradient;
        //cout << "Gradient: " << Gradient << endl;
        axpy_prod(F,gamma,phi); //phi = F * gamma
        phi *= -1.;
        cout << "F: " << F << endl;
}

void CVariableMetric::DoIteration()
{
        double NewRMS = 1e90;
        const unsigned int maxsearch = 20;
        unsigned int searchiterations = 0;
        double Fdenominator;
        double usedstepsize = currentstepsize;
        double mudenom,muenum, mu;
        double newgradsearchprodukt;
        
        oldgamma = gamma;
        oldmodel = CurrentModel;
        
        axpy_prod(PartialDerivs,phi,b,true); // b = part * phi
        muenum = -inner_prod(prod(gamma,ModelCovarInv),phi);
        mudenom = inner_prod(prod(phi,ModelCovarInv),phi) + inner_prod(prod(b,DataCovarInv),b);
        //mu = muenum/mudenom;
        mu = 1.;
        const double gradsearchprodukt =  inner_prod(gamma,phi);
        const double searchangle = acos(gradsearchprodukt/(norm_2(gamma)*norm_2(phi))) *180./M_PI;
        const double Wolfe2 = gradsearchprodukt * Wolfe2Factor;
        cout << endl << "Iteration : " << currentiteration << endl;
        cout << "CMGT: " << CMGT << endl;
        cout << "CurrentMisfit: " << CurrentMisfit<< endl;
        cout << "CDMisfit: " << CDMisfit << endl;
        cout << "Gamma: " << gamma << endl;
        cout << "Phi: " << phi << endl;
        cout << "Angle: " << searchangle <<  endl;
        cout << "CurrentModel: " << CurrentModel << endl;
        do
        {
                cout << "Orig Stepsize: " << mu << endl; 
                mu = LineSearch->SearchStepsize(mu,CurrentRMS,phi,CurrentModel,InputData,gamma);
                cout << "New Stepsize: " << mu << endl;
                CurrentModel = oldmodel + mu * phi;
                CalcPartialDerivs();
                CalcGradient();
                axpy_prod(ModelCovar,trans(PartialDerivs),CMGT,true);
                axpy_prod(DataCovarInv,CurrentMisfit,CDMisfit,true);
                axpy_prod(CMGT,CDMisfit,gamma,true);
                cout << "Gamma: " << gamma << endl;
                //gamma = Gradient;
                //cout << "Gradient: " << Gradient << endl;
                newgradsearchprodukt = inner_prod(gamma,phi);
                
                cout << "NewGradSearchProdukt: " << newgradsearchprodukt << " ";
                cout << "Wolfe2: " << Wolfe2 << endl;
                cout << "CurrentModel: " << CurrentModel << endl;
                cout << "CurrentMisfit: " << CurrentMisfit<< endl;      
                mu *= 10.;
        }
        while (newgradsearchprodukt *1.001 < Wolfe2);
        
        cout << "New Mu: " << mu << endl; 
        
        deltamodel = CurrentModel - oldmodel;
        deltagamma = gamma - oldgamma;
        gplib::rvec v = prod(F,deltagamma);
        gplib::rvec u = deltamodel - v;
        //Update of Hessian using the BFGS formula after Nocedal 1992 p. 17
        double yTs = inner_prod(deltagamma,deltamodel);
        cout << "yTs: " << yTs << endl;
        if (yTs > 0) // we can update the Hessian using the BFGS formula
        {
                cout << "F before: " << F << endl;
                cout << "DeltaModel: " << deltamodel << endl;
                axpy_prod(F,deltamodel,deltamodelF,true);
                cout << "DeltaModelF: " << deltamodelF << endl;
                Fdenominator = inner_prod(deltamodelF,deltamodel);
                deltaF = -outer_prod(deltamodelF,deltamodelF)/Fdenominator;
                cout << "Delta F: " << deltaF << endl;
                //axpy_prod(deltaF,ModelCovarInv,F,false);
                cout << "F after first: " << F << endl;
                cout << "DeltaGamma: " << deltagamma << endl;
                cout << "dgammaTdgamma: " << outer_prod(deltagamma,deltagamma) << endl;
                cout << "dgamma*dgamma: " << inner_prod(deltagamma,deltamodel) << endl;
                cout << "DeltaF: " << outer_prod(deltagamma,deltagamma)/inner_prod(deltagamma,deltamodel) << endl;
                F += outer_prod(u,u)/inner_prod(u,deltagamma); 
                //F += outer_prod(deltagamma,deltagamma)/inner_prod(deltagamma,deltamodel); //needs to incorporate Covariance !!
                cout << "F after second: " << F << endl;
        }
        else // reset hessian to identity matrix
        {
                cout << "yTs: " << yTs << " Reset Hessian ! " << endl; 
                ublas::identity_matrix<double> Iparams(nparams);
                F.assign(Iparams);
        }
        //InvertMatrix(F,Hinv);
        axpy_prod(F,gamma,phi,true); //changed F for Hinv to use inverse of Hessian !!!
        phi *= -1.;
        cout << "|gamma| : " << norm_2(gamma) << endl;
        cout << "Abort: " << epsilon * (1. + CurrentRMS) << endl;
        if (norm_2(gamma) < epsilon * (1. + CurrentRMS))
                throw ConvergenceException();
                
        cout  << endl << endl;
        ++currentiteration;
        FirstIteration = false;
}