C1DRecObjective.cpp

Go to the documentation of this file.
00001 #include "C1DRecObjective.h"
00002 #include "miscfunc.h"
00003 #include <iostream>
00004 #include <numeric>
00005 #include <algorithm>
00006 #include <cstdio>
00007 #include <boost/bind.hpp>
00008 #include <boost/cast.hpp>
00009 #include <gsl/gsl_math.h>
00010 #include "CollapseModel.h"
00011 #include "SeisTools.h"
00012 
00013 using namespace std;
00014 
00015 //constructor without absolute velocities
00016 C1DRecObjective::C1DRecObjective(boost::shared_ptr<const SeismicDataComp> TheRecData,
00017     const int myshift, const double myomega, const double mysigma,
00018     const double myc, const double myslowness, const RecCalc::trfmethod method) :
00019   RecCalculator(myshift, myomega, mysigma, myc, false, method),
00020       ObservedData(TheRecData)
00021   {
00022     slowness = myslowness; //copy parameter value 
00023     poisson = sqrt(3.); //set some default values
00024     errorlevel = 0.01;
00025     starttime = -1; //these don't make sense, will be checked 
00026     endtime = -1; // when calculating performance
00027     RecCalculator.SetNormalize(true);// we normalize the input data and the synthetics to ensure consistency
00028   }
00029 
00030 C1DRecObjective::~C1DRecObjective()
00031   {
00032   }
00033 
00034 C1DRecObjective::C1DRecObjective(const C1DRecObjective &Old) :
00035   PlottableObjective(Old), RecSynthData(Old.RecSynthData),
00036       slowness(Old.slowness), Model(Old.Model),
00037       RecCalculator(Old.RecCalculator), errorlevel(Old.errorlevel),
00038       poisson(Old.poisson), starttime(Old.starttime), endtime(Old.endtime),
00039       ObservedData(Old.ObservedData)
00040   {
00041   }
00042 
00043 C1DRecObjective& C1DRecObjective::operator= (const C1DRecObjective& source)
00044   {
00045     if (this == &source) return *this;
00046     PlottableObjective::operator=(source);
00047     RecSynthData = source.RecSynthData;
00048     slowness = source.slowness;
00049     Model = source.Model;
00050     RecCalculator = source.RecCalculator;
00051     poisson = source.poisson;
00052     errorlevel = source.errorlevel;
00053     starttime = source.starttime;
00054     endtime = source.endtime;
00055     ObservedData = source.ObservedData;
00056     return *this;
00057   }
00058 
00059 void C1DRecObjective::PreParallel(const ttranscribed &member)
00060   {
00061     const unsigned int nfulllayers = member.size()/2; //The model vector contains thickness and SVel, so we have size/2 layers
00062 
00063     Model.SetSlowness(slowness); //copy slowness
00064     ttranscribed thickness(nfulllayers), velocity(nfulllayers);
00065     copy(member.begin(), member.begin()+nfulllayers, thickness.begin()); //copy values
00066     copy(member.begin()+nfulllayers, member.begin()+2*nfulllayers,
00067         velocity.begin());
00068     CollapseModel(thickness, velocity);//reduce to minimum number of different layers
00069     const unsigned int ncolllayers = velocity.size();
00070     Model.Init(ncolllayers); //initialize model with number of reduced layers
00071     copy(thickness.begin(), thickness.end(), Model.SetThickness().begin());
00072     copy(velocity.begin(), velocity.end(), Model.SetSVelocity().begin());
00073     transform(velocity.begin(), velocity.end(), Model.SetDensity().begin(), CalcDensity());
00074     transform(velocity.begin(), velocity.end(), //PVelocity is calculated by multiplication of SVel by poisson's ratio
00075         Model.SetPVelocity().begin(), boost::bind(multiplies<double>(), _1, poisson));
00076     Model.SetDt(ObservedData->GetDt()); // set dt
00077     Model.SetNpts(ObservedData->GetData().size()); //set number of points
00078     std::string temp(tmpnam(NULL));
00079     tempname = temp;
00080     RecCalculator.SynthPreParallel(tempname, Model, RecSynthData, true);
00081   }
00082 
00083 double C1DRecObjective::PostParallel(const ttranscribed &member)
00084   {
00085     RecCalculator.SynthPostParallel(tempname, Model, RecSynthData, true);
00086     return GetRMS(); //return misfit calculated by SafeParallel
00087   }
00088 
00089 void C1DRecObjective::SafeParallel(const ttranscribed &member)
00090   {
00091     RecCalculator.SynthSafeParallel(tempname, Model, RecSynthData, true); // Calculate forward model
00092 
00093     size_t startpoint = 0; //default start is begin
00094     size_t endpoint = ObservedData->GetData().size(); //and default end is end of data
00095     if (starttime > 0) //if starttime was set
00096       {
00097         startpoint = boost::numeric_cast<size_t>((starttime-ObservedData->GetB())/ObservedData->GetDt());
00098         // calculate datapoint corresponding to starttime
00099           }
00100         if (endtime > 0) //and same for endtime
00101 
00102           {
00103             endpoint = min(boost::numeric_cast<size_t>((endtime-ObservedData->GetB())/ObservedData->GetDt()),ObservedData->GetData().size());
00104           }
00105         assert(endpoint > startpoint);
00106         tmisfit LocalMisfit(endpoint-startpoint);//we need vectors of the right size for misfit
00107         tdata LocalData(endpoint-startpoint);//and  data
00108         double returnvalue = 0;//init returnvalue
00109         double currvalue = 0.0;
00110         for (size_t i = startpoint; i < endpoint; ++i)//for the desired range
00111 
00112           {
00113             currvalue = ObservedData->GetData().at(i) - RecSynthData.GetData().at(i); //calculate  difference
00114             currvalue /= errorlevel; //adjust for error
00115             currvalue = gsl_pow_int(currvalue,GetFitExponent());
00116             //cout << ObservedData.GetData().at(i) << " " <<  RecSynthData.GetData().at(i) << " " << currvalue << endl;
00117             assert(gsl_finite(currvalue));
00118             LocalMisfit(i-startpoint)= fabs(currvalue);
00119             LocalData(i-startpoint)= RecSynthData.GetData().at(i); //copy data (for levmar algorithm)
00120             returnvalue += fabs(currvalue); // add misfit to returnvalue
00121           }
00122         SetMisfit(LocalMisfit); // set appropriate members of baseclass
00123         SetSynthData(LocalData);
00124 
00125         SetRMS(pow(returnvalue/LocalMisfit.size(),1.0/GetFitExponent())); //store misfit in local class
00126       }

Generated on Fri Jul 4 15:30:20 2008 for GPLIB++ by  doxygen 1.5.5