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