CorrGauss Class Reference

#include <CorrGauss.h>

List of all members.

Public Member Functions

 CorrGauss ()
 CorrGauss (TMatrixD &covar)
 CorrGauss (Int_t nrows, Int_t ncols, const double *data)
 ~CorrGauss ()
void SetCovarMatrix (TMatrixD &covar)
void SetCovarMatrix (Int_t nrows, Int_t ncols, const double *data)
void DoCholesky (int method=1)
TMatrixD * GetCovarMat () const
double GetCovarElement (int rindx, int cindx) const
TMatrixD * GetDecompMat ()
double GetDecompElement (int rindx, int cindx)
void GetRandomParSet (std::vector< double > &pars, bool correlated=true)
void SetSeed (int seed)

Private Attributes

TRandom * randgen
TMatrixD * V
TMatrixD * C
int randseed

Detailed Description

Definition at line 6 of file CorrGauss.h.


Constructor & Destructor Documentation

CorrGauss::CorrGauss (  ) 

Definition at line 9 of file CorrGauss.cxx.

00009                     :
00010   randgen(0),
00011   V(0),
00012   C(0),
00013   randseed(65539)
00014 {}

CorrGauss::CorrGauss ( TMatrixD &  covar  ) 

Definition at line 16 of file CorrGauss.cxx.

References V.

00016                                    :
00017   randgen(0),
00018   C(0),
00019   randseed(65539)
00020 {
00021   V = new TMatrixD(covar);
00022 }

CorrGauss::CorrGauss ( Int_t  nrows,
Int_t  ncols,
const double *  data 
)

Definition at line 24 of file CorrGauss.cxx.

References V.

00024                                                                 :
00025   randgen(0),
00026   C(0),
00027   randseed(65539)
00028 {
00029   V = new TMatrixD(nrows,ncols,data);
00030 }

CorrGauss::~CorrGauss (  ) 

Definition at line 32 of file CorrGauss.cxx.

References randgen, and V.

00033 {
00034 
00035   if(randgen){
00036     delete randgen;
00037   }
00038   randgen=0;
00039   if(V){
00040     delete V;
00041   }
00042   V=0;
00043   if(C){
00044     delete C;
00045   }
00046   C=0;
00047 
00048 }


Member Function Documentation

void CorrGauss::DoCholesky ( int  method = 1  ) 

Definition at line 77 of file CorrGauss.cxx.

References MuELoss::e, NPAR, and V.

Referenced by GetDecompElement(), and GetDecompMat().

00078 {
00079   //finds a matrix C such that C=C*C_T
00080   //The matrix C can then be used to transform a vector of 
00081   //independent gaussian distruibuted random numbers into a
00082   //vector of correlated random numbers.
00083 
00084   //two methods have been implemented,default method==1 uses root functions
00085   //to decompose the matrix.  Note that the matrix returned by root
00086   //is the transpose of the one that we want
00087 
00088   //method==2 uses code taken from:
00089   //Monte Carolo Theory and Practice, 
00090   //F. James, Rep. Prog. Phys. Vol 43, 1980, p. 1182 
00091   //(except upper limit on sum in last eq. should be j-1, not i-1)
00092 
00093   //method==3 does both and compares the two.
00094 
00095   //root class is probably better protected against strange cases, but 
00096   //I don't know how speedy it is, method==2 is provided first to check against
00097   //previous work, but also in case root method proves to be too slow.
00098 
00099 
00100   if(C){
00101     delete C;
00102   }
00103   C=0;
00104 
00105   if(V==0){
00106     cerr<<"You have not input the covariance matrix"<<endl;
00107     cerr<<"You must provide a covariance matrix before attempting to decompose it"<<endl;
00108     return;
00109   }
00110   if(V->GetNrows()!=V->GetNcols()){
00111     cerr<<"Nrows!=NCols"<<endl;
00112     cerr<<"Covariance matrix must be square"<<endl;
00113     return;
00114   }
00115 
00116   int NPAR = V->GetNrows();
00117   if(NPAR<1){
00118     cerr<<"Nrows==0"<<endl;
00119     cerr<<"Covariance matrix must have some rows"<<endl;
00120     return;
00121   }
00122 
00123   if(method>3||method<1) method=1;
00124 
00125   if(method==1){//use root to do it
00126     TDecompChol *tdc=new TDecompChol(*V);
00127     bool worked=tdc->Decompose();
00128     if(!worked){
00129       cerr<<"Root can not Decompose this covariance matrix"<<endl;
00130       return;
00131     }
00132     TMatrixD *tcheck = new TMatrixD(tdc->GetU());
00133     C=new TMatrixD(tcheck->T());
00134     if(tcheck){
00135       delete tcheck;
00136     }
00137     tcheck=0;
00138     if(tdc){
00139       delete tdc;
00140     }
00141     tdc=0;
00142   }
00143   else if(method==2||method==3){ //use my old code
00144     C=new TMatrixD(V->GetNrows(),V->GetNcols());
00145     for(int i=0;i<NPAR;i++){
00146       double denom=(*V)(0,0);
00147       if(denom<=0){
00148         cerr<<"Can not do Cholesky since v[0][0]="<<denom<<endl;
00149         return;
00150       }
00151       (*C)(i,0)=(*V)(i,0)/sqrt(denom);
00152       for(int j=1;j<NPAR;j++){
00153         if(j<i){
00154           double s=0.;
00155           for(int k=0;k<=j-1;k++){
00156             s+=(*C)(i,k)*(*C)(j,k);
00157           }
00158           if((*C)(j,j)!=0){
00159             (*C)(i,j)=((*V)(i,j)-s)/(*C)(j,j);
00160           }
00161           else{
00162             (*C)(i,j)=0.;
00163           }
00164         }      
00165         if(i==j){
00166           double s=0.;
00167           for(int k=0;k<=i-1;k++){
00168             s+=(*C)(i,k)*(*C)(i,k);
00169           }
00170           if((*V)(i,i)-s>=0){
00171             (*C)(j,j)=sqrt((*V)(i,i)-s);
00172           }
00173           else{
00174             (*C)(j,j)=0;
00175           }
00176         }
00177         if(j>i){
00178           (*C)(i,j)=0.;
00179         }
00180       }
00181     }
00182     //double check, c*c_T=V
00183     for(int i=0;i<NPAR;i++){
00184       for(int j=0;j<NPAR;j++){
00185         double s=0.;
00186         for(int k=0;k<NPAR;k++){
00187           s+=(*C)(i,k)*(*C)(j,k);
00188         }
00189         if(fabs(s-(*V)(i,j))>0.000001){
00190           cout<<"Does not check! "<<s<<" v["<<i<<"]["<<j<<"]="<<(*V)(i,j)<<endl;
00191         }
00192       }
00193     }
00194   }
00195 
00196   if(method==3){//do both, then compare
00197     TDecompChol *tdc=new TDecompChol(*V);
00198     bool worked=tdc->Decompose();
00199     if(!worked){
00200       cerr<<"Root can not Decompose this covariance matrix"<<endl;
00201       return;
00202     }
00203     TMatrixD *tcheck = new TMatrixD(tdc->GetU());
00204     TMatrixD *check=new TMatrixD(tcheck->T());
00205     
00206     if(tdc){
00207       delete tdc;
00208     }
00209     tdc=0;    
00210     bool match=true;
00211     for(int i=0;i<C->GetNrows();i++){
00212       for(int j=0;j<C->GetNrows();j++){
00213         if(fabs((*C)(i,j)-(*check)(i,j))>1e-11){
00214           match=false;
00215           cerr<<"Matrices are not the same!"<<endl;
00216           cerr<<"C["<<i<<","<<j<<"]="<<(*C)(i,j)
00217               <<" check["<<i<<","<<j<<"]="<<(*check)(i,j)
00218               <<"Diff "<<(*C)(i,j)-(*check)(i,j)<<endl;
00219         }       
00220       }
00221     }
00222     if(match){
00223       cout<<"The methods match"<<endl;
00224     }
00225     if(tcheck){
00226       delete tcheck;
00227     }
00228     tcheck=0;
00229     if(check){
00230       delete check;
00231     }
00232     check=0;
00233 
00234   }
00235   
00236 }

double CorrGauss::GetCovarElement ( int  rindx,
int  cindx 
) const

Definition at line 249 of file CorrGauss.cxx.

References V.

00250 { 
00251   //returns the (rindx,cindx) element of the covariance matrix
00252 
00253   if(V){
00254     if(rindx>=V->GetNrows()||cindx>=V->GetNcols()){
00255       cerr<<"Wrong indices, asking for ("<<rindx<<","<<cindx<<")"<<endl;
00256       cerr<<"From a "<<V->GetNrows()<<"x"<<V->GetNcols()<<" matrix"<<endl;
00257       return -999.;
00258     }
00259     return (*V)(rindx,cindx);
00260   }
00261   cerr<<"No covariance matrix"<<endl;
00262   return -999.;
00263  
00264 
00265 }

TMatrixD * CorrGauss::GetCovarMat (  )  const

Definition at line 238 of file CorrGauss.cxx.

References V.

00239 {
00240   //returns the covariance matrix
00241 
00242   if(V){
00243     return V;
00244   }
00245   cerr<<"No covariance matrix"<<endl;
00246   return 0;
00247 }

double CorrGauss::GetDecompElement ( int  rindx,
int  cindx 
)

Definition at line 285 of file CorrGauss.cxx.

References DoCholesky().

00286 {
00287   //returns the (rindx,cindx) element of the decomposition matrix
00288   //if it doesn't exist, trys to find it using DoCholesky
00289   //if is still can't find it, it gives up.
00290 
00291 
00292   if(!C){
00293     DoCholesky();
00294   }
00295   if(C){
00296     if(rindx>=C->GetNrows()||cindx>=C->GetNcols()){
00297       cerr<<"Wrong indices, asking for ("<<rindx<<","<<cindx<<")"<<endl;
00298       cerr<<"From a "<<C->GetNrows()<<"x"<<C->GetNcols()<<" matrix"<<endl;
00299       return -999.;
00300     }
00301     return (*C)(rindx,cindx);
00302   }
00303   cerr<<"No covariance matrix"<<endl;
00304   return -999.;
00305 }

TMatrixD * CorrGauss::GetDecompMat (  ) 

Definition at line 267 of file CorrGauss.cxx.

References DoCholesky().

00268 {
00269   //returns the decomposition matrix
00270   //if it doesn't exist, trys to find it using DoCholesky
00271   //if is still can't find it, it gives up.
00272 
00273   if(C){
00274     return C;
00275   }
00276   DoCholesky();
00277   if(C){
00278     return C;
00279   }
00280   cerr<<"No decomposed matrix"<<endl;
00281   return 0;
00282 }

void CorrGauss::GetRandomParSet ( std::vector< double > &  pars,
bool  correlated = true 
)
void CorrGauss::SetCovarMatrix ( Int_t  nrows,
Int_t  ncols,
const double *  data 
)

Definition at line 64 of file CorrGauss.cxx.

References V.

00065 {
00066   //set a new covariance matrix
00067   if(V){
00068     delete V;
00069   }
00070   V = new TMatrixD(nrows,ncols,data);
00071   if(C){
00072     delete C;
00073   }
00074   C=0;
00075 }

void CorrGauss::SetCovarMatrix ( TMatrixD &  covar  ) 

Definition at line 51 of file CorrGauss.cxx.

References V.

00052 {
00053   //set a covariance matrix
00054   if(V){
00055     delete V;
00056   }
00057   V = new TMatrixD(covar);
00058   if(C){
00059     delete C;
00060   }
00061   C=0;
00062 }

void CorrGauss::SetSeed ( int  seed  )  [inline]

Definition at line 26 of file CorrGauss.h.

References randseed.

00026 { randseed = seed; }


Member Data Documentation

TMatrixD* CorrGauss::C [private]

Definition at line 32 of file CorrGauss.h.

TRandom* CorrGauss::randgen [private]

Definition at line 29 of file CorrGauss.h.

Referenced by ~CorrGauss().

int CorrGauss::randseed [private]

Definition at line 34 of file CorrGauss.h.

Referenced by SetSeed().

TMatrixD* CorrGauss::V [private]

The documentation for this class was generated from the following files:

Generated on 14 Dec 2017 for loon by  doxygen 1.6.1