/*               Minuit example with many local minima.
                         September 24, 2002
             Joel Heinrich  --- University of Pennsylvania

The function to be minimized is:

    f(x,y) =  sin^2(pi*x) + sin^2(pi*y) + (x-sqrt(2)*y)^2

which has a unique global minimum at (0,0), but many local minima.
Designed to be compiled with standard C or C++; it calls the Minuit
fortran routines directly.  The user must supply starting values
on the command line.

It has been compiled and executed on (assuming the appropriate setup):

Linux:
cc minexam.c $CERNLIB/libpacklib.a /usr/lib/libg2c.a -lm -o minexam
g++ minexam.cc $CERNLIB/libpacklib.a /usr/lib/libg2c.a -lm -o minexam
KCC minexam.cc $CERNLIB/libpacklib.a /usr/lib/libg2c.a -lm -o minexam

cdfsga:
cc minexam.c $CRNLIB/libpacklib.a /usr/lib/libftn.so -lm -o minexam
CC minexam.cc $CRNLIB/libpacklib.a /usr/lib/libftn.so -lm -o minexam

fcdfsgi2
cc minexam.c $CRNLIB/libpacklib.a /usr/lib32/libftn.so -lm -o minexam
CC minexam.cc $CRNLIB/libpacklib.a /usr/lib32/libftn.so -lm -o minexam
KCC minexam.cc $CRNLIB/libpacklib.a /usr/lib32/libftn.so -lm -o minexam

SUN Solaris
cc -c minexam.c; f77 minexam.o /cern/2000/lib/libpacklib.a -o minexam
c++ -c minexam.cc; f77 minexam.o /cern/2000/lib/libpacklib.a -o minexam

*/


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

typedef void FCN(const int* npar, double grad[], double* fval,
		 const double xval[],const int* iflag, void*);

#if defined(__cplusplus)
extern "C" {
#endif

/*
      To call fortran directly from C or C++ (Unix conventions)

   1. Append an _ to the fortran-routine name.
   2. All fortran-visible arguments are passed through a pointer
      (or, alternatively, by reference from C++).
   3. Each character string argument requires an additional "invisible"
      argument (integer passed by value), containing the number of characters,
      appended to the end of the argument list.
   4. From C++ (only), must declare the fortran routines extern "C".
*/
  
void mninit_(const int* ird,const int* iwr,const int* isav);
void mnparm_(const int* nun, const char chnam[], const double* stval,
	     const double* step, const double* bnd1, const double* bnd2,
	     int* ierflg, int nchnam);
void mnseti_(const char ctitle[], int nctitle);
void mnexcm_(FCN fcn, const char chcom[], const double arglis[],
	     const int* narg, int* ierflg, void*, int nchcom);
void mnpout_(const int* num, char chnam[], double* val, double* error,
	     double* bnd1, double* bnd2, int* ivarbl, int nchnam);
int intrac_(void* dummy);

#if defined(__cplusplus)
}
#endif



FCN testfun;

int main(int argc, char* argv[]){
  const char *parlist[] = { "first", "second" };
  const int nparms = sizeof(parlist)/sizeof(parlist[0]);

  if(argc < nparms+1) {
    fprintf(stderr,"Must supply %d starting "
	    "values as command line arguments\n",nparms);
    return 1;
  }

  {
    const int ird=5, iwr=6, isav=7;
    mninit_(&ird, &iwr, &isav);
  }

  {
    int num=1, ierflg=0;
    const double step=0.1, bnd1=0.0, bnd2=0.0;

    for(;num<=nparms;++num) {
      const double stval = strtod(argv[num],NULL);
      mnparm_(&num,parlist[num-1],&stval,&step,&bnd1,&bnd2,&ierflg,
	      strlen(parlist[num-1]));
      if(ierflg) {
	fprintf(stderr, "Error on call to mnparm for parameter %s\n",
		parlist[num-1]);
	return ierflg;
      }
    }
  }

  {
    const char title[] = "Minuit example with many local minima";
     mnseti_(title,strlen(title));
  }

  {
    const char chcom[]="SET GRAD";
    const int narg=0;
    int ierflg=0;
    mnexcm_(testfun, chcom, NULL, &narg, &ierflg, NULL, strlen(chcom));
    if(ierflg) {
      fprintf(stderr, "Error on call to mnexcm for command %s\n",chcom);
      return ierflg;
    }
  }  

  {
    const char chcom[]="SET STRAT";
    const int narg=1;
    int ierflg=0;
    const double level=2;
    mnexcm_(testfun, chcom, &level, &narg, &ierflg, NULL, strlen(chcom));
    if(ierflg) {
      fprintf(stderr, "Error on call to mnexcm for command %s\n",chcom);
      return ierflg;
    }
  }

  {
    const char chcom[]="MIGRAD";
    const int narg=0;
    int ierflg=0;
    mnexcm_(testfun, chcom, NULL, &narg, &ierflg, NULL, strlen(chcom));
    if(ierflg) {
      fprintf(stderr, "Error on call to mnexcm for command %s\n",chcom);
      return ierflg;
    }
  }

  {
    int num=1,ivarbl;
    char chnam[11];
    double val,error,bnd1,bnd2;
    printf("\n=====================================================\n");
    for(;num<=nparms;++num) {
      memset(chnam,0,sizeof(chnam));
      mnpout_(&num,chnam,&val,&error,&bnd1,&bnd2,&ivarbl,sizeof(chnam)-1);
      printf("Parameter %s = %f +/- %f\n",chnam,val,error);
    }
  }

  return 0;
}

void testfun(const int* npar, double grad[], double* fval,
	     const double xval[],const int* iflag, void* ignore) {
  /*       sin^2(pi*x) + sin^2(pi*y) + (x-sqrt(2)*y)^2         */
  const double d_pi   = 3.14159265358979323846;
  const double d_sqrt2  = 1.41421356237309504880;
  const double x = xval[0], y = xval[1];
  const double fax = fmod(fabs(x),1.0), fay = fmod(fabs(y),1.0);
  const double sx = sin(fax*d_pi), sy = sin(fay*d_pi);
  const double d = y*d_sqrt2 - x;
  *fval = sx*sx + sy*sy + d*d;
  if ( *iflag == 2 ) {
    grad[0] = d_pi*sin(fax*2.0*d_pi) - 2.0*d;
    grad[1] = d_pi*sin(fay*2.0*d_pi) + d_sqrt2*2.0*d;
  }
  return;
}

int intrac_(void* dummy) {
  return 0;
}