GPLIB++
RecCalc.cpp
Go to the documentation of this file.
1 #include "RecCalc.h"
2 #include "types.h"
3 #include "convert.h"
4 #include "miscfunc.h"
5 #include "IterDecon.h"
6 #include <iostream>
7 #include <fstream>
8 #include <algorithm>
9 #include <boost/bind.hpp>
10 #include <boost/cast.hpp>
11 #include <boost/filesystem.hpp>
12 #include <numeric>
13 #include <functional>
14 #include <cassert>
15 #include <sys/types.h>
16 #include <unistd.h>
17 #include "statutils.h"
18 #include "WFunc.h"
19 #include "CRFSpecElement.h"
20 #include "SeisTools.h"
21 #include "Adaptors.h"
22 #include "Util.h"
23 
24 using namespace std;
25 namespace gplib
26  {
27  RecCalc::RecCalc(int myshift, double mysigma, double myc, bool multicalc,
28  trfmethod themethod) :
29  method(themethod), Spectrum(multicalc), c(myc), sigma(mysigma), shift(myshift), normalize(
30  false)
31  {
32 
33  }
34 
35  RecCalc::RecCalc(const RecCalc &Old) :
36  method(Old.method), Spectrum(Old.Spectrum), c(Old.c), sigma(Old.sigma), shift(
37  Old.shift), normalize(Old.normalize), RadComp(Old.RadComp), VerComp(
38  Old.VerComp)
39  {
40 
41  }
42 
44  {
45  if (this == &source)
46  return *this;
47  method = source.method;
48  Spectrum = source.Spectrum;
49  c = source.c;
50  sigma = source.sigma;
51  shift = source.shift;
52  normalize = source.normalize;
53  RadComp = source.RadComp;
54  VerComp = source.VerComp;
55  return *this;
56  }
57 
59  {
60  }
61 
62  void RecCalc::SpectralDivision(const SeismicDataComp &RComp,
63  const SeismicDataComp &VComp, SeismicDataComp &Receiver)
64  {
65  const double omegastep = 2. * PI / (RComp.GetDt() * RComp.GetData().size());
66  // Make local copies for trimming and mean removal
67  ttsdata LocalRad(RComp.GetData());
68  ttsdata LocalVer(VComp.GetData());
69  // Make both components the same length
70  if (LocalRad.size() > LocalVer.size())
71  LocalRad.erase(LocalRad.end() - (LocalRad.size() - LocalVer.size()),
72  LocalRad.end());
73  else if (LocalVer.size() > LocalRad.size())
74  LocalVer.erase(LocalVer.end() - (LocalVer.size() - LocalRad.size()),
75  LocalVer.end());
76  //remove the mean
77  SubMean(LocalVer.begin(), LocalVer.end());
78  SubMean(LocalRad.begin(), LocalRad.end());
79  //apply a window function for spectral calculation
80  ApplyWindow(LocalVer.begin(), LocalVer.end(), LocalVer.begin(), Steep());
81  ApplyWindow(LocalRad.begin(), LocalRad.end(), LocalRad.begin(), Steep());
82  //assign space for spectra
83  tcompdata RadSpec(LocalRad.size() / 2 + 1);
84  tcompdata VerSpec(LocalVer.size() / 2 + 1);
85  //calculate spetrca
86  Spectrum.CalcSpectrum(LocalRad.begin(), LocalRad.end(), RadSpec.begin(),
87  RadSpec.end());
88  Spectrum.CalcSpectrum(LocalVer.begin(), LocalVer.end(), VerSpec.begin(),
89  VerSpec.end());
90  tcompdata RecSpec(RadSpec.size(), 0);
91 
92  Receiver.GetData().assign(LocalRad.size(), 0);
93  // determine the waterlevel threshold from the maximum vertical spectrum and the waterlevel parameter c
94  tcomp cmaxelem = *max_element(VerSpec.begin(), VerSpec.end(),
95  gplib::absLess<tcomp, tcomp>());
96  cmaxelem *= conj(cmaxelem) * c;
97 
98  // the denominator for the gaussian filter
99  const double denom = 1. / (4.0 * sigma * sigma);
100  //calculate the spectrum of the receiver function
101  transform(RadSpec.begin(), RadSpec.end(), VerSpec.begin(), RecSpec.begin(),
102  CalcSpectralElement(cmaxelem));
103  //determine the phaseshift for the desired shift in points
104  const dcomp phaseshift(0, shift);
105  // apply the gaussian filter
106  /*for (auto recit = RecSpec.begin(); recit < RecSpec.end(); ++recit)
107  {
108  double freq = omegastep * distance(RecSpec.begin(), recit);
109  *recit *= exp(-(freq * freq) * denom - phaseshift * freq);
110  }*/
111  for (size_t i = 0; i < RecSpec.size(); ++i)
112  {
113  double freq = omegastep * i;
114  RecSpec[i] *= exp(-(freq * freq) * denom - phaseshift * freq);
115  }
116  //transform receiver function spectrum back to time domain
117  Spectrum.CalcTimeSeries(RecSpec.begin(), RecSpec.end(),
118  Receiver.GetData().begin(), Receiver.GetData().end());
119  }
120 
121  void RecCalc::IterativeDeconvolution(const SeismicDataComp &RComp,
122  const SeismicDataComp &VComp, SeismicDataComp &Receiver)
123  {
124  gplib::rvec Filter(RComp.GetData().size());
125  const double dt = RComp.GetDt();
126  double normfactor = 1.; // /(sigma * sqrt(2.* acos(-1.)));
127  const int shiftpoints = boost::numeric_cast<int>(shift / dt);
128  const double variancesq = sigma * dt * sigma * dt;
129  const unsigned int maxit = 100;
130  const double minimprove = 1e-6; // we want to stop when the improvement is less than 0.1%
131  const int maxpoints = Filter.size() / 2; //we store it in a signed variable to avoid compiler warning for the loop below
132  //the loop index HAS TO BE SIGNED !
133  for (int i = 0; i <= maxpoints; ++i)
134  {
135  Filter(i) = normfactor * exp(-std::pow((i - shiftpoints), 2.) * variancesq);
136  Filter(Filter.size() - i - 1) = normfactor
137  * exp(-std::pow((-i - 1 - shiftpoints), 2.) * variancesq);
138  }
139 
140  gplib::rvec RadFilter(RComp.GetData().size()), VerFilter(VComp.GetData().size());
141  copy(RComp.GetData().begin(), RComp.GetData().end(), RadFilter.begin());
142  copy(VComp.GetData().begin(), VComp.GetData().end(), VerFilter.begin());
143  SubMean(RadFilter.begin(), RadFilter.end());
144  SubMean(VerFilter.begin(), VerFilter.end());
145  Convolve(RadFilter, Filter, RadFilter, Spectrum);
146  Convolve(VerFilter, Filter, VerFilter, Spectrum);
147  gplib::rvec Current(RadFilter), Output(RadFilter);
148 
149  IterDecon Decon(VerFilter.size(), Spectrum);
150  const double power = ublas::prec_inner_prod(RadFilter, RadFilter);
151  double error = power;
152  double lasterror = 1.;
153  double improvement = 1.;
154  unsigned int iter = 0;
155  while (improvement > minimprove && iter < maxit)
156  {
157  Decon.AdaptFilter(VerFilter, Current);
158  Decon.CalcOutput(VerFilter, Output);
159 
160  Current = RadFilter - Output;
161  error = ublas::prec_inner_prod(Current, Current) / power;
162  improvement = lasterror - error;
163  lasterror = error;
164  ++iter;
165 
166  }
167  std::cout << "Iterations: " << iter << " Error: " << error << std::endl;
168  Convolve(Decon.GetWeights(), Filter, Current, Spectrum); //we re-use current to store the filtered weights, to save another object
169  Receiver.GetData().assign(Current.size(), 0);
170  copy(Current.begin(), Current.end(), Receiver.GetData().begin());
171  }
172 
174  const SeismicDataComp &VerComp, SeismicDataComp &Receiver)
175  {
176  Receiver.CopyHeader(RadComp);
177  switch (method)
178  {
179  case iterdecon:
180  IterativeDeconvolution(RadComp, VerComp, Receiver);
181  break;
182  case specdiv:
183  SpectralDivision(RadComp, VerComp, Receiver);
184  break;
185  }
186  // set parameters for the timeseries object
187  Receiver.SetB(-shift);
188  Receiver.SetDt(RadComp.GetDt());
189  if (normalize)
190  {
191  Normalize(Receiver.GetData());
192  }
193  }
194 
195  void RecCalc::CalcRecSynth(const std::string &filename, ResPkModel &Model,
196  SeismicDataComp &Receiver, const bool cleanfiles)
197  {
198  SynthPreParallel(filename, Model, Receiver, cleanfiles);
199  SynthSafeParallel(filename, Model, Receiver, cleanfiles);
200  SynthPostParallel(filename, Model, Receiver, cleanfiles);
201  }
202 
203  void RecCalc::SynthPreParallel(const std::string &filename, ResPkModel &Model,
204  SeismicDataComp &Receiver, const bool cleanfiles)
205  {
206  Model.WriteRunFile(filename);
207  Model.WriteModel(filename + ".mod");
208  }
209 
210  void RecCalc::SynthSafeParallel(const std::string &filename, ResPkModel &Model,
211  SeismicDataComp &Receiver, const bool cleanfiles)
212  {
213  int result = std::system(("bash ./" + filename).c_str());
214  if (result != 0)
215  throw gplib::FatalException(
216  "Problem in synthetic receiver function calculation !");
217  }
218 
219  void RecCalc::SynthPostParallel(const std::string &filename, ResPkModel &Model,
220  SeismicDataComp &Receiver, const bool cleanfiles)
221  {
222  RadComp.ReadData(filename + ".mod_sp.r");
223  VerComp.ReadData(filename + ".mod_sp.z");
224  Receiver = RadComp;
225  CalcRecData(RadComp, VerComp, Receiver);
226  if (cleanfiles)
227  {
228  CleanFiles(filename);
229  }
230  }
231 
232  }
trfmethod
There are several ways to calculate receiver functions.
Definition: RecCalc.h:18
virtual ~RecCalc()
Definition: RecCalc.cpp:58
int ReadData(const std::string &filename, tseismicdataformat format=sac)
Read in data from a file, as we cannot determine the type from the ending we have to provide it...
void CalcSpectrum(_InputIterator tsbegin, _InputIterator tsend, _OutputIterator freqbegin, _OutputIterator freqend)
The member function to calculate a spectrum from a time series.
Definition: TsSpectrum.h:80
This class is used to calculate receiver functions from seismic data.
Definition: RecCalc.h:14
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 void WriteRunFile(const std::string &filename)
Write out a script that performs a forward calculation for the current model.
Definition: ResPkModel.cpp:122
void SubMean(InputIterator begin, InputIterator end, typename std::iterator_traits< InputIterator >::value_type mean)
Substract the mean from each element in the container, mean is passed as a parameter.
Definition: statutils.h:85
This functor rises steeply at the edges and then has a wide range where it is unity.
Definition: WFunc.h:51
void Normalize(std::vector< double > &Trace)
Definition: SeisTools.h:39
void CalcRecSynth(const std::string &filename, ResPkModel &Model, SeismicDataComp &Receiver, const bool cleanfiles=false)
Calculate synthetic receiver funtions from a model.
Definition: RecCalc.cpp:195
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.
void CalcRecData(const SeismicDataComp &RadComp, const SeismicDataComp &VerComp, SeismicDataComp &Receiver)
Calculate Receiver functions from two data components.
Definition: RecCalc.cpp:173
void ApplyWindow(InputIterator inbegin, InputIterator inend, OutputIterator outbegin, WindowFunctype WFunc, double relshift=0.0)
Apply one of the above window functions to a range.
Definition: WFunc.h:104
void SetB(const double theb)
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
void SetDt(const double dt)
Set delta t in s.
RecCalc & operator=(const RecCalc &source)
Definition: RecCalc.cpp:43
RecCalc(int myshift, double mysigma, double myc, bool multicalc=false, trfmethod themethod=iterdecon)
The constructor takes the essential parameters that shouldn't change during different calculations...
Definition: RecCalc.cpp:27
virtual void WriteModel(const std::string filename)
Write the model in its native format to a file.
Definition: ResPkModel.cpp:55
void CalcTimeSeries(_InputIterator freqbegin, _InputIterator freqend, _OutputIterator tsbegin, _OutputIterator tsend)
The member function to calculate a time series from (complex) spectra.
Definition: TsSpectrum.h:106
double GetDt() const
Return dt in s.
This class calculates one spectral element of the receiver function from the two input spectral eleme...
This class stores and writes model for the respktn 1D seismic code that we use for receiver function ...
Definition: ResPkModel.h:18
The basic exception class for all errors that arise in gplib.