NI3.java

/*
 * Class:        NI3
 * Description:
 * Environment:  Java
 * Software:     SSJ
 * Copyright (C) 2001  Pierre L'Ecuyer and Université de Montréal
 * Organization: DIRO, Université de Montréal
 * @author       Nabil Channouf
 * @since
 *
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *
 */
package umontreal.ssj.probdistmulti.norta;
 
import umontreal.ssj.probdist.*;

public class NI3 extends NortaInitDisc {
   private double tolerance; /*
                              * Desired accuracy for the root-finder algorithm (epsilon in paragraph
                              * "Method NI3" of section 3 in paper).
                              */

   public NI3(double rX, DiscreteDistributionInt dist1, DiscreteDistributionInt dist2, double tr, double tolerance) {
      super(rX, dist1, dist2, tr);
      this.tolerance = tolerance;
      computeParams();
   }

   public double computeCorr() {
      final double ITMAX = 100; // Maximum number of iterations.
      double xl, xh; /*
                      * Left and right endpoints of the root bracket at all iterations.
                      */
      double b = 0.0; // The returned solution.
      double f, df; // Function and its derivative evaluations.
      double dx; // Correction term.
      double dxold; // The root correction at one iteration before the last one.
      double temp;// The root at one iteration before the last one.
      double ccc = rX * sd1 * sd2 + mu1 * mu2; // Precompute constant.
 
      if (rX == 0.0)
         return 0.0;
      if (rX > 0.0) { // Orient the search
         xl = 0.0;
         xh = 1.0;
      } else {
         xl = -1.0;
         xh = 0.0;
      }
 
      b = 2 * Math.sin(Math.PI * rX / 6); // Initial guess
      dxold = xh - xl;
      dx = dxold;
      f = integ(b) - ccc;
      df = deriv(b);
 
      for (int i = 1; i <= ITMAX; i++) { // Begin the search
         if ((((b - xh) * df - f) * ((b - xl) * df - f) > 0.0) || (Math.abs(2.0 * f) > Math.abs(dxold * df))) {
            // Do bisection if solution is out of range
            // or not decreasing fast enough
            dxold = dx;
            dx = 0.5 * (xh - xl);
            b = xl + dx;
            if (xl == b) // Change in root is negligible.
               return b; // Accept this root
         } else {
            dxold = dx;
            dx = f / df;
            temp = b;
            b -= dx;
            if (temp == b)
               return b;
         }
         if (Math.abs(dx) < tolerance)
            return b; // Convergence check
         f = integ(b) - ccc;
         df = deriv(b);
         if (f < 0.0)
            xl = b; // Maintain the brackets on the root
         else
            xh = b;
      }
      return b;
   }
 
   public String toString() {
      // To display the inputs.
      String desc = super.toString();
      desc += "tolerance : " + tolerance + "\n";
      return desc;
   }
 
}