C1DRecObjective.cpp

Go to the documentation of this file.

#include "C1DRecObjective.h"
#include "miscfunc.h"
#include <iostream>
#include <numeric>
#include <algorithm>
#include <cstdio>
#include <boost/bind.hpp>
#include <boost/cast.hpp>
#include <gsl/gsl_math.h>
#include "CollapseModel.h"
#include "SeisTools.h"

using namespace std;

//constructor without absolute velocities
C1DRecObjective::C1DRecObjective(boost::shared_ptr<const SeismicDataComp> TheRecData,
    const int myshift, const double myomega, const double mysigma,
    const double myc, const double myslowness, const RecCalc::trfmethod method) :
  RecCalculator(myshift, myomega, mysigma, myc, false, method),
      ObservedData(TheRecData)
  {
    slowness = myslowness; //copy parameter value 
    poisson = sqrt(3.); //set some default values
    errorlevel = 0.01;
    starttime = -1; //these don't make sense, will be checked 
    endtime = -1; // when calculating performance
    RecCalculator.SetNormalize(true);// we normalize the input data and the synthetics to ensure consistency
  }

C1DRecObjective::~C1DRecObjective()
  {
  }

C1DRecObjective::C1DRecObjective(const C1DRecObjective &Old) :
  PlottableObjective(Old), RecSynthData(Old.RecSynthData),
      slowness(Old.slowness), Model(Old.Model),
      RecCalculator(Old.RecCalculator), errorlevel(Old.errorlevel),
      poisson(Old.poisson), starttime(Old.starttime), endtime(Old.endtime),
      ObservedData(Old.ObservedData)
  {
  }

C1DRecObjective& C1DRecObjective::operator= (const C1DRecObjective& source)
  {
    if (this == &source) return *this;
    PlottableObjective::operator=(source);
    RecSynthData = source.RecSynthData;
    slowness = source.slowness;
    Model = source.Model;
    RecCalculator = source.RecCalculator;
    poisson = source.poisson;
    errorlevel = source.errorlevel;
    starttime = source.starttime;
    endtime = source.endtime;
    ObservedData = source.ObservedData;
    return *this;
  }

void C1DRecObjective::PreParallel(const ttranscribed &member)
  {
    const unsigned int nfulllayers = member.size()/2; //The model vector contains thickness and SVel, so we have size/2 layers

    Model.SetSlowness(slowness); //copy slowness
    ttranscribed thickness(nfulllayers), velocity(nfulllayers);
    copy(member.begin(), member.begin()+nfulllayers, thickness.begin()); //copy values
    copy(member.begin()+nfulllayers, member.begin()+2*nfulllayers,
        velocity.begin());
    CollapseModel(thickness, velocity);//reduce to minimum number of different layers
    const unsigned int ncolllayers = velocity.size();
    Model.Init(ncolllayers); //initialize model with number of reduced layers
    copy(thickness.begin(), thickness.end(), Model.SetThickness().begin());
    copy(velocity.begin(), velocity.end(), Model.SetSVelocity().begin());
    transform(velocity.begin(), velocity.end(), Model.SetDensity().begin(), CalcDensity());
    transform(velocity.begin(), velocity.end(), //PVelocity is calculated by multiplication of SVel by poisson's ratio
        Model.SetPVelocity().begin(), boost::bind(multiplies<double>(), _1, poisson));
    Model.SetDt(ObservedData->GetDt()); // set dt
    Model.SetNpts(ObservedData->GetData().size()); //set number of points
    std::string temp(tmpnam(NULL));
    tempname = temp;
    RecCalculator.SynthPreParallel(tempname, Model, RecSynthData, true);
  }

double C1DRecObjective::PostParallel(const ttranscribed &member)
  {
    RecCalculator.SynthPostParallel(tempname, Model, RecSynthData, true);
    return GetRMS(); //return misfit calculated by SafeParallel
  }

void C1DRecObjective::SafeParallel(const ttranscribed &member)
  {
    RecCalculator.SynthSafeParallel(tempname, Model, RecSynthData, true); // Calculate forward model

    size_t startpoint = 0; //default start is begin
    size_t endpoint = ObservedData->GetData().size(); //and default end is end of data
    if (starttime > 0) //if starttime was set
      {
        startpoint = boost::numeric_cast<size_t>((starttime-ObservedData->GetB())/ObservedData->GetDt());
        // calculate datapoint corresponding to starttime
          }
        if (endtime > 0) //and same for endtime

          {
            endpoint = min(boost::numeric_cast<size_t>((endtime-ObservedData->GetB())/ObservedData->GetDt()),ObservedData->GetData().size());
          }
        assert(endpoint > startpoint);
        tmisfit LocalMisfit(endpoint-startpoint);//we need vectors of the right size for misfit
        tdata LocalData(endpoint-startpoint);//and  data
        double returnvalue = 0;//init returnvalue
        double currvalue = 0.0;
        for (size_t i = startpoint; i < endpoint; ++i)//for the desired range

          {
            currvalue = ObservedData->GetData().at(i) - RecSynthData.GetData().at(i); //calculate  difference
            currvalue /= errorlevel; //adjust for error
            currvalue = gsl_pow_int(currvalue,GetFitExponent());
            //cout << ObservedData.GetData().at(i) << " " <<  RecSynthData.GetData().at(i) << " " << currvalue << endl;
            assert(gsl_finite(currvalue));
            LocalMisfit(i-startpoint)= fabs(currvalue);
            LocalData(i-startpoint)= RecSynthData.GetData().at(i); //copy data (for levmar algorithm)
            returnvalue += fabs(currvalue); // add misfit to returnvalue
          }
        SetMisfit(LocalMisfit); // set appropriate members of baseclass
        SetSynthData(LocalData);

        SetRMS(pow(returnvalue/LocalMisfit.size(),1.0/GetFitExponent())); //store misfit in local class
      }

---