C1DRecObjective.h

Go to the documentation of this file.

#ifndef C1DRECOBJECTIVE_H
#define C1DRECOBJECTIVE_H
#include "PlottableObjective.h"
#include "RecCalc.h"
#include "SeismicDataComp.h"
#include <boost/shared_ptr.hpp>
#include <numeric>
#include <algorithm>
namespace gplib
  {
    //! Calculate the misfit between observed receiver function for a given 1D model by calculating a synthetic receiver function from that model
    /*! The constructor takes a few essential parameters that are expected not to change during the
     * program (in most cases inversion), e.g. the data to fit. See also #C1DRecObjective for the necessary
     * constructor parameters.
     *
     * There are two ways to calculate the Misfit for a given model, either calling GeneralObjective#CalcPerformance,
     * or, for parallel programs, calling the three function PreParallel, SafeParallel and PostParallel in that order.
     * In both cases you have to give the parameter #member, which is a ::ttranscribed of size 2n, where
     * n is the number of layers of the model. The first n entries are the layer thicknesses in km and the last n entries
     * are the S-wave velocities in km/s. Densities and P-Wave velocities are calculated from S-Velocities by hard-coded relationships.
     */
    class C1DRecObjective: public PlottableObjective
      {
    private:
      //! The object holding the calculated synthetic data
      SeismicDataComp RecSynthData;
      //! The slowness used in the calculation
      double slowness;
      //! The object holding the 1D model for forward calculation
      ResPkModel Model;
      //! The object used to calculate the synthetics
      RecCalc RecCalculator;
      //! The relative errorfloor
      double errorlevel;
      //! The absolute errorfloor is calculated from the relative errorfloor
      double errorvalue;
      //! Poisson's ratio vp/vs, this might become part of the model
      double poisson;
      //! We might not want to use the whole receiver function for misfit calculation, the startpoint is calculated from the starttime in SetTimeWindow
      int startpoint;
      //! We might not want to use the whole receiver function for misfit calculation, the endpoint is calculated from the endtime in SetTimeWindow
      int endpoint;
      //! The measured data is stored in a shared pointer.
      /*! For parallel implementations we often work with several instances of this
       * class that all work with the same ObservedData. To save memory we use a shared_pointer
       * and to avoid collisions SeismicDataComp is const */
      boost::shared_ptr<const SeismicDataComp> ObservedData;
    protected:
      const SeismicDataComp &GetObservedData()
        {
          return *ObservedData;
        }
    public:
      //! return a pointer to a copy of the current object
      virtual C1DRecObjective *clone() const
        {
          return new C1DRecObjective(*this);
        }
      //! Set the time window used for misfit calculations, start and end are in seconds
      void SetTimeWindow(const double start, const double end);
      //! Set the errorlevel for fit, this is relative to the maximum amplitude, not for each individual data point
      void SetErrorLevel(const double level)
        {
          errorlevel = level;
          const double recmaxamp = *std::max_element(ObservedData->GetData().begin(),
              ObservedData->GetData().end());
          // we specify the error relative to the maximum amplitude
          //this is usually the initial correlation peak
          // here we calculate its value
          errorvalue = errorlevel * recmaxamp;
        }
      //! Set poisson's ratio, at the moment the same for all layers, used for calculating P-velocity
      void SetPoisson(const double ratio)
        {
          poisson = ratio;
        }
      //! Write the synthetic data to a sac file with name filename, makes only sense after calculating the misfit
      virtual void WriteData(const std::string &filename)
        {
          RecSynthData.WriteAsSac(filename);
        }
      //! Write the current model to ascii file for calculations
      void WriteModel(const std::string &filename)
        {
          Model.WriteModel(filename);
        }
      //! Write the current model to ascii file for plotting
      void WritePlot(const std::string &filename)
        {
          Model.WritePlot(filename);
        }
      //! We have to write runfiles before parallel execution
      virtual void PreParallel(const ttranscribed &member);
      //! We also clean up files serially
      virtual double PostParallel(const ttranscribed &member);
      //! Calculate the misfit between the data calculated from model vector member and measured data given in the constructor.
      virtual void SafeParallel(const ttranscribed &member);
      C1DRecObjective(const C1DRecObjective &Old);
      C1DRecObjective& operator=(const C1DRecObjective& source);
      //! The constructor needs a few essential parameters
      /*! @param TheRecData   Shared pointer to the measured receiver function
       *  @param myshift    the time shift used for calculating the measured receiver function
       *  @param mysigma sigma used for calculating the measured receiver function
       *  @param myc water level used for calculating the measured receiver function
       *  @param myslowness slowness used for calculating the measured receiver function
       * @param method The method used to calculate the observed data. Can be specdiv or iterdecon
       * @param normalized Is the measured data normalized to an initial correlation peak of 1
       * @param InWave Is the incoming plane wave a p-wave or a s-wave
       */
      C1DRecObjective(boost::shared_ptr<const SeismicDataComp> TheRecData, int myshift,
          double mysigma, double myc, double myslowness, RecCalc::trfmethod method =
              RecCalc::specdiv, bool normalized = true, ResPkModel::WaveType InWave =
              ResPkModel::PWave);

      virtual ~C1DRecObjective();
      friend class AbsVelRecObjective;
      };
  }
#endif // C1DRECOBJECTIVE_H