source: anac2020/agentKT/src/main/java/geniusweb/exampleparties/simpleshaop/Cobyla.java

Last change on this file was 1, checked in by wouter, 4 years ago

#1910 added anac2020 parties
File size: 56.6 KB
RevLine 
[1]1/*
2 * jcobyla
3 *
4 * The MIT License
5 *
6 * Copyright (c) 2012 Anders Gustafsson, Cureos AB.
7 *
8 * Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files
9 * (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge,
10 * publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so,
11 * subject to the following conditions:
12 *
13 * The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
14 *
15 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
16 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE
17 * FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
18 * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
19 *
20 * Remarks:
21 *
22 * The original Fortran 77 version of this code was by Michael Powell (M.J.D.Powell @ damtp.cam.ac.uk)
23 * The Fortran 90 version was by Alan Miller (Alan.Miller @ vic.cmis.csiro.au). Latest revision - 30 October 1998
24 */
25package geniusweb.exampleparties.simpleshaop;
26
27/**
28 * Constrained Optimization BY Linear Approximation in Java.
29 *
30 * COBYLA2 is an implementation of Powell’s nonlinear derivative–free constrained optimization that uses
31 * a linear approximation approach. The algorithm is a sequential trust–region algorithm that employs linear
32 * approximations to the objective and constraint functions, where the approximations are formed by linear
33 * interpolation at n + 1 points in the space of the variables and tries to maintain a regular–shaped simplex
34 * over iterations.
35 *
36 * It solves nonsmooth NLP with a moderate number of variables (about 100). Inequality constraints only.
37 *
38 * The initial point X is taken as one vertex of the initial simplex with zero being another, so, X should
39 * not be entered as the zero vector.
40 *
41 * @author Anders Gustafsson, Cureos AB.
42 */
43public class Cobyla
44{
45 /**
46 * Minimizes the objective function F with respect to a set of inequality constraints CON,
47 * and returns the optimal variable array. F and CON may be non-linear, and should preferably be smooth.
48 *
49 * @param calcfc Interface implementation for calculating objective function and constraints.
50 * @param n Number of variables.
51 * @param m Number of constraints.
52 * @param x On input initial values of the variables (zero-based array). On output
53 * optimal values of the variables obtained in the COBYLA minimization.
54 * @param rhobeg Initial size of the simplex.
55 * @param rhoend Final value of the simplex.
56 * @param iprint Print level, 0 <= iprint <= 3, where 0 provides no output and
57 * 3 provides full output to the console.
58 * @param maxfun Maximum number of function evaluations before terminating.
59 * @return Exit status of the COBYLA2 optimization.
60 */
61 public static CobylaExitStatus FindMinimum(final Calcfc calcfc, int n, int m, double[] x, double rhobeg, double rhoend, int iprint, int maxfun)
62 {
63 // This subroutine minimizes an objective function F(X) subject to M
64 // inequality constraints on X, where X is a vector of variables that has
65 // N components. The algorithm employs linear approximations to the
66 // objective and constraint functions, the approximations being formed by
67 // linear interpolation at N+1 points in the space of the variables.
68 // We regard these interpolation points as vertices of a simplex. The
69 // parameter RHO controls the size of the simplex and it is reduced
70 // automatically from RHOBEG to RHOEND. For each RHO the subroutine tries
71 // to achieve a good vector of variables for the current size, and then
72 // RHO is reduced until the value RHOEND is reached. Therefore RHOBEG and
73 // RHOEND should be set to reasonable initial changes to and the required
74 // accuracy in the variables respectively, but this accuracy should be
75 // viewed as a subject for experimentation because it is not guaranteed.
76 // The subroutine has an advantage over many of its competitors, however,
77 // which is that it treats each constraint individually when calculating
78 // a change to the variables, instead of lumping the constraints together
79 // into a single penalty function. The name of the subroutine is derived
80 // from the phrase Constrained Optimization BY Linear Approximations.
81
82 // The user must set the values of N, M, RHOBEG and RHOEND, and must
83 // provide an initial vector of variables in X. Further, the value of
84 // IPRINT should be set to 0, 1, 2 or 3, which controls the amount of
85 // printing during the calculation. Specifically, there is no output if
86 // IPRINT=0 and there is output only at the end of the calculation if
87 // IPRINT=1. Otherwise each new value of RHO and SIGMA is printed.
88 // Further, the vector of variables and some function information are
89 // given either when RHO is reduced or when each new value of F(X) is
90 // computed in the cases IPRINT=2 or IPRINT=3 respectively. Here SIGMA
91 // is a penalty parameter, it being assumed that a change to X is an
92 // improvement if it reduces the merit function
93 // F(X)+SIGMA*MAX(0.0, - C1(X), - C2(X),..., - CM(X)),
94 // where C1,C2,...,CM denote the constraint functions that should become
95 // nonnegative eventually, at least to the precision of RHOEND. In the
96 // printed output the displayed term that is multiplied by SIGMA is
97 // called MAXCV, which stands for 'MAXimum Constraint Violation'. The
98 // argument ITERS is an integer variable that must be set by the user to a
99 // limit on the number of calls of CALCFC, the purpose of this routine being
100 // given below. The value of ITERS will be altered to the number of calls
101 // of CALCFC that are made.
102
103 // In order to define the objective and constraint functions, we require
104 // a subroutine that has the name and arguments
105 // SUBROUTINE CALCFC (N,M,X,F,CON)
106 // DIMENSION X(:),CON(:) .
107 // The values of N and M are fixed and have been defined already, while
108 // X is now the current vector of variables. The subroutine should return
109 // the objective and constraint functions at X in F and CON(1),CON(2),
110 // ...,CON(M). Note that we are trying to adjust X so that F(X) is as
111 // small as possible subject to the constraint functions being nonnegative.
112
113 // Local variables
114 int mpp = m + 2;
115
116 // Internal base-1 X array
117 double[] iox = new double[n + 1];
118 System.arraycopy(x, 0, iox, 1, n);
119
120 // Internal representation of the objective and constraints calculation method,
121 // accounting for that X and CON arrays in the cobylb method are base-1 arrays.
122 Calcfc fcalcfc = new Calcfc()
123 {
124 /**
125 *
126 * @param n the value of n
127 * @param m the value of m
128 * @param x the value of x
129 * @param con the value of con
130 * @return the double
131 */
132 @Override
133 public double Compute(int n, int m, double[] x, double[] con)
134 {
135 double[] ix = new double[n];
136 System.arraycopy(x, 1, ix, 0, n);
137 double[] ocon = new double[m];
138 double f = calcfc.Compute(n, m, ix, ocon);
139 System.arraycopy(ocon, 0, con, 1, m);
140 return f;
141 }
142 };
143
144 CobylaExitStatus status = cobylb(fcalcfc, n, m, mpp, iox, rhobeg, rhoend, iprint, maxfun);
145 System.arraycopy(iox, 1, x, 0, n);
146
147 return status;
148 }
149
150 private static CobylaExitStatus cobylb(Calcfc calcfc, int n, int m, int mpp, double[] x,
151 double rhobeg, double rhoend, int iprint, int maxfun)
152 {
153 // N.B. Arguments CON, SIM, SIMI, DATMAT, A, VSIG, VETA, SIGBAR, DX, W & IACT
154 // have been removed.
155
156 // Set the initial values of some parameters. The last column of SIM holds
157 // the optimal vertex of the current simplex, and the preceding N columns
158 // hold the displacements from the optimal vertex to the other vertices.
159 // Further, SIMI holds the inverse of the matrix that is contained in the
160 // first N columns of SIM.
161
162 // Local variables
163
164 CobylaExitStatus status;
165
166 double alpha = 0.25;
167 double beta = 2.1;
168 double gamma = 0.5;
169 double delta = 1.1;
170
171 double f = 0.0;
172 double resmax = 0.0;
173 double total;
174
175 int np = n + 1;
176 int mp = m + 1;
177 double rho = rhobeg;
178 double parmu = 0.0;
179
180 boolean iflag = false;
181 boolean ifull = false;
182 double parsig = 0.0;
183 double prerec = 0.0;
184 double prerem = 0.0;
185
186 double[] con = new double[1 + mpp];
187 double[][] sim = new double[1 + n][1 + np];
188 double[][] simi = new double[1 + n][1 + n];
189 double[][] datmat = new double[1 + mpp][1 + np];
190 double[][] a = new double[1 + n][1 + mp];
191 double[] vsig = new double[1 + n];
192 double[] veta = new double[1 + n];
193 double[] sigbar = new double[1 + n];
194 double[] dx = new double[1 + n];
195 double[] w = new double[1 + n];
196
197 if (iprint >= 2) System.out.format("%nThe initial value of RHO is %13.6f and PARMU is set to zero.%n", rho);
198
199 int nfvals = 0;
200 double temp = 1.0 / rho;
201
202 for (int i = 1; i <= n; ++i)
203 {
204 sim[i][np] = x[i];
205 sim[i][i] = rho;
206 simi[i][i] = temp;
207 }
208
209 int jdrop = np;
210 boolean ibrnch = false;
211
212 // Make the next call of the user-supplied subroutine CALCFC. These
213 // instructions are also used for calling CALCFC during the iterations of
214 // the algorithm.
215
216 L_40:
217 do
218 {
219 if (nfvals >= maxfun && nfvals > 0)
220 {
221 status = CobylaExitStatus.MaxIterationsReached;
222 break L_40;
223 }
224
225 ++nfvals;
226
227 f = calcfc.Compute(n, m, x, con);
228 resmax = 0.0; for (int k = 1; k <= m; ++k) resmax = Math.max(resmax, -con[k]);
229
230 if (nfvals == iprint - 1 || iprint == 3)
231 {
232 PrintIterationResult(nfvals, f, resmax, x, n);
233 }
234
235 con[mp] = f;
236 con[mpp] = resmax;
237
238 // Set the recently calculated function values in a column of DATMAT. This
239 // array has a column for each vertex of the current simplex, the entries of
240 // each column being the values of the constraint functions (if any)
241 // followed by the objective function and the greatest constraint violation
242 // at the vertex.
243
244 boolean skipVertexIdent = true;
245 if (!ibrnch)
246 {
247 skipVertexIdent = false;
248
249 for (int i = 1; i <= mpp; ++i) datmat[i][jdrop] = con[i];
250
251 if (nfvals <= np)
252 {
253 // Exchange the new vertex of the initial simplex with the optimal vertex if
254 // necessary. Then, if the initial simplex is not complete, pick its next
255 // vertex and calculate the function values there.
256
257 if (jdrop <= n)
258 {
259 if (datmat[mp][np] <= f)
260 {
261 x[jdrop] = sim[jdrop][np];
262 }
263 else
264 {
265 sim[jdrop][np] = x[jdrop];
266 for (int k = 1; k <= mpp; ++k)
267 {
268 datmat[k][jdrop] = datmat[k][np];
269 datmat[k][np] = con[k];
270 }
271 for (int k = 1; k <= jdrop; ++k)
272 {
273 sim[jdrop][k] = -rho;
274 temp = 0.0; for (int i = k; i <= jdrop; ++i) temp -= simi[i][k];
275 simi[jdrop][k] = temp;
276 }
277 }
278 }
279 if (nfvals <= n)
280 {
281 jdrop = nfvals;
282 x[jdrop] += rho;
283 continue L_40;
284 }
285 }
286
287 ibrnch = true;
288 }
289
290 L_140:
291 do
292 {
293 L_550:
294 do
295 {
296 if (!skipVertexIdent)
297 {
298 // Identify the optimal vertex of the current simplex.
299
300 double phimin = datmat[mp][np] + parmu * datmat[mpp][np];
301 int nbest = np;
302
303 for (int j = 1; j <= n; ++j)
304 {
305 temp = datmat[mp][j] + parmu * datmat[mpp][j];
306 if (temp < phimin)
307 {
308 nbest = j;
309 phimin = temp;
310 }
311 else if (temp == phimin && parmu == 0.0 && datmat[mpp][j] < datmat[mpp][nbest])
312 {
313 nbest = j;
314 }
315 }
316
317 // Switch the best vertex into pole position if it is not there already,
318 // and also update SIM, SIMI and DATMAT.
319
320 if (nbest <= n)
321 {
322 for (int i = 1; i <= mpp; ++i)
323 {
324 temp = datmat[i][np];
325 datmat[i][np] = datmat[i][nbest];
326 datmat[i][nbest] = temp;
327 }
328 for (int i = 1; i <= n; ++i)
329 {
330 temp = sim[i][nbest];
331 sim[i][nbest] = 0.0;
332 sim[i][np] += temp;
333
334 double tempa = 0.0;
335 for (int k = 1; k <= n; ++k)
336 {
337 sim[i][k] -= temp;
338 tempa -= simi[k][i];
339 }
340 simi[nbest][i] = tempa;
341 }
342 }
343
344 // Make an error return if SIGI is a poor approximation to the inverse of
345 // the leading N by N submatrix of SIG.
346
347 double error = 0.0;
348 for (int i = 1; i <= n; ++i)
349 {
350 for (int j = 1; j <= n; ++j)
351 {
352 temp = DOT_PRODUCT(PART(ROW(simi, i), 1, n), PART(COL(sim, j), 1, n)) - (i == j ? 1.0 : 0.0);
353 error = Math.max(error, Math.abs(temp));
354 }
355 }
356 if (error > 0.1)
357 {
358 status = CobylaExitStatus.DivergingRoundingErrors;
359 break L_40;
360 }
361
362 // Calculate the coefficients of the linear approximations to the objective
363 // and constraint functions, placing minus the objective function gradient
364 // after the constraint gradients in the array A. The vector W is used for
365 // working space.
366
367 for (int k = 1; k <= mp; ++k)
368 {
369 con[k] = -datmat[k][np];
370 for (int j = 1; j <= n; ++j) w[j] = datmat[k][j] + con[k];
371
372 for (int i = 1; i <= n; ++i)
373 {
374 a[i][k] = (k == mp ? -1.0 : 1.0) * DOT_PRODUCT(PART(w, 1, n), PART(COL(simi, i), 1, n));
375 }
376 }
377
378 // Calculate the values of sigma and eta, and set IFLAG = 0 if the current
379 // simplex is not acceptable.
380
381 iflag = true;
382 parsig = alpha * rho;
383 double pareta = beta * rho;
384
385 for (int j = 1; j <= n; ++j)
386 {
387 double wsig = 0.0; for (int k = 1; k <= n; ++k) wsig += simi[j][k] * simi[j][k];
388 double weta = 0.0; for (int k = 1; k <= n; ++k) weta += sim[k][j] * sim[k][j];
389 vsig[j] = 1.0 / Math.sqrt(wsig);
390 veta[j] = Math.sqrt(weta);
391 if (vsig[j] < parsig || veta[j] > pareta) iflag = false;
392 }
393
394 // If a new vertex is needed to improve acceptability, then decide which
395 // vertex to drop from the simplex.
396
397 if (!ibrnch && !iflag)
398 {
399 jdrop = 0;
400 temp = pareta;
401 for (int j = 1; j <= n; ++j)
402 {
403 if (veta[j] > temp)
404 {
405 jdrop = j;
406 temp = veta[j];
407 }
408 }
409 if (jdrop == 0)
410 {
411 for (int j = 1; j <= n; ++j)
412 {
413 if (vsig[j] < temp)
414 {
415 jdrop = j;
416 temp = vsig[j];
417 }
418 }
419 }
420
421 // Calculate the step to the new vertex and its sign.
422
423 temp = gamma * rho * vsig[jdrop];
424 for (int k = 1; k <= n; ++k) dx[k] = temp * simi[jdrop][k];
425 double cvmaxp = 0.0;
426 double cvmaxm = 0.0;
427
428 total = 0.0;
429 for (int k = 1; k <= mp; ++k)
430 {
431 total = DOT_PRODUCT(PART(COL(a, k), 1, n), PART(dx, 1, n));
432 if (k < mp)
433 {
434 temp = datmat[k][np];
435 cvmaxp = Math.max(cvmaxp, -total - temp);
436 cvmaxm = Math.max(cvmaxm, total - temp);
437 }
438 }
439 double dxsign = parmu * (cvmaxp - cvmaxm) > 2.0 * total ? -1.0 : 1.0;
440
441 // Update the elements of SIM and SIMI, and set the next X.
442
443 temp = 0.0;
444 for (int i = 1; i <= n; ++i)
445 {
446 dx[i] = dxsign * dx[i];
447 sim[i][jdrop] = dx[i];
448 temp += simi[jdrop][i] * dx[i];
449 }
450 for (int k = 1; k <= n; ++k) simi[jdrop][k] /= temp;
451
452 for (int j = 1; j <= n; ++j)
453 {
454 if (j != jdrop)
455 {
456 temp = DOT_PRODUCT(PART(ROW(simi, j), 1, n), PART(dx, 1, n));
457 for (int k = 1; k <= n; ++k) simi[j][k] -= temp * simi[jdrop][k];
458 }
459 x[j] = sim[j][np] + dx[j];
460 }
461 continue L_40;
462 }
463
464 // Calculate DX = x(*)-x(0).
465 // Branch if the length of DX is less than 0.5*RHO.
466
467 ifull = trstlp(n, m, a, con, rho, dx);
468 if (!ifull)
469 {
470 temp = 0.0; for (int k = 1; k <= n; ++k) temp += dx[k] * dx[k];
471 if (temp < 0.25 * rho * rho)
472 {
473 ibrnch = true;
474 break L_550;
475 }
476 }
477
478 // Predict the change to F and the new maximum constraint violation if the
479 // variables are altered from x(0) to x(0) + DX.
480
481 total = 0.0;
482 double resnew = 0.0;
483 con[mp] = 0.0;
484 for (int k = 1; k <= mp; ++k)
485 {
486 total = con[k] - DOT_PRODUCT(PART(COL(a, k), 1, n), PART(dx, 1, n));
487 if (k < mp) resnew = Math.max(resnew, total);
488 }
489
490 // Increase PARMU if necessary and branch back if this change alters the
491 // optimal vertex. Otherwise PREREM and PREREC will be set to the predicted
492 // reductions in the merit function and the maximum constraint violation
493 // respectively.
494
495 prerec = datmat[mpp][np] - resnew;
496 double barmu = prerec > 0.0 ? total / prerec : 0.0;
497 if (parmu < 1.5 * barmu)
498 {
499 parmu = 2.0 * barmu;
500 if (iprint >= 2) System.out.format("%nIncrease in PARMU to %13.6f%n", parmu);
501 double phi = datmat[mp][np] + parmu * datmat[mpp][np];
502 for (int j = 1; j <= n; ++j)
503 {
504 temp = datmat[mp][j] + parmu * datmat[mpp][j];
505 if (temp < phi || (temp == phi && parmu == 0.0 && datmat[mpp][j] < datmat[mpp][np])) continue L_140;
506 }
507 }
508 prerem = parmu * prerec - total;
509
510 // Calculate the constraint and objective functions at x(*).
511 // Then find the actual reduction in the merit function.
512
513 for (int k = 1; k <= n; ++k) x[k] = sim[k][np] + dx[k];
514 ibrnch = true;
515 continue L_40;
516 }
517
518 skipVertexIdent = false;
519 double vmold = datmat[mp][np] + parmu * datmat[mpp][np];
520 double vmnew = f + parmu * resmax;
521 double trured = vmold - vmnew;
522 if (parmu == 0.0 && f == datmat[mp][np])
523 {
524 prerem = prerec;
525 trured = datmat[mpp][np] - resmax;
526 }
527
528 // Begin the operations that decide whether x(*) should replace one of the
529 // vertices of the current simplex, the change being mandatory if TRURED is
530 // positive. Firstly, JDROP is set to the index of the vertex that is to be
531 // replaced.
532
533 double ratio = trured <= 0.0 ? 1.0 : 0.0;
534 jdrop = 0;
535 for (int j = 1; j <= n; ++j)
536 {
537 temp = Math.abs(DOT_PRODUCT(PART(ROW(simi, j), 1, n), PART(dx, 1, n)));
538 if (temp > ratio)
539 {
540 jdrop = j;
541 ratio = temp;
542 }
543 sigbar[j] = temp * vsig[j];
544 }
545
546 // Calculate the value of ell.
547
548 double edgmax = delta * rho;
549 int l = 0;
550 for (int j = 1; j <= n; ++j)
551 {
552 if (sigbar[j] >= parsig || sigbar[j] >= vsig[j])
553 {
554 temp = veta[j];
555 if (trured > 0.0)
556 {
557 temp = 0.0; for (int k = 1; k <= n; ++k) temp += Math.pow(dx[k] - sim[k][j], 2.0);
558 temp = Math.sqrt(temp);
559 }
560 if (temp > edgmax)
561 {
562 l = j;
563 edgmax = temp;
564 }
565 }
566 }
567 if (l > 0) jdrop = l;
568
569 if (jdrop != 0)
570 {
571 // Revise the simplex by updating the elements of SIM, SIMI and DATMAT.
572
573 temp = 0.0;
574 for (int i = 1; i <= n; ++i)
575 {
576 sim[i][jdrop] = dx[i];
577 temp += simi[jdrop][i] * dx[i];
578 }
579 for (int k = 1; k <= n; ++k) simi[jdrop][k] /= temp;
580 for (int j = 1; j <= n; ++j)
581 {
582 if (j != jdrop)
583 {
584 temp = DOT_PRODUCT(PART(ROW(simi, j), 1, n), PART(dx, 1, n));
585 for (int k = 1; k <= n; ++k) simi[j][k] -= temp * simi[jdrop][k];
586 }
587 }
588 for (int k = 1; k <= mpp; ++k) datmat[k][jdrop] = con[k];
589
590 // Branch back for further iterations with the current RHO.
591
592 if (trured > 0.0 && trured >= 0.1 * prerem) continue L_140;
593 }
594 } while (false);
595
596 if (!iflag)
597 {
598 ibrnch = false;
599 continue L_140;
600 }
601
602 if (rho <= rhoend)
603 {
604 status = CobylaExitStatus.Normal;
605 break L_40;
606 }
607
608 // Otherwise reduce RHO if it is not at its least value and reset PARMU.
609
610 double cmin = 0.0, cmax = 0.0;
611
612 rho *= 0.5;
613 if (rho <= 1.5 * rhoend) rho = rhoend;
614 if (parmu > 0.0)
615 {
616 double denom = 0.0;
617 for (int k = 1; k <= mp; ++k)
618 {
619 cmin = datmat[k][np];
620 cmax = cmin;
621 for (int i = 1; i <= n; ++i)
622 {
623 cmin = Math.min(cmin, datmat[k][i]);
624 cmax = Math.max(cmax, datmat[k][i]);
625 }
626 if (k <= m && cmin < 0.5 * cmax)
627 {
628 temp = Math.max(cmax, 0.0) - cmin;
629 denom = denom <= 0.0 ? temp : Math.min(denom, temp);
630 }
631 }
632 if (denom == 0.0)
633 {
634 parmu = 0.0;
635 }
636 else if (cmax - cmin < parmu * denom)
637 {
638 parmu = (cmax - cmin) / denom;
639 }
640 }
641 if (iprint >= 2)
642 System.out.format("%nReduction in RHO to %1$13.6f and PARMU = %2$13.6f%n", rho, parmu);
643 if (iprint == 2)
644 PrintIterationResult(nfvals, datmat[mp][np], datmat[mpp][np], COL(sim, np), n);
645
646 } while (true);
647 } while (true);
648
649 switch (status)
650 {
651 case Normal:
652 if (iprint >= 1) System.out.format("%nNormal return from subroutine COBYLA%n");
653 if (ifull)
654 {
655 if (iprint >= 1) PrintIterationResult(nfvals, f, resmax, x, n);
656 return status;
657 }
658 break;
659 case MaxIterationsReached:
660 if (iprint >= 1)
661 System.out.format("%nReturn from subroutine COBYLA because the MAXFUN limit has been reached.%n");
662 break;
663 case DivergingRoundingErrors:
664 if (iprint >= 1)
665 System.out.format("%nReturn from subroutine COBYLA because rounding errors are becoming damaging.%n");
666 break;
667 }
668
669 for (int k = 1; k <= n; ++k) x[k] = sim[k][np];
670 f = datmat[mp][np];
671 resmax = datmat[mpp][np];
672 if (iprint >= 1) PrintIterationResult(nfvals, f, resmax, x, n);
673
674 return status;
675 }
676
677 private static boolean trstlp(int n, int m, double[][] a, double[] b, double rho, double[] dx)
678 {
679 // N.B. Arguments Z, ZDOTA, VMULTC, SDIRN, DXNEW, VMULTD & IACT have been removed.
680
681 // This subroutine calculates an N-component vector DX by applying the
682 // following two stages. In the first stage, DX is set to the shortest
683 // vector that minimizes the greatest violation of the constraints
684 // A(1,K)*DX(1)+A(2,K)*DX(2)+...+A(N,K)*DX(N) .GE. B(K), K = 2,3,...,M,
685 // subject to the Euclidean length of DX being at most RHO. If its length is
686 // strictly less than RHO, then we use the resultant freedom in DX to
687 // minimize the objective function
688 // -A(1,M+1)*DX(1) - A(2,M+1)*DX(2) - ... - A(N,M+1)*DX(N)
689 // subject to no increase in any greatest constraint violation. This
690 // notation allows the gradient of the objective function to be regarded as
691 // the gradient of a constraint. Therefore the two stages are distinguished
692 // by MCON .EQ. M and MCON .GT. M respectively. It is possible that a
693 // degeneracy may prevent DX from attaining the target length RHO. Then the
694 // value IFULL = 0 would be set, but usually IFULL = 1 on return.
695
696 // In general NACT is the number of constraints in the active set and
697 // IACT(1),...,IACT(NACT) are their indices, while the remainder of IACT
698 // contains a permutation of the remaining constraint indices. Further, Z
699 // is an orthogonal matrix whose first NACT columns can be regarded as the
700 // result of Gram-Schmidt applied to the active constraint gradients. For
701 // J = 1,2,...,NACT, the number ZDOTA(J) is the scalar product of the J-th
702 // column of Z with the gradient of the J-th active constraint. DX is the
703 // current vector of variables and here the residuals of the active
704 // constraints should be zero. Further, the active constraints have
705 // nonnegative Lagrange multipliers that are held at the beginning of
706 // VMULTC. The remainder of this vector holds the residuals of the inactive
707 // constraints at DX, the ordering of the components of VMULTC being in
708 // agreement with the permutation of the indices of the constraints that is
709 // in IACT. All these residuals are nonnegative, which is achieved by the
710 // shift RESMAX that makes the least residual zero.
711
712 // Initialize Z and some other variables. The value of RESMAX will be
713 // appropriate to DX = 0, while ICON will be the index of a most violated
714 // constraint if RESMAX is positive. Usually during the first stage the
715 // vector SDIRN gives a search direction that reduces all the active
716 // constraint violations by one simultaneously.
717
718 // Local variables
719
720 double temp;
721
722 int nactx = 0;
723 double resold = 0.0;
724
725 double[][] z = new double[1 + n][1 + n];
726 double[] zdota = new double[2 + m];
727 double[] vmultc = new double[2 + m];
728 double[] sdirn = new double[1 + n];
729 double[] dxnew = new double[1 + n];
730 double[] vmultd = new double[2 + m];
731 int[] iact = new int[2 + m];
732
733 int mcon = m;
734 int nact = 0;
735 for (int i = 1; i <= n; ++i)
736 {
737 z[i][i] = 1.0;
738 dx[i] = 0.0;
739 }
740
741 int icon = 0;
742 double resmax = 0.0;
743 if (m >= 1)
744 {
745 for (int k = 1; k <= m; ++k)
746 {
747 if (b[k] > resmax)
748 {
749 resmax = b[k];
750 icon = k;
751 }
752 }
753 for (int k = 1; k <= m; ++k)
754 {
755 iact[k] = k;
756 vmultc[k] = resmax - b[k];
757 }
758 }
759
760 // End the current stage of the calculation if 3 consecutive iterations
761 // have either failed to reduce the best calculated value of the objective
762 // function or to increase the number of active constraints since the best
763 // value was calculated. This strategy prevents cycling, but there is a
764 // remote possibility that it will cause premature termination.
765
766 boolean first = true;
767 do
768 {
769 L_60:
770 do
771 {
772 if (!first || (first && resmax == 0.0))
773 {
774 mcon = m + 1;
775 icon = mcon;
776 iact[mcon] = mcon;
777 vmultc[mcon] = 0.0;
778 }
779 first = false;
780
781 double optold = 0.0;
782 int icount = 0;
783
784 double step, stpful;
785
786 L_70:
787 do
788 {
789 double optnew = mcon == m ? resmax : -DOT_PRODUCT(PART(dx, 1, n), PART(COL(a, mcon), 1, n));
790
791 if (icount == 0 || optnew < optold)
792 {
793 optold = optnew;
794 nactx = nact;
795 icount = 3;
796 }
797 else if (nact > nactx)
798 {
799 nactx = nact;
800 icount = 3;
801 }
802 else
803 {
804 --icount;
805 }
806 if (icount == 0) break L_60;
807
808 // If ICON exceeds NACT, then we add the constraint with index IACT(ICON) to
809 // the active set. Apply Givens rotations so that the last N-NACT-1 columns
810 // of Z are orthogonal to the gradient of the new constraint, a scalar
811 // product being set to zero if its nonzero value could be due to computer
812 // rounding errors. The array DXNEW is used for working space.
813
814 double ratio;
815 if (icon <= nact)
816 {
817 if (icon < nact)
818 {
819 // Delete the constraint that has the index IACT(ICON) from the active set.
820
821 int isave = iact[icon];
822 double vsave = vmultc[icon];
823 int k = icon;
824 do
825 {
826 int kp = k + 1;
827 int kk = iact[kp];
828 double sp = DOT_PRODUCT(PART(COL(z, k), 1, n), PART(COL(a, kk), 1, n));
829 temp = Math.sqrt(sp * sp + zdota[kp] * zdota[kp]);
830 double alpha = zdota[kp] / temp;
831 double beta = sp / temp;
832 zdota[kp] = alpha * zdota[k];
833 zdota[k] = temp;
834 for (int i = 1; i <= n; ++i)
835 {
836 temp = alpha * z[i][kp] + beta * z[i][k];
837 z[i][kp] = alpha * z[i][k] - beta * z[i][kp];
838 z[i][k] = temp;
839 }
840 iact[k] = kk;
841 vmultc[k] = vmultc[kp];
842 k = kp;
843 } while (k < nact);
844
845 iact[k] = isave;
846 vmultc[k] = vsave;
847 }
848 --nact;
849
850 // If stage one is in progress, then set SDIRN to the direction of the next
851 // change to the current vector of variables.
852
853 if (mcon > m)
854 {
855 // Pick the next search direction of stage two.
856
857 temp = 1.0 / zdota[nact];
858 for (int k = 1; k <= n; ++k) sdirn[k] = temp * z[k][nact];
859 }
860 else
861 {
862 temp = DOT_PRODUCT(PART(sdirn, 1, n), PART(COL(z, nact + 1), 1, n));
863 for (int k = 1; k <= n; ++k) sdirn[k] -= temp * z[k][nact + 1];
864 }
865 }
866 else
867 {
868 int kk = iact[icon];
869 for (int k = 1; k <= n; ++k) dxnew[k] = a[k][kk];
870 double tot = 0.0;
871
872 {
873 int k = n;
874 while (k > nact)
875 {
876 double sp = 0.0;
877 double spabs = 0.0;
878 for (int i = 1; i <= n; ++i)
879 {
880 temp = z[i][k] * dxnew[i];
881 sp += temp;
882 spabs += Math.abs(temp);
883 }
884 double acca = spabs + 0.1 * Math.abs(sp);
885 double accb = spabs + 0.2 * Math.abs(sp);
886 if (spabs >= acca || acca >= accb) sp = 0.0;
887 if (tot == 0.0)
888 {
889 tot = sp;
890 }
891 else
892 {
893 int kp = k + 1;
894 temp = Math.sqrt(sp * sp + tot * tot);
895 double alpha = sp / temp;
896 double beta = tot / temp;
897 tot = temp;
898 for (int i = 1; i <= n; ++i)
899 {
900 temp = alpha * z[i][k] + beta * z[i][kp];
901 z[i][kp] = alpha * z[i][kp] - beta * z[i][k];
902 z[i][k] = temp;
903 }
904 }
905 --k;
906 }
907 }
908
909 if (tot == 0.0)
910 {
911 // The next instruction is reached if a deletion has to be made from the
912 // active set in order to make room for the new active constraint, because
913 // the new constraint gradient is a linear combination of the gradients of
914 // the old active constraints. Set the elements of VMULTD to the multipliers
915 // of the linear combination. Further, set IOUT to the index of the
916 // constraint to be deleted, but branch if no suitable index can be found.
917
918 ratio = -1.0;
919 {
920 int k = nact;
921 do
922 {
923 double zdotv = 0.0;
924 double zdvabs = 0.0;
925
926 for (int i = 1; i <= n; ++i)
927 {
928 temp = z[i][k] * dxnew[i];
929 zdotv += temp;
930 zdvabs += Math.abs(temp);
931 }
932 double acca = zdvabs + 0.1 * Math.abs(zdotv);
933 double accb = zdvabs + 0.2 * Math.abs(zdotv);
934 if (zdvabs < acca && acca < accb)
935 {
936 temp = zdotv / zdota[k];
937 if (temp > 0.0 && iact[k] <= m)
938 {
939 double tempa = vmultc[k] / temp;
940 if (ratio < 0.0 || tempa < ratio) ratio = tempa;
941 }
942
943 if (k >= 2)
944 {
945 int kw = iact[k];
946 for (int i = 1; i <= n; ++i) dxnew[i] -= temp * a[i][kw];
947 }
948 vmultd[k] = temp;
949 }
950 else
951 {
952 vmultd[k] = 0.0;
953 }
954 } while (--k > 0);
955 }
956 if (ratio < 0.0) break L_60;
957
958 // Revise the Lagrange multipliers and reorder the active constraints so
959 // that the one to be replaced is at the end of the list. Also calculate the
960 // new value of ZDOTA(NACT) and branch if it is not acceptable.
961
962 for (int k = 1; k <= nact; ++k)
963 vmultc[k] = Math.max(0.0, vmultc[k] - ratio * vmultd[k]);
964 if (icon < nact)
965 {
966 int isave = iact[icon];
967 double vsave = vmultc[icon];
968 int k = icon;
969 do
970 {
971 int kp = k + 1;
972 int kw = iact[kp];
973 double sp = DOT_PRODUCT(PART(COL(z, k), 1, n), PART(COL(a, kw), 1, n));
974 temp = Math.sqrt(sp * sp + zdota[kp] * zdota[kp]);
975 double alpha = zdota[kp] / temp;
976 double beta = sp / temp;
977 zdota[kp] = alpha * zdota[k];
978 zdota[k] = temp;
979 for (int i = 1; i <= n; ++i)
980 {
981 temp = alpha * z[i][kp] + beta * z[i][k];
982 z[i][kp] = alpha * z[i][k] - beta * z[i][kp];
983 z[i][k] = temp;
984 }
985 iact[k] = kw;
986 vmultc[k] = vmultc[kp];
987 k = kp;
988 } while (k < nact);
989 iact[k] = isave;
990 vmultc[k] = vsave;
991 }
992 temp = DOT_PRODUCT(PART(COL(z, nact), 1, n), PART(COL(a, kk), 1, n));
993 if (temp == 0.0) break L_60;
994 zdota[nact] = temp;
995 vmultc[icon] = 0.0;
996 vmultc[nact] = ratio;
997 }
998 else
999 {
1000 // Add the new constraint if this can be done without a deletion from the
1001 // active set.
1002
1003 ++nact;
1004 zdota[nact] = tot;
1005 vmultc[icon] = vmultc[nact];
1006 vmultc[nact] = 0.0;
1007 }
1008
1009 // Update IACT and ensure that the objective function continues to be
1010 // treated as the last active constraint when MCON>M.
1011
1012 iact[icon] = iact[nact];
1013 iact[nact] = kk;
1014 if (mcon > m && kk != mcon)
1015 {
1016 int k = nact - 1;
1017 double sp = DOT_PRODUCT(PART(COL(z, k), 1, n), PART(COL(a, kk), 1, n));
1018 temp = Math.sqrt(sp * sp + zdota[nact] * zdota[nact]);
1019 double alpha = zdota[nact] / temp;
1020 double beta = sp / temp;
1021 zdota[nact] = alpha * zdota[k];
1022 zdota[k] = temp;
1023 for (int i = 1; i <= n; ++i)
1024 {
1025 temp = alpha * z[i][nact] + beta * z[i][k];
1026 z[i][nact] = alpha * z[i][k] - beta * z[i][nact];
1027 z[i][k] = temp;
1028 }
1029 iact[nact] = iact[k];
1030 iact[k] = kk;
1031 temp = vmultc[k];
1032 vmultc[k] = vmultc[nact];
1033 vmultc[nact] = temp;
1034 }
1035
1036 // If stage one is in progress, then set SDIRN to the direction of the next
1037 // change to the current vector of variables.
1038
1039 if (mcon > m)
1040 {
1041 // Pick the next search direction of stage two.
1042
1043 temp = 1.0 / zdota[nact];
1044 for (int k = 1; k <= n; ++k) sdirn[k] = temp * z[k][nact];
1045 }
1046 else
1047 {
1048 kk = iact[nact];
1049 temp = (DOT_PRODUCT(PART(sdirn, 1, n), PART(COL(a, kk), 1, n)) - 1.0) / zdota[nact];
1050 for (int k = 1; k <= n; ++k) sdirn[k] -= temp * z[k][nact];
1051 }
1052 }
1053
1054 // Calculate the step to the boundary of the trust region or take the step
1055 // that reduces RESMAX to zero. The two statements below that include the
1056 // factor 1.0E-6 prevent some harmless underflows that occurred in a test
1057 // calculation. Further, we skip the step if it could be zero within a
1058 // reasonable tolerance for computer rounding errors.
1059
1060 double dd = rho * rho;
1061 double sd = 0.0;
1062 double ss = 0.0;
1063 for (int i = 1; i <= n; ++i)
1064 {
1065 if (Math.abs(dx[i]) >= 1.0E-6 * rho) dd -= dx[i] * dx[i];
1066 sd += dx[i] * sdirn[i];
1067 ss += sdirn[i] * sdirn[i];
1068 }
1069 if (dd <= 0.0) break L_60;
1070 temp = Math.sqrt(ss * dd);
1071 if (Math.abs(sd) >= 1.0E-6 * temp) temp = Math.sqrt(ss * dd + sd * sd);
1072 stpful = dd / (temp + sd);
1073 step = stpful;
1074 if (mcon == m)
1075 {
1076 double acca = step + 0.1 * resmax;
1077 double accb = step + 0.2 * resmax;
1078 if (step >= acca || acca >= accb) break L_70;
1079 step = Math.min(step, resmax);
1080 }
1081
1082 // Set DXNEW to the new variables if STEP is the steplength, and reduce
1083 // RESMAX to the corresponding maximum residual if stage one is being done.
1084 // Because DXNEW will be changed during the calculation of some Lagrange
1085 // multipliers, it will be restored to the following value later.
1086
1087 for (int k = 1; k <= n; ++k) dxnew[k] = dx[k] + step * sdirn[k];
1088 if (mcon == m)
1089 {
1090 resold = resmax;
1091 resmax = 0.0;
1092 for (int k = 1; k <= nact; ++k)
1093 {
1094 int kk = iact[k];
1095 temp = b[kk] - DOT_PRODUCT(PART(COL(a, kk), 1, n), PART(dxnew, 1, n));
1096 resmax = Math.max(resmax, temp);
1097 }
1098 }
1099
1100 // Set VMULTD to the VMULTC vector that would occur if DX became DXNEW. A
1101 // device is included to force VMULTD(K) = 0.0 if deviations from this value
1102 // can be attributed to computer rounding errors. First calculate the new
1103 // Lagrange multipliers.
1104
1105 {
1106 int k = nact;
1107 do
1108 {
1109 double zdotw = 0.0;
1110 double zdwabs = 0.0;
1111 for (int i = 1; i <= n; ++i)
1112 {
1113 temp = z[i][k] * dxnew[i];
1114 zdotw += temp;
1115 zdwabs += Math.abs(temp);
1116 }
1117 double acca = zdwabs + 0.1 * Math.abs(zdotw);
1118 double accb = zdwabs + 0.2 * Math.abs(zdotw);
1119 if (zdwabs >= acca || acca >= accb) zdotw = 0.0;
1120 vmultd[k] = zdotw / zdota[k];
1121 if (k >= 2)
1122 {
1123 int kk = iact[k];
1124 for (int i = 1; i <= n; ++i) dxnew[i] -= vmultd[k] * a[i][kk];
1125 }
1126 } while (k-- >= 2);
1127 if (mcon > m) vmultd[nact] = Math.max(0.0, vmultd[nact]);
1128 }
1129
1130 // Complete VMULTC by finding the new constraint residuals.
1131
1132 for (int k = 1; k <= n; ++k) dxnew[k] = dx[k] + step * sdirn[k];
1133 if (mcon > nact)
1134 {
1135 int kl = nact + 1;
1136 for (int k = kl; k <= mcon; ++k)
1137 {
1138 int kk = iact[k];
1139 double total = resmax - b[kk];
1140 double sumabs = resmax + Math.abs(b[kk]);
1141 for (int i = 1; i <= n; ++i)
1142 {
1143 temp = a[i][kk] * dxnew[i];
1144 total += temp;
1145 sumabs += Math.abs(temp);
1146 }
1147 double acca = sumabs + 0.1 * Math.abs(total);
1148 double accb = sumabs + 0.2 * Math.abs(total);
1149 if (sumabs >= acca || acca >= accb) total = 0.0;
1150 vmultd[k] = total;
1151 }
1152 }
1153
1154 // Calculate the fraction of the step from DX to DXNEW that will be taken.
1155
1156 ratio = 1.0;
1157 icon = 0;
1158 for (int k = 1; k <= mcon; ++k)
1159 {
1160 if (vmultd[k] < 0.0)
1161 {
1162 temp = vmultc[k] / (vmultc[k] - vmultd[k]);
1163 if (temp < ratio)
1164 {
1165 ratio = temp;
1166 icon = k;
1167 }
1168 }
1169 }
1170
1171 // Update DX, VMULTC and RESMAX.
1172
1173 temp = 1.0 - ratio;
1174 for (int k = 1; k <= n; ++k) dx[k] = temp * dx[k] + ratio * dxnew[k];
1175 for (int k = 1; k <= mcon; ++k)
1176 vmultc[k] = Math.max(0.0, temp * vmultc[k] + ratio * vmultd[k]);
1177 if (mcon == m) resmax = resold + ratio * (resmax - resold);
1178
1179 // If the full step is not acceptable then begin another iteration.
1180 // Otherwise switch to stage two or end the calculation.
1181
1182 } while (icon > 0);
1183
1184 if (step == stpful) return true;
1185
1186 } while (true);
1187
1188 // We employ any freedom that may be available to reduce the objective
1189 // function before returning a DX whose length is less than RHO.
1190
1191 } while (mcon == m);
1192
1193 return false;
1194 }
1195
1196 private static void PrintIterationResult(int nfvals, double f, double resmax, double[] x, int n)
1197 {
1198 System.out.format("%nNFVALS = %1$5d F = %2$13.6f MAXCV = %3$13.6e%n", nfvals, f, resmax);
1199 System.out.format("X = %s%n", FORMAT(PART(x, 1, n)));
1200 }
1201
1202 private static double[] ROW(double[][] src, int rowidx)
1203 {
1204 int cols = src[0].length;
1205 double[] dest = new double[cols];
1206 for (int col = 0; col < cols; ++col) dest[col] = src[rowidx][col];
1207 return dest;
1208 }
1209
1210 private static double[] COL(double[][] src, int colidx)
1211 {
1212 int rows = src.length;
1213 double[] dest = new double[rows];
1214 for (int row = 0; row < rows; ++row) dest[row] = src[row][colidx];
1215 return dest;
1216 }
1217
1218 private static double[] PART(double[] src, int from, int to)
1219 {
1220 double[] dest = new double[to - from + 1];
1221 int destidx = 0;
1222 for (int srcidx = from; srcidx <= to; ++srcidx, ++destidx) dest[destidx] = src[srcidx];
1223 return dest;
1224 }
1225
1226 private static String FORMAT(double[] x)
1227 {
1228 String fmt = "";
1229 for (int i = 0; i < x.length; ++i) fmt = fmt + String.format("%13.6f", x[i]);
1230 return fmt;
1231 }
1232
1233 private static double DOT_PRODUCT(double[] lhs, double[] rhs)
1234 {
1235 double sum = 0.0; for (int i = 0; i < lhs.length; ++i) sum += lhs[i] * rhs[i];
1236 return sum;
1237 }
1238}
Note: See TracBrowser for help on using the repository browser.