LArSoft  v06_85_00
Liquid Argon Software toolkit - http://larsoft.org/
GaussianEliminationAlg.cxx
Go to the documentation of this file.
1 
11 #include "GaussianEliminationAlg.h"
12 #include <stdexcept>
13 #include <cmath>
14 #include <iostream>
15 #include <limits>
16 
18 {
19  fDistanceStepSize = step;
20  fDistanceMax = max;
21 
23  throw std::runtime_error("Error in GaussianEliminationAlg: Cannot construct with negative step or max.");
24 
26 }
27 
29 {
30 
31  fDistanceLookupTable.clear();
33 
34  //let's be extra safe, and make sure we have more sampling points than we will need
35  double x_val=0.0;
36  while(x_val <= fDistanceMax+fDistanceStepSize){
37  fDistanceLookupTable.push_back( std::exp(x_val*x_val*0.5*-1) );
38  x_val += fDistanceStepSize;
39  }
40  //do one more to be sure to push beyond...
41  fDistanceLookupTable.push_back( std::exp(x_val*x_val*0.5*-1) );
42 }
43 
45 {
46  double d_abs = std::abs(d);
47  if(d_abs > fDistanceMax)
48  return 0.0;
49 
50  size_t low_bin = std::floor(d_abs/fDistanceStepSize);
51  return fDistanceLookupTable[low_bin] -
52  (d_abs/fDistanceStepSize-(double)low_bin)*(fDistanceLookupTable[low_bin]-fDistanceLookupTable[low_bin+1]);
53 
54 }
55 
56 const std::vector<float>& util::GaussianEliminationAlg::SolveEquations(const std::vector<float>& meanVector,
57  const std::vector<float>& sigmaVector,
58  const std::vector<float>& heightVector)
59 {
60 
61  if(meanVector.size()!=sigmaVector.size() ||
62  meanVector.size()!=heightVector.size())
63  throw std::runtime_error("Error in GaussianEliminationAlg: Vector sizes not equal.");
64 
65  FillAugmentedMatrix(meanVector,sigmaVector,heightVector);
67 
68  return fSolutions;
69 
70 }
71 
72 void util::GaussianEliminationAlg::FillAugmentedMatrix(const std::vector<float>& meanVector,
73  const std::vector<float>& sigmaVector,
74  const std::vector<float>& heightVector)
75 {
76 
77  fMatrix.resize(meanVector.size());
78  for(size_t i=0; i<meanVector.size(); i++){
79 
80  fMatrix[i].resize(meanVector.size()+1);
81  for(size_t j=0; j<meanVector.size(); j++){
82  if(sigmaVector[j] < std::numeric_limits<float>::epsilon()){
83  if(i==j)
84  fMatrix[i][j] = 1.0;
85  else
86  fMatrix[i][j] = 0.0;
87  }
88  else
89  fMatrix[i][j] = GetDistance( (meanVector[i]-meanVector[j])/sigmaVector[j] );
90  }
91  fMatrix[i][meanVector.size()] = heightVector[i];
92 
93  }
94 
95 }
96 
98 {
99 
100  fSolutions.resize(fMatrix.size(),0.0);
101 
102  for(size_t i=0; i<fMatrix.size(); i++){
103 
104  for(size_t j=i+1; j<(fMatrix[i].size()-1); j++){
105  float scale_value = fMatrix[j][i] / fMatrix[i][i];
106  for(size_t k=i; k<fMatrix[i].size(); k++)
107  fMatrix[j][k] -= fMatrix[i][k]*scale_value;
108  }//end column loop
109 
110  }//end row loop
111 
112  for(int i=fMatrix.size()-1; i>=0; i--){
113  fSolutions[i] = fMatrix[i].back();
114 
115  for(size_t j=i+1; j<fMatrix.size(); j++)
116  fSolutions[i] -= fMatrix[i][j]*fSolutions[j];
117 
118  fSolutions[i] /= fMatrix[i][i];
119  }
120 
121 }
122 
124 {
125  std::cout << "GaussianEliminationAlg." << std::endl;
126 
127  std::cout << "\tLookup table (step=" << fDistanceStepSize << ", max=" << fDistanceMax << ")" << std::endl;
128  for(size_t i=0; i<fDistanceLookupTable.size(); i++)
129  std::cout << "\t\tGaussian(" << fDistanceStepSize*i << ") = " << fDistanceLookupTable[i] << std::endl;
130 
131 
132  std::cout << "\tAugmented matrix " << std::endl;
133  for(size_t i=0; i<fMatrix.size(); i++){
134  std::cout << "\t\t | ";
135  for(size_t j=0; j<fMatrix[i].size()-1; j++)
136  std::cout << fMatrix[i][j] << " ";
137  std::cout << " | " << fMatrix[i][fMatrix[i].size()-1] << " |" << std::endl;
138  }
139 
140  std::cout << "\tSolutions" << std::endl;
141  for(size_t i=0; i<fSolutions.size(); i++)
142  std::cout << "\t\t" << i << " " << fSolutions[i] << std::endl;
143 }
void FillAugmentedMatrix(const std::vector< float > &meanVector, const std::vector< float > &sigmaVector, const std::vector< float > &heightVector)
Int_t max
Definition: plot.C:27
std::vector< std::vector< double > > fMatrix
Float_t d
Definition: plot.C:237
const std::vector< float > & SolveEquations(const std::vector< float > &meanVector, const std::vector< float > &sigmaVector, const std::vector< float > &heightVector)
std::vector< double > fDistanceLookupTable