Go to the documentation of this file.
1 #include <iostream>
2 #include <numeric>
3 #include <algorithm>
4 #include <boost/bind.hpp>
5 #include <boost/cast.hpp>
6 #include "NumUtil.h"
7 #include "CollapseModel.h"
8 #include "SeisTools.h"
9 #include "convert.h"
10 #include "C1DRecObjective.h"
11 #include "miscfunc.h"
12 #include "FatalException.h"
14 using namespace std;
16 namespace gplib
17  {
18  C1DRecObjective::C1DRecObjective(boost::shared_ptr<const SeismicDataComp> TheRecData,
19  int myshift, double mysigma, double myc, double myslowness,
20  RecCalc::trfmethod method, bool normalized, ResPkModel::WaveType InWave) :
21  RecCalculator(myshift, mysigma, myc, true, method), ObservedData(TheRecData)
22  {
23  Model.SetInputWave(InWave);
24  slowness = myslowness; //copy parameter value
25  poisson = sqrt(3.); //set some default values
26  //we specify a default errorlevel of 1%
27  SetErrorLevel(0.01);
28  //by default we fit the complete receiver function
29  startpoint = 0;
30  endpoint = ObservedData->GetData().size();
31  // we determine the maximum amplitude of the receiver function
33  const double recmaxamp = *max_element(ObservedData->GetData().begin(),
34  ObservedData->GetData().end());
35  const bool isone = fcmp(recmaxamp, 1.0, 1e-4) == 0;
36  if (isone && !normalized)
37  throw FatalException(
38  "Input data has maximum amplitude of 1, but no normalization used");
39  if (!isone && normalized)
40  throw FatalException(
41  "Input data does not have an amplitude of 1, but normalization is used");
42  RecCalculator.SetNormalize(normalized); // we can choose we want to model normalized data
43  }
46  {
47  }
50  PlottableObjective(Old), RecSynthData(Old.RecSynthData), slowness(Old.slowness), Model(
51  Old.Model), RecCalculator(Old.RecCalculator), errorlevel(Old.errorlevel), errorvalue(
52  Old.errorvalue), poisson(Old.poisson), startpoint(Old.startpoint), endpoint(
53  Old.endpoint), ObservedData(Old.ObservedData)
54  {
55  }
57  void C1DRecObjective::SetTimeWindow(const double start, const double end)
58  {
59  //calculate the startpoint corresponding to the starttime
60  //if the values do not make sense, we use the the minimum
61  // and maximum possible indices respectively
62  startpoint = max(boost::numeric_cast<size_t>(0),
63  boost::numeric_cast<size_t>(
64  (start - ObservedData->GetB()) / ObservedData->GetDt()));
65  // calculate datapoint corresponding to the endtime
66  endpoint = min(
67  boost::numeric_cast<size_t>(
68  (end - ObservedData->GetB()) / ObservedData->GetDt()),
69  ObservedData->GetData().size());
70  //the only thing that remains to be checked is whether the start is before the end
71  if (endpoint < startpoint)
72  throw FatalException(
73  "Start index " + stringify(startpoint) + " and end index "
74  + stringify(endpoint) + " do not make sense !");
75  }
78  {
79  if (this == &source)
80  return *this;
82  RecSynthData = source.RecSynthData;
83  slowness = source.slowness;
84  Model = source.Model;
85  RecCalculator = source.RecCalculator;
86  poisson = source.poisson;
87  errorlevel = source.errorlevel;
88  errorvalue = source.errorvalue;
89  startpoint = source.startpoint;
90  endpoint = source.endpoint;
91  ObservedData = source.ObservedData;
92  return *this;
93  }
96  {
97  const unsigned int nfulllayers = member.size() / 2; //The model vector contains thickness and SVel, so we have size/2 layers
99  Model.SetSlowness(slowness); //copy slowness
100  ttranscribed thickness(nfulllayers), velocity(nfulllayers);
101  copy(member.begin(), member.begin() + nfulllayers, thickness.begin()); //copy values
102  copy(member.begin() + nfulllayers, member.begin() + 2 * nfulllayers,
103  velocity.begin());
104  CollapseModel(thickness, velocity); //reduce to minimum number of different layers
105  const unsigned int ncolllayers = velocity.size();
106  Model.Init(ncolllayers); //initialize model with number of reduced layers
107  copy(thickness.begin(), thickness.end(), Model.SetThickness().begin());
108  copy(velocity.begin(), velocity.end(), Model.SetSVelocity().begin());
109  transform(velocity.begin(), velocity.end(), Model.SetDensity().begin(),
110  CalcDensity());
111  transform(velocity.begin(), velocity.end(), //PVelocity is calculated by multiplication of SVel by poisson's ratio
112  Model.SetPVelocity().begin(), boost::bind(multiplies<double>(), _1, poisson));
113  Model.SetDt(ObservedData->GetDt()); // set dt
114  Model.SetNpts(ObservedData->GetData().size()); //set number of points
115  RecCalculator.SynthPreParallel(GetParallelID(), Model, RecSynthData, true);
116  }
119  {
120  RecCalculator.SynthPostParallel(GetParallelID(), Model, RecSynthData, true);
121  SetMisfit().resize(endpoint - startpoint); //we need vectors of the right size for misfit
122  SetSynthData().resize(endpoint - startpoint); //and data
123  double returnvalue = 0.0; //init returnvalue
124  //go through the receiver function from the starttime to the endtime point by point
125  for (int i = startpoint; i < endpoint; ++i)
127  {
128  //calculate the misfit for each point and add up
129  //the relative errorlevel is 0. as it does not make sense for receiver functions
130  //therefore always errorvalue will be used as an absolue error
131  returnvalue += CalcMisfit(ObservedData->GetData().at(i),
132  RecSynthData.GetData().at(i), errorvalue, 0.0, i - startpoint);
133  }
134  //copy the head information from the observed receiver function
135  RecSynthData.CopyHeader(*ObservedData.get());
136  SetRMS(std::pow(returnvalue / (endpoint - startpoint), 1.0 / GetFitExponent())); //store misfit in local class
137  return GetRMS(); //return misfit calculated by SafeParallel
138  }
141  {
142  RecCalculator.SynthSafeParallel(GetParallelID(), Model, RecSynthData, true); // Calculate forward model
143  }
144  }
There are several ways to calculate receiver functions.
Definition: RecCalc.h:18
ublas::vector< double > ttranscribed
Definition: gentypes.h:21
std::vector< double > & GetData()
Access for data vector, for ease of use and efficiency we return a reference.
void SynthSafeParallel(const std::string &filename, ResPkModel &Model, SeismicDataComp &Receiver, const bool cleanfiles=false)
All operations that are safe to execute in parallel.
Definition: RecCalc.cpp:210
virtual double PostParallel(const ttranscribed &member)
We also clean up files serially.
trealdata & SetDensity()
Read-write access to the vector of densities in g/cm^3 in each layer.
Definition: SeismicModel.h:96
Calculate density from a given S-velocity, the formula is taken from Owen et al. JGR 89...
Definition: SeisTools.h:15
This only adds a few plotting functions to GeneralObjective to define a common interface.
C1DRecObjective(const C1DRecObjective &Old)
void SetNormalize(const bool what)
Change whether the output receiver function is normalized to a maximum amplitude of 1...
Definition: RecCalc.h:55
void SynthPostParallel(const std::string &filename, ResPkModel &Model, SeismicDataComp &Receiver, const bool cleanfiles=false)
Operations of the synthetic receiver function calculation that are not safe in parallel and hafe to b...
Definition: RecCalc.cpp:219
void CopyHeader(const SeismicDataComp &source)
Copy the information in the header from another object.
int GetFitExponent()
Get the Fit exponent.
void SetErrorLevel(const double level)
Set the errorlevel for fit, this is relative to the maximum amplitude, not for each individual data p...
trealdata & SetThickness()
Read-write access to the vector of thicknesses in km in each layer.
Definition: SeismicModel.h:106
double CalcMisfit(const double measured, const double predicted, const double measerror, const double errorlevel, const int index)
void SetNpts(const unsigned int s)
Set the number of points for synthetic seismogram calculation.
Definition: SeismicModel.h:46
trealdata & SetPVelocity()
Read-write access to the vector of P-velocities in km/s in each layer.
Definition: SeismicModel.h:76
void SynthPreParallel(const std::string &filename, ResPkModel &Model, SeismicDataComp &Receiver, const bool cleanfiles=false)
The three Synth*Parallel methods provide alternative acces to the steps in CalcRecSynth for safe para...
Definition: RecCalc.cpp:203
const std::string & GetParallelID()
Derived classes need to read the ParallelId for their forward calculations.
double GetRMS()
Get the current RMS.
void Init(const int nlayers)
Init provides a convenient way to allocate memory in all structures for a given number of layers...
trealdata & SetSVelocity()
Read-write access to the vector of S-velocities in km/s in each layer.
Definition: SeismicModel.h:86
tdata & SetSynthData()
Only derived classes can write access the Synthetic data.
virtual void SafeParallel(const ttranscribed &member)
Calculate the misfit between the data calculated from model vector member and measured data given in ...
tmisfit & SetMisfit()
Only derived classes can write access the Misfit.
void SetInputWave(WaveType w)
Definition: ResPkModel.h:31
virtual void PreParallel(const ttranscribed &member)
We have to write runfiles before parallel execution.
PlottableObjective & operator=(const PlottableObjective &source)
Calculate the misfit between observed receiver function for a given 1D model by calculating a synthet...
C1DRecObjective & operator=(const C1DRecObjective &source)
void CollapseModel(ttranscribed &Thickness, ttranscribed &ParmValue)
Remove layers with identical parameters from the model and collapse them into a single layer each...
void SetTimeWindow(const double start, const double end)
Set the time window used for misfit calculations, start and end are in seconds.
void SetRMS(const double x)
void SetDt(const double s)
Set the time between two samples in s, this is for synthetic forward calculations.
Definition: SeismicModel.h:66
The basic exception class for all errors that arise in gplib.
void SetSlowness(const double s)
Set the slowness in s/km for the synthetic forward calculation.
Definition: ResPkModel.h:41