/* Distributions de masse pour n=1,2,3,... tirages */
#include "sopnamsp.h"
#include "machdefs.h"
// Pour faire les histogrammes des distributions de masses
//  pour ntirages:
// ex: on a un pixel qui a N galaxies, quelle est la distribution
//     statistique de masse dans un tel pixel ?
// > cmvschdist -a -v -m 1e+7,1e+14 -n 140 -r 200 -N 100000
#include <iostream>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <unistd.h>
#include "timing.h"
#include "histos.h"
#include "ntuple.h"
#include "srandgen.h"
#include "perandom.h"

#include "geneutils.h"
#include "cosmocalc.h"
#include "schechter.h"

#include <vector>

void usage(void);
void usage(void)
{
 cout<<"cmvschdist -m xmin,xmax -a -n nbin -r ngmax,ngmin -N nalea"<<endl
     <<"  xmin,xmax : limites en masse"<<endl
     <<"  nbin : nombre de bins"<<endl
     <<"  ngmax,ngmin : distribution pour tirage de ngmin a ngmax galaxies"<<endl
     <<"  nalea : nombre de tirages aleatoires"<<endl
     <<"  -v : back verification of random trials"<<endl
     <<"  -a : init auto de l\'aleatoire"<<endl;
}

int main(int narg,char *arg[])
{
 double h75 = 0.71 / 0.75;
 double nstar = 0.006*pow(h75,3.);  //  
 double mstar = pow(10.,9.8/(h75*h75));  // MSol
 double alpha = -1.37;
 cout<<"h75= "<<h75<<" nstar= "<<nstar<<"  mstar="<<mstar<<"  alpha="<<alpha<<endl;

 double xmin=1e8, xmax=1e13;
 int npt = 100;
 int ngmax = 200, ngmin = 1;
 unsigned long long nalea = 10000;
 bool verif = false;

 char c;
 while((c = getopt(narg,arg,"havm:N:n:r:")) != -1) {
   switch (c) {
   case 'm' :
     sscanf(optarg,"%lf,%lf",&xmin,&xmax);
     if(xmin<=0.) xmin=1e6;
     if(xmax<=0.) xmax=xmin;
     if(xmin>xmax) xmax=10.*xmin;
     break;
   case 'r' :
     sscanf(optarg,"%d,%d",&ngmax,&ngmin);
     if(ngmin<=0) ngmin=1;
     if(ngmax<=0) ngmax=ngmin;
     if(ngmin>ngmax) ngmin=ngmax;
     break;
   case 'n' :
     sscanf(optarg,"%d",&npt);
     if(npt<=0) npt = 100;
     break;
   case 'N' :
     sscanf(optarg,"%llu",&nalea);
     if(nalea<=0) nalea=10000;
     break;
   case 'v' :
     verif = true;
     break;
   case 'a' :
     Auto_Ini_Ranf(5);
     break;
   case 'h' :
   default :
     usage(); return -1;
   }
 }

 double lnx1=log10(xmin), lnx2=log10(xmax), dlnx = (lnx2-lnx1)/npt;
 cout<<"xmin="<<xmin<<" ("<<lnx1<<") xmax="<<xmax<<" ("<<lnx2<<"), dlnx="<<dlnx<<endl;
 cout<<"npt="<<npt<<",  nalea="<<nalea<<endl;
 cout<<"ngmin="<<ngmin<<" ngmax="<<ngmax<<endl;
 cout<<"verif="<<verif<<endl;


 //-------m*dn/dm
 cout<<"> Schechter m*dn/dm nstar="<<nstar<<" mstar="<<mstar<<" alpha="<<alpha<<endl;
 Schechter sch(nstar,mstar,alpha);
 sch.SetOutValue(1);  // on veut m*dN/dm
 Histo hmdndm(lnx1,lnx2,npt); hmdndm.ReCenterBin();
 FuncToHisto(sch,hmdndm,true);
 FunRan tirhmdndm = FunRan(hmdndm,true);

 //------- Construct histo
 int nbhist = ngmax-ngmin+1;
 cout<<"> Creating "<<nbhist<<" histos"<<endl;
 vector<Histo> vhisto;  vector<string> vhname;
 for(int i=ngmin;i<=ngmax;i++) {
   Histo h(hmdndm); h.Zero();
   vhisto.push_back(h);
   char str[32]; sprintf(str,"h%d",i);
   vhname.push_back(string(str));
 }

 //------- Random
 InitTim();
 cout<<"> Random: "<<nalea<<" trials"<<endl;
 int lpmod = nalea/25; if(lpmod<=0) lpmod=1;
 double bigsum = 0.; long nbigsum=0;
 for(unsigned long long ia=0;ia<nalea;ia++) {
   if(ia%lpmod==0) cout<<"... "<<ia<<endl;
   double sum = 0.;
   for(int i=1;i<=ngmax;i++) {
     //double l10m = tirhmdndm.Random();
     double l10m = tirhmdndm.RandomInterp();
     double m = pow(10.,l10m);
     sum += m;
     bigsum += m; nbigsum++;
     int ipo = i-ngmin;
     if(ipo<0) continue;
     double v = log10(sum/(double)i);
     vhisto[ipo].Add(v);
   }
   if(ia%lpmod==0) PrtTim(" ");
 }
 PrtTim("End Random loop");
 if(nbigsum>0) {
   bigsum /= (double)nbigsum;
   cout<<"Mean mass : "<<bigsum<<" ("<<log10(fabs(bigsum))<<")"<<endl;
 }

 //------- Generation des classes de tirage aleatoire et des histos de verif
 vector<FunRan> vtir;  vector<string> vtname;
 vector<Histo> vthisto;  vector<string> vthname;
 if(verif) {
   cout<<"> Creating "<<nbhist<<" FunRan and Histos for back-verif"<<endl;
   for(int i=ngmin;i<=ngmax;i++) {
     FunRan t(vhisto[i-ngmin],true);
     vtir.push_back(t);
     char str[32]; sprintf(str,"t%d",i);
     vtname.push_back(string(str));
     Histo h(vhisto[i-ngmin]); h.Zero();
     vthisto.push_back(h);
     sprintf(str,"th%d",i);
     vthname.push_back(string(str));
   }

   cout<<"> Checking tirage"<<endl;
   for(unsigned long long ia=0;ia<nalea;ia++) {
     for(unsigned int i=0;i<vtir.size();i++) {
       //double v = vtir[i].Random();
       double v = vtir[i].RandomInterp();
       vthisto[i].Add(v);
     }
   }
   PrtTim("End Random loop back-tirage");
 }

 //------- Ecriture ppf
 cout<<"Ecriture"<<endl;
 string tag = "cmvschdist.ppf";
 POutPersist pos(tag);
 tag = "hmdndm"; pos.PutObject(hmdndm,tag);
 Histo hdum2(tirhmdndm);
 tag = "tirhmdndm"; pos.PutObject(hdum2,tag);
 if(vhisto.size()>0)
   for(unsigned int i=0;i<vhisto.size();i++)
     if(vhisto[i].NEntries()>0) pos.PutObject(vhisto[i],vhname[i]);
 if(vthisto.size()>0)
   for(unsigned int i=0;i<vthisto.size();i++)
     if(vthisto[i].NEntries()>0) pos.PutObject(vthisto[i],vthname[i]);

 return 0;
}

/*
delobjs *
openppf cmvschdist.ppf

set n 70
echo ${h70.vmax}

set h h$n
set th th$n

disp $h
disp $th "same red"

n/plot $h.val%x ! ! "connectpoints"
n/plot $th.val%x ! ! "nsta connectpoints same red"

n/plot $h.log10(val)%x val>0. ! "connectpoints"
n/plot $th.log10(val)%x val>0. ! "nsta connectpoints same red"

c++exec \
for(int i=0;i<$h.NBins();i++) \
  if($h(i)>0.) { \
  cout<<i<<" m="<<$h.BinCenter(i)<<" ("<<pow(10.,$h.BinCenter(i))<<") h="<<$h(i)<<endl; \
  break; \
} \
for(int i=$h.NBins()-1;i>=0;i--) \
  if($h(i)>0.) { \
  cout<<i<<" m="<<$h.BinCenter(i)<<" ("<<pow(10.,$h.BinCenter(i))<<") h="<<$h(i)<<endl; \
  break; \
}

# Compare evolution
disp h200                
disp h100 "same red"
disp h50 "same blue"
disp h10 "same green"
disp h1 "same yellow"
*/