AbsVelRecObjective.cpp

Go to the documentation of this file.

#include "AbsVelRecObjective.h"
#include <gsl/gsl_math.h>

AbsVelRecObjective::~AbsVelRecObjective()
  {
  }

//constructor including absolute velocities
AbsVelRecObjective::AbsVelRecObjective(boost::shared_ptr<const SeismicDataComp> TheRecData,
    SurfaceWaveData &AbsVel, const int myshift, const double myomega,
    const double mysigma, const double myc, const double myslowness,
    const RecCalc::trfmethod method) :
  C1DRecObjective(TheRecData, myshift, myomega, mysigma, myc, myslowness,
      method), CalcAbsVel(myomega,mysigma,myc, method), MeasuredAbsVel(AbsVel)
  {
    recweight = 1;
    absvelweight = 1;
  }

AbsVelRecObjective::AbsVelRecObjective(const AbsVelRecObjective &Old) :
  C1DRecObjective(Old), CalcAbsVel(Old.CalcAbsVel),
      MeasuredAbsVel(Old.MeasuredAbsVel), SynthAbsVel(Old.SynthAbsVel)
  {
    recweight = Old.recweight;
    absvelweight = Old.absvelweight;
  }

AbsVelRecObjective& AbsVelRecObjective::operator=(const AbsVelRecObjective& source)
  {
    if (this == &source) return *this;
    C1DRecObjective::operator=(source);
    CalcAbsVel = source.CalcAbsVel;
    MeasuredAbsVel = source.MeasuredAbsVel;
    SynthAbsVel = source.SynthAbsVel;
    recweight = source.recweight;
    absvelweight = source.absvelweight;
    return *this;
  }

void AbsVelRecObjective::SafeParallel(const ttranscribed &member)
  {
    C1DRecObjective::SafeParallel(member); //call inherited version to calculate receiver function misfit
    double returnvalue = pow(GetRMS(), GetFitExponent()) * GetMisfit().size(); // calculate raw misfit from RMS
    returnvalue *= recweight; // and multiply by weight
    if (absvelweight > 0.0 && !MeasuredAbsVel.GetPhaseVelocities().empty())
      {
        SynthAbsVel.SetPeriods() = MeasuredAbsVel.GetPeriods();
        SeismicDataComp Rad(RecCalculator.GetRadComp()),
            Ver(RecCalculator.GetVerComp());
        CalcAbsVel.CalcRFVel(slowness, Rad, Ver,
            SynthAbsVel.SetPhaseVelocities());

        assert(SynthAbsVel.GetPhaseVelocities().size() == MeasuredAbsVel.GetPhaseVelocities().size());
        const unsigned int nvel = SynthAbsVel.GetPhaseVelocities().size();
        const unsigned int nrec =GetMisfit().size(); // the number of previous receiver function points 


        tmisfit LocalMisfit(nrec+nvel); // we need vectors of the right size for misfit, this is for levmarq
        tdata LocalData(nrec+nvel); // and  data, levmarq needs to know the size before, so we do not consider weights for the size
        copy(GetMisfit().begin(), GetMisfit().end(), LocalMisfit.begin());
        copy(GetSynthData().begin(), GetSynthData().end(), LocalData.begin());

        double currvalue = 0.0;
        for (unsigned int i = 0; i < nvel; ++i)
          {
            currvalue = SynthAbsVel.GetPhaseVelocities().at(i) - MeasuredAbsVel.GetPhaseVelocities().at(i);
            currvalue /= MeasuredAbsVel.GetPhaseVelocities().at(i) * errorlevel;
            currvalue = gsl_pow_int(currvalue, GetFitExponent());
            LocalMisfit(i+nrec)= fabs(currvalue);
            LocalData(i+nrec)= SynthAbsVel.GetPhaseVelocities().at(i); //copy data (for levmar algorithm)
            assert(gsl_finite(currvalue));
            returnvalue += absvelweight * fabs(currvalue);
          }
        SetMisfit(LocalMisfit); // set appropriate members of baseclass
        SetSynthData(LocalData);
        unsigned int ndata;
        double divisor = 1.0;
        if (recweight > 0) // we have to divide by the right amount of datapoints
          {
            ndata = nrec+nvel;
            divisor = nrec * recweight + nvel*absvelweight;
          }
        else
          {
            ndata = nvel;
            divisor = ndata * absvelweight;
          }
        SetRMS(pow(returnvalue/divisor, 1.0/GetFitExponent())); //store new misfit in local class
      }

  }

---