001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017package org.apache.commons.math4.legacy.optim.nonlinear.scalar.noderiv;
018
019import java.util.Comparator;
020import java.util.List;
021import java.util.ArrayList;
022import java.util.Collections;
023import java.util.function.UnaryOperator;
024import java.util.function.DoublePredicate;
025import org.apache.commons.rng.UniformRandomProvider;
026import org.apache.commons.rng.simple.RandomSource;
027import org.apache.commons.rng.sampling.distribution.ContinuousSampler;
028import org.apache.commons.rng.sampling.distribution.ContinuousUniformSampler;
029import org.apache.commons.math4.legacy.analysis.MultivariateFunction;
030import org.apache.commons.math4.legacy.optim.PointValuePair;
031
032/**
033 * DSSA algorithm.
034 *
035 * Described in
036 * <blockquote>
037 *  <em>Abdel-Rahman Hedar and Masao Fukushima (2002)</em>,
038 *  <b>
039 *   Hybrid simulated annealing and direct search method
040 *   for nonlinear unconstrained global optimization
041 *  </b>,
042 *  Optimization Methods and Software, 17:5, 891-912,
043 *  DOI: 10.1080/1055678021000030084
044 * </blockquote>
045 *
046 * <p>
047 * A note about the {@link #HedarFukushimaTransform(double) "shrink" factor}:
048 * Per DSSA's description, the simplex must keep its size during the simulated
049 * annealing (SA) phase to avoid premature convergence.  This assumes that the
050 * best candidates from the SA phase will each subsequently serve as starting
051 * point for another optimization to hone in on the local optimum.
052 * Values lower than 1 and no subsequent "best list" search correspond to the
053 * "SSA" algorithm in the above paper.
054 */
055public class HedarFukushimaTransform
056    implements Simplex.TransformFactory {
057    /** Shrinkage coefficient. */
058    private final double sigma;
059    /** Sampler for reflection coefficient. */
060    private final ContinuousSampler alphaSampler;
061    /** No shrink indicator. */
062    private final boolean noShrink;
063
064    /**
065     * @param sigma Shrink factor.
066     * @param rng Random generator.
067     * @throws IllegalArgumentException if {@code sigma <= 0} or
068     * {@code sigma > 1}.
069     */
070    public HedarFukushimaTransform(double sigma,
071                                   UniformRandomProvider rng) {
072        if (sigma <= 0 ||
073            sigma > 1) {
074            throw new IllegalArgumentException("Shrink factor out of range: " +
075                                               sigma);
076        }
077
078        this.sigma = sigma;
079        alphaSampler = ContinuousUniformSampler.of(rng, 0.9, 1.1);
080        noShrink = sigma == 1d;
081    }
082
083    /**
084     * @param sigma Shrink factor.
085     * @throws IllegalArgumentException if {@code sigma <= 0} or
086     * {@code sigma > 1}.
087     */
088    public HedarFukushimaTransform(double sigma) {
089        this(sigma, RandomSource.KISS.create());
090    }
091
092    /**
093     * Disable shrinking of the simplex (as mandated by DSSA).
094     */
095    public HedarFukushimaTransform() {
096        this(1d);
097    }
098
099    /** {@inheritDoc} */
100    @Override
101    public UnaryOperator<Simplex> create(final MultivariateFunction evaluationFunction,
102                                         final Comparator<PointValuePair> comparator,
103                                         final DoublePredicate saAcceptance) {
104        if (saAcceptance == null) {
105            throw new IllegalArgumentException("Missing SA acceptance test");
106        }
107
108        return original -> transform(original,
109                                     saAcceptance,
110                                     evaluationFunction,
111                                     comparator);
112    }
113
114    /**
115     * Simulated annealing step (at fixed temperature).
116     *
117     * @param original Current simplex.  Points must be sorted from best to worst.
118     * @param sa Simulated annealing acceptance test.
119     * @param eval Evaluation function.
120     * @param comp Objective function comparator.
121     * @return a new instance.
122     */
123    private Simplex transform(Simplex original,
124                              DoublePredicate sa,
125                              MultivariateFunction eval,
126                              Comparator<PointValuePair> comp) {
127        final int size = original.getSize();
128        // Current best point.
129        final PointValuePair best = original.get(0);
130        final double bestValue = best.getValue();
131
132        for (int k = 1; k < size; k++) {
133            // Perform reflections of the "k" worst points.
134            final List<PointValuePair> reflected = reflectPoints(original, k, eval);
135            Collections.sort(reflected, comp);
136
137            // Check whether the best of the reflected points is better than the
138            // current overall best.
139            final PointValuePair candidate = reflected.get(0);
140            final boolean candidateIsBetter = comp.compare(candidate, best) < 0;
141            final boolean candidateIsAccepted = candidateIsBetter ||
142                sa.test(candidate.getValue() - bestValue);
143
144            if (candidateIsAccepted) {
145                // Replace worst points with the reflected points.
146                return original.replaceLast(reflected);
147            }
148        }
149
150        // No direction provided a better point.
151        return noShrink ?
152            original :
153            original.shrink(sigma, eval);
154    }
155
156    /**
157     * @param simplex Current simplex (whose points must be sorted, from best
158     * to worst).
159     * @param nPoints Number of points to reflect.
160     * The {@code nPoints} worst points will be reflected through the centroid
161     * of the {@code n + 1 - nPoints} best points.
162     * @param eval Evaluation function.
163     * @return the (unsorted) list of reflected points.
164     * @throws IllegalArgumentException if {@code nPoints < 1} or
165     * {@code nPoints > n}.
166     */
167    private List<PointValuePair> reflectPoints(Simplex simplex,
168                                               int nPoints,
169                                               MultivariateFunction eval) {
170        final int size = simplex.getSize();
171        if (nPoints < 1 ||
172            nPoints >= size) {
173            throw new IllegalArgumentException("Out of range: " + nPoints);
174        }
175
176        final int nCentroid = size - nPoints;
177        final List<PointValuePair> centroidList = simplex.asList(0, nCentroid);
178        final List<PointValuePair> reflectList = simplex.asList(nCentroid, size);
179
180        final double[] centroid = Simplex.centroid(centroidList);
181
182        final List<PointValuePair> reflected = new ArrayList<>(nPoints);
183        for (int i = 0; i < reflectList.size(); i++) {
184            reflected.add(newReflectedPoint(reflectList.get(i),
185                                            centroid,
186                                            eval));
187        }
188
189        return reflected;
190    }
191
192    /**
193     * @param point Current point.
194     * @param centroid Coordinates through which reflection must be performed.
195     * @param eval Evaluation function.
196     * @return a new point with Cartesian coordinates set to the reflection
197     * of {@code point} through {@code centroid}.
198     */
199    private PointValuePair newReflectedPoint(PointValuePair point,
200                                             double[] centroid,
201                                             MultivariateFunction eval) {
202        final double alpha = alphaSampler.sample();
203        return Simplex.newPoint(centroid,
204                                -alpha,
205                                point.getPoint(),
206                                eval);
207    }
208
209    /** {@inheritDoc} */
210    @Override
211    public String toString() {
212        return "Hedar-Fukushima [s=" + sigma + "]";
213    }
214}