SurfaceWaveObjective.cpp

Go to the documentation of this file.

#include "SurfaceWaveObjective.h"
#include <boost/bind.hpp>
#include "SeisTools.h"
#include <gsl/gsl_math.h>
#include <numeric>
#include <boost/shared_ptr.hpp>
#include "Sdisp96Model.h"

using namespace std;

SurfaceWaveObjective::SurfaceWaveObjective(const SurfaceWaveData &Data) :
  MeasuredData(Data)
  {
    poisson=sqrt(3.);
    errorlevel = 0.02;
  }

SurfaceWaveObjective& SurfaceWaveObjective::operator=(const SurfaceWaveObjective& source)
  {
    if (this == &source) return *this;
    MeasuredData = source.MeasuredData;
    SynthData = source.SynthData;
    Synthetic = source.Synthetic;
    errorlevel = source.errorlevel;
    poisson = source.poisson;
    return *this;
  }

SurfaceWaveObjective::SurfaceWaveObjective(const SurfaceWaveObjective &Old) :
  MeasuredData(Old.MeasuredData), SynthData(Old.SynthData),
      Synthetic(Old.Synthetic)
  {
    errorlevel = Old.errorlevel;
    poisson = Old.poisson;
  }

SurfaceWaveObjective::~SurfaceWaveObjective()
  {
  }

void SurfaceWaveObjective::PreParallel(const ttranscribed &member)
  {
    std::string temp(tmpnam(NULL));
    tempname = temp;
    const unsigned int nfulllayers = member.size()/2;
    boost::shared_ptr<Sdisp96Model> Model(new Sdisp96Model);
    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());

    Model->SetThicknesses().assign(nfulllayers, 0.0);
    Model->SetShVelocities().assign(nfulllayers, 0.0);
    Model->SetSvVelocities().assign(nfulllayers, 0.0);
    Model->SetPhVelocities().assign(nfulllayers, 0.0);
    Model->SetPvVelocities().assign(nfulllayers, 0.0);
    Model->SetDensities().assign(nfulllayers, 0.0);
    Model->SetEta().assign(nfulllayers, 0.0);
    Model->SetQmu().assign(nfulllayers, 0.0);
    Model->SetQkappa().assign(nfulllayers, 0.0);

    copy(thickness.begin(), thickness.end(), Model->SetThicknesses().begin());
    copy(velocity.begin(), velocity.end(), Model->SetShVelocities().begin());
    copy(velocity.begin(), velocity.end(), Model->SetSvVelocities().begin());
    transform(velocity.begin(), velocity.end(), Model->SetDensities().begin(), CalcDensity());
    transform(velocity.begin(), velocity.end(), //PVelocity is calculated by multiplication of SVel by poisson's ratio
        Model->SetPhVelocities().begin(), boost::bind(multiplies<double>(), _1, poisson));
    copy(Model->GetPhVelocities().begin(), Model->GetPhVelocities().end(), Model->SetPvVelocities().begin());
    fill_n(Model->SetEta().begin(), nfulllayers, 1.0);
    fill_n(Model->SetQmu().begin(), nfulllayers, 100.0);
    fill_n(Model->SetQkappa().begin(), nfulllayers, 99999.00);
    //  for (int i = 0; i < Model->GetThicknesses().size(); ++i) // vector size changes, so we have to use the current size 
    //  {
    //          double currdepth = accumulate(Model->GetThicknesses().begin(),Model->GetThicknesses().begin()+i,0.0);
    //          double maxdepth = Model->GetMaxDepth(currdepth);
    //          i = Model->SplitLayer(i,maxdepth); //this function manipulates the value of i
    //  }
    //cout << "After: " << Model.GetThicknesses().size() << endl;
    //Model.MergeModel(Background);
    Synthetic.SetModel(Model);
    if (Synthetic.GetCalculationPeriods().empty())
      {
        Synthetic.SetCalculationPeriods(MeasuredData.GetPeriods());
      }
    Synthetic.PreParallel(tempname);
  }

double SurfaceWaveObjective::PostParallel(const ttranscribed &member)
  {
    return GetRMS();

  }

void SurfaceWaveObjective::SafeParallel(const ttranscribed &member)
  {

    SynthData = Synthetic.SafeParallel(tempname);
    assert(errorlevel> 1e-4);
    // an errorlevel < 0.01% is absolutely unrealistic
    const unsigned int ndata = min(SynthData.GetPhaseVelocities().size(), MeasuredData.GetPhaseVelocities().size());
    tmisfit LocalMisfit(ndata); // we need vectors of the right size for misfit
    tdata LocalData(ndata); // and  data
    double returnvalue = 0; // init returnvalue
    double currvalue = 0.0;
    for (unsigned int i = 0; i < ndata; ++i)
      {
        currvalue = SynthData.GetPhaseVelocities().at(i) - MeasuredData.GetPhaseVelocities().at(i);
        currvalue /= MeasuredData.GetPhaseVelocities().at(i) * errorlevel;
        currvalue = gsl_pow_int(currvalue, GetFitExponent());
        LocalMisfit(i)= fabs(currvalue);
        LocalData(i)= SynthData.GetPhaseVelocities().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())); //return misfit
  }

---