cadianiso.cpp

Go to the documentation of this file.

#include "cadijoint.h"
#include <iostream>
#include "Aniso1DMTObjective.h"
#include "MTStation.h"
#include "gentypes.h"
#include "CLevanisoConf.h"
#include "MTFitSetup.h"
#include <boost/bind.hpp>
#include "PTensor1DMTObjective.h"
#include <boost/shared_ptr.hpp>
#include "PTensorMTStation.h"

using namespace std;
using namespace gplib;

boost::shared_ptr<PlottableObjective> MTObjective;
CLevanisoConf Configuration;

void misfit(float *p, float *misfit, int *m)
  {
    const int nlayers = abs(*m) / 4;
    const int nparam = nlayers * 4;
    ttranscribed mtmember(nparam);
    for (int i = 0; i < nparam; ++i)
      {
        mtmember(i) = p[i];
      }
    *misfit = 0.0;
    if (Configuration.weights.at(0) > 0)
      {
        //cout << "Calculating MT misfit" << endl;
        *misfit += Configuration.weights.at(0) * MTObjective->CalcPerformance(
            mtmember);
      }
  }

void init(int* nd, float* ranges)
  {
    MTStation MTData;
    PTensorMTStation PTData;
    Configuration.GetData("levaniso.conf");
    try
      {
        MTData.GetData(Configuration.mtinputdata);
      } catch (FatalException &e)
      {
        cerr << e.what() << endl;
      }

    if (Configuration.mtfit == "ptensor")
      {

        PTData.GetData(Configuration.ptensordata);
        boost::shared_ptr<PTensor1DMTObjective> PTObjective(
            new PTensor1DMTObjective(PTData));
        MTObjective = PTObjective;
      }
    else
      {

        MTData.GetData(Configuration.mtinputdata);
        boost::shared_ptr<Aniso1DMTObjective> AnisoObjective(
            new Aniso1DMTObjective(MTData));
        SetupMTFitParameters(Configuration, *AnisoObjective);
        MTObjective = AnisoObjective;
      }

    const int nlayers = Configuration.minres.size();
    const int nparams = nlayers * 4;
    *nd = nparams;
    ttranscribed mtmember(nparams);
    copy(Configuration.startres.begin(), Configuration.startres.end(),
        mtmember.begin());
    copy(Configuration.startthick.begin(), Configuration.startthick.end(),
        mtmember.begin() + nlayers);
    copy(Configuration.startaniso.begin(), Configuration.startaniso.end(),
        mtmember.begin() + 2 * nlayers);
    copy(Configuration.startstrike.begin(), Configuration.startstrike.end(),
        mtmember.begin() + 3 * nlayers);
    MTObjective->CalcPerformance(mtmember);

    for (int i = 0; i < nlayers; ++i)
      {
        ranges[2 * i] = Configuration.minres.at(i);
        ranges[2 * i + 1] = Configuration.maxres.at(i);
      }
    for (int i = 0; i < nlayers; ++i)
      {
        ranges[2 * i + 2 * nlayers] = Configuration.minthick.at(i);
        ranges[2 * i + 2 * nlayers + 1] = Configuration.maxthick.at(i);
      }
    for (int i = 0; i < nlayers; ++i)
      {
        ranges[2 * i + 4 * nlayers] = Configuration.minaniso.at(i);
        ranges[2 * i + 4 * nlayers + 1] = Configuration.maxaniso.at(i);
      }
    for (int i = 0; i < nlayers; ++i)
      {
        ranges[2 * i + 6 * nlayers] = Configuration.minstrike.at(i);
        ranges[2 * i + 6 * nlayers + 1] = Configuration.maxstrike.at(i);
      }
    cout << "Init finished ! " << endl;
    ;
  }