AnisoSurfaceWaveObjective.cpp

Go to the documentation of this file.

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

using namespace std;

namespace gplib
  {
    AnisoSurfaceWaveObjective::AnisoSurfaceWaveObjective(
        const ParkSurfaceWaveData &Data, const double ba) :
      MeasuredData(Data)
      {
        //Define here the ratio between E and B
        avelratio = 1.0;
        backazimuth = ba;
        poisson = sqrt(3.);
        errorlevel = 0.02;
      }

    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;
        AnisoSurfaceWaveModel Model;
        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(), //PVelocity 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 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; ++i)
          {
            returnvalue += CalcMisfit(MeasuredData.GetPhaseVelocities().at(i),
                SynthData.GetPhaseVelocities().at(i), 0.0, errorlevel, i);
          }

        SetRMS(pow(returnvalue / ndata, 1.0 / GetFitExponent())); //return misfit
      }
  }