# Homework Help: Help with MINUIT Fortran Code

1. Oct 1, 2016

### RJLiberator

1. The problem statement, all variables and given/known data
For my research we need to utilize a MINUIT (fortran based) code to compute error values.
We are working on an 'example' minuit library to further our understanding of the code. I am not well-versed in FORTRAN, but we are currently having a problem:

1. Understanding where parameters are defined.

2. Relevant equations

3. The attempt at a solution

Our example code (this code is within a library openly accessible online for Minuit Fortran) is here:

Code (Text):

CC Main program
CC Calls configuration MINTIO
CC Calls MINUIT main function. Calls
CC   1) MNINIT
CC   2) MNCLER
CC   5) MNIEX
CC   6) FCN (zeroes)
CC   7) FCN (initial values of parameters).
CC The first parameter is the name of the likelihood estimator (Chi2)
CC Second parameter FUTIL
PROGRAM TESTMIN
EXTERNAL FCN
CC      OPEN (UNIT=5,FILE='DSDQ.DAT',STATUS='OLD')
CC      OPEN (UNIT=6,FILE='DSDQ.OUT',STATUS='NEW',FORM='FORMATTED')
cc      CALL MINTIO(5,6,7)   ! Not needed, default values
CALL MNINIT (5,6,7)
CALL MNPARM(1,"PARAM1", 1.0, 0.01, 1.0, 10.0, OUTFLG)
CALL MNPARM(2,"PARAM2", 1.0, 0.01, 1.0, 10.0, OUTFLG)
WRITE(*,*)OUTFLG
CALL MNCOMD (FCN, "MINIMIZE", dummy, NULL)
cc      CALL MINUIT(FCN,0)   ! User routine is called FCNK0
STOP
END

cc  The User's FCN

CC Likelihood estimator function.
CC NPAR This is the number of parameters to map.
CC GIN Input
CC Value of Chi2 (real number)

SUBROUTINE FCN(NPAR,GIN,F,X,IFLAG,FUTIL)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION X(*),GIN(*)
C   this subroutine does not use FUTIL
PARAMETER (MXBIN=50)
C

C                      calculate chisquare

CHISQ = 0.
CHI2 = 1.0
CHISQ = CHISQ + CHI2

F = CHISQ
RETURN
END

Now, we are able to compile the code and library. We are able to run this code. However, as there is no real defined paramters and data input, the code doesn't spit out anything of value.

So we are working on understanding the code. We have more examples. However, this example code is in the C language and not Fortran.

Code (Text):

/* Simple optimization example. See README for details. */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#define f2cFortran

#include "cfortran.h"
#include "minuitfcn.h"
#include "minuit.h"

#define MAXLINE 256
#define MAXPARAMS 21
#define RATEFILE "n15ag.out"
#define PARAMFILEIN "par.in"
#define PARAMFILEOUT "par.out"

double ZERO=1.0E-61;
double ROI_low=0.07;
double ROI_high=0.35;
double a[MAXPARAMS];
double val[MAXPARAMS];
double err[MAXPARAMS];
double min[MAXPARAMS];
double max[MAXPARAMS];

// Number of data pairs in the input file
int NDATA;

// Number of parameters in fit
int NPARAM;

// These arrays contain the input file
double X[1000];
double Y[1000];

void readRateFile( void )
{
FILE *fp;
double x;
double y;
double dummy1;
double dummy2;

if( ( fp = fopen( RATEFILE, "r") ) == NULL )
{
printf("Could not open input file");
exit( 1 );
}
else
{
printf( "Reading %s\n", RATEFILE );

NDATA = 0;
while( !feof( fp ) )
{
//          fscanf( fp,"%lf %lf %lf %lf ", &x, &y, &dummy1, &dummy2 );
//     printf( "%e %e %e %e \n", x, y, dummy1, dummy2 );
fscanf( fp,"%lf %lf ", &x, &y );
printf( "%e %e \n", x, y );

X[NDATA] =(double) x;
Y[NDATA] =(double) y;
NDATA++;
}
fclose( fp );
}
printf( "I found %d data points\n", NDATA );
return;
}

void readParamFile( void )
{
FILE   *fp;
double p1;
double p2;
double p3;
double p4;
int    dummy;
int    i;
char   name[] = "a";

if( ( fp = fopen( PARAMFILEIN, "r") ) == NULL )
{
printf("Could not open parameter file");
exit( 1 );
}
else
{
printf( "Reading %s\n", PARAMFILEIN );

NPARAM = 0;
while( !feof( fp ) )
{
fscanf( fp,"%lf %lf %lf %lf ", &p1, &p2, &p3, &p4 );
printf( "%e %e %e %e \n", p1, p2, p3, p4 );
val[NPARAM] = p1;
err[NPARAM] = p2;
min[NPARAM] = p3;
max[NPARAM] = p4;

NPARAM++;
}
fclose( fp );
}
printf( "I found %d parameters\n", NPARAM );

for( i = 0; i < NPARAM; i++ )
MNPARM ( i+1,
name,
val[i],
err[i],
min[i],
max[i],
dummy );

return;
}

void fcn (int npar, double* grad, double* fcnval, double* xval, int iflag, void* futil)
{
double rate;
double Chi2;
double t9;
int i;

Chi2 = 0.0;

for( i = 0; i < NDATA; i++ )
{
t9 = X[i];

rate = exp( xval[0]
+ (xval[1]/t9)
+ (xval[2]/pow(t9,1./3.))
+ (xval[3]*pow(t9,1./3.))
+ (xval[4]*t9)
+ (xval[5]*pow(t9,5./3.))
+ (xval[6]*log(t9)) )
+ exp( xval[7]
+ (xval[8]/t9)
+ (xval[9]/pow(t9,1./3.))
+ (xval[10]*pow(t9,1./3.))
+ (xval[11]*t9)
+ (xval[12]*pow(t9,5./3.))
+ (xval[13]*log(t9)) )
+ exp( xval[14]
+ (xval[15]/t9)
+ (xval[16]/pow(t9,1./3.))
+ (xval[17]*pow(t9,1./3.))
+ (xval[18]*t9)
+ (xval[19]*pow(t9,5./3.))
+ (xval[20]*log(t9)) );

//     printf("rate: %e  Y[i]: %e  \n", rate, Y[i] );

if( ( t9 > ROI_low ) && ( t9 < ROI_high ) )
{
if(!(( rate < ZERO )||( t9 < 5.0E-2  )))
//          Chi2 = Chi2 + 1.0E2 * pow( log10( (Y[i]/rate) ) ,2.0);
Chi2 = Chi2 + 1.0E2 * pow( log10( (rate/Y[i]) ) ,2.0);
}
else
{
if(!(( rate < ZERO )||( t9 < 5.0E-2  )))
//      Chi2 = Chi2 + pow( log10( (Y[i]/rate) ) ,2.0);
Chi2 = Chi2 + 1.0E1 *pow( log10( (rate/Y[i]) ) ,2.0);
}

}

printf( "Chi2=%e\n", Chi2 );

*fcnval = (double)Chi2;

}

void formatR7param( double valuefort, char *dummyStr )
{
char strValue[80];
char dummyChar;
int  dummyInt;

// convert number to string
sprintf( strValue,"% 6.5e", valuefort );

//exchange decimail point and first digit
dummyChar = strValue[1];
strValue[1]= strValue[2];
strValue[2]= dummyChar;

// extract exponent
dummyStr[0]= strValue[9];
dummyStr[1]= strValue[10];
dummyStr[2]= strValue[11];
dummyStr[3]= 0;

// convert exponent from string to number
dummyInt = atoi( dummyStr );

// shift exponent to higher value by one unit
dummyInt = dummyInt + 1;

// Convert back to string
sprintf( dummyStr,"%+.2d", dummyInt );

// replace exponent including sign
dummyStr[10] = dummyStr[0];
dummyStr[11] = dummyStr[1];
dummyStr[12] = dummyStr[2];
dummyStr[13] = dummyStr[3];

// shift characters to make space for zero
dummyStr[9]  = strValue[8];
dummyStr[8]  = strValue[7];
dummyStr[7]  = strValue[6];
dummyStr[6]  = strValue[5];
dummyStr[5]  = strValue[4];
dummyStr[4]  = strValue[3];
dummyStr[3]  = strValue[2];
dummyStr[2]  = strValue[1];

// insert zero
dummyStr[1]  = 48;

// insert sign of mantissa
dummyStr[0]  = strValue[0];

return ;

}

int displayR7Params ( int argNUM, char *fitFlags )
{
FILE   *fp;
int    i;
char   name[80];
double finalval;
char   dummyStr[80];
char   strFinalVal[80];
double dummyf1;
double dummyf2;
double dummyf3;
int    dummy;

//    printf("fitFlags  %s \n", fitFlags);
// Display final values
printf("\n\n**************** Reaclib 7 parameters *************\n");
printf( "\n" );
for (i = 0; i < 4; i++)
{
MNPOUT (i+1, name, finalval, dummyf1, dummyf2, dummyf3, dummy);

formatR7param( finalval, strFinalVal );

printf( "%s", strFinalVal );
a[i] = finalval;
}
printf( "\n" );
for (i = 4; i < 7; i++)
{
MNPOUT (i+1, name, finalval, dummyf1, dummyf2, dummyf3, dummy);

formatR7param( finalval, strFinalVal );

printf( "%s", strFinalVal );
a[i] = finalval;
}
printf( "\n" );
printf( "\n" );
for (i = 7; i < 11; i++)
{
MNPOUT (i+1, name, finalval, dummyf1, dummyf2, dummyf3, dummy);

formatR7param( finalval, strFinalVal );

printf( "%s", strFinalVal );
a[i] = finalval;
}
printf( "\n" );
for (i = 11; i < 14; i++)
{
MNPOUT (i+1, name, finalval, dummyf1, dummyf2, dummyf3, dummy);

formatR7param( finalval, strFinalVal );

printf( "%s", strFinalVal );
a[i] = finalval;
}
printf( "\n" );
printf( "\n" );
for (i = 14; i < 18; i++)
{
MNPOUT (i+1, name, finalval, dummyf1, dummyf2, dummyf3, dummy);

formatR7param( finalval, strFinalVal );

printf( "%s", strFinalVal );
a[i] = finalval;
}
printf( "\n" );
for (i = 18; i < 21; i++)
{
MNPOUT (i+1, name, finalval, dummyf1, dummyf2, dummyf3, dummy);

formatR7param( finalval, strFinalVal );

printf( "%s", strFinalVal );
a[i] = finalval;
}
printf( "\n" );
printf( "\n" );

printf("***************************************************\n\n");

if( ( fp = fopen( PARAMFILEOUT, "w") ) == NULL )
{
printf("Could not create par.out");
exit( 1 );
}
else
{
printf ("argNUM: %d\n", argNUM);
if( argNUM < 3 )
{
for( i = 0; i < NPARAM; i++ )
{

fprintf( fp,
"% e    %.1f    %e    %e\n",
a[i],
err[i],
min[i],
max[i]);

}

}
else
{
for( i = 0; i < 7; i++ )
{
if( fitFlags[0] == '1' )
err[i] = 1.0;
if( fitFlags[0] == '0' )
err[i] = 0.0;

fprintf( fp,
"% e    %.1f    %e    %e\n",
a[i],
err[i],
min[i],
max[i]);

}

for( i = 7; i < 14; i++ )
{
if( fitFlags[1] == '1' )
err[i] = 1.0;
if( fitFlags[1] == '0' )
err[i] = 0.0;

fprintf( fp,
"% e    %.1f    %e    %e\n",
a[i],
err[i],
min[i],
max[i]);

}

for( i = 14; i < 21; i++ )
{
if( fitFlags[2] == '1' )
err[i] = 1.0;
if( fitFlags[2] == '0' )
err[i] = 0.0;

fprintf( fp,
"% e    %.1f    %e    %e\n",
a[i],
err[i],
min[i],
max[i]);

}
}
printf("The new parameter set was written out to par.out\n\n");
fclose( fp );
}

return;
}

void writeRateFile( double *a )
{
FILE   *fp;
double rate;
double t9;
int    i;

if( ( fp = fopen( "ratefit.log", "w") ) == NULL )
{
printf("Could not create ratefit.log");
exit( 1 );
}
else
{

for( i = 0; i < NDATA; i++ )
{
t9 = X[i];

rate = exp( a[0]
+ (a[1]/t9)
+ (a[2]/pow(t9,1./3.))
+ (a[3]*pow(t9,1./3.))
+ (a[4]*t9)
+ (a[5]*pow(t9,5./3.))
+ (a[6]*log(t9)) )
+ exp( a[7]
+ (a[8]/t9)
+ (a[9]/pow(t9,1./3.))
+ (a[10]*pow(t9,1./3.))
+ (a[11]*t9)
+ (a[12]*pow(t9,5./3.))
+ (a[13]*log(t9)) )
+ exp( a[14]
+ (a[15]/t9)
+ (a[16]/pow(t9,1./3.))
+ (a[17]*pow(t9,1./3.))
+ (a[18]*t9)
+ (a[19]*pow(t9,5./3.))
+ (a[20]*log(t9)) );

fprintf( fp,
"%e %e %e %f\n",
t9,
Y[i],
( double )rate,
100.0*fabs(1.0-(Y[i]/(( double )rate))) );

}
printf("For your convenience I wrote your data and the fit curve to ratefit.log\n\n");
fclose( fp );
}
return;
}

int main (int argc, char *argv[])
{
int i, j, dummy;
char *command = NULL;
char name[80];

// Allocate memory for input info to minuit
command = ( char * ) malloc( MAXLINE * sizeof (char) );

// Read file with rate points to fit

// Initialize minuit
MNINIT( 5, 6, 7 );

// Fill in parameter values

// Select minimization strategy
switch( atoi(argv[1]) )
{
case 1:
snprintf (command, MAXLINE, "MIGRAD");
break;
case 2:
snprintf (command, MAXLINE, "MINIMIZE");
break;
case 3:
snprintf (command, MAXLINE, "SCAN");
break;
case 4:
snprintf (command, MAXLINE, "SEEK");
break;
case 5:
snprintf (command, MAXLINE, "SIMPLEX");
break;
default:
snprintf (command, MAXLINE, "EXIT");
}

MNCOMD( minuitfcn, command, dummy, NULL );

displayR7Params( argc, argv[2] );
writeRateFile( a );
exit (0);
}

So, we need to understand what's going on, where the parameters are being defined and how the data is being interpreted/inputted.

The parameter that gave us problems is NPAR.

From the FORTRAN MINUIT manual we have:

input parameters:
NPAR = # of currently variable parameters
X = vector of parameters
IFLAG = indicates what is to be calculated
FUTIL = Name of Utilitary routine (if needed)

Output Parameters:
F = The calculated function value
(GRAD) GIN = The optional vector of first derivatives.

Thank you for any help.

2. Oct 6, 2016