C1DRecObjective.cpp

Go to the documentation of this file.

#include <iostream>
#include <numeric>
#include <algorithm>
#include <boost/bind.hpp>
#include <boost/cast.hpp>
#include "NumUtil.h"
#include "CollapseModel.h"
#include "SeisTools.h"
#include "convert.h"
#include "C1DRecObjective.h"
#include "miscfunc.h"
#include "FatalException.h"

using namespace std;

namespace gplib
  {
    C1DRecObjective::C1DRecObjective(boost::shared_ptr<const SeismicDataComp> TheRecData,
        int myshift, double mysigma, double myc, double myslowness,
        RecCalc::trfmethod method, bool normalized, ResPkModel::WaveType InWave) :
        RecCalculator(myshift, mysigma, myc, true, method), ObservedData(TheRecData)
      {
        Model.SetInputWave(InWave);
        slowness = myslowness; //copy parameter value
        poisson = sqrt(3.); //set some default values
        //we specify a default errorlevel of 1%
        SetErrorLevel(0.01);
        //by default we fit the complete receiver function
        startpoint = 0;
        endpoint = ObservedData->GetData().size();
        // we determine the maximum amplitude of the receiver function

        const double recmaxamp = *max_element(ObservedData->GetData().begin(),
            ObservedData->GetData().end());
        const bool isone = fcmp(recmaxamp, 1.0, 1e-4) == 0;
        if (isone && !normalized)
          throw FatalException(
              "Input data has maximum amplitude of 1, but no normalization used");
        if (!isone && normalized)
          throw FatalException(
              "Input data does not have an amplitude of 1, but normalization is used");
        RecCalculator.SetNormalize(normalized); // we can choose we want to model normalized data
      }

    C1DRecObjective::~C1DRecObjective()
      {
      }

    C1DRecObjective::C1DRecObjective(const C1DRecObjective &Old) :
        PlottableObjective(Old), RecSynthData(Old.RecSynthData), slowness(Old.slowness), Model(
            Old.Model), RecCalculator(Old.RecCalculator), errorlevel(Old.errorlevel), errorvalue(
            Old.errorvalue), poisson(Old.poisson), startpoint(Old.startpoint), endpoint(
            Old.endpoint), ObservedData(Old.ObservedData)
      {
      }

    void C1DRecObjective::SetTimeWindow(const double start, const double end)
      {
        //calculate the startpoint corresponding to the starttime
        //if the values do not make sense, we use the the minimum
        // and maximum possible indices respectively
        startpoint = max(boost::numeric_cast<size_t>(0),
            boost::numeric_cast<size_t>(
                (start - ObservedData->GetB()) / ObservedData->GetDt()));
        // calculate datapoint corresponding to the endtime
        endpoint = min(
            boost::numeric_cast<size_t>(
                (end - ObservedData->GetB()) / ObservedData->GetDt()),
            ObservedData->GetData().size());
        //the only thing that remains to be checked is whether the start is before the end
        if (endpoint < startpoint)
          throw FatalException(
              "Start index " + stringify(startpoint) + " and end index "
                  + stringify(endpoint) + " do not make sense !");
      }

    C1DRecObjective& C1DRecObjective::operator=(const C1DRecObjective& source)
      {
        if (this == &source)
          return *this;
        PlottableObjective::operator=(source);
        RecSynthData = source.RecSynthData;
        slowness = source.slowness;
        Model = source.Model;
        RecCalculator = source.RecCalculator;
        poisson = source.poisson;
        errorlevel = source.errorlevel;
        errorvalue = source.errorvalue;
        startpoint = source.startpoint;
        endpoint = source.endpoint;
        ObservedData = source.ObservedData;
        return *this;
      }

    void C1DRecObjective::PreParallel(const ttranscribed &member)
      {
        const unsigned int nfulllayers = member.size() / 2; //The model vector contains thickness and SVel, so we have size/2 layers

        Model.SetSlowness(slowness); //copy slowness
        ttranscribed thickness(nfulllayers), velocity(nfulllayers);
        copy(member.begin(), member.begin() + nfulllayers, thickness.begin()); //copy values
        copy(member.begin() + nfulllayers, member.begin() + 2 * nfulllayers,
            velocity.begin());
        CollapseModel(thickness, velocity); //reduce to minimum number of different layers
        const unsigned int ncolllayers = velocity.size();
        Model.Init(ncolllayers); //initialize model with number of reduced layers
        copy(thickness.begin(), thickness.end(), Model.SetThickness().begin());
        copy(velocity.begin(), velocity.end(), Model.SetSVelocity().begin());
        transform(velocity.begin(), velocity.end(), Model.SetDensity().begin(),
            CalcDensity());
        transform(velocity.begin(), velocity.end(), //PVelocity is calculated by multiplication of SVel by poisson's ratio
            Model.SetPVelocity().begin(), boost::bind(multiplies<double>(), _1, poisson));
        Model.SetDt(ObservedData->GetDt()); // set dt
        Model.SetNpts(ObservedData->GetData().size()); //set number of points
        RecCalculator.SynthPreParallel(GetParallelID(), Model, RecSynthData, true);
      }

    double C1DRecObjective::PostParallel(const ttranscribed &member)
      {
        RecCalculator.SynthPostParallel(GetParallelID(), Model, RecSynthData, true);
        SetMisfit().resize(endpoint - startpoint); //we need vectors of the right size for misfit
        SetSynthData().resize(endpoint - startpoint); //and  data
        double returnvalue = 0.0; //init returnvalue
        //go through the receiver function from the starttime to the endtime point by point
        for (int i = startpoint; i < endpoint; ++i)

          {
            //calculate the misfit for each point and add up
            //the relative errorlevel is 0. as it does not make sense for receiver functions
            //therefore always errorvalue will be used as an absolue error
            returnvalue += CalcMisfit(ObservedData->GetData().at(i),
                RecSynthData.GetData().at(i), errorvalue, 0.0, i - startpoint);
          }
        //copy the head information from the observed receiver function
        RecSynthData.CopyHeader(*ObservedData.get());
        SetRMS(std::pow(returnvalue / (endpoint - startpoint), 1.0 / GetFitExponent())); //store misfit in local class
        return GetRMS(); //return misfit calculated by SafeParallel
      }

    void C1DRecObjective::SafeParallel(const ttranscribed &member)
      {
        RecCalculator.SynthSafeParallel(GetParallelID(), Model, RecSynthData, true); // Calculate forward model
      }
  }