PTensor1DMTObjective.cpp

Go to the documentation of this file.

#include "PTensor1DMTObjective.h"
#include "Member2Aniso.h"
#include <gsl/gsl_math.h>
#include <cassert>
#include <numeric>
using namespace std;

namespace gplib
  {
    PTensor1DMTObjective::PTensor1DMTObjective(
        const PTensorMTStation &LocalMTData) :
      errorlevel(0.01), // we set a standard error level of 1%
          MeasuredData(LocalMTData)
      {
      }

    PTensor1DMTObjective::~PTensor1DMTObjective()
      {
      }

    PTensor1DMTObjective::PTensor1DMTObjective(const PTensor1DMTObjective &Old) :
      PlottableObjective(Old), errorlevel(Old.errorlevel), MeasuredData(
          Old.MeasuredData), AnisoMTSynth(Old.AnisoMTSynth)
      {
      }

    PTensor1DMTObjective& PTensor1DMTObjective::operator=(
        const PTensor1DMTObjective& source)
      {
        if (this == &source)
          return *this;
        PlottableObjective::operator=(source);
        errorlevel = source.errorlevel;
        MeasuredData = source.MeasuredData;
        AnisoMTSynth = source.AnisoMTSynth;
        return *this;
      }

    void PTensor1DMTObjective::SafeParallel(const ttranscribed &member)
      {
        if (AnisoMTSynth.GetFrequencies().empty()) //if we didn't set the frequencies before
          {
            AnisoMTSynth.SetFrequencies(MeasuredData.GetFrequencies());
          }
        Member2Aniso(member, AnisoMTSynth);
        AnisoMTSynth.GetData();

        double returnvalue = 0; // init misfit value

        const unsigned int nelements = 4;
        const unsigned int ndata = nelements * MeasuredData.GetTensor().size(); //we fit 4 real phase tensor elements
        SetMisfit().resize(ndata); //make sure Misfit in base class can hold enough elements
        SetSynthData().resize(ndata); // and same for data

        for (unsigned int i = 0; i < ndata; i += nelements) //we do the 4 elements in one block
          {
            returnvalue += CalcMisfit(
                MeasuredData.at(i / nelements).GetPhi11(), AnisoMTSynth.at(i
                    / nelements).GetPhi11(),
                MeasuredData.at(i / nelements).GetdPhi11(), errorlevel, i);
            returnvalue += CalcMisfit(
                MeasuredData.at(i / nelements).GetPhi12(), AnisoMTSynth.at(i
                    / nelements).GetPhi12(),
                MeasuredData.at(i / nelements).GetdPhi12(), errorlevel, i + 1);
            returnvalue += CalcMisfit(
                MeasuredData.at(i / nelements).GetPhi21(), AnisoMTSynth.at(i
                    / nelements).GetPhi21(),
                MeasuredData.at(i / nelements).GetdPhi21(), errorlevel, i + 2);
            returnvalue += CalcMisfit(
                MeasuredData.at(i / nelements).GetPhi22(), AnisoMTSynth.at(i
                    / nelements).GetPhi22(),
                MeasuredData.at(i / nelements).GetdPhi22(), errorlevel, i + 3);
          }
        SetRMS(pow(returnvalue / ndata, 1.0 / GetFitExponent()));
      }

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