Java Code Examples for org.apache.commons.math3.util.FastMath#pow()

/** Evaluate Taylor expansion of a derivative structure.
 * @param ds array holding the derivative structure
 * @param dsOffset offset of the derivative structure in its array
 * @param delta parameters offsets (Δx, Δy, ...)
 * @return value of the Taylor expansion at x + Δx, y + Δy, ...
 * @throws MathArithmeticException if factorials becomes too large
 */
public double taylor(final double[] ds, final int dsOffset, final double ... delta)
   throws MathArithmeticException {
    double value = 0;
    for (int i = getSize() - 1; i >= 0; --i) {
        final int[] orders = getPartialDerivativeOrders(i);
        double term = ds[dsOffset + i];
        for (int k = 0; k < orders.length; ++k) {
            if (orders[k] > 0) {
                try {
                    term *= FastMath.pow(delta[k], orders[k]) /
                    CombinatoricsUtils.factorial(orders[k]);
                } catch (NotPositiveException e) {
                    // this cannot happen
                    throw new MathInternalError(e);
                }
            }
        }
        value += term;
    }
    return value;
}
 

/**
 * Update the cluster centers.
 */
private void updateClusterCenters() {
    int j = 0;
    final List<CentroidCluster<T>> newClusters = new ArrayList<CentroidCluster<T>>(k);
    for (final CentroidCluster<T> cluster : clusters) {
        final Clusterable center = cluster.getCenter();
        int i = 0;
        double[] arr = new double[center.getPoint().length];
        double sum = 0.0;
        for (final T point : points) {
            final double u = FastMath.pow(membershipMatrix[i][j], fuzziness);
            final double[] pointArr = point.getPoint();
            for (int idx = 0; idx < arr.length; idx++) {
                arr[idx] += u * pointArr[idx];
            }
            sum += u;
            i++;
        }
        MathArrays.scaleInPlace(1.0 / sum, arr);
        newClusters.add(new CentroidCluster<T>(new DoublePoint(arr)));
        j++;
    }
    clusters.clear();
    clusters = newClusters;
}
 

public static void applyAdamUpdater(INDArray gradient, INDArray m, INDArray v, double learningRate, double beta1, double beta2,
                                                     double epsilon, int iteration){

    INDArray oneMinusBeta1Grad = gradient.mul(1.0 - beta1);
    m.muli(beta1).addi(oneMinusBeta1Grad);

    INDArray oneMinusBeta2GradSquared = gradient.mul(gradient).muli(1 - beta2);
    v.muli(beta2).addi(oneMinusBeta2GradSquared);

    double beta1t = FastMath.pow(beta1, iteration + 1);
    double beta2t = FastMath.pow(beta2, iteration + 1);

    double alphat = learningRate * FastMath.sqrt(1 - beta2t) / (1 - beta1t);
    if (Double.isNaN(alphat) || alphat == 0.0)
        alphat = epsilon;
    INDArray sqrtV = Transforms.sqrt(v.dup('c'), false).addi(epsilon);

    gradient.assign(m).muli(alphat).divi(sqrtV);
}
 

private static double yjInvTransSingle(double x, double lam) {
	/*
	 * This is where it gets messy, but we can theorize that
   	 * if the x is < 0 and the lambda meets the appropriate conditions,
   	 * that the x was sub-zero to begin with
	 */
	if(x >= 0) {
		// Case 1: x >= 0 and lambda is not 0
		if(!nearZero(lam)) {
			x *= lam;
			x += 1;
			x = FastMath.pow(x, 1.0 / lam);
			return x - 1;
		}
		
		// Case 2: x >= 0 and lambda is 0
		return FastMath.exp(x) - 1;
	} else {
		// Case 3: lambda does not equal 2
		if(lam != 2.0) {
			x *= -(2.0 - lam);
            x += 1;
            x = FastMath.pow(x, 1.0 / (2.0 - lam));
            x -= 1;
            return -x;
		}
		
		// Case 4: lambda equals 2
		return -(FastMath.exp(-x) - 1);
	}
}
 

@Test
public void testIncreasingTolerance()
  {

  int previousCalls = Integer.MAX_VALUE;
  AdaptiveStepsizeIntegrator integ =
      new DormandPrince853Integrator(0, Double.POSITIVE_INFINITY,
                                     Double.NaN, Double.NaN);
  for (int i = -12; i < -2; ++i) {
    TestProblem1 pb = new TestProblem1();
    double minStep = 0;
    double maxStep = pb.getFinalTime() - pb.getInitialTime();
    double scalAbsoluteTolerance = FastMath.pow(10.0, i);
    double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
    integ.setStepSizeControl(minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);

    TestProblemHandler handler = new TestProblemHandler(pb, integ);
    integ.addStepHandler(handler);
    integ.integrate(pb,
                    pb.getInitialTime(), pb.getInitialState(),
                    pb.getFinalTime(), new double[pb.getDimension()]);

    // the 1.3 factor is only valid for this test
    // and has been obtained from trial and error
    // there is no general relation between local and global errors
    Assert.assertTrue(handler.getMaximalValueError() < (1.3 * scalAbsoluteTolerance));
    Assert.assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);

    int calls = pb.getCalls();
    Assert.assertEquals(integ.getEvaluations(), calls);
    Assert.assertTrue(calls <= previousCalls);
    previousCalls = calls;

  }

}
 

@Test
public void testIncreasingTolerance()
  {

  int previousCalls = Integer.MAX_VALUE;
  for (int i = -12; i < -4; ++i) {
    TestProblem1 pb     = new TestProblem1();
    double minStep      = 0;
    double maxStep      = pb.getFinalTime() - pb.getInitialTime();
    double absTolerance = FastMath.pow(10.0, i);
    double relTolerance = absTolerance;

    FirstOrderIntegrator integ =
      new GraggBulirschStoerIntegrator(minStep, maxStep,
                                       absTolerance, relTolerance);
    TestProblemHandler handler = new TestProblemHandler(pb, integ);
    integ.addStepHandler(handler);
    integ.integrate(pb,
                    pb.getInitialTime(), pb.getInitialState(),
                    pb.getFinalTime(), new double[pb.getDimension()]);

    // the coefficients are only valid for this test
    // and have been obtained from trial and error
    // there is no general relation between local and global errors
    double ratio =  handler.getMaximalValueError() / absTolerance;
    Assert.assertTrue(ratio < 2.4);
    Assert.assertTrue(ratio > 0.02);
    Assert.assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);

    int calls = pb.getCalls();
    Assert.assertEquals(integ.getEvaluations(), calls);
    Assert.assertTrue(calls <= previousCalls);
    previousCalls = calls;

  }

}
 

/** {@inheritDoc} */
public ConfidenceInterval createInterval(int numberOfTrials, int numberOfSuccesses, double confidenceLevel) {
    IntervalUtils.checkParameters(numberOfTrials, numberOfSuccesses, confidenceLevel);
    final double alpha = (1.0 - confidenceLevel) / 2;
    final NormalDistribution normalDistribution = new NormalDistribution();
    final double z = normalDistribution.inverseCumulativeProbability(1 - alpha);
    final double zSquared = FastMath.pow(z, 2);
    final double modifiedNumberOfTrials = numberOfTrials + zSquared;
    final double modifiedSuccessesRatio = (1.0 / modifiedNumberOfTrials) * (numberOfSuccesses + 0.5 * zSquared);
    final double difference = z *
                              FastMath.sqrt(1.0 / modifiedNumberOfTrials * modifiedSuccessesRatio *
                                            (1 - modifiedSuccessesRatio));
    return new ConfidenceInterval(modifiedSuccessesRatio - difference, modifiedSuccessesRatio + difference,
                                  confidenceLevel);
}
 

public double value(double[] x) {
    double f = 0;
    double fac;
    for (int i = 0; i < x.length; ++i) {
        fac = FastMath.pow(axisratio, (i - 1.) / (x.length - 1.));
        if (i == 0 && x[i] < 0)
            fac *= 1.;
        f += fac * fac * x[i] * x[i] + amplitude
        * (1. - FastMath.cos(2. * FastMath.PI * fac * x[i]));
    }
    return f;
}
 

public double value(double[] x) {
    double f = 0;
    double res2 = 0;
    double fac = 0;
    for (int i = 0; i < x.length; ++i) {
        fac = FastMath.pow(axisratio, (i - 1.) / (x.length - 1.));
        f += fac * fac * x[i] * x[i];
        res2 += FastMath.cos(2. * FastMath.PI * fac * x[i]);
    }
    f = (20. - 20. * FastMath.exp(-0.2 * FastMath.sqrt(f / x.length))
            + FastMath.exp(1.) - FastMath.exp(res2 / x.length));
    return f;
}
 

/**
 * Unbiased skewness.
 * See  http://commons.apache.org/proper/commons-math/apidocs/org/apache/commons/math4/stat/descriptive/moment/Skewness.html
 * @return unbiased skewness
 */
@Override
public double getSkewness() {
  //  skewness = [n / (n -1) (n - 2)] sum[(x_i - mean)^3] / std^3
  if(n < 3) {
    return Double.NaN;
  }
  double t1 = (1.0*n)/((n - 1)*(n-2));
  double std = getStandardDeviation();
  return t1*M3/FastMath.pow(std, 3);
}
 

public static double density(final double x, final double shape,
                             final double scale) {
    /*
     * This is a copy of
     * double GammaDistribution.density(double)
     * prior to MATH-753.
     */
    if (x < 0) {
        return 0;
    }
    return FastMath.pow(x / scale, shape - 1) / scale *
           FastMath.exp(-x / scale) / FastMath.exp(logGamma(shape));
}
 

public double value(double[] x) {
    double f = 0;
    double res2 = 0;
    double fac = 0;
    for (int i = 0; i < x.length; ++i) {
        fac = FastMath.pow(axisratio, (i - 1.) / (x.length - 1.));
        f += fac * fac * x[i] * x[i];
        res2 += FastMath.cos(2. * FastMath.PI * fac * x[i]);
    }
    f = (20. - 20. * FastMath.exp(-0.2 * FastMath.sqrt(f / x.length))
            + FastMath.exp(1.) - FastMath.exp(res2 / x.length));
    return f;
}
 

/** {@inheritDoc} */
public UnivariateFunction derivative() {
    return new UnivariateFunction() {
        /** {@inheritDoc} */
        public double value(double x) {
            return p * FastMath.pow(x, p - 1);
        }
    };
}
 

@Test
public void testOnDistortedSine() {
    int numPoints = 100;
    double[] xval = new double[numPoints];
    double[] yval = new double[numPoints];
    double xnoise = 0.1;
    double ynoise = 0.2;

    generateSineData(xval, yval, xnoise, ynoise);

    LoessInterpolator li = new LoessInterpolator(0.3, 4, 1e-12);

    double[] res = li.smooth(xval, yval);

    // Check that the resulting curve differs from
    // the "real" sine less than the jittered one

    double noisyResidualSum = 0;
    double fitResidualSum = 0;

    for(int i = 0; i < numPoints; ++i) {
        double expected = FastMath.sin(xval[i]);
        double noisy = yval[i];
        double fit = res[i];

        noisyResidualSum += FastMath.pow(noisy - expected, 2);
        fitResidualSum += FastMath.pow(fit - expected, 2);
    }

    Assert.assertTrue(fitResidualSum < noisyResidualSum);
}
 

/** {@inheritDoc} */
public double cumulativeProbability(int x) {
    double ret;
    if (x < 0) {
        ret = 0.0;
    } else {
        final double p = probabilityOfSuccess;
        ret = 1.0 - FastMath.pow(1 - p, x + 1);
    }
    return ret;
}
 

/** {@inheritDoc} */
public double value(double x) {
    return FastMath.pow(x, p);
}
 

/** Initialize the integration step.
 * @param forward forward integration indicator
 * @param order order of the method
 * @param scale scaling vector for the state vector (can be shorter than state vector)
 * @param t0 start time
 * @param y0 state vector at t0
 * @param yDot0 first time derivative of y0
 * @param y1 work array for a state vector
 * @param yDot1 work array for the first time derivative of y1
 * @return first integration step
 * @exception MaxCountExceededException if the number of functions evaluations is exceeded
 * @exception DimensionMismatchException if arrays dimensions do not match equations settings
 */
public double initializeStep(final boolean forward, final int order, final double[] scale,
                             final double t0, final double[] y0, final double[] yDot0,
                             final double[] y1, final double[] yDot1)
    throws MaxCountExceededException, DimensionMismatchException {

  if (initialStep > 0) {
    // use the user provided value
    return forward ? initialStep : -initialStep;
  }

  // very rough first guess : h = 0.01 * ||y/scale|| / ||y'/scale||
  // this guess will be used to perform an Euler step
  double ratio;
  double yOnScale2 = 0;
  double yDotOnScale2 = 0;
  for (int j = 0; j < scale.length; ++j) {
    ratio         = y0[j] / scale[j];
    yOnScale2    += ratio * ratio;
    ratio         = yDot0[j] / scale[j];
    yDotOnScale2 += ratio * ratio;
  }

  double h = ((yOnScale2 < 1.0e-10) || (yDotOnScale2 < 1.0e-10)) ?
             1.0e-6 : (0.01 * FastMath.sqrt(yOnScale2 / yDotOnScale2));
  if (! forward) {
    h = -h;
  }

  // perform an Euler step using the preceding rough guess
  for (int j = 0; j < y0.length; ++j) {
    y1[j] = y0[j] + h * yDot0[j];
  }
  computeDerivatives(t0 + h, y1, yDot1);

  // estimate the second derivative of the solution
  double yDDotOnScale = 0;
  for (int j = 0; j < scale.length; ++j) {
    ratio         = (yDot1[j] - yDot0[j]) / scale[j];
    yDDotOnScale += ratio * ratio;
  }
  yDDotOnScale = FastMath.sqrt(yDDotOnScale) / h;

  // step size is computed such that
  // h^order * max (||y'/tol||, ||y''/tol||) = 0.01
  final double maxInv2 = FastMath.max(FastMath.sqrt(yDotOnScale2), yDDotOnScale);
  final double h1 = (maxInv2 < 1.0e-15) ?
                    FastMath.max(1.0e-6, 0.001 * FastMath.abs(h)) :
                    FastMath.pow(0.01 / maxInv2, 1.0 / order);
  h = FastMath.min(100.0 * FastMath.abs(h), h1);
  h = FastMath.max(h, 1.0e-12 * FastMath.abs(t0));  // avoids cancellation when computing t1 - t0
  if (h < getMinStep()) {
    h = getMinStep();
  }
  if (h > getMaxStep()) {
    h = getMaxStep();
  }
  if (! forward) {
    h = -h;
  }

  return h;

}
 

public double value(double[] x) {
    double f = 0;
    for (int i = 0; i < x.length; ++i)
        f += FastMath.pow(factor, i / (x.length - 1.)) * x[i] * x[i];
    return f;
}
 

@Test
public void testDecreasingSteps()
    throws DimensionMismatchException, NumberIsTooSmallException,
           MaxCountExceededException, NoBracketingException {

  TestProblemAbstract[] problems = TestProblemFactory.getProblems();
  for (int k = 0; k < problems.length; ++k) {

    double previousValueError = Double.NaN;
    double previousTimeError = Double.NaN;
    for (int i = 4; i < 10; ++i) {

      TestProblemAbstract pb = problems[k].copy();
      double step = (pb.getFinalTime() - pb.getInitialTime()) * FastMath.pow(2.0, -i);
      FirstOrderIntegrator integ = new MidpointIntegrator(step);
      TestProblemHandler handler = new TestProblemHandler(pb, integ);
      integ.addStepHandler(handler);
      EventHandler[] functions = pb.getEventsHandlers();
      for (int l = 0; l < functions.length; ++l) {
        integ.addEventHandler(functions[l],
                                   Double.POSITIVE_INFINITY, 1.0e-6 * step, 1000);
      }
      double stopTime = integ.integrate(pb,
                                        pb.getInitialTime(), pb.getInitialState(),
                                        pb.getFinalTime(), new double[pb.getDimension()]);
      if (functions.length == 0) {
          Assert.assertEquals(pb.getFinalTime(), stopTime, 1.0e-10);
      }

      double valueError = handler.getMaximalValueError();
      if (i > 4) {
        Assert.assertTrue(valueError < FastMath.abs(previousValueError));
      }
      previousValueError = valueError;

      double timeError = handler.getMaximalTimeError();
      if (i > 4) {
        Assert.assertTrue(timeError <= FastMath.abs(previousTimeError));
      }
      previousTimeError = timeError;

    }

  }

}
 

/**
 * <p>
 * Returns a <code>double</code> whose value is
 * <tt>(this<sup>exponent</sup>)</tt>, returning the result in reduced form.
 * </p>
 *
 * @param exponent
 *            exponent to which this <code>BigFraction</code> is to be raised.
 * @return <tt>this<sup>exponent</sup></tt>.
 */
public double pow(final double exponent) {
    return FastMath.pow(numerator.doubleValue(),   exponent) /
           FastMath.pow(denominator.doubleValue(), exponent);
}