7 #include "CLevmarConf.h"
11 #include <boost/bind.hpp>
12 #include <boost/shared_ptr.hpp>
15 using namespace gplib;
26 void misfit(
double *p,
double *x,
int m,
int n,
void *data)
28 const unsigned int nlayers = m / 3;
29 const unsigned int nparam = nlayers * 2;
31 for (
unsigned int i = 0; i < nparam; ++i)
35 for (
unsigned int i = 0; i < nparam; ++i)
36 recmember(i) = p[i + nlayers];
61 string version =
"$Id: levjoint.cpp 1816 2009-09-07 11:28:35Z mmoorkamp $";
63 cout <<
"Program " << version << endl;
64 cout <<
" Levenberg Marquardt to jointly invert seismic and MT data "
66 cout <<
" look at levmar.conf for configuration and setup " << endl;
67 cout <<
" output files will be of the form " << endl;
68 cout <<
" best_lev + ending " << endl;
78 cerr << e.what() << endl;
88 rfmethod = RecCalc::iterdecon;
100 const int nparams = nlayers * 3;
101 ttranscribed mtmember(nlayers * 2), recmember(nlayers * 2);
105 mtmember.begin() + nlayers);
109 recmember.begin() + nlayers);
114 double *p =
new double[nparams];
115 double *lb =
new double[nparams];
116 double *ub =
new double[nparams];
127 cerr <<
"One weight has to be > 0 !";
131 double opts[LM_OPTS_SZ], info[LM_INFO_SZ];
133 opts[0] = LM_INIT_MU;
137 opts[4] = LM_DIFF_DELTA;
139 for (
int i = 0; i < nlayers; ++i)
145 for (
int i = 0; i < nlayers; ++i)
151 for (
int i = 0; i < nlayers; ++i)
158 for (
int i = 0; i < n; i++)
161 double ret = dlevmar_bc_dif(
misfit, p, x, nparams, n, lb, ub,
163 cout <<
"Levenberg-Marquardt returned " << ret <<
" in " << info[5]
164 <<
"iter, reason " << info[6] << endl;
166 cout << endl <<
" Minimization info:" << endl;
167 for (
int i = 0; i < LM_INFO_SZ; ++i)
168 cout << info[i] <<
" ";
172 for (
int i = 0; i < nlayers * 2; ++i)
175 recbest(i) = p[i + nlayers];
181 cout <<
"MT-Fit: " << mtdiff <<
" Rec-Fit: " << recfit << endl;
182 ostringstream filename;
183 filename <<
"best_lev";
185 cout <<
" Saved as : " << filename.str() << endl;
void SetRecWeight(const double w)
Set the relative weight for the pure receiver function.
trfmethod
There are several ways to calculate receiver functions.
void SetupMTFitParameters(const tConfObject &Configuration, C1DMTObjective &Objective)
ublas::vector< double > ttranscribed
double CalcPerformance(const ttranscribed &member)
For serial execution CalcPerformance calls the three Parallel functions for more convenient use...
void misfit(double *p, double *x, int m, int n, void *data)
virtual void WriteData(const std::string &filename)
Write current data to a file.
void WritePlot(const std::string &filename)
Write the current model to ascii file for plotting.
virtual void WriteData(const std::string &filename)
Write out the receiver function and absolute velocity data, the absolute velocity data gets ...
boost::shared_ptr< const SeismicDataComp > RecData(new SeismicDataComp(Configuration.recinputdata, SeismicDataComp::sac))
virtual void WritePlot(const std::string &filename)
write the current model for plotting to a file
virtual void GetData(const std::string filename)
read in data from file, determines format by ending
A class to read, write and store fundamental mode surface wave dispersion data.
CLevmarConf Configuration("levmar.conf")
const tdata & GetSynthData()
Return the current synthetic data.
void SetErrorLevel(const double level)
Set the errorlevel for fit, this is relative to the maximum amplitude, not for each individual data p...
The class MTStation is used to store the transfer functions and related information for a MT-site...
void SetAbsVelWeight(const double w)
Set the relative weight for the absolute velocity information.
virtual void ReadAscii(const std::string &filename)
Read a file in general ascii format, i.e lines with period velocity each.
const tmisfit & GetMisfit()
Return the misfit vector.
void WriteModel(const std::string &filename)
Write the current model to ascii file for calculations.
This objective function calculates the weighted misfit for a receiver function and the corresponding ...
void SetPoisson(const double ratio)
Set poisson's ratio, at the moment the same for all layers, used for calculating P-velocity.
AbsVelRecObjective RecObjective(RecData, RFAbsVel, Configuration.shift, Configuration.sigma, Configuration.wlevel, Configuration.slowness)
This class implements the forward calculation for a 1D isotropic MT model.
void SetTimeWindow(const double start, const double end)
Set the time window used for misfit calculations, start and end are in seconds.
virtual void WriteModel(const std::string &filename)
write the current model to a file
The basic exception class for all errors that arise in gplib.
Iso1DMTObjective MTObjective(MTData)