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