source: Sophya/trunk/Cosmo/SimLSS/cmvconcherr.cc@ 3282

#include "sopnamsp.h"
#include "machdefs.h"
#include <iostream>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <unistd.h>
#include "timing.h"
#include "histerr.h"
#include "hist2err.h"
#include "fitshisterr.h"

int ObjectType(string nameh,string nameppf);
int ConcatHistoErr(string nameh,vector<string> ppfname,int wrtyp=2,bool do_cov=false);
int ConcatHisto2DErr(string nameh,vector<string> ppfname,int wrtyp=2,int ibinkt=-1,int ibinkz=-1);
void usage(void);
void usage(void)
{
cout
<<"cmvconcherr -w wtyp -n name_histoerr file1.ppf file2.ppf ..."<<endl
<<"  name_histoerr: object name"<<endl
<<"  wtyp: bit0 write ASCII, bit1 write PPF, bit2 write FITS (all=7)"<<endl
<<"  -C : compute covariance for 1D spectrum"<<endl
<<"  -Z ibinKt : compute covariance of line Pk(ibinKt,.) along Kz for 2D"<<endl
<<"  -T ibinKz : compute covariance of line Pk(.,ibinKz) along Kt for 2D"<<endl
<<endl;
}


//---------------------------------------------------------------
int main(int narg,char *arg[])
{
 // "hpkgen" "hpkgenf" "hpkrec" "hpkgen2" "hpkgenf2" "hpkrec2"
 string nameh = "";
 vector<string> ppfname;
 int Write_Type = 2;
 bool Do_Cov = false;
 int ibinKt=-1, ibinKz=-1;

 // --- Decodage arguments
 char c;
 while((c = getopt(narg,arg,"hCn:w:Z:T:")) != -1) {
  switch (c) {
  case 'n' :
    nameh = optarg;
    break;
  case 'w' :
    sscanf(optarg,"%d",&Write_Type);
    break;
  case 'C' :
    Do_Cov = true;
    break;
  case 'Z' :
    sscanf(optarg,"%d",&ibinKt);
    break;
  case 'T' :
    sscanf(optarg,"%d",&ibinKz);
    break;
  case 'h' :
  default :
    usage(); return -1;
  }
 }
 if(nameh.size()<=0 || optind>=narg) {usage(); return -1;}
 for (int i=optind;i<narg;i++) ppfname.push_back(arg[i]);
 cout<<"nameh = "<<nameh
     <<" , number of file to be read = "<<ppfname.size()
     <<" , write_type="<<Write_Type<<" , cov="<<Do_Cov<<endl;

 // --- Quel type d'objet ?
 int obtyp = ObjectType(nameh,ppfname[0]);

 // --- lecture et moyenne des HistoErr ou Histo2DErr
 int nread=0;
 if(obtyp==1) nread = ConcatHistoErr(nameh,ppfname,Write_Type,Do_Cov);
 else if(obtyp==2) nread = ConcatHisto2DErr(nameh,ppfname,Write_Type,ibinKt,ibinKz);
 else {
   cout<<"Type of object "<<nameh<<" not decoded"<<endl;
   return -1;
 }
 cout<<"Number of file read: "<<nread<<endl;
 return 0;
}

//---------------------------------------------------------------
int ObjectType(string nameh,string nameppf)
// retourne: 1 si HistoErr
//           2 si Histo2DErr,
//           0 si autre objet
//          -1 si il n'y a pas d'objet de nom "nameh"
{
 PInPersist pis(nameppf.c_str());

 bool found_tag = pis.GotoNameTag(nameh.c_str());
 if(!found_tag) return -1;

 PPersist *pps = pis.ReadObject();
 AnyDataObj *obj = pps->DataObj();
 cout<<"Object type read from input PPF file : "<<typeid(*obj).name()<< endl;

 HistoErr *herr = dynamic_cast<HistoErr *>(obj);
 if(herr) return 1;

 Histo2DErr *herr2 = dynamic_cast<Histo2DErr *>(obj);
 if(herr2) return 2;

 return 0;
}

//---------------------------------------------------------------
int ConcatHistoErr(string nameh,vector<string> ppfname,int wrtyp,bool do_cov)
{
 if(ppfname.size()<=0) return 0;

 HistoErr *herrconc = NULL;
 TVector<r_8> Tsum; TMatrix<r_8> Tsum2;
 int nherr=0, nread=0, itest=0;
 double sum=0., sum2=0., nsum=0;
 for (unsigned int ifile=0;ifile<ppfname.size();ifile++) {
   PInPersist pis(ppfname[ifile].c_str());
   HistoErr herr;
   pis.GetObject(herr,nameh);
   cout<<ifile<<" : "<<ppfname[ifile]<<"  nbin="<<herr.NBins()<<endl;

   if(ifile==0) { 
     herrconc = new HistoErr(herr.XMin(),herr.XMax(),herr.NBins());
     nherr = herr.NBins();
     itest = int(herrconc->NBins()/5.+0.5);
     if(do_cov) {Tsum.ReSize(nherr); Tsum=0.; Tsum2.ReSize(nherr,nherr); Tsum2=0.;}
   } else if(nherr!=herr.NBins()) {
     cout<<"BAD NUMBER OF BINS"<<endl;
     continue;
   }

   if(herr.NMean()!=1) cout<<"ATTENTION: file="<<ppfname[ifile]
                           <<" nmean="<<herr.NMean()<<endl;
   nread++;
   for(int i=0;i<nherr;i++) herrconc->AddBin(i,herr(i));
   sum+=herr(itest); sum2+=herr(itest)*herr(itest); nsum++;

   if(do_cov) {
     for(int i=0;i<nherr;i++) {
       Tsum(i) += herr(i);
       for(int j=0;j<nherr;j++) Tsum2(i,j) += herr(i)*herr(j);
     }
   }

 }
 herrconc->ToVariance();
 if(nsum>0) {sum/=nsum; sum2/=nsum; sum2-=sum*sum;}
 cout<<"Test bin "<<itest<<" mean="<<sum<<" sigma^2="<<sum2<<"  nsum="<<nsum<<endl;
 cout<<"         "<<" mean="<<(*herrconc)(itest)
                  <<" sigma^2="<<herrconc->Error2(itest)<<endl;

 // --- Covariance
 if(do_cov && nread>1) {
   Tsum /= nread;
   Tsum2 /= nread;
   for(int i=0;i<nherr;i++)
     for(int j=0;j<nherr;j++) Tsum2(i,j) -= Tsum(i)*Tsum(j);
 }

 // --- ecriture
 if(wrtyp&1) {   // ecriture ASCII
   string asname = "cmvconcherr.data";
   herrconc->WriteASCII(asname);
 }
 if(wrtyp&2 || do_cov) {   // ecriture PPF
   string tagobs = "cmvconcherr.ppf";
   POutPersist posobs(tagobs);
   tagobs = "herrconc"; posobs.PutObject(*herrconc,tagobs);
   if(do_cov) {
     tagobs = "mean"; posobs.PutObject(Tsum,tagobs);
     tagobs = "cov"; posobs.PutObject(Tsum2,tagobs);
   }
 }
 if(wrtyp&4) {   // ecriture FITS
   FitsInOutFile fio("!cmvconcherr.fits",FitsInOutFile::Fits_Create);
   fio << *herrconc;
 }

 delete herrconc;

 return nread;
}

//---------------------------------------------------------------
int ConcatHisto2DErr(string nameh,vector<string> ppfname,int wrtyp,int ibinkt,int ibinkz)
{
 if(ppfname.size()<=0) return 0;

 Histo2DErr *herrconc = NULL;
 TVector<r_8> Tsum_kz; TMatrix<r_8> Tsum2_kz;
 TVector<r_8> Tsum_kt; TMatrix<r_8> Tsum2_kt;
 int nherrx=0, nherry=0, nread=0, itestx=0, itesty=0;
 double sum=0., sum2=0., nsum=0;
 for (unsigned int ifile=0;ifile<ppfname.size();ifile++) {
   PInPersist pis(ppfname[ifile].c_str());
   Histo2DErr herr;
   pis.GetObject(herr,nameh);
   cout<<ifile<<" : "<<ppfname[ifile]<<"  nbin="
       <<herr.NBinX()<<","<<herr.NBinY()<<endl;

   if(ifile==0) { 
     herrconc = new Histo2DErr(herr.XMin(),herr.XMax(),herr.NBinX()
                              ,herr.YMin(),herr.YMax(),herr.NBinY());
     nherrx = herr.NBinX(); nherry = herr.NBinY();
     itestx = int(herrconc->NBinX()/5.+0.5);
     itesty = int(herrconc->NBinY()/5.+0.5);
     if(ibinkt>=0) { // matrix de covariance pour une ligne Kz
       if(ibinkt>=nherrx) ibinkt = nherrx-1;
       Tsum_kz.ReSize(nherry); Tsum_kz=0.;
       Tsum2_kz.ReSize(nherry,nherry); Tsum2_kz=0.;
     }
     if(ibinkz>=0) { // matrix de covariance pour une ligne Kt
       if(ibinkz>=nherry) ibinkz = nherry-1;
       Tsum_kt.ReSize(nherrx); Tsum_kt=0.;
       Tsum2_kt.ReSize(nherrx,nherrx); Tsum2_kt=0.;
     }
   } else if(nherrx!=herr.NBinX() || nherry!=herr.NBinY()) {
     cout<<"BAD NUMBER OF BINS"<<endl;
     continue;
   }

   if(herr.NMean()!=1) cout<<"ATTENTION: file="<<ppfname[ifile]
                           <<" nmean="<<herr.NMean()<<endl;
   nread++;
   for(int i=0;i<nherrx;i++)
     for(int j=0;j<nherry;j++) herrconc->AddBin(i,j,herr(i,j));
   sum+=herr(itestx,itesty); sum2+=herr(itestx,itesty)*herr(itestx,itesty); nsum++;

   if(ibinkt>=0) {
     for(int i=0;i<Tsum_kz.Size();i++) {
       Tsum_kz(i) += herr(ibinkt,i);
       for(int j=0;j<Tsum_kz.Size();j++) Tsum2_kz(i,j) += herr(ibinkt,i)*herr(ibinkt,j);
     }
   }

   if(ibinkz>=0) {
     for(int i=0;i<Tsum_kt.Size();i++) {
       Tsum_kt(i) += herr(i,ibinkz);
       for(int j=0;j<Tsum_kt.Size();j++) Tsum2_kt(i,j) += herr(i,ibinkz)*herr(j,ibinkz);
     }
   }

 }
 herrconc->ToVariance();
 if(nsum>0) {sum/=nsum; sum2/=nsum; sum2-=sum*sum;}
 cout<<"Test bin "<<itestx<<","<<itesty
     <<" mean="<<sum<<" sigma^2="<<sum2<<"  nsum="<<nsum<<endl;
 cout<<"         "<<" mean="<<(*herrconc)(itestx,itesty)
                  <<" sigma^2="<<herrconc->Error2(itestx,itesty)<<endl;

 // --- Covariance
 if(ibinkt>=0 && nread>1) {
   Tsum_kz /= nread;
   Tsum2_kz /= nread;
   for(int i=0;i<Tsum_kz.Size();i++)
     for(int j=0;j<Tsum_kz.Size();j++) Tsum2_kz(i,j) -= Tsum_kz(i)*Tsum_kz(j);
 }
 if(ibinkz>=0 && nread>1) {
   Tsum_kt /= nread;
   Tsum2_kt /= nread;
   for(int i=0;i<Tsum_kt.Size();i++)
     for(int j=0;j<Tsum_kt.Size();j++) Tsum2_kt(i,j) -= Tsum_kt(i)*Tsum_kt(j);
 }

 // --- ecriture
 if(wrtyp&1) {   // ecriture ASCII
   string asname = "cmvconcherr2.data";
   herrconc->WriteASCII(asname);
 }
 if(wrtyp&2) {   // ecriture PPF
   string tagobs = "cmvconcherr2.ppf";
   POutPersist posobs(tagobs);
   tagobs = "herrconc2"; posobs.PutObject(*herrconc,tagobs);
   if(ibinkt>=0) {
     char str[32];
     sprintf(str,"kzm%d",ibinkt);
     posobs.PutObject(Tsum_kz,str);
     sprintf(str,"kzc%d",ibinkt);
     posobs.PutObject(Tsum2_kz,str);
   }
   if(ibinkz>=0) {
     char str[32];
     sprintf(str,"ktm%d",ibinkz);
     posobs.PutObject(Tsum_kt,str);
     sprintf(str,"ktc%d",ibinkz);
     posobs.PutObject(Tsum2_kt,str);
   }
 }
 if(wrtyp&4) {   // ecriture FITS
   FitsInOutFile fio("!cmvconcherr2.fits",FitsInOutFile::Fits_Create);
   fio << *herrconc;
 }

 delete herrconc;

 return nread;
}

/*
####
#### HistoErr 1D
####
openppf cmvconcherr.ppf

zone 2 2
disp herrconc "hbincont err"
disp herrconc "hbinerr"
disp herrconc "hbinent"

zone
n/plot herrconc.val%log10(x) x>0 ! "connectpoints"
n/plot herrconc.sqrt(err2)%log10(x) x>0&&err2>0 ! "connectpoints same red"

n/plot herrconc.sqrt(err2)/val%log10(x) x>0&&err2>0&&val>0 ! "connectpoints"

disp mean
imag cov

del cor
c++exec \
TMatrix<r_8> cor(cov,false); cor = 0.; \
for(int i=0;i<cor.NRows();i++) { \
  for(int j=0;j<cor.NCols();j++) { \
    double v = cov(i,i)*cov(j,j); \
    if(v<=0.) continue; \
    cor(i,j) = cov(i,j)/sqrt(v); \
} } \
KeepObj(cor); \
cout<<"Matrice cor cree "<<endl;

imag cor

####
#### Histo2DErr 2D
####
openppf cmvconcherr2.ppf

zone 2 2
imag herrconc2 "hbincont"
imag herrconc2 "hbinerr"
imag herrconc2 "hbinent"

zone

set mean kzm140
set cov kzc140

set mean ktm32
set cov ktc32

n/plot $mean.val%n ! ! "connectpoints"
disp $cov
n/plot $cov.val%c r==32 ! "connectpoints"

set cor ${cov}cor
del $cor
c++exec \
TMatrix<r_8> $cor($cov,false); $cor = 0.; \
for(int i=0;i<$cor.NRows();i++) { \
  for(int j=0;j<$cor.NCols();j++) { \
    double v = $cov(i,i)*$cov(j,j); \
    if(v<=0.) continue; \
    $cor(i,j) = $cov(i,j)/sqrt(v); \
} } \
KeepObj($cor); \
cout<<"Matrice cor cree "<<endl;

disp $cor
n/plot $cor.val%c r==32 ! "connectpoints"
*/