|
|
|
#include <GenericHouseholder.h>
Public Member Functions | |
| GenericHouseholder (bool normalise=false) | |
| std::vector< float > | iterate (const std::vector< std::vector< float > > &eventMatrix, const std::vector< float > &energyVector) |
| run the Householder Algorithm. Returns the vector of calibration coefficients. | |
| std::vector< float > | iterate (const std::vector< std::vector< float > > &eventMatrix, const std::vector< float > &energyVector, const int nIter) |
| run the Householder Algorithm several times (nIter). Returns the final vector of calibration coefficients. | |
| ~GenericHouseholder () | |
| Destructor. | |
Private Member Functions | |
| bool | decompose (const int m, const int n, std::vector< std::vector< float > > &qr, std::vector< float > &alpha, std::vector< int > &pivot) |
| void | solve (int m, int n, const std::vector< std::vector< float > > &qr, const std::vector< float > &alpha, const std::vector< int > &pivot, std::vector< float > &r, std::vector< float > &y) |
Private Attributes | |
| bool | normaliseFlag |
Generic implementation of the QR decomposition of a matrix using Householder transformation
Definition at line 15 of file GenericHouseholder.h.
| GenericHouseholder::GenericHouseholder | ( | bool | normalise = false | ) |
Default constructor CAVEAT: use normalise = true only if you know what you're doing!
Definition at line 14 of file GenericHouseholder.cc.
:normaliseFlag(normalise) { }
| GenericHouseholder::~GenericHouseholder | ( | ) |
| bool GenericHouseholder::decompose | ( | const int | m, |
| const int | n, | ||
| std::vector< std::vector< float > > & | qr, | ||
| std::vector< float > & | alpha, | ||
| std::vector< int > & | pivot | ||
| ) | [private] |
make decomposition input: m=number of events, n=number of channels, qr=event matrix output: qr = new event matrix, alpha, pivot returns a boolean value, true if decomposition worked, false if it didn't
Definition at line 207 of file GenericHouseholder.cc.
References beta, gather_cfg::cout, i, j, gen::k, m, n, mathSSE::sqrt(), and detailsBasic3DVector::y.
Referenced by iterate().
{
int i,j,jbar,k;
float beta,sigma,alphak,qrkk;
std::vector<float> y(n);
std::vector<float> sum(n);
std::cout << "Householder::decompose() started" << std::endl;
for (j=0;j<n;j++) {
// jth column sum
sum[j]=0.;
for (i=0;i<m;i++)
// std::cout << "0: qr[i][j]" << qr[i][j] << " i = " << i << " j = " << j << std::endl;
sum[j]+=qr[i][j]*qr[i][j];
pivot[j] = j;
}
for (k=0;k<n;k++) {
// kth Householder transformation
sigma = sum[k];
jbar = k;
for (j=k+1;j<n;j++) {
if (sigma < sum[j]) {
sigma = sum[j];
jbar=j;
}
}
if (jbar != k) {
// column interchange
i = pivot[k];
pivot[k]=pivot[jbar];
pivot[jbar]=i;
sum[jbar]=sum[k];
sum[k]=sigma;
for (i=0;i<m;i++) {
sigma=qr[i][k];
qr[i][k]=qr[i][jbar];
// std::cout << "A: qr[i][k]" << qr[i][k] << " i = " << i << " k = " << k << std::endl;
qr[i][jbar]=sigma;
// std::cout << "B: qr[i][jbar]" << qr[i][k] << " i = " << i << " jbar = " << jbar << std::endl;
}
} // end column interchange
sigma=0.;
for (i=k;i<m;i++){
sigma+=qr[i][k]*qr[i][k];
// std::cout << "C: qr[i][k]" << qr[i][k] << " i = " << i << " k = " << k << std::endl;
}
if (sigma == 0.) {
std::cout << "Householder::decompose() failed" << std::endl;
return false;
}
qrkk = qr[k][k];
if (qrkk < 0.)
alpha[k] = sqrt(sigma);
else
alpha[k] = sqrt(sigma) * (-1.);
alphak = alpha[k];
beta = 1/(sigma-qrkk*alphak);
qr[k][k]=qrkk-alphak;
for (j=k+1;j<n;j++) {
y[j]=0.;
for (i=k;i<m;i++)
y[j]+=qr[i][k]*qr[i][j];
y[j]*=beta;
}
for (j=k+1;j<n;j++) {
for (i=k;i<m;i++) {
qr[i][j]-=qr[i][k]*y[j];
sum[j]-=qr[k][j]*qr[k][j];
}
}
} // end of kth householder transformation
std::cout << "Householder::decompose() finished" << std::endl;
return true;
}
| std::vector< float > GenericHouseholder::iterate | ( | const std::vector< std::vector< float > > & | eventMatrix, |
| const std::vector< float > & | energyVector | ||
| ) |
run the Householder Algorithm. Returns the vector of calibration coefficients.
Definition at line 59 of file GenericHouseholder.cc.
References funct::A, alpha, gather_cfg::cout, decompose(), alignCSCRings::e, i, j, m, n, normaliseFlag, alignCSCRings::r, and solve().
{
// An implementation of the Householder in-situ calibration algorithm
// (Least squares minimisation of residual R=b-Ax, with QR decomposition of A)
// A: matrix of channel response for all calib events
// x: vector of channel calibration coefficients
// b: vector of energies
// adapted from the original code by Matt Probert 9/08/01.
std::vector<float> solution;
unsigned int m=eventMatrix.size(); // Number of events to calibrate with
unsigned int n=eventMatrix[0].size(); // Number of channel coefficients to optimize
std::cout << "Householder::runIter(): starting calibration optimization:" << std::endl;
std::cout << " Events:" << m << ", channels: " << n << std::endl;
// Sanity check
if (m != energyVector.size())
{
std::cout << "Householder::runIter(): matrix dimensions non-conformant. " << std::endl;
std::cout << " energyVector.size()=" << energyVector.size() << std::endl;
std::cout << " eventMatrix[0].size()=" << eventMatrix[0].size() << std::endl;
std::cout << " ****************** ERROR *********************" << std::endl;
return solution; // empty vector
}
// Reserve workspace
float e25p;
unsigned int i,j;
std::vector<std::vector<float> > A(eventMatrix);
std::vector<float> energies(energyVector);
float normalisation = 0.;
// Normalise if normaliseFlag is set
if (normaliseFlag)
{
std::cout << "Householder::iterate(): Normalising event data" << std::endl;
std::cout << " WARNING: assuming 5x5 filtering has already been done" << std::endl;
for (i=0; i<m; i++)
{
e25p = 0.;
for (j=0;j<n;j++){
e25p += eventMatrix[i][j]; // lorenzo -> trying to use ESetup which already performs calibs on rechits
}
e25p /= energyVector[i];
normalisation += e25p; // SUM e25p for all events
}
normalisation/=m;
std::cout << " Normalisation = " << normalisation << std::endl;
for (i=0;i<energies.size();++i)
energies[i]*=normalisation;
}
// This is where the work goes on...
// matrix decomposition
std::vector<std::vector<float> > Acopy(A);
std::vector<float> alpha(n);
std::vector<int> pivot(n);
if( !decompose(m, n, A, alpha, pivot)) {
std::cout << "Householder::runIter(): Failed: Singular condition in decomposition."
<< std::endl;
std::cout << "***************** PROBLEM in DECOMPOSITION *************************"<<std::endl;
return solution; // empty vector
}
/* DBL_EPSILON: Difference between 1.0 and the minimum float greater than 1.0 */
float etasqr = DBL_EPSILON*DBL_EPSILON;
std::cout<<"LOOK at DBL_EPSILON :"<<DBL_EPSILON<<std::endl;
std::vector<float> r(energies); // copy energies vector
std::vector<float> e(n);
// apply transformations to rhs - find solution vector
solution.assign(n,0.);
solve(m,n,A,alpha,pivot,r,solution);
// compute residual vector r
for (i=0;i<m;i++) {
r[i]=energies[i];
for (j=0;j<n;j++)
r[i]-=Acopy[i][j]*solution[j];
}
// compute first correction vector e
solve(m,n,A,alpha,pivot,r,e);
float normy0=0.;
float norme1=0.;
float norme0;
for (i=0;i<n;i++) {
normy0+=solution[i]*solution[i];
norme1+=e[i]*e[i];
}
std::cout << "Householder::runIter(): applying first correction" << std::endl;
std::cout << " normy0 = " << normy0 << std::endl;
std::cout << " norme1 = " << norme1 << std::endl;
// not attempt at obtaining the solution is made unless the norm of the first
// correction is significantly smaller than the norm of the initial solution
if (norme1>(0.0625*normy0)) {
std::cout << "Householder::runIter(): first correction is too large. Failed."
<< std::endl;
}
// improve the solution
for (i=0;i<n;i++)
solution[i]+=e[i];
std::cout << "Householder::runIter(): improving solution...." << std::endl;
// only continue iteration if the correction was significant
while (norme1>(etasqr*normy0)) {
std::cout << "Householder::runIter(): norme1 = " << norme1 << std::endl;
for (i=0;i<m;i++) {
r[i] = energies[i];
for (j=0;j<n;j++)
r[i]-=Acopy[i][j]*solution[j];
}
// compute next correction vector
solve(m,n,A,alpha,pivot,r,e);
norme0=norme1;
norme1=0.;
for (i=0;i<n;i++)
norme1+=e[i]*e[i];
// terminate iteration if the norm of the new correction failed to decrease
// significantly compared to the norm of the previous correction
if (norme1>(0.0625*norme0))
break;
// apply correction vector
for (i=0;i<n;i++)
solution[i]+=e[i];
}
return solution;
}
| std::vector< float > GenericHouseholder::iterate | ( | const std::vector< std::vector< float > > & | eventMatrix, |
| const std::vector< float > & | energyVector, | ||
| const int | nIter | ||
| ) |
run the Householder Algorithm several times (nIter). Returns the final vector of calibration coefficients.
Definition at line 25 of file GenericHouseholder.cc.
References i.
{
std::vector<float> solution;
std::vector<float> theCalibVector(energyVector.size(),1.);
std::vector<std::vector<float> > myEventMatrix(eventMatrix);
int Nevents = eventMatrix.size(); // Number of events to calibrate with
int Nchannels = eventMatrix[0].size(); // Number of channel coefficients
// Iterate the correction
for (int iter=1;iter<=nIter;iter++)
{
// make one iteration
solution = iterate(myEventMatrix, energyVector);
if (solution.empty()) return solution;
// R.O.: or throw an exception, what's the standard CMS way ?
// re-calibrate eventMatrix with solution
for (int i=0; i<Nchannels; i++)
{
for (int ievent = 0; ievent<Nevents; ievent++)
{
myEventMatrix[ievent][i] *= solution[i];
}
// save solution into theCalibVector
theCalibVector[i] *= solution[i];
}
} // end iterate the correction
return theCalibVector;
}
| void GenericHouseholder::solve | ( | int | m, |
| int | n, | ||
| const std::vector< std::vector< float > > & | qr, | ||
| const std::vector< float > & | alpha, | ||
| const std::vector< int > & | pivot, | ||
| std::vector< float > & | r, | ||
| std::vector< float > & | y | ||
| ) | [private] |
Apply transformations to rhs output: r = ?, y = solution
Definition at line 301 of file GenericHouseholder.cc.
References gather_cfg::cout, i, j, m, n, and z.
Referenced by iterate().
{
std::vector<float> z(n,0.);
float gamma;
int i,j;
std::cout << "Householder::solve() begin" << std::endl;
for (j=0;j<n;j++) {
// apply jth transformation to the right hand side
gamma=0.;
for (i=j;i<m;i++)
gamma+=qr[i][j]*r[i];
gamma/=(alpha[j]*qr[j][j]);
for (i=j;i<m;i++)
r[i]+=gamma*qr[i][j];
}
// std::cout<<"OK1:"<<std::endl;
z[n-1]=r[n-1]/alpha[n-1];
// std::cout<<"OK2:"<<std::endl;
for (i=n-2;i>=0;i--) {
z[i]= r[i];
for (j=i+1;j<n;j++)
z[i]-=qr[i][j]*z[j];
z[i]/=alpha[i];
}
// std::cout<<"OK3:"<<std::endl;
for (i=0;i<n;i++)
y[pivot[i]]=z[i];
std::cout << "Householder::solve() finished." << std::endl;
}
bool GenericHouseholder::normaliseFlag [private] |
Definition at line 42 of file GenericHouseholder.h.
Referenced by iterate().
1.7.3