Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members | Related Pages

NC::SpectrumInterpolatorFancy Class Reference

Given a set of reference histograms, interpolate between them. More...

#include <NCSpectrumInterpolator.h>

Inheritance diagram for NC::SpectrumInterpolatorFancy:

NC::ISpectrumInterpolator List of all members.

Protected Member Functions

int intintpow (int a, int b) const
 Private utility function. Calculates a^b using integer arithmetic.
double doubleintpow (double a, int b) const
 Private utility function. Calculates a^b where b is an integer.
virtual TMatrixD FindCoeffs (const std::vector< std::vector< double > > &xs, const std::vector< double > &x) const
 Do the bulk of the work.

Detailed Description

Given a set of reference histograms, interpolate between them.

Consider a curve $y(x)$ fit to a set of $N$ points ${(y_n, x_n)}$ by a chi-squared statistic

\[ \chi^2(x)=\sum_n{(y(x_n)-y_n)^2\over\sigma(x, x_n)^2} \]

where the weight accorded to any point varies with $x$.

Let $y$ be a linear combination of some basis functions $L_j(x)$

\[ y(x)=\sum_ja_jL_j(x) \]

We want to solve for the $ \vec a $ that gives the best fit ie $ {{\rm d}\chi^2\over{\rm d}a_i}=0 $ for all $i$

\[ {1\over2}{{\rm d}\chi^2\over{\rm d}a_i}= \sum_n{(\sum_ka_kL_k(x_n)-y_n)L_i(x_n)\over\sigma_n^2}=0 \]

\[ \to \sum_n{\sum_ka_kL_k(x_n)L_i(x_n)\over\sigma_n^2}= \sum_n{L_i(x_n)y_n\over\sigma_n^2} \]

Naming the terms we have

\[ \sum_kM_{ik}a_k=\sum_nV_{in}y_n \]

and solving for $a$ gives

\[ a_k=(M^{-1})_{ki}V_{in}y_n \]

The final result $y(x) = \sum_k a_kL_k(x)$ can now be computed.

\[ y=(M^{-1})_{ki}V_{in}y_nL_k(x)=L^T(x)M^{-1}V\vec y \]

So the coefficients of all the $y_n$'s ($L^T(x)M^{-1}V$) can be precomputed as they are functions only of $x,x_n,y_n$.

The calculation is general and works for multidimensional $x$. Because the final combination of the $y_n$'s is linear and the coefficients can be precomputed we can operate directly over entire histograms.

The $L$ functions are chosen to be one-term polynomials, making up the full $y(x)$ polynomial. So, we need to compute the full set of possible combinations of powers to raise things to so we know what terms are possible.

\[ L_k(x) = \prod_d(x^d)^{P_{kd}} \]

Where $x$ is the input vector in D-dimensional space $\vec x = [x^1,...,x^D]$ and $P$ is the matrix of powers. above, pows in the code.

The calculation of the coefficients is achieved in ISpectrumInterpolator::FindCoeffs. $L$ is called basis, $M$ is mat, $V$ is vecPart, $\sigma^2$ is denom and $P$ is pows.

The function chosen for $\sigma$ is $\sigma(x,x_n)=|x-x_n|$

This has the feature that $1\over\sigma^2$ diverges whenever $x$ is any of the $x_n$'s

This is good in that it forces the fit function to go through that point by giving it infinite weight, but bad in that it causes divide-by-zeros. We detect when $x$ is very close to an $x_n$ and in that case construct the coefficient vector by hand so that it picks out the corresponding $y_n$.

Todo:
Try multiplying up by all the weights, so that, instead of the weight going infinite, it goes to zero on all the other terms. This should prevent divide-by-zero but we may still have problems with the matrix being non-invertible.

Definition at line 183 of file NCSpectrumInterpolator.h.

Member Function Documentation

double NC::SpectrumInterpolatorFancy::doubleintpow (  double  a,
int  b
)  const [protected]
 

Private utility function. Calculates a^b where b is an integer.

Definition at line 124 of file NCSpectrumInterpolator.cxx.

Referenced by FindCoeffs().

00125 {
00126   double ret = 1;
00127   while(b--) ret *= a;
00128   return ret;
00129 }

TMatrixD NC::SpectrumInterpolatorFancy::FindCoeffs (  const std::vector< std::vector< double > > &  xs,
const std::vector< double > &  x
)  const [protected, virtual]
 

Do the bulk of the work.

Implements NC::ISpectrumInterpolator.

Definition at line 321 of file NCSpectrumInterpolator.cxx.

References doubleintpow(), intintpow(), and SQR().

00322 {
00323   assert(xs[0].size() == x.size());
00324 
00325   // 2 to fit a general quadratic, 1 to fit a general linear function, etc.
00326   const int fitOrder = 1;
00327 
00328   // The dimensionality of the space we are interpolating in
00329   const int D = x.size();
00330 
00331   // The size of the matrices, which is the number of coefficients we
00332   // need to determine, which is the number of different products of the
00333   // components of the x coordinate each raised to a power up to fitOrder
00334   const int K = intintpow(fitOrder+1, D);
00335 
00336   assert(xs.size() >= (unsigned int)K);
00337 
00338 
00339   // What powers to raise the xs to for each row in the matrix
00340   vector<vector<int> > pows(K);
00341   for(int k = 0; k < K; ++k) pows[k] = vector<int>(xs.size());
00342 
00343   // Calculate every combination of powers
00344   // This produces a matrix like |012012012012012012012012012|
00345   // (where k is the col number, |000111222000111222000111222|
00346   // and d the row number)       |000000000111111111222222222|
00347   for(int k = 0; k < K; ++k)
00348     for(int d = 0; d < D; ++d)
00349       pows[k][d] = (k%intintpow(fitOrder+1, d+1))/intintpow(fitOrder+1, d);
00350 
00351 
00352   // Calculate the value of the different polynomial parts at x
00353   TMatrixD basis(1, K);
00354   for(int k = 0; k < K; ++k){
00355     double basisAdd = 1;
00356     for(int d = 0; d < D; ++d){
00357       basisAdd *= doubleintpow(x[d], pows[k][d]);
00358     }
00359     basis(0, k) = basisAdd;
00360   } // end for k
00361 
00362 
00363   // The number of data points we are fitting
00364   const int N = xs.size();
00365 
00366   // How close is really close?
00367   const double eps = 1e-12;
00368   // If evaluating really close to a known point, just use that point
00369   // since when  x-x_n -> 0 we get a divide-by-zero
00370   for(int n = 0; n < N; ++n){
00371     double dist = 0;
00372     for(int d = 0; d < D; ++d){
00373       dist += SQR(xs[n][d]-x[d]);
00374     }
00375     if(dist < eps){
00376       TMatrixD coeffs(1, N);
00377       coeffs(0, n) = 1;
00378       return coeffs;
00379     }
00380   } // end for n
00381 
00382   TMatrixD mat(K, K);
00383 
00384   TMatrixD vecPart(K, N);
00385 
00386   // Loop over all the data points
00387   for(int n = 0; n < N; ++n){
00388     // The weight by which this point is discounted in the fit
00389     // denom=1 corresponds to the usual unweighted fit
00390     // TODO - only taking weight in one dimension
00391     double denom = 0;
00392     for(int d = 0; d < D; ++d) denom += SQR(x[d]-xs[n][d]);
00393 
00394     // Loop over the rows in the vector
00395     for(int j = 0; j < K; ++j){
00396       double vecadd = 1/denom;
00397       for(int d = 0; d < D; ++d)
00398         vecadd *= doubleintpow(xs[n][d], pows[j][d]);
00399       vecPart(j, n) = vecadd;
00400     }
00401 
00402     // Loop over the cells of the matrix
00403     for(int j = 0; j < K; ++j){
00404       for(int i = 0; i < K; ++i){
00405         double matadd = 1/denom;
00406         for(int d = 0; d < D; ++d){
00407           matadd *= doubleintpow(xs[n][d], pows[j][d]);
00408           matadd *= doubleintpow(xs[n][d], pows[i][d]);
00409         }
00410         mat(j, i) += matadd;
00411       } // end for i
00412     } // end for j
00413   } // end for n
00414   mat.Invert();
00415 
00416   return basis*mat*vecPart;
00417 }

int NC::SpectrumInterpolatorFancy::intintpow (  int  a,
int  b
)  const [protected]
 

Private utility function. Calculates a^b using integer arithmetic.

Definition at line 116 of file NCSpectrumInterpolator.cxx.

Referenced by FindCoeffs().

00117 {
00118   int ret = 1;
00119   while(b--) ret *= a;
00120   return ret;
00121 }

The documentation for this class was generated from the following files:
Generated on Mon Feb 15 11:10:38 2010 for loon by  doxygen 1.3.9.1