levjoint.cpp

Go to the documentation of this file.

#include <iostream>
#include "Iso1DMTObjective.h"
#include "MTStation.h"
#include "gentypes.h"
#include "SeismicDataComp.h"
#include "AbsVelRecObjective.h"
#include "CLevmarConf.h"
#include "lm.h"
#include "MTFitSetup.h"
#include "SurfaceWaveData.h"
#include <boost/bind.hpp>
#include <boost/shared_ptr.hpp>

using namespace std;
using namespace gplib;

CLevmarConf Configuration("levmar.conf");
MTStation MTData, Best;
boost::shared_ptr<const SeismicDataComp> RecData(new SeismicDataComp(
    Configuration.recinputdata, SeismicDataComp::sac));
Iso1DMTObjective MTObjective(MTData);
SurfaceWaveData RFAbsVel;
AbsVelRecObjective RecObjective(RecData, RFAbsVel, Configuration.shift,
    Configuration.sigma, Configuration.wlevel, Configuration.slowness);

void misfit(double *p, double *x, int m, int n, void *data)
  {
    const unsigned int nlayers = m / 3;
    const unsigned int nparam = nlayers * 2;
    ttranscribed mtmember(nparam), recmember(nparam);
    for (unsigned int i = 0; i < nparam; ++i)
      {
        mtmember(i) = p[i];
      }
    for (unsigned int i = 0; i < nparam; ++i)
      recmember(i) = p[i + nlayers];
    int offset = 0;
    if (Configuration.weights.at(0) > 0)
      {
        MTObjective.CalcPerformance(mtmember);
        for (size_t i = 0; i < MTObjective.GetMisfit().size(); ++i)
          {
            x[i] = MTObjective.GetMisfit()(i) * Configuration.weights.at(0);
          }
        offset = MTObjective.GetMisfit().size();
      }
    if (Configuration.weights.at(1) > 0)
      {
        RecObjective.CalcPerformance(recmember);
        for (size_t i = 0; i < RecObjective.GetMisfit().size(); ++i)
          {
            x[i + offset] = RecObjective.GetMisfit()(i)
                * Configuration.weights.at(1);
          }
      }
  }

int main()
  {

    string version = "$Id: levjoint.cpp 1816 2009-09-07 11:28:35Z mmoorkamp $";
    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);
        RFAbsVel.ReadAscii(Configuration.rfabsveldata);
      } catch (FatalException &e)
      {
        cerr << e.what() << endl;
        return -1;
      }

    //we have to do this again becuase the global objects were initialized with the wrong inputdata
    MTObjective = Iso1DMTObjective(MTData);
    SetupMTFitParameters(Configuration, MTObjective);
    RecCalc::trfmethod rfmethod = RecCalc::specdiv;
    if (Configuration.recmethod == "iterdecon")
      {
        rfmethod = RecCalc::iterdecon;
      }
    RecObjective = AbsVelRecObjective(RecData, RFAbsVel, Configuration.shift,
        Configuration.sigma, Configuration.wlevel, Configuration.slowness,
        rfmethod);
    RecObjective.SetPoisson(Configuration.poisson);
    RecObjective.SetTimeWindow(Configuration.starttime, Configuration.endtime);
    RecObjective.SetErrorLevel(Configuration.recerror);
    RecObjective.SetRecWeight(Configuration.recweight);
    RecObjective.SetAbsVelWeight(Configuration.absvelweight);

    const int nlayers = Configuration.minres.size();
    const int nparams = nlayers * 3;
    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);

    //double p[nparams],lb[nparams],ub[nparams];
    double *p = new double[nparams];
    double *lb = new double[nparams];
    double *ub = new double[nparams];
    double *x;

    int n = 0;
    if (Configuration.weights.at(0) > 0)
      n += MTObjective.GetSynthData().size();
    if (Configuration.weights.at(1) > 0)
      n += RecObjective.GetSynthData().size();

    if (n == 0)
      {
        cerr << "One weight has to be > 0 !";
        return 100;
      }
    x = new double[n];
    double opts[LM_OPTS_SZ], info[LM_INFO_SZ];

    opts[0] = LM_INIT_MU;
    opts[1] = 1E-15;
    opts[2] = 1E-15;
    opts[3] = 1E-20;
    opts[4] = LM_DIFF_DELTA; // relevant only if the finite difference jacobian version is used

    for (int i = 0; i < nlayers; ++i)
      {
        p[i] = Configuration.startres.at(i);
        lb[i] = Configuration.minres.at(i);
        ub[i] = Configuration.maxres.at(i);
      }
    for (int i = 0; i < nlayers; ++i)
      {
        p[i + nlayers] = Configuration.startthick.at(i);
        lb[i + nlayers] = Configuration.minthick.at(i);
        ub[i + nlayers] = Configuration.maxthick.at(i);
      }
    for (int i = 0; i < nlayers; ++i)
      {
        p[i + 2 * nlayers] = Configuration.startsvel.at(i);
        lb[i + 2 * nlayers] = Configuration.minsvel.at(i);
        ub[i + 2 * nlayers] = Configuration.maxsvel.at(i);
      }

    for (int i = 0; i < n; i++)
      x[i] = 0.0;

    double ret = dlevmar_bc_dif(misfit, p, x, nparams, n, lb, ub,
        Configuration.maxiter, opts, info, NULL, NULL, NULL); // no jacobian
    cout << "Levenberg-Marquardt returned " << ret << " in " << info[5]
        << "iter, reason " << info[6] << endl;

    cout << endl << " Minimization info:" << endl;
    for (int i = 0; i < LM_INFO_SZ; ++i)
      cout << info[i] << " ";
    cout << endl;

    ttranscribed mtbest(nlayers * 2), recbest(nlayers * 2);
    for (int i = 0; i < nlayers * 2; ++i)
      {
        mtbest(i) = p[i];
        recbest(i) = p[i + nlayers];
      }
    cout << endl;
    double mtdiff = MTObjective.CalcPerformance(mtbest);
    double recfit = RecObjective.CalcPerformance(recbest);

    cout << "MT-Fit: " << mtdiff << " Rec-Fit: " << recfit << endl;
    ostringstream filename;
    filename << "best_lev";

    cout << " Saved as : " << filename.str() << endl;
    MTObjective.WriteData(filename.str());
    MTObjective.WriteModel(filename.str() + "_mt.mod");
    MTObjective.WritePlot(filename.str() + "_mt.plot");
    RecObjective.WriteData(filename.str() + ".rec");
    RecObjective.WriteModel(filename.str() + "_rec.mod");
    RecObjective.WritePlot(filename.str() + "_rec");
    delete[] x;
    delete[] p;
    delete[] lb;
    delete[] ub;
    return 0;
  }