source: Sophya/trunk/Cosmo/SimLSS/cmvfitpk.cc@ 3891

#include "sopnamsp.h"
#include "machdefs.h"
#include <iostream>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <unistd.h>

#include "ntuple.h"
#include "histerr.h"

#include "constcosmo.h"
#include "pkspectrum.h"

#include "Minuit2/FunctionMinimum.h"
#include "Minuit2/MnMigrad.h"
#include "Minuit2/MnUserParameters.h"
#include "Minuit2/MnPrint.h"
#include "Minuit2/MnMinos.h"
#include "Minuit2/MinosError.h"
#include "Minuit2/FCNBase.h"
using namespace ROOT::Minuit2;

namespace ROOT { namespace Minuit2 {
class MyFCN : public FCNBase {
public:
  MyFCN(HistoErr& herr,PkEisenstein& pk);
  void Init(const vector<double>& par) const;
  double Func(double x,const vector<double>& par) const;
  double operator()(const vector<double>& par) const;
  double Up(void) const {return 1.;}
private:
  HistoErr& herr_;
  PkEisenstein& pk_;
};
} }  // namespace ROOT + Minuit2

void usage(void);
void usage(void)
{
cout<<"usage: cmvfitpk [-k kmin,kmax] -e conchpkrec.ppf observ3d_1.ppf observ3d_2.ppf ..."<<endl;
}

int main(int narg,char *arg[])
{
 // --- Cosmography definition (WMAP)
 double zref = 0.175;
 double ob0 = 0.0444356;
 double h100=0.71, om0=0.267804, ol0=0.73;

 // --- Spectrum and variance definition
 double ns = 1., as = 1.;
 double R=8./h100, sigmaR = 1.;
 double kmin=1e-5,kmax=1000.;
 int npt = 10000;
 double eps=1.e-3;

 // --- les fichiers
 double klim[2]={0.,-1.};
 string fconc = "";

 // --- les options de print/plot
 int nprint = 3, nplot = -1;  // -1 = all

 // --- Reading arguments
 char c;
 while((c = getopt(narg,arg,"hk:e:")) != -1) {
  switch (c) {
  case 'k' :
    sscanf(optarg,"%lf,%lf",&klim[0],&klim[1]);
    break;
  case 'e' :
    fconc = optarg;
    break;
  case 'h' :
  default :
    usage(); return -1;
  }
 }

 if(optind>=narg) return -2;

 // --- Create spectrum
 cout<<endl<<"\n--- Create Spectrum"<<endl;
 InitialPowerLaw pkini(ns,as);
 bool nobaryon = false;
 double ocdm0 = om0-ob0;
 TransfertEisenstein tf(h100,ocdm0,ob0,T_CMB_Par,nobaryon);
 GrowthEisenstein growth(om0,ol0); cout<<"...Growth="<<growth(zref)<<endl;
 PkEisenstein pkz(pkini,tf,growth,zref);

 // --- Compute variance and normalize spectrum
 pkz.SetZ(0.);
 cout<<endl<<"\n--- Compute variance for top-hat R="<<R
     <<" at z="<<pkz.GetZ()<<endl;
 VarianceSpectrum varpk_th(pkz,R,VarianceSpectrum::TOPHAT);
 double kfind_th = varpk_th.FindMaximum(kmin,kmax,eps);
 double pkmax_th = varpk_th(kfind_th);
 cout<<"kfind_th = "<<kfind_th<<" ("<<log10(kfind_th)<<"), integrand="<<pkmax_th<<endl;
 double k1=kmin, k2=kmax;
 int rc = varpk_th.FindLimits(pkmax_th/1.e4,k1,k2,eps);
 cout<<"limit_th: rc="<<rc<<" : "<<k1<<" ("<<log10(k1)<<") , "
     <<k2<<" ("<<log10(k2)<<")"<<endl;

 double ldlk = (log10(k2)-log10(k1))/npt;
 varpk_th.SetInteg(0.01,ldlk,-1.,4);
 double sr2 = varpk_th.Variance(k1,k2);
 cout<<"varpk_th="<<sr2<<"  ->  sigma="<<sqrt(sr2)<<endl;

 double normpkz = sigmaR*sigmaR/sr2;
 pkz.SetScale(normpkz);
 cout<<"Spectrum normalisation = "<<pkz.GetScale()<<endl;
 pkz.SetZ(zref);

 // --- Reading reference mean spectrum
 cout<<endl<<"\n--- Reading spectrum form file "<<fconc<<endl;
 HistoErr herrconc;
 {
 PInPersist pis(fconc.c_str());
 pis >> PPFNameTag("herrconc") >> herrconc;
 }
 cout<<herrconc;
 if(herrconc.NBins()<1) return -3;
 double b_ini = herrconc(2*herrconc.NBins()/3);
 cout<<"b_ini = "<<b_ini<<endl;
 if(klim[1]<=klim[0]) klim[1] = herrconc.XMax();

 // --- NTuple et ppf
 const int npar = 7;  // Number of parameters in the fit
 TMatrix<r_8> Cov(npar,npar); Cov = 0.;
 POutPersist pos("cmvfitpk.ppf");
 const int nvar = 2*npar+1;
 const char *vname[nvar] = {
       "xi2"
       ,"b","a","oc","ob","ol","h","n"
       ,"eb","ea","eoc","eob","eol","eh","en"
       };
 NTuple nt(nvar,vname);
 double xnt[nvar];

 // ***
 // --- Boucle sur les fichiers a fitter
 // ***
 int nfile=0, ncov=0;
 for(int ifile=optind;ifile<narg;ifile++) {

   // --- lecture du fichier
   cout<<">>> Reading["<<nfile<<","<<ifile<<"]: "<<arg[ifile]<<endl;
   HistoErr herr;
   PInPersist pis(arg[ifile]);
   pis >> PPFNameTag("hpkrec") >> herr;
   if(herr.NBins()!=herrconc.NBins()) {
     cout<<"... bad number of bins: "<<herr.NBins()<<endl;
     continue;
   }

   // --- Mise en forme
   unsigned int nbinok = 0;
   for(int i=0;i<herr.NBins();i++) {
     double k = herr.BinCenter(i);
     if(i==0 || k<klim[0] || k>klim[1]) {herr.SetErr2(i,-1.); continue;}
     herr.SetErr2(i,herrconc.Error2(i));
     nbinok++;
   }

   // --- Initialisation de minuit
    MnUserParameters upar;
    upar.Add("B",b_ini,b_ini/100.);
    upar.Add("A",1.,0.01);
    upar.Add("Ocdm",ocdm0,0.001,ocdm0/2.,2.*ocdm0);
    //upar.Fix("Ocdm");
    upar.Add("Ob",ob0,0.001,ob0/10.,10.*ob0);
    //upar.Fix("Ob");
    upar.Add("Ol",ol0,0.001,0.1,2.);
    upar.Fix("Ol");
    upar.Add("h100",h100,0.01,10.,150.);
    upar.Fix("h100");
    upar.Add("ns",ns,0.001,0.7,1.5);
    upar.Fix("ns");
    MyFCN fcn(herr,pkz);
    if(nbinok<upar.VariableParameters()) {
      cout<<"... nombre de degres de liberte <0: nbinok="<<nbinok<<endl;
      continue;
    }

   // --- Fit minuit
   unsigned int nfcall = 1000; double tolerance = 0.05;
   MnMigrad migrad(fcn,upar,MnStrategy(1));
   FunctionMinimum min = migrad(nfcall,tolerance);
   if(!min.IsValid()) {
     cout<<"!!!!!!!!! Echec Minuit"<<endl;
     //cout<<"minimum: "<<min<<endl;
     continue;
   }

   MnUserParameterState ustate = min.UserState();
   double xi2r = ustate.Fval()/(nbinok-ustate.VariableParameters());
   vector<double> parfit  = ustate.Parameters().Params();
   vector<double> eparfit = ustate.Parameters().Errors();
   if(ustate.HasCovariance()) {
     MnUserCovariance ucov = ustate.Covariance();
     for(unsigned int i=0;i<upar.VariableParameters();i++)
       for(unsigned int j=0;j<upar.VariableParameters();j++)
         Cov(ustate.ExtOfInt(i),ustate.ExtOfInt(j)) += ucov(i,j);
     ncov++;
   }

   // --- Recuperation des informations et remplissage NTuple
   for(int i=0;i<nvar;i++) xnt[i] = 0.;
   xnt[0] = xi2r;
   for(int i=0;i<npar;i++) xnt[1+i] = parfit[i];
   for(int i=0;i<npar;i++) xnt[npar+1+i] = eparfit[i];
   nt.Fill(xnt);

   // --- stoquage des histos de plots
   if(nfile<nplot || nplot<0) {
     HistoErr herrf(herr); herrf.Zero();
     HistoErr herrd(herr); herrd.Zero();
     fcn.Init(parfit);
     for(int i=1;i<herr.NBins();i++) {
       double f = fcn.Func(herr.BinCenter(i),parfit);
       herrf.SetBin(i,f,herr.Error(i),1.);
       herrd.SetBin(i,herr(i)-f,herr.Error(i),1.);
     }
     char str[64];
     sprintf(str,"h%d",nfile); pos.PutObject(herr,str);
     sprintf(str,"hf%d",nfile); pos.PutObject(herrf,str);
     sprintf(str,"hd%d",nfile); pos.PutObject(herrd,str);
   }

   // --- Print de debug
   if(nfile<nprint) {
     cout<<"minimum: "<<min<<endl;
     for(int i=0;i<nvar;i++) cout<<"["<<i<<"] = "<<xnt[i]<<endl;

     MnMinos Minos(fcn,min);
     pair<double,double> eminos[npar];
     for(int ip=0;ip<npar;ip++) {
       eminos[ip].first = eminos[ip].second = 0.;
       if(!ustate.Parameter(ip).IsFixed()) eminos[ip] = Minos(ip);
       cout<<"Minos: "<<ip<<" err=["<<eminos[ip].first<<","<<eminos[ip].second<<"]"<<endl;
     }
   }

   nfile++;
 }

 // --- Fin du job
 if(nfile>0) pos.PutObject(nt,"nt");
 if(ncov>0) {Cov /= (double)ncov; pos.PutObject(Cov,"cov");}

 return 0;
}

//--------------------------------------------------
MyFCN::MyFCN(HistoErr& herr,PkEisenstein& pk)
  : herr_(herr) , pk_(pk)
{
}

void MyFCN::Init(const vector<double>& par) const
{
 double A=par[1], Ocdm0=par[2], Ob0=par[3], Ol0=par[4], h100=par[5], ns=par[6];
 pk_.GetPkIni().SetSlope(ns);
 pk_.GetPkIni().SetNorm(A);
 pk_.GetTransfert().SetParTo(h100,Ocdm0,Ob0);
 pk_.GetGrowthFactor().SetParTo(Ocdm0+Ob0,Ol0);
}

double MyFCN::Func(double x,const vector<double>& par) const
// par: [0]=B,  [1]=A,  [2]=Ocdm0,  [3]=Ob0,  [4]=Ol0, [5]=h100, [6]=npow
{
 Init(par);
 return pk_(x) + par[0];
}

double MyFCN::operator()(const vector<double>& par) const
{
 Init(par);
 double xi2 = 0.;
 for(int i=0;i<herr_.NBins();i++) {
   double e2 = herr_.Error2(i);
   if(e2<=0.) continue;
   double x = herr_.BinCenter(i);
   double v = herr_(i);
   double f = Func(x,par);
   xi2 += (v-f)*(v-f)/e2;
 }
 return xi2;
}

/*
openppf cmvfitpk.ppf

print nt

set i 0
zone
disp h$i "nsta hbincont err"
disp hf$i "nsta same red"

zone
disp hd$i "nsta hbincont err red"
disp hd$i "nsta same"

zone
n/plot nt.xi2

zone
n/plot nt.a
n/plot nt.b
n/plot nt.b%a ! ! "crossmarker5"

n/plot nt.oc
n/plot nt.ob
n/plot nt.oc+ob
n/plot nt.ob%oc ! ! "crossmarker5"

n/plot nt.ea
n/plot nt.eb
n/plot nt.ea%eb ! ! "crossmarker5"

disp cov "zoomx30"
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 "zoomx30"

 */