C1DMTObjective.cpp

Go to the documentation of this file.

#include "C1DMTObjective.h"
#include "Adaptors.h"
#include <iostream>
#include <algorithm>
#include <functional>
#include <numeric>
#include <cassert>
#include <boost/numeric/ublas/io.hpp>
#include <boost/function.hpp>
#include <boost/bind.hpp>
#include <gsl/gsl_math.h>
#include "miscfunc.h"

namespace ublas = boost::numeric::ublas;
using namespace std;



C1DMTObjective::C1DMTObjective(const MTStation  &LocalMTData): 
        // we always need a reference to the data to fit and an object to calculate synthetic data
MTData(LocalMTData)
{ 
        DataFunctions.push_back(&MTTensor::GetPhixy); // if nothing will be set in SetFitParameters
        ErrorFunctions.push_back(&MTTensor::GetdPhixy); // we use these values by default
        ErrorLevels.push_back(0.01);
        FitparametersSet = false; //but we want to know whether the fit parameters have been set to issue a warning
}

C1DMTObjective::~C1DMTObjective()
{}

C1DMTObjective::C1DMTObjective(const C1DMTObjective &Old):
PlottableObjective(Old),
DataFunctions(Old.DataFunctions),
ErrorFunctions(Old.ErrorFunctions),
ErrorLevels(Old.ErrorLevels),
FitparametersSet(Old.FitparametersSet),
MTData(Old.MTData)
{
}

C1DMTObjective& C1DMTObjective::operator= (const C1DMTObjective& source)
{
        if (this == &source) return *this;
        PlottableObjective::operator=(source);
        DataFunctions.resize(source.DataFunctions.size());
        ErrorFunctions.resize(source.ErrorFunctions.size());
        ErrorLevels.resize(source.ErrorLevels.size());
        copy(source.DataFunctions.begin(),source.DataFunctions.end(),DataFunctions.begin());
        copy(source.ErrorFunctions.begin(),source.ErrorFunctions.end(),ErrorFunctions.begin());
        copy(source.ErrorLevels.begin(),source.ErrorLevels.end(),ErrorLevels.begin());
        MTData = source.MTData;  
        FitparametersSet = source.FitparametersSet;
        return *this;
}

/*! Instead of seperate access functions, we force the user to set all of them if he doesn't want defaults 
 *
 *  This functions takes three parameters
 * @param TheDataV    A vector containing pointers to functions that return the data to fit as a double when given a *MTTensor
 * @param TheErrorV    A similar vector of pointers to functions that return the corresponding error
 * @param TheErrLevel The relative errorlevel that corresponds to the data function
 * 
 * All three vectors have to be of the same size and have at least one entry to work
 */
void C1DMTObjective::SetFitParameters(const datafuncvector_t TheDataV,const datafuncvector_t TheErrorV
        ,const std::vector<double> TheErrLevel)
{
        if (TheDataV.size() != TheErrorV.size() || TheDataV.size() != TheErrLevel.size())
                throw CFatalException("Vectors of fit parameters have incompatible size !");
        if (TheDataV.empty())
                throw CFatalException("Cannot work with empty vectors !");
        DataFunctions.resize(TheDataV.size());
        ErrorFunctions.resize(TheErrorV.size());
        ErrorLevels.resize(TheErrLevel.size());
        copy(TheDataV.begin(),TheDataV.end(),DataFunctions.begin());
        copy(TheErrorV.begin(),TheErrorV.end(),ErrorFunctions.begin());
        copy(TheErrLevel.begin(),TheErrLevel.end(),ErrorLevels.begin());
        FitparametersSet = true; //signal that the parameters have been selected by the user
}

void C1DMTObjective::AppendFitParameters(const datafunc_t TheDataFunc, const datafunc_t TheErrorFunc, const double TheErrLevel)
{
        if (!FitparametersSet) // if we didn't set anything before
        {
                DataFunctions.clear();
                ErrorFunctions.clear();
                ErrorLevels.clear();
                DataFunctions.push_back(TheDataFunc); //we overwrite the defaults
                ErrorFunctions.push_back(TheErrorFunc);
                ErrorLevels.push_back(TheErrLevel);
                FitparametersSet = true;
        }
        else //if we have set something before
        {
                DataFunctions.push_back(TheDataFunc); // we append the new values
                ErrorFunctions.push_back(TheErrorFunc);
                ErrorLevels.push_back(TheErrLevel);
        }
}

/*! Calculate synthetics from the model given by member and compare with Data given
 * the fitparameters
 */
void C1DMTObjective::SafeParallel(const ttranscribed &member)
{
        
        double returnvalue = 0; // init misfit value 
        
        
        if (!FitparametersSet)
        {
                FitparametersSet = true; // we only issue this warning once
                cerr << "Warning errorlevel and Fitparameters not set, using hardcoded default ! " << endl;
                cerr << "You should change your program !" << endl;
        }
        if (GetMTSynth().GetFrequencies().empty()) //if we didn't set the frequencies before
        { 
                GetMTSynth().SetFrequencies(MTData.GetFrequencies()); //make sure the modelled frequencies match the data frequencies
        }
        CalcSynthData(member);
        
        double currvalue = 0; //some declarations for the loop below
        unsigned int datacount = 0;
        double currdata = 0.0;
        typedef vector<MTTensor>::const_iterator tensorit_t; // local convenience typedef
        typedef datafuncvector_t::iterator funcit_t;
        funcit_t dfuncit = DataFunctions.begin();
        funcit_t efuncit = ErrorFunctions.begin();
        std::vector<double>::const_iterator elevelit = ErrorLevels.begin();
        const unsigned int ndata = DataFunctions.size() * MTData.GetMTData().size();
        SetMisfit().resize(ndata); //make sure Misfit in base class can hold enough elements
        SetSynthData().resize(ndata); // and same for data
        const double absolutemin = 1e-6;
        double currabsoluteerror = 0.0;
        for (/* init done above*/ ; dfuncit != DataFunctions.end(); ++dfuncit, ++efuncit, ++ elevelit) //for all types of parameters to calculate misfit for
        {
                for (tensorit_t datait = MTData.GetMTData().begin(), synthit = GetMTSynth().GetMTData().begin(); datait < MTData.GetMTData().end();
                                ++datait,++synthit) // go through the frequencies
                        {
                                currdata = (*dfuncit)(&*datait); //calculate the current parameter using the function pointer in the vector DataFunctions
                                currvalue = currdata - (*dfuncit)(&*synthit); //calculate difference to corresponding synthetic data
                                assert(gsl_finite(currvalue));
                                currabsoluteerror = max(fabs(currdata) * (*elevelit), (*efuncit)(&*datait) ); //determine maximum of dataerror and absolute errorlevel
                                currvalue /= max(absolutemin,currabsoluteerror); //synthetic data can have 0 error, so we have to have an absolute minimum
                                currvalue = gsl_pow_int(currvalue,GetFitExponent());
                                assert(gsl_finite(currvalue));
                                SetMisfit()(datacount) = currvalue; //set the misfit vector
                                SetSynthData()(datacount) = (*dfuncit)(&*synthit); // and synthetic data vector
                                returnvalue += fabs(currvalue); // add up for xi^2 value that will be returned
                                ++datacount; 
                        }
        }
        SetRMS( pow(returnvalue/datacount,1.0/GetFitExponent())); 
}

//! All calculation has been done in SafeParallel we only return the stored result
double C1DMTObjective::PostParallel(const ttranscribed &member)
{
        return GetRMS();
}

---