GradFunc.cpp

Go to the documentation of this file.

#include "GradFunc.h"
#include "NumUtil.h"
#include <iostream>

//! Calculate the Gradient dO(m)/dm of the objective function Objective at the point given by Model, return as a real vector
gplib::rvec CalcGradient(GeneralObjective &Objective, const gplib::rvec &Model)
  {
    const double delta = 0.001;
    const double poschange = 1. + delta;
    const double negchange = 1. - delta;

    const unsigned int nparams = Model.size();
    double CurrentRMS = Objective.CalcPerformance(Model);

    gplib::rvec posmodel(Model);
    gplib::rvec negmodel(Model);
    gplib::rvec Gradient(nparams);
    for (unsigned int i = 0; i < nparams; ++i)
      {
        if (i > 0) //revert the change of the last iteration
          {
            posmodel(i - 1) = Model(i - 1);
            negmodel(i - 1) = Model(i - 1);
          }
        posmodel(i) = Model(i) * poschange;
        negmodel(i) = Model(i) * negchange;
        double posrms = Objective.CalcPerformance(posmodel);
        double negrms = Objective.CalcPerformance(negmodel);

        //cout << "Pos Model: " << posmodel << " Data:" << posdata << endl;
        //cout << "Neg Model: " << negmodel << " Data:" << negdata << endl;
        double denominator = 1. / (2. * delta * Model(i));
        if (fcmp(denominator, 0.0, 1e-10) != 0)
          Gradient(i) = (posrms - negrms) * denominator;
        else
          Gradient(i) = 0.0;
      }
    return Gradient;
  }

//! Calculate the matrix of partial derivatives (or sensitivity matrix) dD/dm, for objective function Objective at point Model
gplib::rmat CalcPartialDerivs(GeneralObjective &Objective,
    const gplib::rvec &Model, const gplib::rvec varfactor)
  {
    unsigned int i, j;

    const double delta = 0.001;
    const double poschange = 1. + delta;
    const double negchange = 1. - delta;
    double denominator;
    const int nparams = Model.size();
    Objective.CalcPerformance(Model);
    const int ndata = Objective.GetSynthData().size();
    gplib::rvec posdata(Objective.GetSynthData());
    gplib::rvec negdata(Objective.GetSynthData());
    gplib::rvec posmodel(Model);
    gplib::rvec negmodel(Model);
    gplib::rmat PartialDerivs(ndata, nparams);

    for (i = 0; i < nparams; ++i)
      {
        if (i > 0)
          {
            posmodel(i - 1) = Model(i - 1);
            negmodel(i - 1) = Model(i - 1);
          }
        posmodel(i) = Model(i) * poschange;
        negmodel(i) = Model(i) * negchange;
        Objective.CalcPerformance(posmodel);
        posdata = Objective.GetSynthData();
        Objective.CalcPerformance(negmodel);
        negdata = Objective.GetSynthData();
        if (fcmp(Model(i), 0.0, 1e-8) != 0)
          denominator = 1. / (2. * delta * Model(i));
        else
          denominator = 0.0;
        for (j = 0; j < ndata; ++j)
          {
            PartialDerivs(j, i) = (posdata(j) - negdata(j)) * denominator
                * varfactor(i);
          }

      }
    return PartialDerivs;
  }

void ModelAnalysis::SetModel(const gplib::rvec &TheModel)
  {
    Model = TheModel;
    haveSVD = false;
    s.resize(Model.size());
    sinv.resize(Model.size(), Model.size());
    vt.resize(Model.size(), Model.size());
    vnull.resize(Model.size(), Model.size());
    varfactor.resize(Model.size());
    for (size_t i = 0; i < varfactor.size(); ++i)
      varfactor(i) = 1.0;
  }

void ModelAnalysis::CalcSensitivityMatrix()
  {
    //SensitivityMatrix.resize(Objective.GetSynthData().size(),Model.size(),false);
    SensitivityMatrix = CalcPartialDerivs(Objective, Model, varfactor);
  }

void ModelAnalysis::CalcSVD()
  {
    if (!haveSVD)
      {
        CalcSensitivityMatrix();
        //std::cout << "ModelAnalyis: " << SensitivityMatrix << std::endl;
        OldSensitivityMatrix = SensitivityMatrix;
      }
    else
      {
        SensitivityMatrix = OldSensitivityMatrix;
      }
    if (changedThreshold)
      {
        u.resize(SensitivityMatrix.size1(), Model.size());
        boost::numeric::bindings::lapack::gesvd(SensitivityMatrix, s, u, vt);
        ApplyThreshold();

      }
    haveSVD = true;
    changedThreshold = false;
  }

void ModelAnalysis::ApplyThreshold()
  {
    const double maxeval = *max_element(s.begin(), s.end());
    std::cout << "Maximum Eval: " << maxeval << std::endl;
    const double AbsoluteThreshold = maxeval * SVDThreshold;
    std::cout << "Eval Threshold: " << AbsoluteThreshold << std::endl;
    for (int i = 0; i < s.size(); ++i)
      {
        column(sinv, i) *= 0.0;
        if (fcmp(s(i), AbsoluteThreshold, std::numeric_limits<double>::epsilon()) == -1
            || (fcmp(s(i), 0.0, std::numeric_limits<double>::epsilon()) == 0))
          {
            //std::cout << "Eval: " << s(i) << " set to zero !" << std::endl;
            column(vnull, i) = row(vt, i);
            column(u, i) *= 0.0;
            row(vt, i) *= 0.0;
            sinv(i, i) = 0.0;
          }
        else
          {
            column(vnull, i) *= 0.0;
            sinv(i, i) = 1. / s(i);
          }
      }
  }

gplib::rmat ModelAnalysis::GetModelResolution()
  {
    CalcSVD();
    return prec_prod(trans(vt), vt);
  }
gplib::rmat ModelAnalysis::GetDataResolution()
  {
    CalcSVD();
    return prec_prod(u, trans(u));
  }
gplib::rmat ModelAnalysis::GetModelNullspace()
  {
    CalcSVD();
    return prec_prod(vnull, trans(vnull));
  }

//! The calculation of the Model Covariance is based on Menke, p. 122 formula 7.37
gplib::rmat ModelAnalysis::GetModelCovariance(const gplib::rmat &DataCovar)
  {
    CalcSVD();
    gplib::rmat temp1 = prec_prod(sinv, trans(u));
    SensitivityInv = prec_prod(trans(vt), temp1);
    //gplib::rmat temp2 = prec_prod(DataCovar,trans(SensitivityInv));
    //return prec_prod(SensitivityInv,temp2);
    return prec_prod(SensitivityInv, trans(SensitivityInv));
  }
gplib::rmat ModelAnalysis::GetModelCorrelation(const gplib::rmat &DataCovar)
  {
    gplib::rmat temp1 = GetModelCovariance(DataCovar);
    gplib::rmat temp2 = temp1;
    for (size_t i = 0; i < temp1.size1(); ++i)
      for (size_t j = 0; j < temp1.size2(); ++j)
        {
          if (std::abs(temp1(i, i)) < 1e-38 || std::abs(temp1(j, j))
              < 1e-38)
            temp2(i, j) = 0.0;
          else
            {
              temp2(i, j) = temp1(i, j) / (sqrt(temp1(i, i))
                  * sqrt(temp1(j, j)));
              //std::cout << i << " " << j << " " << temp2(i,j) << " " << temp1(i,j) << " " << temp1(i,i) << " " << temp1(j,j) << std::endl;
            }
        }
    return temp2;
  }

gplib::rvec ModelAnalysis::GetSingularValues()
  {
    CalcSVD();
    return s;
  }