GblData.cc

Go to the documentation of this file.

/*
 * GblData.cpp
 *
 *  Created on: Aug 18, 2011
 *      Author: kleinwrt
 */

#include "Alignment/ReferenceTrajectories/interface/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),
    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) {

    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 aDiff = theValue - thePrediction;
    return aDiff * aDiff * thePrecision * theDownWeight;
  }

  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];
    }
  }
}