ReneTest.java

package org.freehep.math.minuit.example.sim;

import java.util.List;
import org.freehep.math.minuit.FCNBase;
import org.freehep.math.minuit.FunctionMinimum;
import org.freehep.math.minuit.MnMigrad;
import org.freehep.math.minuit.MnPlot;
import org.freehep.math.minuit.MnScan;
import org.freehep.math.minuit.MnStrategy;
import org.freehep.math.minuit.MnUserParameters;
import org.freehep.math.minuit.Point;
/**
 *
 * @version $Id: ReneTest.java 8584 2006-08-10 23:06:37Z duns $
 */
public class ReneTest
{
   public static void main(String[] args)
   {
      double[] tmp =
      {38.,36.,46.,52.,54.,52.,61.,52.,64.,77.,
       60.,56.,78.,71.,81.,83.,89.,96.,118.,96.,
       109.,111.,107.,107.,135.,156.,196.,137.,
       160.,153.,185.,222.,251.,270.,329.,422.,
       543.,832.,1390.,2835.,3462.,2030.,1130.,
       657.,469.,411.,375.,295.,281.,281.,289.,
       273.,297.,256.,274.,287.,280.,274.,286.,
       279.,293.,314.,285.,322.,307.,313.,324.,
       351.,314.,314.,301.,361.,332.,342.,338.,
       396.,356.,344.,395.,416.,406.,411.,422.,
       393.,393.,409.,455.,427.,448.,459.,403.,
       441.,510.,501.,502.,482.,487.,506.,506.,
       526.,517.,534.,509.,482.,591.,569.,518.,
       609.,569.,598.,627.,617.,610.,662.,666.,
       652.,671.,647.,650.,701.};
       
       double[] measurements = tmp.clone();
       
       ReneFcn theFCN = new ReneFcn(measurements);
       
       MnUserParameters upar = new MnUserParameters();
       upar.add("p0", 100., 10.);
       upar.add("p1", 100., 10.);
       upar.add("p2", 100., 10.);
       upar.add("p3", 100., 10.);
       upar.add("p4", 1., 0.3);
       upar.add("p5", 1., 0.3);
       
       System.out.println("Initial parameters: "+upar);
       
       System.out.println("start migrad");
       MnMigrad migrad = new MnMigrad(theFCN, upar);
       FunctionMinimum min = migrad.minimize();
       if(!min.isValid())
       {
          //try with higher strategy
          System.out.println("FM is invalid, try with strategy = 2.");
          MnMigrad migrad2 = new MnMigrad(theFCN, min.userState(), new MnStrategy(2));
          min = migrad2.minimize();
       }
       System.out.println("minimum: "+min);
       {
          double[] params = {1,1,1,1,1,1};
          double[] error = {1,1,1,1,1,1};
          MnScan scan = new MnScan(theFCN, params, error);
          System.out.println("scan parameters: "+scan.parameters());
          MnPlot plot = new MnPlot();
          for(int i = 0; i < upar.variableParameters(); i++)
          {
             List<Point> xy = scan.scan(i);
             plot.plot(xy);
          }
          System.out.println("scan parameters: "+scan.parameters());
       }
       
       {
          double[] params = {1,1,1,1,1,1};
          double[] error = {1,1,1,1,1,1};
          MnScan scan = new MnScan(theFCN, params, error);
          System.out.println("scan parameters: "+scan.parameters());
          FunctionMinimum min2 = scan.minimize();
          //     std::cout<<min<<std::endl;
          System.out.println("scan parameters: "+scan.parameters());
       }
       
   }
   static class ReneFcn implements FCNBase
   {
      ReneFcn(double[] meas)
      {
         theMeasurements = meas;
      }
      public double errorDef()
      {
         return 1;
      }
      
      public double valueOf(double[] par)
      {
         double a = par[2];
         double b = par[1];
         double c = par[0];
         double p0 = par[3];
         double p1 = par[4];
         double p2 = par[5];
         double fval = 0.;
         for( int i = 0; i < theMeasurements.length; i++)
         {
            double ni = theMeasurements[i];
            if(ni < 1.e-10) continue;
            double xi = (i+1.)/40. - 1./80.; //xi=0-3
            double ei = ni;
            double nexp = a*xi*xi + b*xi + c + (0.5*p0*p1/Math.PI)/Math.max(1.e-10, (xi-p2)*(xi-p2) + 0.25*p1*p1);
            fval += (ni-nexp)*(ni-nexp)/ei;
         }
         return fval;
      }
      private double[] theMeasurements;
   }
}