MultiGaussianStateCombiner1D Class Reference

#include <MultiGaussianStateCombiner1D.h>

Public Member Functions
SingleGaussianState1D	combine (const MultiGaussianState1D &theState) const
SingleGaussianState1D	combine (const VSC &theComponents) const

Private Types
typedef std::vector< SingleGaussianState1D >	VSC

Detailed Description

Class to collapse (combine) a Gaussian mixture of states into one. (c.f. R. Fruewirth et.al., Comp.Phys.Comm 100 (1997) 1

Definition at line 13 of file MultiGaussianStateCombiner1D.h.

Member Typedef Documentation

typedef std::vector<SingleGaussianState1D> MultiGaussianStateCombiner1D::VSC

private

Definition at line 15 of file MultiGaussianStateCombiner1D.h.

Member Function Documentation

SingleGaussianState1D MultiGaussianStateCombiner1D::combine

(

const MultiGaussianState1D &

theState

)

const

Definition at line 7 of file MultiGaussianStateCombiner1D.cc.

References MultiGaussianState1D::components().

Referenced by MultiGaussianState1D::checkCombinedState().

                                                                                                      {
  return combine(theState.components());
}

SingleGaussianState1D MultiGaussianStateCombiner1D::combine

(

const VSC &

theComponents

)

const

Definition at line 11 of file MultiGaussianStateCombiner1D.cc.

References gather_cfg::cout, Exception, and mps_merge::weight.

                                                                                          {
  if (theComponents.empty()) {
    throw cms::Exception("LogicError") << "MultiGaussianStateCombiner1D:: state container to collapse is empty";
    //     return SingleState(SingleState::Vector(),SingleState::Matrix(),0.);
    //     return SingleState(SingleState::Vector(),
    //             SingleState::Matrix(),0.);
    return SingleGaussianState1D();
  }

  const SingleGaussianState1D firstState(theComponents.front());
  if (theComponents.size() == 1)
    return firstState;

  //   int size = firstState.mean().num_row();
  double meanMean(0.);
  double weightSum(0.);
  //   double weight;
  double measCovar1(0.);
  double measCovar2(0.);
  for (VSC::const_iterator mixtureIter1 = theComponents.begin(); mixtureIter1 != theComponents.end(); mixtureIter1++) {
    double weight = mixtureIter1->weight();
    weightSum += weight;

    double mean1 = mixtureIter1->mean();
    meanMean += weight * mean1;
    measCovar1 += weight * mixtureIter1->variance();

    for (VSC::const_iterator mixtureIter2 = mixtureIter1 + 1; mixtureIter2 != theComponents.end(); mixtureIter2++) {
      double posDiff = mean1 - mixtureIter2->mean();
      //       SingleState::Matrix s(1,1); //stupid trick to make CLHEP work decently
      //       measCovar2 +=weight * (*mixtureIter2).weight() *
      //                    s.similarity(posDiff.T().T());
      //       SingleState::Matrix mat;
      //       for ( unsigned int i1=0; i1<N; i1++ ) {
      //    for ( unsigned int i2=0; i2<=i1; i2++ )  mat(i1,i2) = posDiff(i1)*posDiff(i2);
      //       }
      //       measCovar2 += weight * (*mixtureIter2).weight() * mat;
      //
      // TensorProd yields a general matrix - need to convert to a symm. matrix
      double covGen = posDiff * posDiff;
      //       double covSym(covGen.LowerBlock());
      measCovar2 += weight * mixtureIter2->weight() * covGen;
    }
  }

  double measCovar;
  if (weightSum < DBL_MIN) {
    std::cout << "MultiGaussianStateCombiner1D:: New state has total weight of 0." << std::endl;
    //     meanMean = SingleState::Vector(size,0);
    meanMean = 0.;
    measCovar = 0.;
    weightSum = 0.;
  } else {
    weightSum = 1. / weightSum;
    meanMean *= weightSum;
    measCovar1 *= weightSum;
    measCovar2 *= weightSum * weightSum;
    measCovar = measCovar1 + measCovar2;
    //     measCovar = measCovar1/weightSum + measCovar2/weightSum/weightSum;
  }

  return SingleGaussianState1D(meanMean, measCovar, weightSum);
}