commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From er...@apache.org
Subject svn commit: r1420684 [4/15] - in /commons/proper/math/trunk/src: main/java/org/apache/commons/math3/exception/ main/java/org/apache/commons/math3/exception/util/ main/java/org/apache/commons/math3/fitting/ main/java/org/apache/commons/math3/optim/ main...
Date Wed, 12 Dec 2012 14:11:04 GMT
Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optim/nonlinear/scalar/MultivariateFunctionPenaltyAdapter.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optim/nonlinear/scalar/MultivariateFunctionPenaltyAdapter.java?rev=1420684&view=auto
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optim/nonlinear/scalar/MultivariateFunctionPenaltyAdapter.java (added)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optim/nonlinear/scalar/MultivariateFunctionPenaltyAdapter.java Wed Dec 12 14:10:38 2012
@@ -0,0 +1,187 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You 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 org.apache.commons.math3.optim.nonlinear.scalar;
+
+import org.apache.commons.math3.analysis.MultivariateFunction;
+import org.apache.commons.math3.exception.DimensionMismatchException;
+import org.apache.commons.math3.exception.NumberIsTooSmallException;
+import org.apache.commons.math3.util.FastMath;
+import org.apache.commons.math3.util.MathUtils;
+
+/**
+ * <p>Adapter extending bounded {@link MultivariateFunction} to an unbouded
+ * domain using a penalty function.</p>
+ *
+ * <p>
+ * This adapter can be used to wrap functions subject to simple bounds on
+ * parameters so they can be used by optimizers that do <em>not</em> directly
+ * support simple bounds.
+ * </p>
+ * <p>
+ * The principle is that the user function that will be wrapped will see its
+ * parameters bounded as required, i.e when its {@code value} method is called
+ * with argument array {@code point}, the elements array will fulfill requirement
+ * {@code lower[i] <= point[i] <= upper[i]} for all i. Some of the components
+ * may be unbounded or bounded only on one side if the corresponding bound is
+ * set to an infinite value. The optimizer will not manage the user function by
+ * itself, but it will handle this adapter and it is this adapter that will take
+ * care the bounds are fulfilled. The adapter {@link #value(double[])} method will
+ * be called by the optimizer with unbound parameters, and the adapter will check
+ * if the parameters is within range or not. If it is in range, then the underlying
+ * user function will be called, and if it is not the value of a penalty function
+ * will be returned instead.
+ * </p>
+ * <p>
+ * This adapter is only a poor-man's solution to simple bounds optimization
+ * constraints that can be used with simple optimizers like
+ * {@link org.apache.commons.math3.optim.nonlinear.scalar.noderiv.SimplexOptimizer
+ * SimplexOptimizer}.
+ * A better solution is to use an optimizer that directly supports simple bounds like
+ * {@link org.apache.commons.math3.optim.nonlinear.scalar.noderiv.CMAESOptimizer
+ * CMAESOptimizer} or
+ * {@link org.apache.commons.math3.optim.nonlinear.scalar.noderiv.BOBYQAOptimizer
+ * BOBYQAOptimizer}.
+ * One caveat of this poor-man's solution is that if start point or start simplex
+ * is completely outside of the allowed range, only the penalty function is used,
+ * and the optimizer may converge without ever entering the range.
+ * </p>
+ *
+ * @see MultivariateFunctionMappingAdapter
+ *
+ * @version $Id: MultivariateFunctionPenaltyAdapter.java 1416643 2012-12-03 19:37:14Z tn $
+ * @since 3.0
+ */
+public class MultivariateFunctionPenaltyAdapter
+    implements MultivariateFunction {
+    /** Underlying bounded function. */
+    private final MultivariateFunction bounded;
+    /** Lower bounds. */
+    private final double[] lower;
+    /** Upper bounds. */
+    private final double[] upper;
+    /** Penalty offset. */
+    private final double offset;
+    /** Penalty scales. */
+    private final double[] scale;
+
+    /**
+     * Simple constructor.
+     * <p>
+     * When the optimizer provided points are out of range, the value of the
+     * penalty function will be used instead of the value of the underlying
+     * function. In order for this penalty to be effective in rejecting this
+     * point during the optimization process, the penalty function value should
+     * be defined with care. This value is computed as:
+     * <pre>
+     *   penalty(point) = offset + &sum;<sub>i</sub>[scale[i] * &radic;|point[i]-boundary[i]|]
+     * </pre>
+     * where indices i correspond to all the components that violates their boundaries.
+     * </p>
+     * <p>
+     * So when attempting a function minimization, offset should be larger than
+     * the maximum expected value of the underlying function and scale components
+     * should all be positive. When attempting a function maximization, offset
+     * should be lesser than the minimum expected value of the underlying function
+     * and scale components should all be negative.
+     * minimization, and lesser than the minimum expected value of the underlying
+     * function when attempting maximization.
+     * </p>
+     * <p>
+     * These choices for the penalty function have two properties. First, all out
+     * of range points will return a function value that is worse than the value
+     * returned by any in range point. Second, the penalty is worse for large
+     * boundaries violation than for small violations, so the optimizer has an hint
+     * about the direction in which it should search for acceptable points.
+     * </p>
+     * @param bounded bounded function
+     * @param lower lower bounds for each element of the input parameters array
+     * (some elements may be set to {@code Double.NEGATIVE_INFINITY} for
+     * unbounded values)
+     * @param upper upper bounds for each element of the input parameters array
+     * (some elements may be set to {@code Double.POSITIVE_INFINITY} for
+     * unbounded values)
+     * @param offset base offset of the penalty function
+     * @param scale scale of the penalty function
+     * @exception DimensionMismatchException if lower bounds, upper bounds and
+     * scales are not consistent, either according to dimension or to bounadary
+     * values
+     */
+    public MultivariateFunctionPenaltyAdapter(final MultivariateFunction bounded,
+                                              final double[] lower, final double[] upper,
+                                              final double offset, final double[] scale) {
+
+        // safety checks
+        MathUtils.checkNotNull(lower);
+        MathUtils.checkNotNull(upper);
+        MathUtils.checkNotNull(scale);
+        if (lower.length != upper.length) {
+            throw new DimensionMismatchException(lower.length, upper.length);
+        }
+        if (lower.length != scale.length) {
+            throw new DimensionMismatchException(lower.length, scale.length);
+        }
+        for (int i = 0; i < lower.length; ++i) {
+            // note the following test is written in such a way it also fails for NaN
+            if (!(upper[i] >= lower[i])) {
+                throw new NumberIsTooSmallException(upper[i], lower[i], true);
+            }
+        }
+
+        this.bounded = bounded;
+        this.lower   = lower.clone();
+        this.upper   = upper.clone();
+        this.offset  = offset;
+        this.scale   = scale.clone();
+    }
+
+    /**
+     * Computes the underlying function value from an unbounded point.
+     * <p>
+     * This method simply returns the value of the underlying function
+     * if the unbounded point already fulfills the bounds, and compute
+     * a replacement value using the offset and scale if bounds are
+     * violated, without calling the function at all.
+     * </p>
+     * @param point unbounded point
+     * @return either underlying function value or penalty function value
+     */
+    public double value(double[] point) {
+
+        for (int i = 0; i < scale.length; ++i) {
+            if ((point[i] < lower[i]) || (point[i] > upper[i])) {
+                // bound violation starting at this component
+                double sum = 0;
+                for (int j = i; j < scale.length; ++j) {
+                    final double overshoot;
+                    if (point[j] < lower[j]) {
+                        overshoot = scale[j] * (lower[j] - point[j]);
+                    } else if (point[j] > upper[j]) {
+                        overshoot = scale[j] * (point[j] - upper[j]);
+                    } else {
+                        overshoot = 0;
+                    }
+                    sum += FastMath.sqrt(overshoot);
+                }
+                return offset + sum;
+            }
+        }
+
+        // all boundaries are fulfilled, we are in the expected
+        // domain of the underlying function
+        return bounded.value(point);
+    }
+}

Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optim/nonlinear/scalar/MultivariateFunctionPenaltyAdapter.java
------------------------------------------------------------------------------
    svn:eol-style = native

Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optim/nonlinear/scalar/MultivariateOptimizer.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optim/nonlinear/scalar/MultivariateOptimizer.java?rev=1420684&view=auto
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optim/nonlinear/scalar/MultivariateOptimizer.java (added)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optim/nonlinear/scalar/MultivariateOptimizer.java Wed Dec 12 14:10:38 2012
@@ -0,0 +1,119 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You 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 org.apache.commons.math3.optim.nonlinear.scalar;
+
+import org.apache.commons.math3.analysis.MultivariateFunction;
+import org.apache.commons.math3.optim.BaseMultivariateOptimizer;
+import org.apache.commons.math3.optim.OptimizationData;
+import org.apache.commons.math3.optim.GoalType;
+import org.apache.commons.math3.optim.ConvergenceChecker;
+import org.apache.commons.math3.optim.PointValuePair;
+import org.apache.commons.math3.optim.ObjectiveFunction;
+import org.apache.commons.math3.exception.TooManyEvaluationsException;
+
+/**
+ * Base class for a multivariate scalar function optimizer.
+ *
+ * @version $Id$
+ * @since 3.1
+ */
+public abstract class MultivariateOptimizer
+    extends BaseMultivariateOptimizer<PointValuePair> {
+    /** Objective function. */
+    private MultivariateFunction function;
+    /** Type of optimization. */
+    private GoalType goal;
+
+    /**
+     * @param checker Convergence checker.
+     */
+    protected MultivariateOptimizer(ConvergenceChecker<PointValuePair> checker) {
+        super(checker);
+    }
+
+    /**
+     * {@inheritDoc}
+     *
+     * @param optData Optimization data. The following data will be looked for:
+     * <ul>
+     *  <li>{@link org.apache.commons.math3.optim.MaxEval}</li>
+     *  <li>{@link org.apache.commons.math3.optim.InitialGuess}</li>
+     *  <li>{@link org.apache.commons.math3.optim.SimpleBounds}</li>
+     *  <li>{@link ObjectiveFunction}</li>
+     *  <li>{@link GoalType}</li>
+     * </ul>
+     * @return {@inheritDoc}
+     * @throws TooManyEvaluationsException if the maximal number of
+     * evaluations is exceeded.
+     */
+    @Override
+    public PointValuePair optimize(OptimizationData... optData)
+        throws TooManyEvaluationsException {
+         // Retrieve settings.
+        parseOptimizationData(optData);
+        // Set up base class and perform computation.
+        return super.optimize(optData);
+    }
+
+    /**
+     * Scans the list of (required and optional) optimization data that
+     * characterize the problem.
+     *
+     * @param optData Optimization data.
+     * The following data will be looked for:
+     * <ul>
+     *  <li>{@link ObjectiveFunction}</li>
+     *  <li>{@link GoalType}</li>
+     * </ul>
+     */
+    private void parseOptimizationData(OptimizationData... optData) {
+        // The existing values (as set by the previous call) are reused if
+        // not provided in the argument list.
+        for (OptimizationData data : optData) {
+            if (data instanceof GoalType) {
+                goal = (GoalType) data;
+                continue;
+            }
+            if  (data instanceof ObjectiveFunction) {
+                function = ((ObjectiveFunction) data).getObjectiveFunction();
+                continue;
+            }
+        }
+    }
+
+    /**
+     * @return the optimization type.
+     */
+    public GoalType getGoalType() {
+        return goal;
+    }
+
+    /**
+     * Computes the objective function value.
+     * This method <em>must</em> be called by subclasses to enforce the
+     * evaluation counter limit.
+     *
+     * @param params Point at which the objective function must be evaluated.
+     * @return the objective function value at the specified point.
+     * @throws TooManyEvaluationsException if the maximal number of
+     * evaluations is exceeded.
+     */
+    protected double computeObjectiveValue(double[] params) {
+        super.incrementEvaluationCount();
+        return function.value(params);
+    }
+}

Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optim/nonlinear/scalar/MultivariateOptimizer.java
------------------------------------------------------------------------------
    svn:eol-style = native

Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optim/nonlinear/scalar/ObjectiveFunctionGradient.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optim/nonlinear/scalar/ObjectiveFunctionGradient.java?rev=1420684&view=auto
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optim/nonlinear/scalar/ObjectiveFunctionGradient.java (added)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optim/nonlinear/scalar/ObjectiveFunctionGradient.java Wed Dec 12 14:10:38 2012
@@ -0,0 +1,47 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You 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 org.apache.commons.math3.optim.nonlinear.scalar;
+
+import org.apache.commons.math3.analysis.MultivariateVectorFunction;
+import org.apache.commons.math3.optim.OptimizationData;
+
+/**
+ * Gradient of the scalar function to be optimized.
+ *
+ * @version $Id$
+ * @since 3.1
+ */
+public class ObjectiveFunctionGradient implements OptimizationData {
+    /** Function to be optimized. */
+    private final MultivariateVectorFunction gradient;
+
+    /**
+     * @param g Gradient of the function to be optimized.
+     */
+    public ObjectiveFunctionGradient(MultivariateVectorFunction g) {
+        gradient = g;
+    }
+
+    /**
+     * Gets the gradient of the function to be optimized.
+     *
+     * @return the objective function gradient.
+     */
+    public MultivariateVectorFunction getObjectiveFunctionGradient() {
+        return gradient;
+    }
+}

Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optim/nonlinear/scalar/ObjectiveFunctionGradient.java
------------------------------------------------------------------------------
    svn:eol-style = native

Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optim/nonlinear/scalar/gradient/NonLinearConjugateGradientOptimizer.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optim/nonlinear/scalar/gradient/NonLinearConjugateGradientOptimizer.java?rev=1420684&view=auto
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optim/nonlinear/scalar/gradient/NonLinearConjugateGradientOptimizer.java (added)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optim/nonlinear/scalar/gradient/NonLinearConjugateGradientOptimizer.java Wed Dec 12 14:10:38 2012
@@ -0,0 +1,393 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You 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 org.apache.commons.math3.optim.nonlinear.scalar.gradient;
+
+import org.apache.commons.math3.analysis.UnivariateFunction;
+import org.apache.commons.math3.analysis.solvers.BrentSolver;
+import org.apache.commons.math3.analysis.solvers.UnivariateSolver;
+import org.apache.commons.math3.exception.MathInternalError;
+import org.apache.commons.math3.exception.MathIllegalStateException;
+import org.apache.commons.math3.exception.TooManyEvaluationsException;
+import org.apache.commons.math3.exception.util.LocalizedFormats;
+import org.apache.commons.math3.optim.OptimizationData;
+import org.apache.commons.math3.optim.GoalType;
+import org.apache.commons.math3.optim.PointValuePair;
+import org.apache.commons.math3.optim.ConvergenceChecker;
+import org.apache.commons.math3.optim.nonlinear.scalar.GradientMultivariateOptimizer;
+import org.apache.commons.math3.util.FastMath;
+
+/**
+ * Non-linear conjugate gradient optimizer.
+ * <p>
+ * This class supports both the Fletcher-Reeves and the Polak-Ribière
+ * update formulas for the conjugate search directions.
+ * It also supports optional preconditioning.
+ * </p>
+ *
+ * @version $Id: NonLinearConjugateGradientOptimizer.java 1416643 2012-12-03 19:37:14Z tn $
+ * @since 2.0
+ */
+public class NonLinearConjugateGradientOptimizer
+    extends GradientMultivariateOptimizer {
+    /** Update formula for the beta parameter. */
+    private final Formula updateFormula;
+    /** Preconditioner (may be null). */
+    private final Preconditioner preconditioner;
+    /** solver to use in the line search (may be null). */
+    private final UnivariateSolver solver;
+    /** Initial step used to bracket the optimum in line search. */
+    private double initialStep = 1;
+
+    /**
+     * Constructor with default {@link BrentSolver line search solver} and
+     * {@link IdentityPreconditioner preconditioner}.
+     *
+     * @param updateFormula formula to use for updating the &beta; parameter,
+     * must be one of {@link Formula#FLETCHER_REEVES} or
+     * {@link Formula#POLAK_RIBIERE}.
+     * @param checker Convergence checker.
+     */
+    public NonLinearConjugateGradientOptimizer(final Formula updateFormula,
+                                               ConvergenceChecker<PointValuePair> checker) {
+        this(updateFormula,
+             checker,
+             new BrentSolver(),
+             new IdentityPreconditioner());
+    }
+
+    /**
+     * Available choices of update formulas for the updating the parameter
+     * that is used to compute the successive conjugate search directions.
+     * For non-linear conjugate gradients, there are
+     * two formulas:
+     * <ul>
+     *   <li>Fletcher-Reeves formula</li>
+     *   <li>Polak-Ribière formula</li>
+     * </ul>
+     *
+     * On the one hand, the Fletcher-Reeves formula is guaranteed to converge
+     * if the start point is close enough of the optimum whether the
+     * Polak-Ribière formula may not converge in rare cases. On the
+     * other hand, the Polak-Ribière formula is often faster when it
+     * does converge. Polak-Ribière is often used.
+     *
+     * @since 2.0
+     */
+    public static enum Formula {
+        /** Fletcher-Reeves formula. */
+        FLETCHER_REEVES,
+        /** Polak-Ribière formula. */
+        POLAK_RIBIERE
+    }
+
+    /**
+     * The initial step is a factor with respect to the search direction
+     * (which itself is roughly related to the gradient of the function).
+     * <br/>
+     * It is used to find an interval that brackets the optimum in line
+     * search.
+     *
+     * @since 3.1
+     */
+    public static class BracketingStep implements OptimizationData {
+        /** Initial step. */
+        private final double initialStep;
+
+        /**
+         * @param step Initial step for the bracket search.
+         */
+        public BracketingStep(double step) {
+            initialStep = step;
+        }
+
+        /**
+         * Gets the initial step.
+         *
+         * @return the initial step.
+         */
+        public double getBracketingStep() {
+            return initialStep;
+        }
+    }
+
+    /**
+     * Constructor with default {@link IdentityPreconditioner preconditioner}.
+     *
+     * @param updateFormula formula to use for updating the &beta; parameter,
+     * must be one of {@link Formula#FLETCHER_REEVES} or
+     * {@link Formula#POLAK_RIBIERE}.
+     * @param checker Convergence checker.
+     * @param lineSearchSolver Solver to use during line search.
+     */
+    public NonLinearConjugateGradientOptimizer(final Formula updateFormula,
+                                               ConvergenceChecker<PointValuePair> checker,
+                                               final UnivariateSolver lineSearchSolver) {
+        this(updateFormula,
+             checker,
+             lineSearchSolver,
+             new IdentityPreconditioner());
+    }
+
+    /**
+     * @param updateFormula formula to use for updating the &beta; parameter,
+     * must be one of {@link Formula#FLETCHER_REEVES} or
+     * {@link Formula#POLAK_RIBIERE}.
+     * @param checker Convergence checker.
+     * @param lineSearchSolver Solver to use during line search.
+     * @param preconditioner Preconditioner.
+     */
+    public NonLinearConjugateGradientOptimizer(final Formula updateFormula,
+                                               ConvergenceChecker<PointValuePair> checker,
+                                               final UnivariateSolver lineSearchSolver,
+                                               final Preconditioner preconditioner) {
+        super(checker);
+
+        this.updateFormula = updateFormula;
+        solver = lineSearchSolver;
+        this.preconditioner = preconditioner;
+        initialStep = 1;
+    }
+
+    /**
+     * {@inheritDoc}
+     *
+     * @param optData Optimization data.
+     * The following data will be looked for:
+     * <ul>
+     *  <li>{@link org.apache.commons.math3.optim.MaxEval}</li>
+     *  <li>{@link org.apache.commons.math3.optim.InitialGuess}</li>
+     *  <li>{@link org.apache.commons.math3.optim.SimpleBounds}</li>
+     *  <li>{@link org.apache.commons.math3.optim.GoalType}</li>
+     *  <li>{@link org.apache.commons.math3.optim.ObjectiveFunction}</li>
+     *  <li>{@link org.apache.commons.math3.optim.nonlinear.scalar.ObjectiveFunctionGradient}</li>
+     *  <li>{@link BracketingStep}</li>
+     * </ul>
+     * @return {@inheritDoc}
+     * @throws TooManyEvaluationsException if the maximal number of
+     * evaluations (of the objective function) is exceeded.
+     */
+    @Override
+    public PointValuePair optimize(OptimizationData... optData)
+        throws TooManyEvaluationsException {
+         // Retrieve settings.
+        parseOptimizationData(optData);
+        // Set up base class and perform computation.
+        return super.optimize(optData);
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    protected PointValuePair doOptimize() {
+        final ConvergenceChecker<PointValuePair> checker = getConvergenceChecker();
+        final double[] point = getStartPoint();
+        final GoalType goal = getGoalType();
+        final int n = point.length;
+        double[] r = computeObjectiveGradient(point);
+        if (goal == GoalType.MINIMIZE) {
+            for (int i = 0; i < n; i++) {
+                r[i] = -r[i];
+            }
+        }
+
+        // Initial search direction.
+        double[] steepestDescent = preconditioner.precondition(point, r);
+        double[] searchDirection = steepestDescent.clone();
+
+        double delta = 0;
+        for (int i = 0; i < n; ++i) {
+            delta += r[i] * searchDirection[i];
+        }
+
+        PointValuePair current = null;
+        int iter = 0;
+        int maxEval = getMaxEvaluations();
+        while (true) {
+            ++iter;
+
+            final double objective = computeObjectiveValue(point);
+            PointValuePair previous = current;
+            current = new PointValuePair(point, objective);
+            if (previous != null) {
+                if (checker.converged(iter, previous, current)) {
+                    // We have found an optimum.
+                    return current;
+                }
+            }
+
+            // Find the optimal step in the search direction.
+            final UnivariateFunction lsf = new LineSearchFunction(point, searchDirection);
+            final double uB = findUpperBound(lsf, 0, initialStep);
+            // XXX Last parameters is set to a value close to zero in order to
+            // work around the divergence problem in the "testCircleFitting"
+            // unit test (see MATH-439).
+            final double step = solver.solve(maxEval, lsf, 0, uB, 1e-15);
+            maxEval -= solver.getEvaluations(); // Subtract used up evaluations.
+
+            // Validate new point.
+            for (int i = 0; i < point.length; ++i) {
+                point[i] += step * searchDirection[i];
+            }
+
+            r = computeObjectiveGradient(point);
+            if (goal == GoalType.MINIMIZE) {
+                for (int i = 0; i < n; ++i) {
+                    r[i] = -r[i];
+                }
+            }
+
+            // Compute beta.
+            final double deltaOld = delta;
+            final double[] newSteepestDescent = preconditioner.precondition(point, r);
+            delta = 0;
+            for (int i = 0; i < n; ++i) {
+                delta += r[i] * newSteepestDescent[i];
+            }
+
+            final double beta;
+            switch (updateFormula) {
+            case FLETCHER_REEVES:
+                beta = delta / deltaOld;
+                break;
+            case POLAK_RIBIERE:
+                double deltaMid = 0;
+                for (int i = 0; i < r.length; ++i) {
+                    deltaMid += r[i] * steepestDescent[i];
+                }
+                beta = (delta - deltaMid) / deltaOld;
+                break;
+            default:
+                // Should never happen.
+                throw new MathInternalError();
+            }
+            steepestDescent = newSteepestDescent;
+
+            // Compute conjugate search direction.
+            if (iter % n == 0 ||
+                beta < 0) {
+                // Break conjugation: reset search direction.
+                searchDirection = steepestDescent.clone();
+            } else {
+                // Compute new conjugate search direction.
+                for (int i = 0; i < n; ++i) {
+                    searchDirection[i] = steepestDescent[i] + beta * searchDirection[i];
+                }
+            }
+        }
+    }
+
+    /**
+     * Scans the list of (required and optional) optimization data that
+     * characterize the problem.
+     *
+     * @param optData Optimization data.
+     * The following data will be looked for:
+     * <ul>
+     *  <li>{@link InitialStep}</li>
+     * </ul>
+     */
+    private void parseOptimizationData(OptimizationData... optData) {
+        // The existing values (as set by the previous call) are reused if
+        // not provided in the argument list.
+        for (OptimizationData data : optData) {
+            if  (data instanceof BracketingStep) {
+                initialStep = ((BracketingStep) data).getBracketingStep();
+                // If more data must be parsed, this statement _must_ be
+                // changed to "continue".
+                break;
+            }
+        }
+    }
+
+    /**
+     * Finds the upper bound b ensuring bracketing of a root between a and b.
+     *
+     * @param f function whose root must be bracketed.
+     * @param a lower bound of the interval.
+     * @param h initial step to try.
+     * @return b such that f(a) and f(b) have opposite signs.
+     * @throws MathIllegalStateException if no bracket can be found.
+     */
+    private double findUpperBound(final UnivariateFunction f,
+                                  final double a, final double h) {
+        final double yA = f.value(a);
+        double yB = yA;
+        for (double step = h; step < Double.MAX_VALUE; step *= FastMath.max(2, yA / yB)) {
+            final double b = a + step;
+            yB = f.value(b);
+            if (yA * yB <= 0) {
+                return b;
+            }
+        }
+        throw new MathIllegalStateException(LocalizedFormats.UNABLE_TO_BRACKET_OPTIMUM_IN_LINE_SEARCH);
+    }
+
+    /** Default identity preconditioner. */
+    public static class IdentityPreconditioner implements Preconditioner {
+        /** {@inheritDoc} */
+        public double[] precondition(double[] variables, double[] r) {
+            return r.clone();
+        }
+    }
+
+    /**
+     * Internal class for line search.
+     * <p>
+     * The function represented by this class is the dot product of
+     * the objective function gradient and the search direction. Its
+     * value is zero when the gradient is orthogonal to the search
+     * direction, i.e. when the objective function value is a local
+     * extremum along the search direction.
+     * </p>
+     */
+    private class LineSearchFunction implements UnivariateFunction {
+        /** Current point. */
+        private final double[] currentPoint;
+        /** Search direction. */
+        private final double[] searchDirection;
+
+        /**
+         * @param point Current point.
+         * @param direction Search direction.
+         */
+        public LineSearchFunction(double[] point,
+                                  double[] direction) {
+            currentPoint = point.clone();
+            searchDirection = direction.clone();
+        }
+
+        /** {@inheritDoc} */
+        public double value(double x) {
+            // current point in the search direction
+            final double[] shiftedPoint = currentPoint.clone();
+            for (int i = 0; i < shiftedPoint.length; ++i) {
+                shiftedPoint[i] += x * searchDirection[i];
+            }
+
+            // gradient of the objective function
+            final double[] gradient = computeObjectiveGradient(shiftedPoint);
+
+            // dot product with the search direction
+            double dotProduct = 0;
+            for (int i = 0; i < gradient.length; ++i) {
+                dotProduct += gradient[i] * searchDirection[i];
+            }
+
+            return dotProduct;
+        }
+    }
+}

Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optim/nonlinear/scalar/gradient/NonLinearConjugateGradientOptimizer.java
------------------------------------------------------------------------------
    svn:eol-style = native

Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optim/nonlinear/scalar/gradient/Preconditioner.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optim/nonlinear/scalar/gradient/Preconditioner.java?rev=1420684&view=auto
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optim/nonlinear/scalar/gradient/Preconditioner.java (added)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optim/nonlinear/scalar/gradient/Preconditioner.java Wed Dec 12 14:10:38 2012
@@ -0,0 +1,45 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You 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 org.apache.commons.math3.optim.nonlinear.scalar.gradient;
+
+/**
+ * This interface represents a preconditioner for differentiable scalar
+ * objective function optimizers.
+ * @version $Id: Preconditioner.java 1416643 2012-12-03 19:37:14Z tn $
+ * @since 2.0
+ */
+public interface Preconditioner {
+    /**
+     * Precondition a search direction.
+     * <p>
+     * The returned preconditioned search direction must be computed fast or
+     * the algorithm performances will drop drastically. A classical approach
+     * is to compute only the diagonal elements of the hessian and to divide
+     * the raw search direction by these elements if they are all positive.
+     * If at least one of them is negative, it is safer to return a clone of
+     * the raw search direction as if the hessian was the identity matrix. The
+     * rationale for this simplified choice is that a negative diagonal element
+     * means the current point is far from the optimum and preconditioning will
+     * not be efficient anyway in this case.
+     * </p>
+     * @param point current point at which the search direction was computed
+     * @param r raw search direction (i.e. opposite of the gradient)
+     * @return approximation of H<sup>-1</sup>r where H is the objective function hessian
+     */
+    double[] precondition(double[] point, double[] r);
+}

Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optim/nonlinear/scalar/gradient/Preconditioner.java
------------------------------------------------------------------------------
    svn:eol-style = native

Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optim/nonlinear/scalar/gradient/package-info.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optim/nonlinear/scalar/gradient/package-info.java?rev=1420684&view=auto
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optim/nonlinear/scalar/gradient/package-info.java (added)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optim/nonlinear/scalar/gradient/package-info.java Wed Dec 12 14:10:38 2012
@@ -0,0 +1,21 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You 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 org.apache.commons.math3.optim.nonlinear.scalar.gradient;
+
+/**
+ * This package provides optimization algorithms that require derivatives.
+ */

Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optim/nonlinear/scalar/gradient/package-info.java
------------------------------------------------------------------------------
    svn:eol-style = native

Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optim/nonlinear/scalar/noderiv/AbstractSimplex.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optim/nonlinear/scalar/noderiv/AbstractSimplex.java?rev=1420684&view=auto
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optim/nonlinear/scalar/noderiv/AbstractSimplex.java (added)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optim/nonlinear/scalar/noderiv/AbstractSimplex.java Wed Dec 12 14:10:38 2012
@@ -0,0 +1,346 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You 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 org.apache.commons.math3.optim.nonlinear.scalar.noderiv;
+
+import java.util.Arrays;
+import java.util.Comparator;
+
+import org.apache.commons.math3.analysis.MultivariateFunction;
+import org.apache.commons.math3.exception.NotStrictlyPositiveException;
+import org.apache.commons.math3.exception.DimensionMismatchException;
+import org.apache.commons.math3.exception.ZeroException;
+import org.apache.commons.math3.exception.OutOfRangeException;
+import org.apache.commons.math3.exception.NullArgumentException;
+import org.apache.commons.math3.exception.MathIllegalArgumentException;
+import org.apache.commons.math3.exception.util.LocalizedFormats;
+import org.apache.commons.math3.optim.PointValuePair;
+import org.apache.commons.math3.optim.OptimizationData;
+
+/**
+ * This class implements the simplex concept.
+ * It is intended to be used in conjunction with {@link SimplexOptimizer}.
+ * <br/>
+ * The initial configuration of the simplex is set by the constructors
+ * {@link #AbstractSimplex(double[])} or {@link #AbstractSimplex(double[][])}.
+ * The other {@link #AbstractSimplex(int) constructor} will set all steps
+ * to 1, thus building a default configuration from a unit hypercube.
+ * <br/>
+ * Users <em>must</em> call the {@link #build(double[]) build} method in order
+ * to create the data structure that will be acted on by the other methods of
+ * this class.
+ *
+ * @see SimplexOptimizer
+ * @version $Id: AbstractSimplex.java 1397759 2012-10-13 01:12:58Z erans $
+ * @since 3.0
+ */
+public abstract class AbstractSimplex implements OptimizationData {
+    /** Simplex. */
+    private PointValuePair[] simplex;
+    /** Start simplex configuration. */
+    private double[][] startConfiguration;
+    /** Simplex dimension (must be equal to {@code simplex.length - 1}). */
+    private final int dimension;
+
+    /**
+     * Build a unit hypercube simplex.
+     *
+     * @param n Dimension of the simplex.
+     */
+    protected AbstractSimplex(int n) {
+        this(n, 1d);
+    }
+
+    /**
+     * Build a hypercube simplex with the given side length.
+     *
+     * @param n Dimension of the simplex.
+     * @param sideLength Length of the sides of the hypercube.
+     */
+    protected AbstractSimplex(int n,
+                              double sideLength) {
+        this(createHypercubeSteps(n, sideLength));
+    }
+
+    /**
+     * The start configuration for simplex is built from a box parallel to
+     * the canonical axes of the space. The simplex is the subset of vertices
+     * of a box parallel to the canonical axes. It is built as the path followed
+     * while traveling from one vertex of the box to the diagonally opposite
+     * vertex moving only along the box edges. The first vertex of the box will
+     * be located at the start point of the optimization.
+     * As an example, in dimension 3 a simplex has 4 vertices. Setting the
+     * steps to (1, 10, 2) and the start point to (1, 1, 1) would imply the
+     * start simplex would be: { (1, 1, 1), (2, 1, 1), (2, 11, 1), (2, 11, 3) }.
+     * The first vertex would be set to the start point at (1, 1, 1) and the
+     * last vertex would be set to the diagonally opposite vertex at (2, 11, 3).
+     *
+     * @param steps Steps along the canonical axes representing box edges. They
+     * may be negative but not zero.
+     * @throws NullArgumentException if {@code steps} is {@code null}.
+     * @throws ZeroException if one of the steps is zero.
+     */
+    protected AbstractSimplex(final double[] steps) {
+        if (steps == null) {
+            throw new NullArgumentException();
+        }
+        if (steps.length == 0) {
+            throw new ZeroException();
+        }
+        dimension = steps.length;
+
+        // Only the relative position of the n final vertices with respect
+        // to the first one are stored.
+        startConfiguration = new double[dimension][dimension];
+        for (int i = 0; i < dimension; i++) {
+            final double[] vertexI = startConfiguration[i];
+            for (int j = 0; j < i + 1; j++) {
+                if (steps[j] == 0) {
+                    throw new ZeroException(LocalizedFormats.EQUAL_VERTICES_IN_SIMPLEX);
+                }
+                System.arraycopy(steps, 0, vertexI, 0, j + 1);
+            }
+        }
+    }
+
+    /**
+     * The real initial simplex will be set up by moving the reference
+     * simplex such that its first point is located at the start point of the
+     * optimization.
+     *
+     * @param referenceSimplex Reference simplex.
+     * @throws NotStrictlyPositiveException if the reference simplex does not
+     * contain at least one point.
+     * @throws DimensionMismatchException if there is a dimension mismatch
+     * in the reference simplex.
+     * @throws IllegalArgumentException if one of its vertices is duplicated.
+     */
+    protected AbstractSimplex(final double[][] referenceSimplex) {
+        if (referenceSimplex.length <= 0) {
+            throw new NotStrictlyPositiveException(LocalizedFormats.SIMPLEX_NEED_ONE_POINT,
+                                                   referenceSimplex.length);
+        }
+        dimension = referenceSimplex.length - 1;
+
+        // Only the relative position of the n final vertices with respect
+        // to the first one are stored.
+        startConfiguration = new double[dimension][dimension];
+        final double[] ref0 = referenceSimplex[0];
+
+        // Loop over vertices.
+        for (int i = 0; i < referenceSimplex.length; i++) {
+            final double[] refI = referenceSimplex[i];
+
+            // Safety checks.
+            if (refI.length != dimension) {
+                throw new DimensionMismatchException(refI.length, dimension);
+            }
+            for (int j = 0; j < i; j++) {
+                final double[] refJ = referenceSimplex[j];
+                boolean allEquals = true;
+                for (int k = 0; k < dimension; k++) {
+                    if (refI[k] != refJ[k]) {
+                        allEquals = false;
+                        break;
+                    }
+                }
+                if (allEquals) {
+                    throw new MathIllegalArgumentException(LocalizedFormats.EQUAL_VERTICES_IN_SIMPLEX,
+                                                           i, j);
+                }
+            }
+
+            // Store vertex i position relative to vertex 0 position.
+            if (i > 0) {
+                final double[] confI = startConfiguration[i - 1];
+                for (int k = 0; k < dimension; k++) {
+                    confI[k] = refI[k] - ref0[k];
+                }
+            }
+        }
+    }
+
+    /**
+     * Get simplex dimension.
+     *
+     * @return the dimension of the simplex.
+     */
+    public int getDimension() {
+        return dimension;
+    }
+
+    /**
+     * Get simplex size.
+     * After calling the {@link #build(double[]) build} method, this method will
+     * will be equivalent to {@code getDimension() + 1}.
+     *
+     * @return the size of the simplex.
+     */
+    public int getSize() {
+        return simplex.length;
+    }
+
+    /**
+     * Compute the next simplex of the algorithm.
+     *
+     * @param evaluationFunction Evaluation function.
+     * @param comparator Comparator to use to sort simplex vertices from best
+     * to worst.
+     * @throws org.apache.commons.math3.exception.TooManyEvaluationsException
+     * if the algorithm fails to converge.
+     */
+    public abstract void iterate(final MultivariateFunction evaluationFunction,
+                                 final Comparator<PointValuePair> comparator);
+
+    /**
+     * Build an initial simplex.
+     *
+     * @param startPoint First point of the simplex.
+     * @throws DimensionMismatchException if the start point does not match
+     * simplex dimension.
+     */
+    public void build(final double[] startPoint) {
+        if (dimension != startPoint.length) {
+            throw new DimensionMismatchException(dimension, startPoint.length);
+        }
+
+        // Set first vertex.
+        simplex = new PointValuePair[dimension + 1];
+        simplex[0] = new PointValuePair(startPoint, Double.NaN);
+
+        // Set remaining vertices.
+        for (int i = 0; i < dimension; i++) {
+            final double[] confI = startConfiguration[i];
+            final double[] vertexI = new double[dimension];
+            for (int k = 0; k < dimension; k++) {
+                vertexI[k] = startPoint[k] + confI[k];
+            }
+            simplex[i + 1] = new PointValuePair(vertexI, Double.NaN);
+        }
+    }
+
+    /**
+     * Evaluate all the non-evaluated points of the simplex.
+     *
+     * @param evaluationFunction Evaluation function.
+     * @param comparator Comparator to use to sort simplex vertices from best to worst.
+     * @throws org.apache.commons.math3.exception.TooManyEvaluationsException
+     * if the maximal number of evaluations is exceeded.
+     */
+    public void evaluate(final MultivariateFunction evaluationFunction,
+                         final Comparator<PointValuePair> comparator) {
+        // Evaluate the objective function at all non-evaluated simplex points.
+        for (int i = 0; i < simplex.length; i++) {
+            final PointValuePair vertex = simplex[i];
+            final double[] point = vertex.getPointRef();
+            if (Double.isNaN(vertex.getValue())) {
+                simplex[i] = new PointValuePair(point, evaluationFunction.value(point), false);
+            }
+        }
+
+        // Sort the simplex from best to worst.
+        Arrays.sort(simplex, comparator);
+    }
+
+    /**
+     * Replace the worst point of the simplex by a new point.
+     *
+     * @param pointValuePair Point to insert.
+     * @param comparator Comparator to use for sorting the simplex vertices
+     * from best to worst.
+     */
+    protected void replaceWorstPoint(PointValuePair pointValuePair,
+                                     final Comparator<PointValuePair> comparator) {
+        for (int i = 0; i < dimension; i++) {
+            if (comparator.compare(simplex[i], pointValuePair) > 0) {
+                PointValuePair tmp = simplex[i];
+                simplex[i] = pointValuePair;
+                pointValuePair = tmp;
+            }
+        }
+        simplex[dimension] = pointValuePair;
+    }
+
+    /**
+     * Get the points of the simplex.
+     *
+     * @return all the simplex points.
+     */
+    public PointValuePair[] getPoints() {
+        final PointValuePair[] copy = new PointValuePair[simplex.length];
+        System.arraycopy(simplex, 0, copy, 0, simplex.length);
+        return copy;
+    }
+
+    /**
+     * Get the simplex point stored at the requested {@code index}.
+     *
+     * @param index Location.
+     * @return the point at location {@code index}.
+     */
+    public PointValuePair getPoint(int index) {
+        if (index < 0 ||
+            index >= simplex.length) {
+            throw new OutOfRangeException(index, 0, simplex.length - 1);
+        }
+        return simplex[index];
+    }
+
+    /**
+     * Store a new point at location {@code index}.
+     * Note that no deep-copy of {@code point} is performed.
+     *
+     * @param index Location.
+     * @param point New value.
+     */
+    protected void setPoint(int index, PointValuePair point) {
+        if (index < 0 ||
+            index >= simplex.length) {
+            throw new OutOfRangeException(index, 0, simplex.length - 1);
+        }
+        simplex[index] = point;
+    }
+
+    /**
+     * Replace all points.
+     * Note that no deep-copy of {@code points} is performed.
+     *
+     * @param points New Points.
+     */
+    protected void setPoints(PointValuePair[] points) {
+        if (points.length != simplex.length) {
+            throw new DimensionMismatchException(points.length, simplex.length);
+        }
+        simplex = points;
+    }
+
+    /**
+     * Create steps for a unit hypercube.
+     *
+     * @param n Dimension of the hypercube.
+     * @param sideLength Length of the sides of the hypercube.
+     * @return the steps.
+     */
+    private static double[] createHypercubeSteps(int n,
+                                                 double sideLength) {
+        final double[] steps = new double[n];
+        for (int i = 0; i < n; i++) {
+            steps[i] = sideLength;
+        }
+        return steps;
+    }
+}

Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optim/nonlinear/scalar/noderiv/AbstractSimplex.java
------------------------------------------------------------------------------
    svn:eol-style = native



Mime
View raw message