1 /*
2 * Copyright (C) 2012 Alberto Irurueta Carro (alberto@irurueta.com)
3 *
4 * Licensed under the Apache License, Version 2.0 (the "License");
5 * you may not use this file except in compliance with the License.
6 * You may obtain a copy of the License at
7 *
8 * http://www.apache.org/licenses/LICENSE-2.0
9 *
10 * Unless required by applicable law or agreed to in writing, software
11 * distributed under the License is distributed on an "AS IS" BASIS,
12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 * See the License for the specific language governing permissions and
14 * limitations under the License.
15 */
16 package com.irurueta.numerical.optimization;
17
18 import com.irurueta.algebra.Matrix;
19 import com.irurueta.algebra.WrongSizeException;
20 import com.irurueta.numerical.EvaluationException;
21 import com.irurueta.numerical.LockedException;
22 import com.irurueta.numerical.MultiDimensionFunctionEvaluatorListener;
23 import com.irurueta.numerical.NotAvailableException;
24 import com.irurueta.numerical.NotReadyException;
25
26 import java.util.Arrays;
27
28 /**
29 * This class searches for a multi dimension function local minimum.
30 * The local minimum is searched by starting the algorithm at a start point
31 * and a given set of deltas to obtain surrounding points to the start point
32 * where the algorithm will be started.
33 * The Simplex algorithm will increase or decrease such deltas and move the
34 * start point around until the minimum is found.
35 * The implementation of this class is based on Numerical Recipes 3rd ed.
36 * Section 10.5 page 502.
37 */
38 public class SimplexMultiOptimizer extends MultiOptimizer {
39 /**
40 * Maximum number of iterations.
41 */
42 public static final int NMAX = 5000;
43
44 /**
45 * Small value considered to be machine precision.
46 */
47 public static final double TINY = 1e-10;
48
49 /**
50 * Constant defining default tolerance or accuracy to be achieved on the
51 * minimum being estimated by this class.
52 */
53 public static final double DEFAULT_TOLERANCE = 3e-8;
54
55 /**
56 * Minimum allowed tolerance value.
57 */
58 public static final double MIN_TOLERANCE = 0.0;
59
60 /**
61 * The fractional tolerance in the function value such that failure to
62 * decrease by more than this amount on one iteration signals doneness.
63 */
64 private double ftol;
65
66 /**
67 * Number of function evaluations.
68 */
69 private int nfunc;
70
71 /**
72 * Number of points in the simplex.
73 */
74 private int mpts;
75
76 /**
77 * Number of dimensions of current function being optimized.
78 */
79 private int ndim;
80
81 /**
82 * Function values at the vertices of the simples.
83 */
84 private double[] y;
85
86 /**
87 * Current simplex.
88 */
89 private Matrix p;
90
91 /**
92 * Empty constructor.
93 */
94 public SimplexMultiOptimizer() {
95 super();
96 nfunc = mpts = ndim = 0;
97 y = null;
98 p = null;
99 ftol = DEFAULT_TOLERANCE;
100 }
101
102 /**
103 * Constructor.
104 *
105 * @param listener Listener to evaluate a multidimensional function.
106 * @param startPoint Point where the algorithm is started.
107 * @param delta Delta to find other points close to the start point so that
108 * an initial simplex is defined.
109 * @param tolerance Tolerance or accuracy to be expected on estimated local
110 * minimum.
111 * @throws IllegalArgumentException Raised if tolerance is negative.
112 */
113 public SimplexMultiOptimizer(
114 final MultiDimensionFunctionEvaluatorListener listener, final double[] startPoint, final double delta,
115 final double tolerance) {
116 super(listener);
117 nfunc = mpts = ndim = 0;
118 y = null;
119 p = null;
120 internalSetTolerance(tolerance);
121 internalSetSimplex(startPoint, delta);
122 }
123
124 /**
125 * Constructor.
126 *
127 * @param listener Listener to evaluate a multidimensional function.
128 * @param startPoint Point where the algorithm is started.
129 * @param deltas Set of deltas to find other points close to the start point
130 * so that an initial simplex is defined.
131 * @param tolerance Tolerance or accuracy to be expected on estimated local
132 * minimum.
133 * @throws IllegalArgumentException Raised if tolerance is negative or if
134 * deltas don't have the same length as provided start point.
135 */
136 public SimplexMultiOptimizer(
137 final MultiDimensionFunctionEvaluatorListener listener, final double[] startPoint, final double[] deltas,
138 final double tolerance) {
139 super(listener);
140 nfunc = mpts = ndim = 0;
141 y = null;
142 p = null;
143 internalSetTolerance(tolerance);
144 internalSetSimplex(startPoint, deltas);
145 }
146
147 /**
148 * Constructor.
149 *
150 * @param listener Listener to evaluate a multidimensional function.
151 * @param simplex Initial simplex to start the algorithm. The simplex is
152 * a set of points around the minimum to be found.
153 * @param tolerance Tolerance or accuracy to be expected on estimated local
154 * minimum.
155 * @throws IllegalArgumentException Raised if tolerance is negative.
156 */
157 public SimplexMultiOptimizer(
158 final MultiDimensionFunctionEvaluatorListener listener, final Matrix simplex, final double tolerance) {
159 super(listener);
160 nfunc = mpts = ndim = 0;
161 y = null;
162 p = null;
163 internalSetTolerance(tolerance);
164 internalSetSimplex(simplex);
165 }
166
167 /**
168 * Returns current simplex.
169 * A simplex is a set of points delimiting an n-dimensional region.
170 * The simplex can be seen as the n-dimensional version of a bracket of
171 * values.
172 * This class must be started at an initial simplex where a minimum will be
173 * searched.
174 * As the algorithm iterates, the simplex will be moved around and its size
175 * will be reduced until a minimum is found.
176 *
177 * @return Current simplex.
178 * @throws NotAvailableException Raised if not provided or computed.
179 */
180 public Matrix getSimplex() throws NotAvailableException {
181 if (!isSimplexAvailable()) {
182 throw new NotAvailableException();
183 }
184 return p;
185 }
186
187 /**
188 * Sets a simplex defined as a central point and a set of surrounding points
189 * at distance delta. The simplex will be made of the computed surrounding
190 * points.
191 *
192 * @param startPoint Central point.
193 * @param delta Distance of surrounding points.
194 * @throws LockedException Raised if this instance is locked.
195 */
196 public void setSimplex(final double[] startPoint, final double delta) throws LockedException {
197 if (isLocked()) {
198 throw new LockedException();
199 }
200 internalSetSimplex(startPoint, delta);
201 }
202
203 /**
204 * Internal method to set a simplex as a central point and a set of
205 * surrounding points at distance delta. The simplex will be made of the
206 * computed surrounding points.
207 * This method does not check whether this instance is locked.
208 *
209 * @param startPoint Central point.
210 * @param delta Distance of surrounding points.
211 */
212 private void internalSetSimplex(final double[] startPoint, final double delta) {
213 final var deltas = new double[startPoint.length];
214 Arrays.fill(deltas, delta);
215 internalSetSimplex(startPoint, deltas);
216 }
217
218 /**
219 * Sets a simplex defined as a central point and a set of surrounding points
220 * at their corresponding distance deltas[i], where "i" corresponds to one
221 * position of provided array of distances. The simplex will be made of the
222 * computed surrounding points.
223 *
224 * @param startPoint Central point.
225 * @param deltas Distances of surrounding points. Each surrounding point can
226 * have a different distance than the others. The number of provided
227 * distances must be equal to the dimension or length of the start point
228 * array.
229 * @throws LockedException Raised if this instance is locked.
230 * @throws IllegalArgumentException Raised if startPoint and deltas don't
231 * have the same length.
232 */
233 public void setSimplex(final double[] startPoint, final double[] deltas) throws LockedException {
234 if (isLocked()) {
235 throw new LockedException();
236 }
237 internalSetSimplex(startPoint, deltas);
238 }
239
240 /**
241 * Internal method to set a simplex defined as a central point and a set of
242 * surrounding points at their corresponding distance deltas[i], where "i"
243 * corresponds to one position of provided array of distances. The simplex
244 * will be made of the computed surrounding points.
245 *
246 * @param startPoint Central point.
247 * @param deltas Distances of surrounding points. Each surrounding point can
248 * have a different distance than the others. The number of provided
249 * distances must be equal to the dimension or length of the start point
250 * array.
251 * @throws IllegalArgumentException Raised if startPoint and deltas don't
252 * have the same length.
253 */
254 private void internalSetSimplex(final double[] startPoint, final double[] deltas) {
255 final var localndim = startPoint.length;
256
257 if (deltas.length != localndim) {
258 throw new IllegalArgumentException();
259 }
260
261 final Matrix pp;
262 try {
263 pp = new Matrix(localndim + 1, localndim);
264 } catch (final WrongSizeException e) {
265 throw new IllegalArgumentException(e);
266 }
267 for (var i = 0; i < localndim + 1; i++) {
268 for (var j = 0; j < localndim; j++) {
269 pp.setElementAt(i, j, startPoint[j]);
270 }
271 if (i != 0) {
272 pp.setElementAt(i, i - 1, pp.getElementAt(i, i - 1) + deltas[i - 1]);
273 }
274 }
275
276 internalSetSimplex(pp);
277 }
278
279 /**
280 * Sets simplex as a matrix containing on each row a point of the simplex.
281 * The number of columns defines the dimension of the points in the
282 * function, if not properly set, then minimize() method will fail when
283 * being called.
284 *
285 * @param simplex Simplex of points.
286 * @throws LockedException Raised if this instance is locked.
287 */
288 public void setSimplex(final Matrix simplex) throws LockedException {
289 if (isLocked()) {
290 throw new LockedException();
291 }
292 internalSetSimplex(simplex);
293 }
294
295 /**
296 * Internal method to Set simplex as a matrix containing on each row a point
297 * of the simplex.
298 * The number of columns defines the dimension of the points in the
299 * function, if not properly set, then minimize() method will fail when
300 * being called.
301 * This method does not check whether this instance is locked.
302 *
303 * @param simplex Simplex of points.
304 */
305 private void internalSetSimplex(final Matrix simplex) {
306 p = simplex;
307 }
308
309 /**
310 * Returns boolean indicating whether a simplex has been provided and is
311 * available for retrieval.
312 *
313 * @return True if available, false otherwise.
314 */
315 public boolean isSimplexAvailable() {
316 return p != null;
317 }
318
319 /**
320 * Returns function evaluations at simplex points or vertices.
321 *
322 * @return Function evaluations at simplex vertices.
323 * @throws NotAvailableException Raised if not available for retrieval.
324 * This parameter will be available one minimization has been computed.
325 */
326 public double[] getEvaluationsAtSimplex() throws NotAvailableException {
327 if (!areFunctionEvaluationsAvailable()) {
328 throw new NotAvailableException();
329 }
330 if (!areFunctionEvaluationsAvailable()) {
331 throw new NotAvailableException();
332 }
333 return y;
334 }
335
336 /**
337 * Returns boolean indicating whether function evaluations at simplex
338 * vertices are available for retrieval.
339 * Function evaluations at simplex vertices will be available one
340 * minimization has been computed.
341 *
342 * @return True if available, false otherwise.
343 */
344 public boolean areFunctionEvaluationsAvailable() {
345 return y != null;
346 }
347
348 /**
349 * Returns tolerance or accuracy to be expected on estimated local minimum.
350 *
351 * @return Tolerance or accuracy to be expected on estimated local minimum.
352 */
353 public double getTolerance() {
354 return ftol;
355 }
356
357 /**
358 * Sets tolerance or accuracy to be expected on estimated local minimum.
359 *
360 * @param tolerance Tolerance or accuracy to be expected on estimated local
361 * minimum.
362 * @throws LockedException Raised if this instance is locked.
363 * @throws IllegalArgumentException Raised if provided tolerance is
364 * negative.
365 */
366 public void setTolerance(final double tolerance) throws LockedException {
367 if (isLocked()) {
368 throw new LockedException();
369 }
370 internalSetTolerance(tolerance);
371 }
372
373 /**
374 * Internal method to set tolerance or accuracy to be expected on estimated
375 * local minimum.
376 * This method does not check whether this instance is locked.
377 *
378 * @param tolerance Tolerance or accuracy to be expected on estimated local
379 * minimum.
380 * @throws IllegalArgumentException Raised if provided tolerance is
381 * negative.
382 */
383 private void internalSetTolerance(final double tolerance) {
384 if (tolerance < MIN_TOLERANCE) {
385 throw new IllegalArgumentException();
386 }
387 ftol = tolerance;
388 }
389
390 /**
391 * This function estimates a function minimum.
392 * Implementations of this class will usually search a local minimum inside
393 * a given simplex of points. The algorithm will iterate moving the simplex
394 * around and reducing its size until a minimum is found.
395 *
396 * @throws LockedException Raised if this instance is locked, because
397 * estimation is being computed.
398 * @throws NotReadyException Raised if this instance is not ready, because
399 * a listener, a gradient listener and a start point haven't been provided.
400 * @throws OptimizationException Raised if the algorithm failed because of
401 * lack of convergence or because function couldn't be evaluated.
402 */
403 @Override
404 public void minimize() throws LockedException, NotReadyException, OptimizationException {
405 if (isLocked()) {
406 throw new LockedException();
407 }
408 if (!isReady()) {
409 throw new NotReadyException();
410 }
411
412 locked = true;
413
414 // clean previous result and previous function evaluations
415 try {
416 int ihi;
417 int ilo;
418 int inhi;
419 mpts = p.getRows();
420 ndim = p.getColumns();
421
422 final var psum = new double[ndim];
423 final var pmin = new double[ndim];
424 final var x = new double[ndim];
425 final var v1 = new double[1];
426 final var v2 = new double[1];
427
428 // instantiate new function evaluations
429 y = new double[mpts];
430
431 for (var i = 0; i < mpts; i++) {
432 for (var j = 0; j < ndim; j++) {
433 x[j] = p.getElementAt(i, j);
434 }
435 y[i] = listener.evaluate(x);
436 }
437
438 nfunc = 0;
439 getPsum(p, psum);
440
441 for (; ; ) {
442 ilo = 0;
443 if (y[0] > y[1]) {
444 ihi = 0;
445 inhi = 1;
446 } else {
447 ihi = 1;
448 inhi = 0;
449 }
450 for (var i = 0; i < mpts; i++) {
451 if (y[i] <= y[ilo]) {
452 ilo = i;
453 }
454 if (y[i] > y[ihi]) {
455 inhi = ihi;
456 ihi = i;
457 } else if (y[i] > y[inhi] && i != ihi) {
458 inhi = i;
459 }
460 }
461 final var rtol = 2.0 * Math.abs(y[ihi] - y[ilo]) / (Math.abs(y[ihi]) + Math.abs(y[ilo]) + TINY);
462
463 if (rtol < ftol) {
464 v1[0] = y[0];
465 v2[0] = y[ilo];
466 swap(v1, v2);
467 y[0] = v1[0];
468 y[ilo] = v2[0];
469 for (var i = 0; i < ndim; i++) {
470 v1[0] = p.getElementAt(0, i);
471 v2[0] = p.getElementAt(ilo, i);
472 swap(v1, v2);
473 p.setElementAt(0, i, v1[0]);
474 p.setElementAt(ilo, i, v2[0]);
475 pmin[i] = p.getElementAt(0, i);
476 }
477 fmin = y[0];
478
479 // save result
480 xmin = pmin;
481 resultAvailable = true;
482
483 return;
484 }
485 if (nfunc >= NMAX) {
486 throw new OptimizationException();
487 }
488
489 nfunc += 2;
490 var ytry = amotry(p, y, psum, ihi, -1.0, listener);
491
492 if (ytry <= y[ilo]) {
493 amotry(p, y, psum, ihi, 2.0, listener);
494 } else if (ytry >= y[inhi]) {
495 var ysave = y[ihi];
496 ytry = amotry(p, y, psum, ihi, 0.5, listener);
497 if (ytry >= ysave) {
498 for (var i = 0; i < mpts; i++) {
499 if (i != ilo) {
500 for (var j = 0; j < ndim; j++) {
501 psum[j] = 0.5 * (p.getElementAt(i, j) + p.getElementAt(ilo, j));
502 p.setElementAt(i, j, psum[j]);
503 }
504 y[i] = listener.evaluate(psum);
505 }
506 }
507 nfunc += ndim;
508 getPsum(p, psum);
509 }
510 } else {
511 --nfunc;
512 }
513
514 if (iterationCompletedListener != null) {
515 iterationCompletedListener.onIterationCompleted(this, nfunc, NMAX);
516 }
517 }
518 } catch (final EvaluationException e) {
519 throw new OptimizationException(e);
520 } finally {
521 locked = false;
522 }
523 }
524
525 /**
526 * Returns boolean indicating whether this instance is ready to start the
527 * estimation of a local minimum.
528 * An instance is ready once a listener and a simplex are provided.
529 *
530 * @return True if this instance is ready, false otherwise.
531 */
532 @Override
533 public boolean isReady() {
534 return isListenerAvailable() && isSimplexAvailable();
535 }
536
537 /**
538 * Computes the sum of the elements of each matrix column.
539 *
540 * @param p Matrix where columns will be summed.
541 * @param psum Array to contain computed sums. Provided array must have a
542 * length equal to the number of matrix columns. This is an output
543 * parameter.
544 */
545 private void getPsum(final Matrix p, final double[] psum) {
546 for (var j = 0; j < ndim; j++) {
547 var sum = 0.0;
548 for (int i = 0; i < mpts; i++) {
549 sum += p.getElementAt(i, j);
550 }
551 psum[j] = sum;
552 }
553 }
554
555 /**
556 * Internal method to move simplex around.
557 *
558 * @param p a matrix.
559 * @param y an array.
560 * @param psum another array.
561 * @param ihi a value.
562 * @param fac another value.
563 * @param listener a listener to evaluate function to optimize.
564 * @return a value returned by the evaluated function.
565 * @throws EvaluationException Raised if function evaluation fails.
566 */
567 private double amotry(final Matrix p, final double[] y, final double[] psum, final int ihi, final double fac,
568 final MultiDimensionFunctionEvaluatorListener listener) throws EvaluationException {
569 final var ptry = new double[ndim];
570 final var fac1 = (1.0 - fac) / ndim;
571 final var fac2 = fac1 - fac;
572 for (var j = 0; j < ndim; j++) {
573 ptry[j] = psum[j] * fac1 - p.getElementAt(ihi, j) * fac2;
574 }
575 final var ytry = listener.evaluate(ptry);
576 if (ytry < y[ihi]) {
577 y[ihi] = ytry;
578 for (var j = 0; j < ndim; j++) {
579 psum[j] += ptry[j] - p.getElementAt(ihi, j);
580 p.setElementAt(ihi, j, ptry[j]);
581 }
582 }
583
584 return ytry;
585 }
586
587 /**
588 * Swaps a and b values.
589 *
590 * @param a value.
591 * @param b value.
592 */
593 private static void swap(final double[] a, final double[] b) {
594 final var tmp = a[0];
595 a[0] = b[0];
596 b[0] = tmp;
597 }
598 }