C1DRecObjective.cpp

Go to the documentation of this file.

#include <iostream>
#include <numeric>
#include <algorithm>
#include <boost/bind.hpp>
#include <boost/cast.hpp>
#include <gsl/gsl_math.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, const int myshift,
        const double mysigma, const double myc, const double myslowness,
        const RecCalc::trfmethod method, const bool normalized) :
      RecCalculator(myshift, mysigma, myc, true, method), ObservedData(
          TheRecData)
      {
        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 = gsl_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(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
      }
  }