source: Sophya/trunk/SophyaExt/MinuitAdapt/minuitadapt.cc@ 2411

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

#include "pexceptions.h"
#include "minuitadapt.h"
#include "localmnpout.h"

#define AVEC_CFORTRAN

#ifdef AVEC_CFORTRAN
#include "Cern/cfortran.h"
#include "Cern/minuit.h"
#else
#include "minuitprot.h"
#endif

/*!
  \class SOPHYA::Minuit Cern Fitter adapter
  \ingroup MinuitAdapt
*/

MinuitAdapt::MinuitAdapt(
  void (*myfcn)(int_4 *,double *,double *,double *,int_4 *,double (*f)(double *)),
  double (*myfutils)(double *x)
  )
  : mNPar(0), iErFlg(0), fcn(myfcn), futils(myfutils)
{
 int_4 i1=5,i2=6,i3=0;
#ifdef AVEC_CFORTRAN
 MNINIT(i1,i2,i3);
#else
 mninit_(&i1,&i2,&i3);
#endif
}

MinuitAdapt::~MinuitAdapt(void)
{
}

void MinuitAdapt::SetTitle(char* title)
{
 int_4 nc=strlen(title);
 if(nc>0) {
#ifdef AVEC_CFORTRAN
   MNSETI(title);
#else
   mnseti_(title,nc);
#endif
 }
}

void MinuitAdapt::Clear(void)
{
 char* str="CLEAR";
 double arglis[1]={0.}; int_4 narg=0;
#ifdef AVEC_CFORTRAN
 MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
#else
 int_4 nc=strlen(str);
 mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
#endif
 mNPar = 0;
}

void MinuitAdapt::Return(void)
{
 char* str="RETURN";
 double arglis[1]={0.}; int_4 narg=0;
#ifdef AVEC_CFORTRAN
 MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
#else
 int_4 nc=strlen(str);
 mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
#endif
}

void MinuitAdapt::PrintLevel(int_4 n)
{
 char* str="SET PRINTOUT";
 double arglis[1]={(double)n}; int_4 narg=1;
#ifdef AVEC_CFORTRAN
 MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
#else
 int_4 nc=strlen(str);
 mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
#endif
}

void MinuitAdapt::SetWidthPage(int_4 n)
{
 char* str="SET WIDTHPAGE";
 double arglis[1]={(double)n}; int_4 narg=1;
#ifdef AVEC_CFORTRAN
 MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
#else
 int_4 nc=strlen(str);
 mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
#endif
}

void MinuitAdapt::SetWarnings(bool w)
{
 char* str="SET WARNINGS";
 if(!w) str="SET NOWARNINGS";
 double arglis[1]={0.}; int_4 narg=0;
#ifdef AVEC_CFORTRAN
 MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
#else
 int_4 nc=strlen(str);
 mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
#endif
}

void MinuitAdapt::SetErrorDef(double err)
// 1=xi2, 0.5=likelyhood
{
 char* str="SET ERRORDEF";
 double arglis[1]={err}; int_4 narg=1;
#ifdef AVEC_CFORTRAN
 MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
#else
 int_4 nc=strlen(str);
 mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
#endif
}

void MinuitAdapt::SetEpsMachine(double eps)
{
 char* str="SET EPSMACHINE";
 double arglis[1]={eps}; int_4 narg=1;
#ifdef AVEC_CFORTRAN
 MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
#else
 int_4 nc=strlen(str);
 mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
#endif
}

void MinuitAdapt::SetStrategy(int_4 n)
{
 char* str="SET STRATEGY";
 double arglis[1]={(double)n}; int_4 narg=1;
#ifdef AVEC_CFORTRAN
 MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
#else
 int_4 nc=strlen(str);
 mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
#endif
}

void MinuitAdapt::SetGradient(int_4 n)
// n<0 : gradient from finite differences computed by Minuit
// n=0 : gradient from user AND finite diff.. then compare
// n>0 : force gradient from user
{
 char* str="SET NOGRADIENT";
 double arglis[1]={1.}; int_4 narg=0;
 if(n==0) str="SET GRADIENT";
 else if(n>0) {str="SET GRADIENT"; narg=1;}
#ifdef AVEC_CFORTRAN
 MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
#else
 int_4 nc=strlen(str);
 mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
#endif
}

void MinuitAdapt::SetRandom(int_4 seed)
// Seed must be integer between 10000 and 900000000
{
 char* str="SET RANDOMGENERATOR";
 double arglis[1]={(double)seed}; int_4 narg=1;
#ifdef AVEC_CFORTRAN
 MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
#else
 int_4 nc=strlen(str);
 mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
#endif
}

void MinuitAdapt::DefineParameter(int_4 n,char* name,double init,double step
                                 ,double vmin,double vmax)
// nom parametre, numero parametre, initial value, initial step size, parameter bounds
// - ATTENTION numero de parametre [1,...]
{
 char str1[16]; sprintf(str1,"P%d",n);
 char* str=str1;
 if(name!=NULL) str=name;
#ifdef AVEC_CFORTRAN
 MNPARM(n,str,init,step,vmin,vmax,iErFlg);
#else
 int_4 nc=strlen(str);
 mnparm_(&n,str,&init,&step,&vmin,&vmax,&iErFlg,nc);
#endif
 mNPar++;
}

void MinuitAdapt::SetParameter(int_4 n,double val)
{
 char* str="SET PARAMETER";
 double arglis[2]={(double)n,val}; int_4 narg=2;
#ifdef AVEC_CFORTRAN
 MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
#else
 int_4 nc=strlen(str);
 mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
#endif
}

void MinuitAdapt::SetLimits(int_4 n,double val1,double val2)
{
 char* str="SET LIMITS";
 double arglis[3]={(double)n,val1,val2}; int_4 narg=3;
#ifdef AVEC_CFORTRAN
 MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
#else
 int_4 nc=strlen(str);
 mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
#endif
}

void MinuitAdapt::SetFix(int_4 n)
{
 char* str="FIX";
 double arglis[1]={(double)n}; int_4 narg=1;
#ifdef AVEC_CFORTRAN
 MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
#else
 int_4 nc=strlen(str);
 mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
#endif
}

void MinuitAdapt::Release(int_4 n)
{
 char* str="RELEASE";
 double arglis[1]={(double)n}; int_4 narg=1;
#ifdef AVEC_CFORTRAN
 MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
#else
 int_4 nc=strlen(str);
 mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
#endif
}

void MinuitAdapt::Restore(bool only_last_fixed)
{
 char* str="RESTORE";
 double arglis[1]={0.}; int_4 narg=0;
 if(only_last_fixed) {arglis[0]=1.; narg=1;}
#ifdef AVEC_CFORTRAN
 MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
#else
 int_4 nc=strlen(str);
 mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
#endif
}

void MinuitAdapt::DoFit(char *method,int_4 maxcall,double tol)
{
 char* str=method;
 double arglis[2]={(double)maxcall,tol}; int_4 narg=2;
#ifdef AVEC_CFORTRAN
 MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
#else
 int_4 nc=strlen(str);
 mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
#endif
}

void MinuitAdapt::Improve(int_4 maxcall)
{
 char* str="IMPROVE";
 double arglis[1]={(double)maxcall}; int_4 narg=1;
#ifdef AVEC_CFORTRAN
 MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
#else
 int_4 nc=strlen(str);
 mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
#endif
}

void MinuitAdapt::Hesse(int_4 maxcall)
{
 char* str="HESSE";
 double arglis[1]={(double)maxcall}; int_4 narg=1;
#ifdef AVEC_CFORTRAN
 MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
#else
 int_4 nc=strlen(str);
 mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
#endif
}

void MinuitAdapt::Minos(int_4 maxcall)
{
 char* str="MINOS";
 double arglis[1]={(double)maxcall}; int_4 narg=1;
#ifdef AVEC_CFORTRAN
 MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
#else
 int_4 nc=strlen(str);
 mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
#endif
}

void MinuitAdapt::Seek(int_4 maxcall,double stdev)
{
 char* str="SEEK";
 double arglis[2]={(double)maxcall,stdev}; int_4 narg=2;
#ifdef AVEC_CFORTRAN
 MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
#else
 int_4 nc=strlen(str);
 mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
#endif
}

void MinuitAdapt::Call(int_4 iflag)
// See meaning in fcn exemple
{
 char* str="CALL";
 double arglis[1]={(double)iflag}; int_4 narg=1;
#ifdef AVEC_CFORTRAN
 MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
#else
 int_4 nc=strlen(str);
 mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
#endif
}

void MinuitAdapt::DrawContour(int_4 n1,int_4 n2,int_4 npts)
{
 char* str="MNCONTOUR";
 double arglis[3]={(double)n1,(double)n2,(double)npts}; int_4 narg=3;
#ifdef AVEC_CFORTRAN
 MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
#else
 int_4 nc=strlen(str);
 mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
#endif
}

int_4 MinuitAdapt::GetContour(int_4 n1,int_4 n2,int_4 npts
                             ,TVector<r_8>& xpt,TVector<r_8>& ypt)
{
 if(npts<5) npts=5;
 int_4 ncontok = 0;
 double* xcont = new double[npts];
 double* ycont = new double[npts];
#ifdef AVEC_CFORTRAN
 MNCONT(fcn,n1,n2,npts,xcont[0],ycont[0],ncontok,futils);
#else
 mncont_(fcn,&n1,&n2,&npts,xcont,ycont,&ncontok,futils);
#endif
 if(ncontok<1) return 0;
 xpt.ReSize(ncontok); ypt.ReSize(ncontok);
 for(int_4 i=0;i<ncontok;i++) {xpt(i)=xcont[i]; ypt(i)=ycont[i];}
 delete [] xcont; delete [] ycont;
 return ncontok;
}

void MinuitAdapt::Scan(int_4 n,double from,double to,int_4 npts)
{
 char* str="SCAN";
 double arglis[4]={(double)n,(double)npts,from,to}; int_4 narg=4;
#ifdef AVEC_CFORTRAN
 MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
#else
 int_4 nc=strlen(str);
 mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
#endif
}

void MinuitAdapt::ShowParameter(void)
{
 char* str="SHOW PARAMETER";
 double arglis[1]={0.}; int_4 narg=0;
#ifdef AVEC_CFORTRAN
 MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
#else
 int_4 nc=strlen(str);
 mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
#endif
}

void MinuitAdapt::ShowRandom(void)
{
 char* str="SHOW RANDOM";
 double arglis[1]={0.}; int_4 narg=0;
#ifdef AVEC_CFORTRAN
 MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
#else
 int_4 nc=strlen(str);
 mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
#endif
}

void MinuitAdapt::ShowFcnValue(void)
{
 char* str="SHOW FCNVALUE";
 double arglis[1]={0.}; int_4 narg=0;
#ifdef AVEC_CFORTRAN
 MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
#else
 int_4 nc=strlen(str);
 mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
#endif
}

void MinuitAdapt::ShowCovariance(void)
{
 char* str="SHOW COVARIANCE";
 double arglis[1]={0.}; int_4 narg=0;
#ifdef AVEC_CFORTRAN
 MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
#else
 int_4 nc=strlen(str);
 mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
#endif
}

void MinuitAdapt::ShowCorrelations(void)
{
 char* str="SHOW CORRELATIONS";
 double arglis[1]={0.}; int_4 narg=0;
#ifdef AVEC_CFORTRAN
 MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
#else
 int_4 nc=strlen(str);
 mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
#endif
}

void MinuitAdapt::ShowEigenValues(void)
{
 char* str="SHOW EIGENVALUES";
 double arglis[1]={0.}; int_4 narg=0;
#ifdef AVEC_CFORTRAN
 MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
#else
 int_4 nc=strlen(str);
 mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
#endif
}

void MinuitAdapt::GetParameter(int_4 n,string& name,double& val,double& err
                              ,double& bound1,double& bound2,int_4& ivarbl)
{
 double p,ep,b1,b2; int_4 i; char str[256]="";
 localmnpout_(&n,str,&p,&ep,&b1,&b2,&i);
 val=p; err=ep; bound1=b1; bound2=b2; ivarbl=i;
 // Protection contre les retours intempestifs du Fortran !
 if(strlen(str)>0)
   for(unsigned int i=0;i<strlen(str);i++)
     if(str[i]==' ') {str[i]='\0'; break;}
 name = str;
}

void MinuitAdapt::GetErrors(int_4 n,double& eplus,double& eminus
                           ,double& eparab,double& globcc)
{
 double ep,em,epa,g;
#ifdef AVEC_CFORTRAN
 MNERRS(n,ep,em,epa,g);
#else
 mnerrs_(&n,&ep,&em,&epa,&g);
#endif
 eplus=ep; eminus=em; eparab=epa; globcc=g;
}

TMatrix<r_8> MinuitAdapt::GetErrorsMatrix(void)
{
 TMatrix<r_8> mat;
 int_4 npar=NPar();
 if(npar<=0) return mat;
 double* mcov = new double[npar*npar];
#ifdef AVEC_CFORTRAN
 MNEMAT(mcov[0],npar);
#else
 mnemat_(mcov,&npar);
#endif
 mat.ReSize(npar,npar);
 for(int i=0;i<npar;i++)
  for(int j=0;j<npar;j++) mat(i,j) = mcov[i*npar+j];
 delete [] mcov;
 return mat;
}

////////////////////////////////////////////////////////////
/*
void fcn(int_4 *npar,double *grad,double *fval,double *xval
        ,int_4 *iflag,double futils(double *))
{
  //cout<<"iflag="<<*iflag<<endl;
  IFLAG[0]++;
  if(*iflag>0 && *iflag<10) IFLAG[*iflag]++;

  // Read input,init,... data values
  //  if(*iflag==1) {...}

  // Instruct Minuit to redefine the problem
  // and forget about previously best fitted values.
  //  if(*iflag==5) {...}

  // Always compute Chi2 or Likelyhood (here iflag==4)
 *fval=0.;
 for(int i=0;i<DIMX;i++) for(int j=0;j<DIMY;j++) {
   double f = Z[i][j]-Gauss2D(X[i],Y[j],xval);
   *fval += f*f/(EZ[i][j]*EZ[i][j]);
 }

 // Compute (optionnal) the first derivative of Chi2 / parameters
 if(*iflag==2) {
   // Return gradient of chi2 (if SET GRA called)
   // C'est DChi2/DPi = -2*sum{(Yi-F(Xi))/EYi^2 * dF/dPi(Xi)}
   double dpar[NPAR];
   for(int j=0;j<NPAR;j++) grad[j]=0.;
   for(int i=0;i<DIMX;i++) for(int j=0;j<DIMY;j++) {
     double f=-2.*(Z[i][j]-Gauss2D(X[i],Y[j],xval))/(EZ[i][j]*EZ[i][j]);
     dGauss2D(X[i],Y[j],xval,dpar);
     for(int k=0;k<NPAR;k++) grad[k]+= f*dpar[k];
   }
 }

 // Called at the end of the fit (on the Miniut RETURN)
 if(*iflag==3) {
   cout<<"Call fcn iflag="<<*iflag<<endl;
   for(int k=0;k<NPAR;k++) cout<<" P"<<k+1<<"="<<xval[k];
   cout<<endl;
   cout<<"Number of fcn calls="<<IFLAG[0]<<endl;
   for(int k=1;k<10;k++)
     cout<<"  iflag="<<k<<" number of calls="<<IFLAG[k]<<endl;
 }
}
*/