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