1. Not finding help here? Sign up for a free 30min tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Help with MINUIT Fortran Code

  1. Oct 1, 2016 #1

    RJLiberator

    User Avatar
    Gold Member

    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   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)
    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
        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.
     
  2. jcsd
  3. Oct 6, 2016 #2
    Thanks for the thread! This is an automated courtesy bump. Sorry you aren't generating responses at the moment. Do you have any further information, come to any new conclusions or is it possible to reword the post? The more details the better.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted



Similar Discussions: Help with MINUIT Fortran Code
  1. Help with fortran (Replies: 1)

  2. Fortran Help! (Replies: 4)

  3. Fortran help (Replies: 16)

  4. Fortran code (Replies: 1)

  5. FORTRAN code question (Replies: 5)

Loading...