AnisoSurfaceWaveObjective.cpp

Go to the documentation of this file.

#include "AnisoSurfaceWaveObjective.h"
#include "AnisoSurfaceWaveModel.h"
#include <boost/bind.hpp>
#include <cmath>
#include "SeisTools.h"

using namespace std;

namespace gplib
  {
AnisoSurfaceWaveObjective::AnisoSurfaceWaveObjective(
    const ParkSurfaceWaveData &Data, const double ba, const double avel) :
  MeasuredData(Data)
  {
    avelratio = avel;
    backazimuth = ba;
    poisson = sqrt(3.);
    errorlevel = 0.01;
  }

AnisoSurfaceWaveObjective::~AnisoSurfaceWaveObjective()
  {
  }

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

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

void AnisoSurfaceWaveObjective::PreParallel(const ttranscribed &member)
  {
    const unsigned int nfulllayers = member.size() / 6;
    boost::shared_ptr<AnisoSurfaceWaveModel> Model(new AnisoSurfaceWaveModel);
    ttranscribed thickness(nfulllayers), velocity(nfulllayers), angle(
        nfulllayers), deltaangle(nfulllayers), anisotropy1(nfulllayers),
        anisotropy2(nfulllayers);
    copy(member.begin(), member.begin() + nfulllayers, thickness.begin()); //copy values
    copy(member.begin() + nfulllayers, member.begin() + 2 * nfulllayers,
        velocity.begin());
    transform(member.begin() + 2 * nfulllayers, member.begin() + 3
        * nfulllayers, angle.begin(), boost::bind(std::minus<double>(), _1,
        backazimuth));
    copy(member.begin() + 3 * nfulllayers, member.begin() + 4 * nfulllayers,
        deltaangle.begin());
    copy(member.begin() + 4 * nfulllayers, member.begin() + 5 * nfulllayers,
        anisotropy1.begin());
    copy(member.begin() + 5 * nfulllayers, member.begin() + 6 * nfulllayers,
        anisotropy2.begin());

    copy(thickness.begin(), thickness.end(), back_inserter(
        Model->SetThicknesses()));
    copy(velocity.begin(), velocity.end(), back_inserter(Model->SetVs()));
    transform(velocity.begin(), velocity.end(), //P Velocity is calculated by multiplication of SVel by poisson's ratio
        back_inserter(Model->SetVp()), boost::bind(multiplies<double> (), _1,
            poisson));
    transform(velocity.begin(), velocity.end(), back_inserter(
        Model->SetDensities()), CalcDensity());
    fill_n(back_inserter(Model->SetTheta()), nfulllayers, 90.0);
    fill_n(back_inserter(Model->SetC()), nfulllayers, 0.0);
    // Set B anisotropy to product of anisotropy1 and anisotropy2
    transform(anisotropy1.begin(), anisotropy1.end(), anisotropy2.begin(),
        back_inserter(Model->SetB()), std::multiplies<double>());
    // Set B = anisotropy2
//    copy(anisotropy2.begin(), anisotropy2.end(), back_inserter(Model->SetB()));
    // Set E to product of B and avelratio
    transform(Model->GetB().begin(), Model->GetB().end(), back_inserter(
        Model->SetE()), boost::bind(multiplies<double> (), _1, avelratio));
    //Set anisotropy angle to sum of angle and deltaangle
    transform(angle.begin(), angle.end(), deltaangle.begin(), back_inserter(
        Model->SetPhi()), std::plus<double>());

    Synthetic.SetModel(Model);
    Synthetic.PreParallel(GetParallelID());
  }

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

void AnisoSurfaceWaveObjective::SafeParallel(const ttranscribed &member)
  {
    SynthData = Synthetic.SafeParallel(GetParallelID());
    assert(errorlevel> 1e-4);
    // an errorlevel < 0.01% is absolutely unrealistic
    const unsigned int ndata = min(SynthData.GetPhaseVelocities().size(),
        MeasuredData.GetPhaseVelocities().size());
    SetMisfit().resize(ndata); // we need vectors of the right size for misfit
    SetSynthData().resize(ndata); // and  data
    double returnvalue = 0.0; // init returnvalue
    for (unsigned int i = 0; i < ndata -2 ; ++i)
      {
        returnvalue += CalcMisfit(MeasuredData.GetPhaseVelocities().at(i),
            SynthData.GetPhaseVelocities().at(i), 0.0, errorlevel, i);
       }

    SetRMS(std::pow(returnvalue / (ndata -2), 1.0 / GetFitExponent())); //return misfit
  }
}