C1DMTObjective.cpp

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

Generated on Fri Jul 4 15:30:20 2008 for GPLIB++ by  doxygen 1.5.5