AbsVelRecObjective.cpp

Go to the documentation of this file.
00001 #include "AbsVelRecObjective.h"
00002 #include <gsl/gsl_math.h>
00003 
00004 namespace gplib
00005   {
00006     AbsVelRecObjective::~AbsVelRecObjective()
00007       {
00008       }
00009 
00010     //constructor including absolute velocities
00011     AbsVelRecObjective::AbsVelRecObjective(boost::shared_ptr<
00012         const SeismicDataComp> TheRecData, SurfaceWaveData &AbsVel,
00013         const int myshift, const double mysigma, const double myc,
00014         const double myslowness, const RecCalc::trfmethod method,
00015         const bool normalized) :
00016       C1DRecObjective(TheRecData, myshift, mysigma, myc, myslowness, method,
00017           normalized), CalcAbsVel(mysigma, myc, method), MeasuredAbsVel(AbsVel)
00018       {
00019         recweight = 1;
00020         absvelweight = 1;
00021       }
00022 
00023     AbsVelRecObjective::AbsVelRecObjective(const AbsVelRecObjective &Old) :
00024       C1DRecObjective(Old), CalcAbsVel(Old.CalcAbsVel), MeasuredAbsVel(
00025           Old.MeasuredAbsVel), SynthAbsVel(Old.SynthAbsVel)
00026       {
00027         recweight = Old.recweight;
00028         absvelweight = Old.absvelweight;
00029       }
00030 
00031     AbsVelRecObjective& AbsVelRecObjective::operator=(
00032         const AbsVelRecObjective& source)
00033       {
00034         if (this == &source)
00035           return *this;
00036         C1DRecObjective::operator=(source);
00037         CalcAbsVel = source.CalcAbsVel;
00038         MeasuredAbsVel = source.MeasuredAbsVel;
00039         SynthAbsVel = source.SynthAbsVel;
00040         recweight = source.recweight;
00041         absvelweight = source.absvelweight;
00042         return *this;
00043       }
00044 
00045     double AbsVelRecObjective::PostParallel(const ttranscribed &member)
00046       {
00047         //call inherited version to calculate receiver function misfit
00048         C1DRecObjective::PostParallel(member);
00049         // calculate raw misfit from RMS
00050         double returnvalue = pow(GetRMS(), GetFitExponent())
00051             * GetMisfit().size();
00052         // and multiply by weight
00053         returnvalue *= recweight;
00054         //if we include absolute velocity information and have data
00055         if (absvelweight > 0.0 && !MeasuredAbsVel.GetPhaseVelocities().empty())
00056           {
00057             //calculate synthetic absolute velocity data
00058             SynthAbsVel.SetPeriods() = MeasuredAbsVel.GetPeriods();
00059             SeismicDataComp Rad(RecCalculator.GetRadComp()), Ver(
00060                 RecCalculator.GetVerComp());
00061             CalcAbsVel.CalcRFVel(slowness, Rad, Ver,
00062                 SynthAbsVel.SetPhaseVelocities());
00063 
00064             assert(SynthAbsVel.GetPhaseVelocities().size() == MeasuredAbsVel.GetPhaseVelocities().size());
00065             const unsigned int nvel = SynthAbsVel.GetPhaseVelocities().size();
00066             const unsigned int nrec = GetMisfit().size(); // the number of previous receiver function points
00067 
00068             //we have to copy the old misfit values from the base class
00069             //so we create local copties
00070             tmisfit LocalMisfit(nrec + nvel); // we need vectors of the right size for misfit, this is for levmarq
00071             tdata LocalData(nrec + nvel); // and  data, levmarq needs to know the size before, so we do not consider weights for the size
00072             copy(GetMisfit().begin(), GetMisfit().end(), LocalMisfit.begin());
00073             copy(GetSynthData().begin(), GetSynthData().end(),
00074                 LocalData.begin());
00075             //after reserving the right amount of memory, we can copy back
00076             SetMisfit(LocalMisfit);
00077             SetSynthData(LocalData);
00078 
00079             for (unsigned int i = 0; i < nvel; ++i)
00080               {
00081                 //we add up the misfit for each velocity estimate
00082                 //we have not stored any measurement error, so we only
00083                 //use the relative error floor, we store the misfit and data after the normal receiver function
00084                 returnvalue += CalcMisfit(
00085                     MeasuredAbsVel.GetPhaseVelocities().at(i),
00086                     SynthAbsVel.GetPhaseVelocities().at(i), 0.0, errorlevel, i
00087                         + nrec);
00088               }
00089 
00090             unsigned int ndata = 0;
00091             double divisor = 1.0;
00092             if (recweight > 0) // we have to divide by the right amount of datapoints
00093               {
00094                 ndata = nrec + nvel;
00095                 divisor = nrec * recweight + nvel * absvelweight;
00096               }
00097             else
00098               {
00099                 ndata = nvel;
00100                 divisor = ndata * absvelweight;
00101               }
00102             SetRMS(pow(returnvalue / divisor, 1.0 / GetFitExponent())); //store new misfit in local class
00103           }
00104         return GetRMS();
00105       }
00106   }

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