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