How Does MINUIT Fortran Define and Utilize Parameters in Optimization?

  • Comp Sci
  • Thread starter RJLiberator
  • Start date
  • Tags
    Code Fortran
In summary, We need to utilize a MINUIT (fortran based) code to compute error values for our research. We are working on an 'example' minuit library to further our understanding of the code. We are currently having a problem understanding where parameters are defined. The example code is within a library accessible online for Minuit Fortran and is written in the C language, not Fortran. The code can be compiled and run, but as there are no defined parameters or data input, it does not produce any useful output. We are working on understanding the code and have more examples available.
  • #1
RJLiberator
Gold Member
1,095
63

Homework Statement


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.

Homework Equations

The Attempt at a Solution



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

Code: 
CC Main program
CC Calls configuration MINTIO
CC Calls MINUIT main function. Calls
CC   1) MNINIT
CC   2) MNCLER
CC   3) MNREAD (reads title)
CC   4) MNREAD (reads parameters)
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)
CC                      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: 
/* 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
    readRateFile( );

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

    // Fill in parameter values
    readParamFile( );

    // 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.
 

1. What is MINUIT Fortran Code?

MINUIT Fortran Code is a computer program used for minimizing a function of several variables. It is commonly used in scientific research, particularly in high energy physics, to find the best fit of data to a theoretical model.

2. How do I install MINUIT Fortran Code?

Installing MINUIT Fortran Code requires a Fortran compiler and an operating system that supports Fortran. Detailed instructions for installation can be found on the MINUIT website, and may vary depending on your specific system.

3. Can MINUIT Fortran Code be used for any type of data analysis?

MINUIT Fortran Code is specifically designed for minimizing a function of several variables. While it is commonly used in high energy physics, it can also be applied to other scientific fields where a function needs to be minimized.

4. How do I use MINUIT Fortran Code to fit my data?

To use MINUIT Fortran Code for data fitting, you will need to write a Fortran program that defines the function you want to minimize and the data you want to fit it to. You will also need to specify the starting values for the fit parameters. Detailed instructions for using MINUIT can be found in the user manual.

5. Are there any alternatives to MINUIT Fortran Code?

Yes, there are several alternatives to MINUIT Fortran Code, such as the Levenberg-Marquardt algorithm and the Nelder-Mead algorithm. These alternatives may have different strengths and weaknesses, so it is important to choose the one that best fits your specific data analysis needs.

Hot Threads

Recent Insights

Back
Top