SeismicModelDiff.cpp

Go to the documentation of this file.

#include "SeismicModelDiff.h"
#include <gsl/gsl_math.h>
#include <iostream>
#include <numeric>
using namespace std;

namespace gplib
  {
    SeismicModelDiff::SeismicModelDiff(const ResPkModel &Seis) :
      Model(Seis)
      {
      }

    SeismicModelDiff::~SeismicModelDiff()
      {
      }

    SeismicModelDiff::SeismicModelDiff(const SeismicModelDiff &Old) :
      GeneralObjective(Old), Model(Old.Model)
      {
      }

    SeismicModelDiff& SeismicModelDiff::operator=(
        const SeismicModelDiff& source)
      {
        if (this == &source)
          return *this;
        GeneralObjective::operator=(source);
        Model = source.Model;
        return *this;
      }

    //! Calculate the roughness of the model given by the parameter member
    void SeismicModelDiff::SafeParallel(const ttranscribed &member)
      {
        const unsigned int length = member.size() / 3; //we have 3 parameters in the model, so size/3 layers
        double roughness = 0.0; // init returnvalue

        vector<double> memberdepth, refdepth;
        partial_sum(member.begin() + length, member.begin() + 2 * length,
            back_inserter(memberdepth));
        partial_sum(Model.GetThickness().begin(), Model.GetThickness().end(),
            back_inserter(refdepth));
        refdepth.back() += 1000.0; // the last layer is a halfspace
        vector<double> layerthick(memberdepth); //vector to contain the depths to interfaces in both models, contruct from inversion model
        copy(refdepth.begin(), refdepth.end(), back_inserter(layerthick)); // and insert reference model interfaces at back
        sort(layerthick.begin(), layerthick.end()); //then put them in right order
        adjacent_difference(layerthick.begin(), layerthick.end(),
            layerthick.begin());
        layerthick.erase(remove(layerthick.begin(), layerthick.end(), 0.0),
            layerthick.end()); //remove layers of 0 thickness
        unsigned int refindex = 0;
        unsigned int memberindex = 0;
        unsigned int thickindex = 0;
        double currthick = 0.0;
        while (refindex < Model.GetThickness().size() && memberindex < length)
          {
            int difference = gsl_fcmp(memberdepth.at(memberindex), refdepth.at(
                refindex), GSL_DBL_EPSILON);
            currthick = layerthick.at(thickindex);
            roughness += gsl_pow_2((member(memberindex + 2 * length)
                - Model.GetSVelocity().at(refindex)) / Model.GetSVelocity().at(
                refindex) * currthick);
            switch (difference)
              {
            case 1: // memberdepth bigger

              ++refindex;
              break;
            case -1:

              ++memberindex;
              break;
            case 0:

              ++refindex;
              ++memberindex;
              break;
              }
            thickindex++;
          }
        SetRMS(roughness);
      }

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