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