View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.apache.commons.math4.legacy.optim.nonlinear.scalar.noderiv;
18  
19  import java.util.Comparator;
20  import java.util.List;
21  import java.util.ArrayList;
22  import java.util.Collections;
23  import java.util.function.UnaryOperator;
24  import java.util.function.DoublePredicate;
25  import org.apache.commons.rng.UniformRandomProvider;
26  import org.apache.commons.rng.simple.RandomSource;
27  import org.apache.commons.rng.sampling.distribution.ContinuousSampler;
28  import org.apache.commons.rng.sampling.distribution.ContinuousUniformSampler;
29  import org.apache.commons.math4.legacy.analysis.MultivariateFunction;
30  import org.apache.commons.math4.legacy.optim.PointValuePair;
31  
32  /**
33   * DSSA algorithm.
34   *
35   * Described in
36   * <blockquote>
37   *  <em>Abdel-Rahman Hedar and Masao Fukushima (2002)</em>,
38   *  <b>
39   *   Hybrid simulated annealing and direct search method
40   *   for nonlinear unconstrained global optimization
41   *  </b>,
42   *  Optimization Methods and Software, 17:5, 891-912,
43   *  DOI: 10.1080/1055678021000030084
44   * </blockquote>
45   *
46   * <p>
47   * A note about the {@link #HedarFukushimaTransform(double) "shrink" factor}:
48   * Per DSSA's description, the simplex must keep its size during the simulated
49   * annealing (SA) phase to avoid premature convergence.  This assumes that the
50   * best candidates from the SA phase will each subsequently serve as starting
51   * point for another optimization to hone in on the local optimum.
52   * Values lower than 1 and no subsequent "best list" search correspond to the
53   * "SSA" algorithm in the above paper.
54   */
55  public class HedarFukushimaTransform
56      implements Simplex.TransformFactory {
57      /** Shrinkage coefficient. */
58      private final double sigma;
59      /** Sampler for reflection coefficient. */
60      private final ContinuousSampler alphaSampler;
61      /** No shrink indicator. */
62      private final boolean noShrink;
63  
64      /**
65       * @param sigma Shrink factor.
66       * @param rng Random generator.
67       * @throws IllegalArgumentException if {@code sigma <= 0} or
68       * {@code sigma > 1}.
69       */
70      public HedarFukushimaTransform(double sigma,
71                                     UniformRandomProvider rng) {
72          if (sigma <= 0 ||
73              sigma > 1) {
74              throw new IllegalArgumentException("Shrink factor out of range: " +
75                                                 sigma);
76          }
77  
78          this.sigma = sigma;
79          alphaSampler = ContinuousUniformSampler.of(rng, 0.9, 1.1);
80          noShrink = sigma == 1d;
81      }
82  
83      /**
84       * @param sigma Shrink factor.
85       * @throws IllegalArgumentException if {@code sigma <= 0} or
86       * {@code sigma > 1}.
87       */
88      public HedarFukushimaTransform(double sigma) {
89          this(sigma, RandomSource.KISS.create());
90      }
91  
92      /**
93       * Disable shrinking of the simplex (as mandated by DSSA).
94       */
95      public HedarFukushimaTransform() {
96          this(1d);
97      }
98  
99      /** {@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 }