org.apache.commons.math3.optim.univariate.SearchInterval Java Examples

/**
 * Given a list of posterior samples, returns an estimate of the posterior mode (using
 * mllib kernel density estimation in {@link KernelDensity} and {@link BrentOptimizer}).
 * Note that estimate may be poor if number of samples is small (resulting in poor kernel density estimation),
 * or if posterior is not unimodal (or is sufficiently pathological otherwise). If the samples contain
 * {@link Double#NaN}, {@link Double#NaN} will be returned.
 * @param samples   posterior samples, cannot be {@code null} and number of samples must be greater than 0
 * @param ctx       {@link JavaSparkContext} used by {@link KernelDensity} for mllib kernel density estimation
 */
public static double calculatePosteriorMode(final List<Double> samples, final JavaSparkContext ctx) {
    Utils.nonNull(samples);
    Utils.validateArg(samples.size() > 0, "Number of samples must be greater than zero.");

    //calculate sample min, max, mean, and standard deviation
    final double sampleMin = Collections.min(samples);
    final double sampleMax = Collections.max(samples);
    final double sampleMean = new Mean().evaluate(Doubles.toArray(samples));
    final double sampleStandardDeviation = new StandardDeviation().evaluate(Doubles.toArray(samples));

    //if samples are all the same or contain NaN, can simply return mean
    if (sampleStandardDeviation == 0. || Double.isNaN(sampleMean)) {
        return sampleMean;
    }

    //use Silverman's rule to set bandwidth for kernel density estimation from sample standard deviation
    //see https://en.wikipedia.org/wiki/Kernel_density_estimation#Practical_estimation_of_the_bandwidth
    final double bandwidth =
            SILVERMANS_RULE_CONSTANT * sampleStandardDeviation * Math.pow(samples.size(), SILVERMANS_RULE_EXPONENT);

    //use kernel density estimation to approximate posterior from samples
    final KernelDensity pdf = new KernelDensity().setSample(ctx.parallelize(samples, 1)).setBandwidth(bandwidth);

    //use Brent optimization to find mode (i.e., maximum) of kernel-density-estimated posterior
    final BrentOptimizer optimizer =
            new BrentOptimizer(RELATIVE_TOLERANCE, RELATIVE_TOLERANCE * (sampleMax - sampleMin));
    final UnivariateObjectiveFunction objective =
            new UnivariateObjectiveFunction(f -> pdf.estimate(new double[] {f})[0]);
    //search for mode within sample range, start near sample mean
    final SearchInterval searchInterval = new SearchInterval(sampleMin, sampleMax, sampleMean);
    return optimizer.optimize(objective, GoalType.MAXIMIZE, searchInterval, BRENT_MAX_EVAL).getPoint();
}
 

/**
 * Find the minimum of the function {@code f(p + alpha * d)}.
 *
 * @param p Starting point.
 * @param d Search direction.
 * @return the optimum.
 * @throws org.apache.commons.math3.exception.TooManyEvaluationsException
 * if the number of evaluations is exceeded.
 */
public UnivariatePointValuePair search(final double[] p, final double[] d) {
    final int n = p.length;
    final UnivariateFunction f = new UnivariateFunction() {
            public double value(double alpha) {
                final double[] x = new double[n];
                for (int i = 0; i < n; i++) {
                    x[i] = p[i] + alpha * d[i];
                }
                final double obj = PowellOptimizer.this.computeObjectiveValue(x);
                return obj;
            }
        };

    final GoalType goal = PowellOptimizer.this.getGoalType();
    bracket.search(f, goal, 0, 1);
    // Passing "MAX_VALUE" as a dummy value because it is the enclosing
    // class that counts the number of evaluations (and will eventually
    // generate the exception).
    return optimize(new MaxEval(Integer.MAX_VALUE),
                    new UnivariateObjectiveFunction(f),
                    goal,
                    new SearchInterval(bracket.getLo(),
                                       bracket.getHi(),
                                       bracket.getMid()));
}
 

/**
 * Find the minimum of the function {@code f(p + alpha * d)}.
 *
 * @param p Starting point.
 * @param d Search direction.
 * @return the optimum.
 * @throws org.apache.commons.math3.exception.TooManyEvaluationsException
 * if the number of evaluations is exceeded.
 */
public UnivariatePointValuePair search(final double[] p, final double[] d) {
    final int n = p.length;
    final UnivariateFunction f = new UnivariateFunction() {
            public double value(double alpha) {
                final double[] x = new double[n];
                for (int i = 0; i < n; i++) {
                    x[i] = p[i] + alpha * d[i];
                }
                final double obj = PowellOptimizer.this.computeObjectiveValue(x);
                return obj;
            }
        };

    final GoalType goal = PowellOptimizer.this.getGoalType();
    bracket.search(f, goal, 0, 1);
    // Passing "MAX_VALUE" as a dummy value because it is the enclosing
    // class that counts the number of evaluations (and will eventually
    // generate the exception).
    return optimize(new MaxEval(Integer.MAX_VALUE),
                    new UnivariateObjectiveFunction(f),
                    goal,
                    new SearchInterval(bracket.getLo(),
                                       bracket.getHi(),
                                       bracket.getMid()));
}
 

/**
 * Find the minimum of the function {@code f(p + alpha * d)}.
 *
 * @param p Starting point.
 * @param d Search direction.
 * @return the optimum.
 * @throws org.apache.commons.math3.exception.TooManyEvaluationsException
 * if the number of evaluations is exceeded.
 */
public UnivariatePointValuePair search(final double[] p, final double[] d) {
    final int n = p.length;
    final UnivariateFunction f = new UnivariateFunction() {
            public double value(double alpha) {
                final double[] x = new double[n];
                for (int i = 0; i < n; i++) {
                    x[i] = p[i] + alpha * d[i];
                }
                final double obj = PowellOptimizer.this.computeObjectiveValue(x);
                return obj;
            }
        };

    final GoalType goal = PowellOptimizer.this.getGoalType();
    bracket.search(f, goal, 0, 1);
    // Passing "MAX_VALUE" as a dummy value because it is the enclosing
    // class that counts the number of evaluations (and will eventually
    // generate the exception).
    return optimize(new MaxEval(Integer.MAX_VALUE),
                    new UnivariateObjectiveFunction(f),
                    goal,
                    new SearchInterval(bracket.getLo(),
                                       bracket.getHi(),
                                       bracket.getMid()));
}
 

/**
 * Find the minimum of the function {@code f(p + alpha * d)}.
 *
 * @param p Starting point.
 * @param d Search direction.
 * @return the optimum.
 * @throws org.apache.commons.math3.exception.TooManyEvaluationsException
 * if the number of evaluations is exceeded.
 */
public UnivariatePointValuePair search(final double[] p, final double[] d) {
    final int n = p.length;
    final UnivariateFunction f = new UnivariateFunction() {
            public double value(double alpha) {
                final double[] x = new double[n];
                for (int i = 0; i < n; i++) {
                    x[i] = p[i] + alpha * d[i];
                }
                final double obj = PowellOptimizer.this.computeObjectiveValue(x);
                return obj;
            }
        };

    final GoalType goal = PowellOptimizer.this.getGoalType();
    bracket.search(f, goal, 0, 1);
    // Passing "MAX_VALUE" as a dummy value because it is the enclosing
    // class that counts the number of evaluations (and will eventually
    // generate the exception).
    return optimize(new MaxEval(Integer.MAX_VALUE),
                    new UnivariateObjectiveFunction(f),
                    goal,
                    new SearchInterval(bracket.getLo(),
                                       bracket.getHi(),
                                       bracket.getMid()));
}
 

/**
 * Find the minimum of the function {@code f(p + alpha * d)}.
 *
 * @param p Starting point.
 * @param d Search direction.
 * @return the optimum.
 * @throws org.apache.commons.math3.exception.TooManyEvaluationsException
 * if the number of evaluations is exceeded.
 */
public UnivariatePointValuePair search(final double[] p, final double[] d) {
    final int n = p.length;
    final UnivariateFunction f = new UnivariateFunction() {
            public double value(double alpha) {
                final double[] x = new double[n];
                for (int i = 0; i < n; i++) {
                    x[i] = p[i] + alpha * d[i];
                }
                final double obj = PowellOptimizer.this.computeObjectiveValue(x);
                return obj;
            }
        };

    final GoalType goal = PowellOptimizer.this.getGoalType();
    bracket.search(f, goal, 0, 1);
    // Passing "MAX_VALUE" as a dummy value because it is the enclosing
    // class that counts the number of evaluations (and will eventually
    // generate the exception).
    return optimize(new MaxEval(Integer.MAX_VALUE),
                    new UnivariateObjectiveFunction(f),
                    goal,
                    new SearchInterval(bracket.getLo(),
                                       bracket.getHi(),
                                       bracket.getMid()));
}
 

private double[] calcNewRhos(final List<ACNVModeledSegment> segments,
                             final List<double[][][]> responsibilitiesBySeg,
                             final double lambda, final double[] rhos, final int[] mVals, final int[] nVals,
                             final JavaSparkContext ctx) {

    // Since, we pass in the entire responsibilities matrix, we need the correct index for each rho.  That, and the
    //  fact that this is a univariate objective function, means we need to create an instance for each rho.  And
    //  then we blast across Spark.

    final List<Pair<? extends Function<Double, Double>, SearchInterval>> objectives = IntStream.range(0, rhos.length)
            .mapToObj(i -> new Pair<>(
                    new Function<Double, Double>() {
                        @Override
                        public Double apply(Double rho) {
                            return calculateESmnObjective(rho, segments, responsibilitiesBySeg, mVals, nVals, lambda, i);
                        }
                    },
                    new SearchInterval(0.0, 1.0, rhos[i])))
            .collect(Collectors.toList());

    final JavaRDD<Pair<? extends Function<Double, Double>, SearchInterval>> objectivesRDD = ctx.parallelize(objectives);

    final List<Double> resultsAsDouble = objectivesRDD
            .map(objective -> optimizeIt(objective.getFirst(), objective.getSecond()))
            .collect();

    return resultsAsDouble.stream().mapToDouble(Double::doubleValue).toArray();
}
 

private double optimizeIt(final Function<Double, Double> objectiveFxn, final SearchInterval searchInterval) {
    final MaxEval BRENT_MAX_EVAL = new MaxEval(1000);
    final double RELATIVE_TOLERANCE = 0.001;
    final double ABSOLUTE_TOLERANCE = 0.001;
    final BrentOptimizer OPTIMIZER = new BrentOptimizer(RELATIVE_TOLERANCE, ABSOLUTE_TOLERANCE);

    final UnivariateObjectiveFunction objective = new UnivariateObjectiveFunction(x -> objectiveFxn.apply(x));
    return OPTIMIZER.optimize(objective, GoalType.MAXIMIZE, searchInterval, BRENT_MAX_EVAL).getPoint();
}
 

/**
 * Finds the number {@code alpha} that optimizes
 * {@code f(startPoint + alpha * direction)}.
 *
 * @param startPoint Starting point.
 * @param direction Search direction.
 * @return the optimum.
 * @throws org.apache.commons.math3.exception.TooManyEvaluationsException
 * if the number of evaluations is exceeded.
 */
public UnivariatePointValuePair search(final double[] startPoint,
                                       final double[] direction) {
    final int n = startPoint.length;
    final UnivariateFunction f = new UnivariateFunction() {
            public double value(double alpha) {
                final double[] x = new double[n];
                for (int i = 0; i < n; i++) {
                    x[i] = startPoint[i] + alpha * direction[i];
                }
                final double obj = mainOptimizer.computeObjectiveValue(x);
                return obj;
            }
        };

    final GoalType goal = mainOptimizer.getGoalType();
    bracket.search(f, goal, 0, initialBracketingRange);
    // Passing "MAX_VALUE" as a dummy value because it is the enclosing
    // class that counts the number of evaluations (and will eventually
    // generate the exception).
    return lineOptimizer.optimize(new MaxEval(Integer.MAX_VALUE),
                                  new UnivariateObjectiveFunction(f),
                                  goal,
                                  new SearchInterval(bracket.getLo(),
                                                     bracket.getHi(),
                                                     bracket.getMid()));
}
 

/**
 * Finds the number {@code alpha} that optimizes
 * {@code f(startPoint + alpha * direction)}.
 *
 * @param startPoint Starting point.
 * @param direction Search direction.
 * @return the optimum.
 * @throws org.apache.commons.math3.exception.TooManyEvaluationsException
 * if the number of evaluations is exceeded.
 */
public UnivariatePointValuePair search(final double[] startPoint,
                                       final double[] direction) {
    final int n = startPoint.length;
    final UnivariateFunction f = new UnivariateFunction() {
            public double value(double alpha) {
                final double[] x = new double[n];
                for (int i = 0; i < n; i++) {
                    x[i] = startPoint[i] + alpha * direction[i];
                }
                final double obj = mainOptimizer.computeObjectiveValue(x);
                return obj;
            }
        };

    final GoalType goal = mainOptimizer.getGoalType();
    bracket.search(f, goal, 0, initialBracketingRange);
    // Passing "MAX_VALUE" as a dummy value because it is the enclosing
    // class that counts the number of evaluations (and will eventually
    // generate the exception).
    return lineOptimizer.optimize(new MaxEval(Integer.MAX_VALUE),
                                  new UnivariateObjectiveFunction(f),
                                  goal,
                                  new SearchInterval(bracket.getLo(),
                                                     bracket.getHi(),
                                                     bracket.getMid()));
}
 

/**
 * Given a list of posterior samples, returns an estimate of the posterior mode (using
 * mllib kernel density estimation in {@link KernelDensity} and {@link BrentOptimizer}).
 * Note that estimate may be poor if number of samples is small (resulting in poor kernel density estimation),
 * or if posterior is not unimodal (or is sufficiently pathological otherwise). If the samples contain
 * {@link Double#NaN}, {@link Double#NaN} will be returned.
 * @param samples   posterior samples, cannot be {@code null} and number of samples must be greater than 0
 * @param ctx       {@link JavaSparkContext} used by {@link KernelDensity} for mllib kernel density estimation
 */
public static double calculatePosteriorMode(final List<Double> samples, final JavaSparkContext ctx) {
    Utils.nonNull(samples);
    Utils.validateArg(samples.size() > 0, "Number of samples must be greater than zero.");

    //calculate sample min, max, mean, and standard deviation
    final double sampleMin = Collections.min(samples);
    final double sampleMax = Collections.max(samples);
    final double sampleMean = new Mean().evaluate(Doubles.toArray(samples));
    final double sampleStandardDeviation = new StandardDeviation().evaluate(Doubles.toArray(samples));

    //if samples are all the same or contain NaN, can simply return mean
    if (sampleStandardDeviation == 0. || Double.isNaN(sampleMean)) {
        return sampleMean;
    }

    //use Silverman's rule to set bandwidth for kernel density estimation from sample standard deviation
    //see https://en.wikipedia.org/wiki/Kernel_density_estimation#Practical_estimation_of_the_bandwidth
    final double bandwidth =
            SILVERMANS_RULE_CONSTANT * sampleStandardDeviation * Math.pow(samples.size(), SILVERMANS_RULE_EXPONENT);

    //use kernel density estimation to approximate posterior from samples
    final KernelDensity pdf = new KernelDensity().setSample(ctx.parallelize(samples, 1)).setBandwidth(bandwidth);

    //use Brent optimization to find mode (i.e., maximum) of kernel-density-estimated posterior
    final BrentOptimizer optimizer =
            new BrentOptimizer(RELATIVE_TOLERANCE, RELATIVE_TOLERANCE * (sampleMax - sampleMin));
    final UnivariateObjectiveFunction objective =
            new UnivariateObjectiveFunction(f -> pdf.estimate(new double[] {f})[0]);
    //search for mode within sample range, start near sample mean
    final SearchInterval searchInterval = new SearchInterval(sampleMin, sampleMax, sampleMean);
    return optimizer.optimize(objective, GoalType.MAXIMIZE, searchInterval, BRENT_MAX_EVAL).getPoint();
}
 

public static double argmax(final Function<Double, Double> function, final double min, final double max, final double guess) {
    final SearchInterval interval = new SearchInterval(min, max, guess);
    return DEFAULT_OPTIMIZER.optimize(new UnivariateObjectiveFunction(function::apply), GoalType.MAXIMIZE, interval, DEFAULT_MAX_EVAL).getPoint();
}
 

public static double argmax(final Function<Double, Double> function, final double min, final double max, final double guess,
                            final double relativeTolerance, final double absoluteTolerance, final int maxEvaluations) {
    final BrentOptimizer optimizer = new BrentOptimizer(relativeTolerance, absoluteTolerance);
    final SearchInterval interval = new SearchInterval(min, max, guess);
    return optimizer.optimize(new UnivariateObjectiveFunction(function::apply), GoalType.MAXIMIZE, interval, new MaxEval(maxEvaluations)).getPoint();
}
 

/**
 * Maximizes a univariate function using a grid search followed by Brent's algorithm.
 *
 * @param fn the likelihood function to minimize
 * @param gridStart the lower bound for the grid search
 * @param gridEnd the upper bound for the grid search
 * @param gridStep step size for the grid search
 * @param relErr relative error tolerance for Brent's algorithm
 * @param absErr absolute error tolerance for Brent's algorithm
 * @param maxIter maximum # of iterations to perform in Brent's algorithm
 * @param maxEval maximum # of Likelihood function evaluations in Brent's algorithm
 *
 * @return the value of the parameter that maximizes the function
 */
public static double maximize(UnivariateFunction fn, double gridStart, double gridEnd,
    double gridStep, double relErr, double absErr, int maxIter, int maxEval) {
  Interval interval = gridSearch(fn, gridStart, gridEnd, gridStep);
  BrentOptimizer bo = new BrentOptimizer(relErr, absErr);
  UnivariatePointValuePair max = bo.optimize(
      new MaxIter(maxIter),
      new MaxEval(maxEval),
      new SearchInterval(interval.getInf(), interval.getSup()),
      new UnivariateObjectiveFunction(fn),
      GoalType.MAXIMIZE);
  return max.getPoint();
}