C1DMTObjective.cpp

Go to the documentation of this file.
00001 #include <iostream>
00002 #include <algorithm>
00003 #include <functional>
00004 #include <numeric>
00005 #include <cassert>
00006 #include <boost/numeric/ublas/io.hpp>
00007 #include <boost/function.hpp>
00008 #include <boost/bind.hpp>
00009 #include <gsl/gsl_math.h>
00010 #include "miscfunc.h"
00011 #include "C1DMTObjective.h"
00012 #include "Adaptors.h"
00013 #include "FatalException.h"
00014 
00015 namespace ublas = boost::numeric::ublas;
00016 using namespace std;
00017 
00018 namespace gplib
00019   {
00020     C1DMTObjective::C1DMTObjective(const MTStation &LocalMTData) :
00021       // we always need a reference to the data to fit and an object to calculate synthetic data
00022           MTData(LocalMTData)
00023       {
00024         // if nothing will be set in SetFitParameters
00025         // we use the xy phase by default
00026         DataFunctions.push_back(&MTTensor::GetPhixy);
00027         ErrorFunctions.push_back(&MTTensor::GetdPhixy);
00028         ErrorLevels.push_back(0.01);
00029         FitparametersSet = false; //but we want to know whether the fit parameters have been set to issue a warning
00030       }
00031 
00032     C1DMTObjective::~C1DMTObjective()
00033       {
00034       }
00035 
00036     C1DMTObjective::C1DMTObjective(const C1DMTObjective &Old) :
00037       PlottableObjective(Old), DataFunctions(Old.DataFunctions),
00038           ErrorFunctions(Old.ErrorFunctions), ErrorLevels(Old.ErrorLevels),
00039           FitparametersSet(Old.FitparametersSet), MTData(Old.MTData)
00040       {
00041       }
00042 
00043     C1DMTObjective& C1DMTObjective::operator=(const C1DMTObjective& source)
00044       {
00045         if (this == &source)
00046           return *this;
00047         PlottableObjective::operator=(source);
00048         DataFunctions.resize(source.DataFunctions.size());
00049         ErrorFunctions.resize(source.ErrorFunctions.size());
00050         ErrorLevels.resize(source.ErrorLevels.size());
00051         copy(source.DataFunctions.begin(), source.DataFunctions.end(),
00052             DataFunctions.begin());
00053         copy(source.ErrorFunctions.begin(), source.ErrorFunctions.end(),
00054             ErrorFunctions.begin());
00055         copy(source.ErrorLevels.begin(), source.ErrorLevels.end(),
00056             ErrorLevels.begin());
00057         MTData = source.MTData;
00058         FitparametersSet = source.FitparametersSet;
00059         return *this;
00060       }
00061 
00062     /*! Instead of seperate access functions, we force the user to set all of them if he doesn't want defaults
00063      *
00064      *  This functions takes three parameters
00065      * @param TheDataV    A vector containing pointers to functions that return the data to fit as a double when given a *MTTensor
00066      * @param TheErrorV    A similar vector of pointers to functions that return the corresponding error
00067      * @param TheErrLevel The relative errorlevel that corresponds to the data function
00068      *
00069      * All three vectors have to be of the same size and have at least one entry to work
00070      */
00071     void C1DMTObjective::SetFitParameters(const datafuncvector_t TheDataV,
00072         const datafuncvector_t TheErrorV, const std::vector<double> TheErrLevel)
00073       {
00074         //we have to make sure all vectors have the same size
00075         if (TheDataV.size() != TheErrorV.size() || TheDataV.size()
00076             != TheErrLevel.size())
00077           throw FatalException(
00078               "Vectors of fit parameters have incompatible size !");
00079         if (TheDataV.empty())
00080           throw FatalException("Cannot work with empty vectors !");
00081         //we resize and copy, so we overwrite anything that was in the vector before
00082         DataFunctions.resize(TheDataV.size());
00083         ErrorFunctions.resize(TheErrorV.size());
00084         ErrorLevels.resize(TheErrLevel.size());
00085         copy(TheDataV.begin(), TheDataV.end(), DataFunctions.begin());
00086         copy(TheErrorV.begin(), TheErrorV.end(), ErrorFunctions.begin());
00087         copy(TheErrLevel.begin(), TheErrLevel.end(), ErrorLevels.begin());
00088         FitparametersSet = true; //signal that the parameters have been selected by the user
00089       }
00090 
00091     /*! Add a new type of MT data to the vector of quantities to fit. The first call overrides
00092      * the default xy phase that is set by the constructor.
00093      * @param TheDataFunc The data to fit. A reference to function that takes a MTTensor* const and returns a double. e.g. &MTTensor::GetPhixy to fit xy phase
00094      * @param TheErrorFunc The error associated with the data. Also a reference to function that takes a MTTensor* const and returns a double. e.g. &MTTensor::GetdPhixy
00095      * @param TheErrLevel The relative error floor for this type of data
00096      */
00097     void C1DMTObjective::AppendFitParameters(const datafunc_t TheDataFunc,
00098         const datafunc_t TheErrorFunc, const double TheErrLevel)
00099       {
00100         if (!FitparametersSet) // if we didn't set anything before
00101           {
00102             DataFunctions.clear();
00103             ErrorFunctions.clear();
00104             ErrorLevels.clear();
00105             DataFunctions.push_back(TheDataFunc); //we overwrite the defaults
00106             ErrorFunctions.push_back(TheErrorFunc);
00107             ErrorLevels.push_back(TheErrLevel);
00108             FitparametersSet = true;
00109           }
00110         else //if we have set something before
00111           {
00112             DataFunctions.push_back(TheDataFunc); // we append the new values
00113             ErrorFunctions.push_back(TheErrorFunc);
00114             ErrorLevels.push_back(TheErrLevel);
00115           }
00116       }
00117 
00118     /*! Calculate synthetics from the model given by member and compare with Data given
00119      * the fitparameters
00120      */
00121     void C1DMTObjective::SafeParallel(const ttranscribed &member)
00122       {
00123 
00124         double returnvalue = 0; // init misfit value
00125 
00126         if (!FitparametersSet)
00127           {
00128             FitparametersSet = true; // we only issue this warning once
00129             cerr
00130                 << "Warning errorlevel and Fitparameters not set, using hardcoded default ! "
00131                 << endl;
00132             cerr << "You should change your program !" << endl;
00133           }
00134         //if we didn't set the frequencies before
00135         if (GetMTSynth().GetFrequencies().empty())
00136           {
00137             //make sure the modelled frequencies match the data frequencies
00138             GetMTSynth().SetFrequencies(MTData.GetFrequencies());
00139           }
00140         //calculate the synthetic data for the model
00141         CalcSynthData(member);
00142 
00143         //the number of data we fit, depends on which aspects of the tensor we fit
00144         unsigned int datacount = 0;
00145         double currdata = 0.0;
00146         typedef vector<MTTensor>::const_iterator tensorit_t; // local convenience typedef
00147         typedef datafuncvector_t::iterator funcit_t;
00148         //we have to step through all types of data to fit
00149         //we achieve this vectors of function that return the necessary quantities and their errors
00150         //for each tensor element we call those two functions
00151         funcit_t dfuncit = DataFunctions.begin();
00152         funcit_t efuncit = ErrorFunctions.begin();
00153         //we can have a different error floor for every type of data
00154         std::vector<double>::const_iterator elevelit = ErrorLevels.begin();
00155         const unsigned int ndata = DataFunctions.size()
00156             * MTData.GetMTData().size();
00157         SetMisfit().resize(ndata); //make sure Misfit in base class can hold enough elements
00158         SetSynthData().resize(ndata); // and same for data
00159 
00160         //for all types of parameters to calculate misfit for
00161         //step through the frequencies and add up the misfit
00162         for (/* init done above*/; dfuncit != DataFunctions.end(); ++dfuncit, ++efuncit, ++elevelit)
00163           {
00164             for (tensorit_t datait = MTData.GetMTData().begin(), synthit =
00165                 GetMTSynth().GetMTData().begin(); datait
00166                 < MTData.GetMTData().end(); ++datait, ++synthit) // go through the frequencies
00167               {
00168                 //calculate the current parameter using the function pointer in the vector DataFunctions
00169                 //the function calls through iterators with elements that are referenced by iterators
00170                 //are cryptic so we store them in temporary variables with better names
00171                 currdata = (*dfuncit)(&*datait);
00172                 const double currpredicted = (*dfuncit)(&*synthit);
00173                 const double currerror = (*efuncit)(&*datait);
00174                 const double errorlevel = *elevelit;
00175                 returnvalue += CalcMisfit(currdata, currpredicted, currerror,
00176                     errorlevel, datacount);
00177                 ++datacount;
00178               }
00179           }
00180         SetRMS(pow(returnvalue / datacount, 1.0 / GetFitExponent()));
00181       }
00182 
00183     //! All calculation has been done in SafeParallel we only return the stored result
00184     double C1DMTObjective::PostParallel(const ttranscribed &member)
00185       {
00186         return GetRMS();
00187       }
00188   }

Generated on Tue Nov 3 13:24:13 2009 for GPLIB++ by  doxygen 1.5.8