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