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