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