C1DRecObjective.cpp

Go to the documentation of this file.
00001 #include <iostream>
00002 #include <numeric>
00003 #include <algorithm>
00004 #include <boost/bind.hpp>
00005 #include <boost/cast.hpp>
00006 #include <gsl/gsl_math.h>
00007 #include "CollapseModel.h"
00008 #include "SeisTools.h"
00009 #include "convert.h"
00010 #include "C1DRecObjective.h"
00011 #include "miscfunc.h"
00012 #include "FatalException.h"
00013 
00014 using namespace std;
00015 
00016 namespace gplib
00017   {
00018     C1DRecObjective::C1DRecObjective(
00019         boost::shared_ptr<const SeismicDataComp> TheRecData, const int myshift,
00020         const double mysigma, const double myc, const double myslowness,
00021         const RecCalc::trfmethod method, const bool normalized) :
00022       RecCalculator(myshift, mysigma, myc, true, method), ObservedData(
00023           TheRecData)
00024       {
00025         slowness = myslowness; //copy parameter value
00026         poisson = sqrt(3.); //set some default values
00027         //we specify a default errorlevel of 1%
00028         SetErrorLevel(0.01);
00029         //by default we fit the complete receiver function
00030         startpoint = 0;
00031         endpoint = ObservedData->GetData().size();
00032         // we determine the maximum amplitude of the receiver function
00033 
00034         const double recmaxamp = *max_element(ObservedData->GetData().begin(),
00035             ObservedData->GetData().end());
00036         const bool isone = gsl_fcmp(recmaxamp, 1.0, 1e-4) == 0;
00037         if (isone && !normalized)
00038           throw FatalException(
00039               "Input data has maximum amplitude of 1, but no normalization used");
00040         if (!isone && normalized)
00041           throw FatalException(
00042               "Input data does not have an amplitude of 1, but normalization is used");
00043         RecCalculator.SetNormalize(normalized);// we can choose we want to model normalized data
00044       }
00045 
00046     C1DRecObjective::~C1DRecObjective()
00047       {
00048       }
00049 
00050     C1DRecObjective::C1DRecObjective(const C1DRecObjective &Old) :
00051       PlottableObjective(Old), RecSynthData(Old.RecSynthData), slowness(
00052           Old.slowness), Model(Old.Model), RecCalculator(Old.RecCalculator),
00053           errorlevel(Old.errorlevel), errorvalue(Old.errorvalue), poisson(
00054               Old.poisson), startpoint(Old.startpoint), endpoint(Old.endpoint),
00055           ObservedData(Old.ObservedData)
00056       {
00057       }
00058 
00059     void C1DRecObjective::SetTimeWindow(const double start, const double end)
00060       {
00061         //calculate the startpoint corresponding to the starttime
00062         //if the values do not make sense, we use the the minimum
00063         // and maximum possible indices respectively
00064         startpoint = max(boost::numeric_cast<size_t>(0), boost::numeric_cast<
00065             size_t>((start - ObservedData->GetB()) / ObservedData->GetDt()));
00066         // calculate datapoint corresponding to the endtime
00067         endpoint = min(boost::numeric_cast<size_t>((end - ObservedData->GetB())
00068             / ObservedData->GetDt()), ObservedData->GetData().size());
00069         //the only thing that remains to be checked is whether the start is before the end
00070         if (endpoint < startpoint)
00071           throw FatalException("Start index " + stringify(startpoint)
00072               + " and end index " + stringify(endpoint)
00073               + " do not make sense !");
00074       }
00075 
00076     C1DRecObjective& C1DRecObjective::operator=(const C1DRecObjective& source)
00077       {
00078         if (this == &source)
00079           return *this;
00080         PlottableObjective::operator=(source);
00081         RecSynthData = source.RecSynthData;
00082         slowness = source.slowness;
00083         Model = source.Model;
00084         RecCalculator = source.RecCalculator;
00085         poisson = source.poisson;
00086         errorlevel = source.errorlevel;
00087         errorvalue = source.errorvalue;
00088         startpoint = source.startpoint;
00089         endpoint = source.endpoint;
00090         ObservedData = source.ObservedData;
00091         return *this;
00092       }
00093 
00094     void C1DRecObjective::PreParallel(const ttranscribed &member)
00095       {
00096         const unsigned int nfulllayers = member.size() / 2; //The model vector contains thickness and SVel, so we have size/2 layers
00097 
00098         Model.SetSlowness(slowness); //copy slowness
00099         ttranscribed thickness(nfulllayers), velocity(nfulllayers);
00100         copy(member.begin(), member.begin() + nfulllayers, thickness.begin()); //copy values
00101         copy(member.begin() + nfulllayers, member.begin() + 2 * nfulllayers,
00102             velocity.begin());
00103         CollapseModel(thickness, velocity);//reduce to minimum number of different layers
00104         const unsigned int ncolllayers = velocity.size();
00105         Model.Init(ncolllayers); //initialize model with number of reduced layers
00106         copy(thickness.begin(), thickness.end(), Model.SetThickness().begin());
00107         copy(velocity.begin(), velocity.end(), Model.SetSVelocity().begin());
00108         transform(velocity.begin(), velocity.end(), Model.SetDensity().begin(),
00109             CalcDensity());
00110         transform(velocity.begin(), velocity.end(), //PVelocity is calculated by multiplication of SVel by poisson's ratio
00111             Model.SetPVelocity().begin(), boost::bind(multiplies<double> (),
00112                 _1, poisson));
00113         Model.SetDt(ObservedData->GetDt()); // set dt
00114         Model.SetNpts(ObservedData->GetData().size()); //set number of points
00115         RecCalculator.SynthPreParallel(GetParallelID(), Model, RecSynthData,
00116             true);
00117       }
00118 
00119     double C1DRecObjective::PostParallel(const ttranscribed &member)
00120       {
00121         RecCalculator.SynthPostParallel(GetParallelID(), Model, RecSynthData,
00122             true);
00123         SetMisfit().resize(endpoint - startpoint);//we need vectors of the right size for misfit
00124         SetSynthData().resize(endpoint - startpoint);//and  data
00125         double returnvalue = 0.0;//init returnvalue
00126         //go through the receiver function from the starttime to the endtime point by point
00127         for (int i = startpoint; i < endpoint; ++i)
00128 
00129           {
00130             //calculate the misfit for each point and add up
00131             //the relative errorlevel is 0. as it does not make sense for receiver functions
00132             //therefore always errorvalue will be used as an absolue error
00133             returnvalue += CalcMisfit(ObservedData->GetData().at(i),
00134                 RecSynthData.GetData().at(i), errorvalue, 0.0, i - startpoint);
00135           }
00136         //copy the head information from the observed receiver function
00137         RecSynthData.CopyHeader(*ObservedData.get());
00138         SetRMS(pow(returnvalue / (endpoint - startpoint), 1.0
00139             / GetFitExponent())); //store misfit in local class
00140         return GetRMS(); //return misfit calculated by SafeParallel
00141       }
00142 
00143     void C1DRecObjective::SafeParallel(const ttranscribed &member)
00144       {
00145         RecCalculator.SynthSafeParallel(GetParallelID(), Model, RecSynthData,
00146             true); // Calculate forward model
00147       }
00148   }

Generated on Tue Nov 3 13:24:13 2009 for GPLIB++ by  doxygen 1.5.8