cadijoint.cpp

Go to the documentation of this file.

#include "cadijoint.h"
#include <iostream>
#include "C1DMTObjective.h"
#include "CMTStation.h"
#include "gentypes.h"
#include "C1DRecObjective.h"
#include "SeismicDataComp.h"
#include "CLevmarConf.h"
#include "MTFitSetup.h"
#include <boost/bind.hpp>

using namespace std;
using namespace gplib;

CMTStation MTData, Best;
SeismicDataComp RecData;
CLevmarConf Configuration;
C1DMTObjective MTObjective(MTData);
C1DRecObjective RecObjective(RecData, Configuration.shift, Configuration.omega,
    Configuration.sigma, Configuration.wlevel, Configuration.slowness);

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

void init(int* nd, float* ranges)
  {

    string version = "$Id: levjoint.cpp 1131 2007-05-15 22:25:06Z max $";
    cout << endl << endl;
    cout << "Program " << version << endl;
    cout << " Levenberg Marquardt to jointly invert seismic and MT data "
        << endl;
    cout << " look at levmar.conf for configuration and setup " << endl;
    cout << " output files will be of the form " << endl;
    cout << " best_lev + ending " << endl;

    Configuration.GetData("levmar.conf");
    try
      {
        MTData.GetData(Configuration.mtinputdata);
        RecData.GetData(Configuration.recinputdata, SeismicDataComp::sac);
      } catch (FatalException &e)
      {
        cerr << e.what() << endl;
      }

    SetupMTFitParameters(Configuration, MTObjective);

    RecObjective = C1DRecObjective(RecData, Configuration.shift,
        Configuration.omega, Configuration.sigma, Configuration.wlevel,
        Configuration.slowness);
    RecObjective.SetPoisson(Configuration.poisson);
    RecObjective.SetTimeWindow(Configuration.starttime, Configuration.endtime);
    RecObjective.SetErrorLevel(Configuration.recerror);

    const int nlayers = Configuration.minres.size();
    const int nparams = nlayers * 3;
    *nd = nparams;
    ttranscribed mtmember(nlayers * 2), recmember(nlayers * 2);
    copy(Configuration.startres.begin(), Configuration.startres.end(),
        mtmember.begin());
    copy(Configuration.startthick.begin(), Configuration.startthick.end(),
        mtmember.begin() + nlayers);
    copy(Configuration.startthick.begin(), Configuration.startthick.end(),
        recmember.begin());
    copy(Configuration.startsvel.begin(), Configuration.startsvel.end(),
        recmember.begin() + nlayers);
    MTObjective.CalcPerformance(mtmember);
    RecObjective.CalcPerformance(recmember);

    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.minsvel.at(i);
        ranges[2 * i + 4 * nlayers + 1] = Configuration.maxsvel.at(i);
      }
    cout << "Init finished ! " << endl;
    ;
  }