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