package org.freehep.math.minuit;

import org.apache.commons.math3.util.FastMath;

/** API class for Minos error analysis (asymmetric errors).
 * Minimization has to be done before and minimum must be valid;
 * possibility to ask only for one side of the Minos error;
 * @version $Id: MnMinos.java 8584 2006-08-10 23:06:37Z duns $
 */
public class MnMinos
{
   /** construct from FCN + minimum */
   public MnMinos(FCNBase fcn, FunctionMinimum min)
   {
      this(fcn,min,MnApplication.DEFAULT_STRATEGY);
   }
   
   /** construct from FCN + minimum + strategy */
   public MnMinos(FCNBase fcn, FunctionMinimum min, int stra)
   {
      this(fcn,min,new MnStrategy(stra));
   }
   
   /** construct from FCN + minimum + strategy */
   public MnMinos(FCNBase fcn, FunctionMinimum min, MnStrategy stra)
   {
      theFCN = fcn;
      theMinimum = min;
      theStrategy = stra;
   }
   public MinosError minos(int par)
   {
      return minos(par,1.);
   }
   public MinosError minos(int par, double errDef)
   {
      return minos(par,errDef, MnApplication.DEFAULT_MAXFCN);
   }   
   /**
    * Causes a MINOS error analysis to be performed on the parameter whose number is
    * specified. MINOS errors may be expensive to calculate, but are very reliable since
    * they take account of non-linearities in the problem as well as parameter correlations,
    * and are in general asymmetric.
    * @param maxcalls Specifies the (approximate) maximum number of function calls per parameter 
    * requested, after which the calculation will be stopped for that parameter.
    */
   public MinosError minos(int par,  double errDef, int maxcalls)
   {
      assert(theMinimum.isValid());
      assert(!theMinimum.userState().parameter(par).isFixed());
      assert(!theMinimum.userState().parameter(par).isConst());
      
      MnCross up = upval(par, errDef, maxcalls);
      MnCross lo = loval(par, errDef, maxcalls);
      
      return new MinosError(par, theMinimum.userState().value(par), lo, up);
   }
   public Point range(int par)
   {
      return range(par,1);
   }
   public Point range(int par, double errDef)
   {
      return range(par,errDef, MnApplication.DEFAULT_MAXFCN);
   }
   /**
    * Causes a MINOS error analysis for external parameter n.
    * @return The lower and upper bounds of parameter 
    */
   public Point range(int par, double errDef, int maxcalls)
   {
      MinosError mnerr = minos(par, errDef, maxcalls);
      return mnerr.range();
   }
   public double lower(int par)
   {
      return lower(par,1);
   }
   public double lower(int par, double errDef)
   {
      return lower(par,errDef, MnApplication.DEFAULT_MAXFCN);
   }  
   /** calculate one side (negative or positive error) of the parameter */
   public double lower(int par, double errDef, int maxcalls)
   {
      MnUserParameterState upar = theMinimum.userState();
      double err = theMinimum.userState().error(par);
      MnCross aopt = loval(par, errDef, maxcalls);
      double lower = aopt.isValid() ? -1.*err*(1.+ aopt.value()) : (aopt.atLimit() ? upar.parameter(par).lowerLimit() : upar.value(par));
      return lower;
   }
   public double upper(int par)
   {
      return upper(par,1);
   }
   public double upper(int par, double errDef)
   {
      return upper(par,errDef, MnApplication.DEFAULT_MAXFCN);
   } 
   public double upper(int par, double errDef, int maxcalls)
   {
      MnUserParameterState upar = theMinimum.userState();
      double err = theMinimum.userState().error(par);
      MnCross aopt = upval(par, errDef, maxcalls);
      double upper = aopt.isValid() ? err*(1.+ aopt.value()) : (aopt.atLimit() ? upar.parameter(par).upperLimit() : upar.value(par));
      return upper;
   }
   public MnCross loval(int par)
   {
      return loval(par,1);
   }
   public MnCross loval(int par, double errDef)
   {
      return loval(par,errDef, MnApplication.DEFAULT_MAXFCN);
   } 
   public MnCross loval(int par, double errDef, int maxcalls)
   {
      errDef *= theMinimum.errorDef();
      assert(theMinimum.isValid());
      assert(!theMinimum.userState().parameter(par).isFixed());
      assert(!theMinimum.userState().parameter(par).isConst());
      if(maxcalls == 0)
      {
         int nvar = theMinimum.userState().variableParameters();
         maxcalls = 2*(nvar+1)*(200 + 100*nvar + 5*nvar*nvar);
      }
      int[] para = {par};
      
      MnUserParameterState upar = theMinimum.userState().clone();
      double err = upar.error(par);
      double val = upar.value(par) - err;
      double[] xmid = {val};
      double[] xdir = {-err};
      
      int ind = upar.intOfExt(par);
      MnAlgebraicSymMatrix m = theMinimum.error().matrix();
      double xunit = FastMath.sqrt(errDef/err);
      for(int i = 0; i < m.nrow(); i++)
      {
         if(i == ind) continue;
         double xdev = xunit*m.get(ind,i);
         int ext = upar.extOfInt(i);
         upar.setValue(ext, upar.value(ext) - xdev);
      }
      
      upar.fix(par);
      upar.setValue(par, val);
      
      double toler = 0.1;
      MnFunctionCross cross = new MnFunctionCross(theFCN, upar, theMinimum.fval(), theStrategy, errDef);
      
      MnCross aopt = cross.cross(para, xmid, xdir, toler, maxcalls);
      
      if(aopt.atLimit())
         System.out.println("MnMinos parameter "+par+" is at lower limit.");
      if(aopt.atMaxFcn())
         System.out.println("MnMinos maximum number of function calls exceeded for parameter "+par);
      if(aopt.newMinimum())
         System.out.println("MnMinos new minimum found while looking for parameter "+par);
      if(!aopt.isValid())
         System.out.println("MnMinos could not find lower value for parameter "+par+".");
      
      return aopt;
      
   }
   public MnCross upval(int par)
   {
      return upval(par,1);
   }
   public MnCross upval(int par, double errDef)
   {
      return upval(par,errDef, MnApplication.DEFAULT_MAXFCN);
   } 
   public MnCross upval(int par, double errDef, int maxcalls)
   {
      errDef *= theMinimum.errorDef();
      assert(theMinimum.isValid());
      assert(!theMinimum.userState().parameter(par).isFixed());
      assert(!theMinimum.userState().parameter(par).isConst());
      if(maxcalls == 0)
      {
         int nvar = theMinimum.userState().variableParameters();
         maxcalls = 2*(nvar+1)*(200 + 100*nvar + 5*nvar*nvar);
      }
      
      int[] para = { par };
      
      MnUserParameterState upar = theMinimum.userState().clone();
      double err = upar.error(par);
      double val = upar.value(par) + err;
      double[] xmid = {val};
      double[] xdir = {err};
      
      int ind = upar.intOfExt(par);
      MnAlgebraicSymMatrix m = theMinimum.error().matrix();
      double xunit = FastMath.sqrt(errDef/err);
      for(int i = 0; i < m.nrow(); i++)
      {
         if(i == ind) continue;
         double xdev = xunit*m.get(ind,i);
         int ext = upar.extOfInt(i);
         upar.setValue(ext, upar.value(ext) + xdev);
      }
      
      upar.fix(par);
      upar.setValue(par, val);
      
      double toler = 0.1;
      MnFunctionCross cross = new MnFunctionCross(theFCN, upar, theMinimum.fval(), theStrategy, errDef);
      MnCross aopt = cross.cross(para, xmid, xdir, toler, maxcalls);
      
      if(aopt.atLimit())
         System.err.println("MnMinos parameter "+par+" is at upper limit.");
      if(aopt.atMaxFcn())
         System.err.println("MnMinos maximum number of function calls exceeded for parameter "+par);
      if(aopt.newMinimum())
         System.err.println("MnMinos new minimum found while looking for parameter "+par);
      if(!aopt.isValid())
         System.err.println("MnMinos could not find upper value for parameter "+par+".");
      
      return aopt;
   }
   
   private FCNBase theFCN;
   private FunctionMinimum theMinimum;
   private MnStrategy theStrategy;
}