levaniso.cpp

Go to the documentation of this file.

#include <iostream>
#include "Aniso1DMTObjective.h"
#include "MTStation.h"
#include "gentypes.h"
#include "CLevanisoConf.h"
#include "lm.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(double *p, double *x, int m, int n, void *data)
  {
    const unsigned int nlayers = m / 4;
    const unsigned int nparam = m;
    ttranscribed mtmember(nparam);
    for (int i = 0; i < nparam; ++i)
      {
        mtmember(i) = p[i];
      }

    if (Configuration.weights.at(0) > 0)
      {
        MTObjective->CalcPerformance(mtmember);
        const unsigned int nelements = MTObjective->GetMisfit().size();
        for (size_t i = 0; i < nelements; ++i)
          {
            x[i] = MTObjective->GetMisfit()(i) * Configuration.weights.at(0);
          }
      }
  }

int main()
  {

    string version = "$Id: levaniso.cpp 1816 2009-09-07 11:28:35Z mmoorkamp $";
    cout << endl << endl;
    cout << "Program " << version << endl;
    cout << " Levenberg Marquardt to  MT data for 1D anisotropic structure"
        << endl;
    cout << " look at levaniso.conf for configuration and setup " << endl;
    cout << " output files will be of the form " << endl;
    cout << " best_lev + ending " << endl;

    //CMTStation MTData,Best;

    Configuration.GetData("levaniso.conf");
    try
      {
        //MTData.GetData(Configuration.mtinputdata);


        //setup MT Objective function


        if (Configuration.mtfit == "ptensor")
          {
            PTensorMTStation PTData;
            PTData.GetData(Configuration.ptensordata);
            boost::shared_ptr<PTensor1DMTObjective> PTObjective(
                new PTensor1DMTObjective(PTData));
            MTObjective = PTObjective;
          }
        else
          {
            MTStation MTData;
            MTData.GetData(Configuration.mtinputdata);
            boost::shared_ptr<Aniso1DMTObjective> AnisoObjective(
                new Aniso1DMTObjective(MTData));
            SetupMTFitParameters(Configuration, *AnisoObjective);
            MTObjective = AnisoObjective;
          }
        MTObjective->SetFitExponent(2);

        const int nlayers = Configuration.minres.size();
        const int nparams = nlayers * 4;
        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);
        double startfit = MTObjective->CalcPerformance(mtmember);

        cout << "Start misfit: " << startfit << endl;

        //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 (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.startaniso.at(i);
            lb[i + 2 * nlayers] = Configuration.minaniso.at(i);
            ub[i + 2 * nlayers] = Configuration.maxaniso.at(i);
          }
        for (int i = 0; i < nlayers; ++i)
          {
            p[i + 3 * nlayers] = Configuration.startstrike.at(i);
            lb[i + 3 * nlayers] = Configuration.minstrike.at(i);
            ub[i + 3 * nlayers] = Configuration.maxstrike.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(nparams);
        for (int i = 0; i < nparams; ++i)
          {
            mtbest(i) = p[i];
          }
        cout << endl;
        double mtdiff = MTObjective->CalcPerformance(mtbest);

        cout << "MT-Fit: " << mtdiff << 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");

        delete[] x;
        delete[] p;
        delete[] lb;
        delete[] ub;
      } catch (FatalException &e)
      {
        cerr << e.what() << endl;
        return -1;
      }
    return 0;
  }