loon: CorrGauss Class Reference

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

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

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

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 }

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 }

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 }

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

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 }

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 }

00308 {
00309   //fills the vector pars with gausian distributed random numbers
00310   //if correlated==true, the elements of pars will be correlated
00311   //according to the covariance matrix provided.
00312   //if correlated==false, the elements will not be correlated
00313   //the gaussians are centered around 0, user should add the nominal value of 
00314   //the parameter to shift the mean outside of this code.
00315 
00316   pars.clear();
00317   vector<double> uncorr;
00318   if(!randgen){
00319     randgen = new TRandom(randseed);
00320   }
00321   for(int i=0;i<C->GetNrows();i++){
00322     //pick uncorrelated random numbers
00323     if(correlated){
00324       uncorr.push_back(randgen->Gaus(0,1));
00325     }
00326     else{
00327       pars.push_back(randgen->Gaus(0,sqrt((*V)(i,i))));
00328     }
00329   }
00330 
00331   if(!correlated){
00332     return;
00333   }
00334 
00335   //should do this with matrix multiplication, but I don't 
00336   //know how fast root is
00337   //transform the vector to get correlated random numbers
00338   for(int j=0;j<C->GetNrows();j++){
00339     double row=0.;
00340     for(int k=0;k<C->GetNcols();k++){
00341       row+=(*C)(j,k)*uncorr[k];
00342     }
00343     //  cout<<"row "<<row<<" beam_syst "<<beam_syst[j]<<endl;
00344     pars.push_back(row);      
00345   }
00346 }

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 }

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 }

00026 { randseed = seed; }


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

CorrGauss Class Reference

Public Member Functions

Private Attributes

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation