GblData.cpp

Go to the documentation of this file.

/*
 * GblData.cpp
 *
 *  Created on: Aug 18, 2011
 *      Author: kleinwrt
 */
 
#include "GblData.h"
using namespace Eigen;
 
namespace gbl {
 
 
GblData::GblData(unsigned int aLabel, dataBlockType aType, double aValue,
                double aPrec, unsigned int aTraj, unsigned int aPoint) :
                theLabel(aLabel), theRow(0), theType(aType), theValue(aValue), thePrecision(
                                aPrec), theTrajectory(aTraj), thePoint(aPoint), theDWMethod(0), theDownWeight(
                                1.), thePrediction(0.), theNumLocal(0), moreParameters(), moreDerivatives() {
 
}
 
GblData::~GblData() {
}
 
 
void GblData::addDerivatives(const std::vector<unsigned int> &index,
                const std::vector<double> &derivatives) {
        for (unsigned int i = 0; i < derivatives.size(); ++i) // any derivatives
                        {
                if (derivatives[i]) {
                        moreParameters.push_back(index[i]);
                        moreDerivatives.push_back(derivatives[i]);
                }
        }
}
 
void GblData::setPrediction(const VVector &aVector) {
 
        thePrediction = 0.;
        if (theNumLocal > 0) {
                for (unsigned int i = 0; i < theNumLocal; ++i) {
                        thePrediction += theDerivatives[i] * aVector(theParameters[i] - 1);
                }
        } else {
                for (unsigned int i = 0; i < moreDerivatives.size(); ++i) {
                        thePrediction += moreDerivatives[i]
                                        * aVector(moreParameters[i] - 1);
                }
        }
}
 
 
double GblData::setDownWeighting(unsigned int aMethod) {
 
        theDWMethod = aMethod;
        double aWeight = 1.;
        double scaledResidual = fabs(theValue - thePrediction) * sqrt(thePrecision);
        if (aMethod == 1) // Tukey
                        {
                if (scaledResidual < 4.6851) {
                        aWeight = (1.0 - 0.045558 * scaledResidual * scaledResidual);
                        aWeight *= aWeight;
                } else {
                        aWeight = 0.;
                }
        } else if (aMethod == 2) //Huber
                        {
                if (scaledResidual >= 1.345) {
                        aWeight = 1.345 / scaledResidual;
                }
        } else if (aMethod == 3) //Cauchy
                        {
                aWeight = 1.0 / (1.0 + (scaledResidual * scaledResidual / 5.6877));
        }
        theDownWeight = aWeight;
        return aWeight;
}
 
 
double GblData::getChi2() const {
        double scaledResidual = fabs(theValue - thePrediction) * sqrt(thePrecision);
        double chi2 = scaledResidual * scaledResidual;
        if (theDWMethod == 1) // Tukey
                        {
                if (scaledResidual < 4.6851) {
                        chi2 = (1.0
                                        - pow(1.0 - 0.045558 * scaledResidual * scaledResidual, 3))
                                        / (3. * 0.045558);
                } else {
                        chi2 = 1.0 / (3. * 0.045558);
                }
        } else if (theDWMethod == 2) //Huber
                        {
                if (scaledResidual >= 1.345) {
                        chi2 = 1.345 * (2. * scaledResidual - 1.345);
                }
        } else if (theDWMethod == 3) //Cauchy
                        {
                chi2 = log(1.0 + (scaledResidual * scaledResidual / 5.6877)) * 5.6877;
        }
        return chi2;
}
 
void GblData::printData() const {
 
        std::cout << " measurement at label " << theLabel << " of type " << theType
                        << " from row " << theRow << ": " << theValue << ", "
                        << thePrecision << std::endl;
        std::cout << "  param " << moreParameters.size() + theNumLocal << ":";
        for (unsigned int i = 0; i < moreParameters.size(); ++i) {
                std::cout << " " << moreParameters[i];
        }
        for (unsigned int i = 0; i < theNumLocal; ++i) {
                std::cout << " " << theParameters[i];
        }
        std::cout << std::endl;
        std::cout << "  deriv " << moreDerivatives.size() + theNumLocal << ":";
        for (unsigned int i = 0; i < moreDerivatives.size(); ++i) {
                std::cout << " " << moreDerivatives[i];
        }
        for (unsigned int i = 0; i < theNumLocal; ++i) {
                std::cout << " " << theDerivatives[i];
        }
        std::cout << std::endl;
}
 
 
unsigned int GblData::getLabel() const {
        return theLabel;
}
 
 
dataBlockType GblData::getType() const {
        return theType;
}
 
 
void GblData::getLocalData(double &aValue, double &aWeight,
                unsigned int &numLocal, unsigned int* &indLocal, double* &derLocal) {
 
        aValue = theValue;
        aWeight = thePrecision * theDownWeight;
        if (theNumLocal > 0) {
                numLocal = theNumLocal;
                indLocal = theParameters;
                derLocal = theDerivatives;
        } else {
                numLocal = moreParameters.size();
                indLocal = &moreParameters[0];
                derLocal = &moreDerivatives[0];
        }
}
 
 
void GblData::getAllData(double &aValue, double &aErr, unsigned int &numLocal,
                unsigned int* &indLocal, double* &derLocal, unsigned int &aTraj,
                unsigned int &aPoint, unsigned int &aRow) {
        aValue = theValue;
        aErr = 1.0 / sqrt(thePrecision);
        if (theNumLocal > 0) {
                numLocal = theNumLocal;
                indLocal = theParameters;
                derLocal = theDerivatives;
        } else {
                numLocal = moreParameters.size();
                indLocal = &moreParameters[0];
                derLocal = &moreDerivatives[0];
        }
        aTraj = theTrajectory;
        aPoint = thePoint;
        aRow = theRow;
}
 
 
void GblData::getResidual(double &aResidual, double &aVariance,
                double &aDownWeight, unsigned int &numLocal, unsigned int* &indLocal,
                double* &derLocal) {
        aResidual = theValue - thePrediction;
        aVariance = 1.0 / thePrecision;
        aDownWeight = theDownWeight;
        if (theNumLocal > 0) {
                numLocal = theNumLocal;
                indLocal = theParameters;
                derLocal = theDerivatives;
        } else {
                numLocal = moreParameters.size();
                indLocal = &moreParameters[0];
                derLocal = &moreDerivatives[0];
        }
}
}