Changeset 2430 for trunk/sources/ALGLIB/leastsquares.cs
- Timestamp:
- 10/15/09 13:22:07 (15 years ago)
- Location:
- trunk/sources/ALGLIB
- Files:
-
- 1 added
- 1 copied
Legend:
- Unmodified
- Added
- Removed
-
trunk/sources/ALGLIB/leastsquares.cs
r2426 r2430 2 2 Copyright (c) 2006-2007, Sergey Bochkanov (ALGLIB project). 3 3 4 Redistribution and use in source and binary forms, with or without 5 modification, are permitted provided that the following conditions are 6 met: 7 8 - Redistributions of source code must retain the above copyright 9 notice, this list of conditions and the following disclaimer. 10 11 - Redistributions in binary form must reproduce the above copyright 12 notice, this list of conditions and the following disclaimer listed 13 in this license in the documentation and/or other materials 14 provided with the distribution. 15 16 - Neither the name of the copyright holders nor the names of its 17 contributors may be used to endorse or promote products derived from 18 this software without specific prior written permission. 19 20 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 21 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 22 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR 23 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT 24 OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, 25 SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 26 LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 27 DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 28 THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 29 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 30 OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 4 >>> SOURCE LICENSE >>> 5 This program is free software; you can redistribute it and/or modify 6 it under the terms of the GNU General Public License as published by 7 the Free Software Foundation (www.fsf.org); either version 2 of the 8 License, or (at your option) any later version. 9 10 This program is distributed in the hope that it will be useful, 11 but WITHOUT ANY WARRANTY; without even the implied warranty of 12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13 GNU General Public License for more details. 14 15 A copy of the GNU General Public License is available at 16 http://www.fsf.org/licensing/licenses 17 18 >>> END OF LICENSE >>> 31 19 *************************************************************************/ 32 20 33 21 using System; 34 22 35 class leastsquares 23 namespace alglib 36 24 { 37 /************************************************************************* 38 Weighted approximation by arbitrary function basis in a space of arbitrary 39 dimension using linear least squares method. 40 41 Input parameters: 42 Y - array[0..N-1] 43 It contains a set of function values in N points. Space 44 dimension and points don't matter. Procedure works with 45 function values in these points and values of basis functions 46 only. 47 48 W - array[0..N-1] 49 It contains weights corresponding to function values. Each 50 summand in square sum of approximation deviations from given 51 values is multiplied by the square of corresponding weight. 52 53 FMatrix-a table of basis functions values, array[0..N-1, 0..M-1]. 54 FMatrix[I, J] - value of J-th basis function in I-th point. 55 56 N - number of points used. N>=1. 57 M - number of basis functions, M>=1. 58 59 Output parameters: 60 C - decomposition coefficients. 61 Array of real numbers whose index goes from 0 to M-1. 62 C[j] - j-th basis function coefficient. 63 64 -- ALGLIB -- 65 Copyright by Bochkanov Sergey 66 *************************************************************************/ 67 public static void buildgeneralleastsquares(ref double[] y, 68 ref double[] w, 69 ref double[,] fmatrix, 70 int n, 71 int m, 72 ref double[] c) 25 public class leastsquares 73 26 { 74 int i = 0; 75 int j = 0; 76 double[,] a = new double[0,0]; 77 double[,] q = new double[0,0]; 78 double[,] vt = new double[0,0]; 79 double[] b = new double[0]; 80 double[] tau = new double[0]; 81 double[,] b2 = new double[0,0]; 82 double[] tauq = new double[0]; 83 double[] taup = new double[0]; 84 double[] d = new double[0]; 85 double[] e = new double[0]; 86 bool isuppera = new bool(); 87 int mi = 0; 88 int ni = 0; 89 double v = 0; 90 int i_ = 0; 91 int i1_ = 0; 92 93 mi = n; 94 ni = m; 95 c = new double[ni-1+1]; 96 97 // 98 // Initialize design matrix. 99 // Here we are making MI>=NI. 100 // 101 a = new double[ni+1, Math.Max(mi, ni)+1]; 102 b = new double[Math.Max(mi, ni)+1]; 103 for(i=1; i<=mi; i++) 27 /************************************************************************* 28 Weighted approximation by arbitrary function basis in a space of arbitrary 29 dimension using linear least squares method. 30 31 Input parameters: 32 Y - array[0..N-1] 33 It contains a set of function values in N points. Space 34 dimension and points don't matter. Procedure works with 35 function values in these points and values of basis functions 36 only. 37 38 W - array[0..N-1] 39 It contains weights corresponding to function values. Each 40 summand in square sum of approximation deviations from given 41 values is multiplied by the square of corresponding weight. 42 43 FMatrix-a table of basis functions values, array[0..N-1, 0..M-1]. 44 FMatrix[I, J] - value of J-th basis function in I-th point. 45 46 N - number of points used. N>=1. 47 M - number of basis functions, M>=1. 48 49 Output parameters: 50 C - decomposition coefficients. 51 Array of real numbers whose index goes from 0 to M-1. 52 C[j] - j-th basis function coefficient. 53 54 -- ALGLIB -- 55 Copyright by Bochkanov Sergey 56 *************************************************************************/ 57 public static void buildgeneralleastsquares(ref double[] y, 58 ref double[] w, 59 ref double[,] fmatrix, 60 int n, 61 int m, 62 ref double[] c) 104 63 { 105 b[i] = w[i-1]*y[i-1]; 64 int i = 0; 65 int j = 0; 66 double[,] a = new double[0,0]; 67 double[,] q = new double[0,0]; 68 double[,] vt = new double[0,0]; 69 double[] b = new double[0]; 70 double[] tau = new double[0]; 71 double[,] b2 = new double[0,0]; 72 double[] tauq = new double[0]; 73 double[] taup = new double[0]; 74 double[] d = new double[0]; 75 double[] e = new double[0]; 76 bool isuppera = new bool(); 77 int mi = 0; 78 int ni = 0; 79 double v = 0; 80 int i_ = 0; 81 int i1_ = 0; 82 83 mi = n; 84 ni = m; 85 c = new double[ni-1+1]; 86 87 // 88 // Initialize design matrix. 89 // Here we are making MI>=NI. 90 // 91 a = new double[ni+1, Math.Max(mi, ni)+1]; 92 b = new double[Math.Max(mi, ni)+1]; 93 for(i=1; i<=mi; i++) 94 { 95 b[i] = w[i-1]*y[i-1]; 96 } 97 for(i=mi+1; i<=ni; i++) 98 { 99 b[i] = 0; 100 } 101 for(j=1; j<=ni; j++) 102 { 103 i1_ = (0) - (1); 104 for(i_=1; i_<=mi;i_++) 105 { 106 a[j,i_] = fmatrix[i_+i1_,j-1]; 107 } 108 } 109 for(j=1; j<=ni; j++) 110 { 111 for(i=mi+1; i<=ni; i++) 112 { 113 a[j,i] = 0; 114 } 115 } 116 for(j=1; j<=ni; j++) 117 { 118 for(i=1; i<=mi; i++) 119 { 120 a[j,i] = a[j,i]*w[i-1]; 121 } 122 } 123 mi = Math.Max(mi, ni); 124 125 // 126 // LQ-decomposition of A' 127 // B2 := Q*B 128 // 129 lq.lqdecomposition(ref a, ni, mi, ref tau); 130 lq.unpackqfromlq(ref a, ni, mi, ref tau, ni, ref q); 131 b2 = new double[1+1, ni+1]; 132 for(j=1; j<=ni; j++) 133 { 134 b2[1,j] = 0; 135 } 136 for(i=1; i<=ni; i++) 137 { 138 v = 0.0; 139 for(i_=1; i_<=mi;i_++) 140 { 141 v += b[i_]*q[i,i_]; 142 } 143 b2[1,i] = v; 144 } 145 146 // 147 // Back from A' to A 148 // Making cols(A)=rows(A) 149 // 150 for(i=1; i<=ni-1; i++) 151 { 152 for(i_=i+1; i_<=ni;i_++) 153 { 154 a[i,i_] = a[i_,i]; 155 } 156 } 157 for(i=2; i<=ni; i++) 158 { 159 for(j=1; j<=i-1; j++) 160 { 161 a[i,j] = 0; 162 } 163 } 164 165 // 166 // Bidiagonal decomposition of A 167 // A = Q * d2 * P' 168 // B2 := (Q'*B2')' 169 // 170 bidiagonal.tobidiagonal(ref a, ni, ni, ref tauq, ref taup); 171 bidiagonal.multiplybyqfrombidiagonal(ref a, ni, ni, ref tauq, ref b2, 1, ni, true, false); 172 bidiagonal.unpackptfrombidiagonal(ref a, ni, ni, ref taup, ni, ref vt); 173 bidiagonal.unpackdiagonalsfrombidiagonal(ref a, ni, ni, ref isuppera, ref d, ref e); 174 175 // 176 // Singular value decomposition of A 177 // A = U * d * V' 178 // B2 := (U'*B2')' 179 // 180 if( !bdsvd.bidiagonalsvddecomposition(ref d, e, ni, isuppera, false, ref b2, 1, ref q, 0, ref vt, ni) ) 181 { 182 for(i=0; i<=ni-1; i++) 183 { 184 c[i] = 0; 185 } 186 return; 187 } 188 189 // 190 // B2 := (d^(-1) * B2')' 191 // 192 if( d[1]!=0 ) 193 { 194 for(i=1; i<=ni; i++) 195 { 196 if( d[i]>AP.Math.MachineEpsilon*10*Math.Sqrt(ni)*d[1] ) 197 { 198 b2[1,i] = b2[1,i]/d[i]; 199 } 200 else 201 { 202 b2[1,i] = 0; 203 } 204 } 205 } 206 207 // 208 // B := (V * B2')' 209 // 210 for(i=1; i<=ni; i++) 211 { 212 b[i] = 0; 213 } 214 for(i=1; i<=ni; i++) 215 { 216 v = b2[1,i]; 217 for(i_=1; i_<=ni;i_++) 218 { 219 b[i_] = b[i_] + v*vt[i,i_]; 220 } 221 } 222 223 // 224 // Out 225 // 226 for(i=1; i<=ni; i++) 227 { 228 c[i-1] = b[i]; 229 } 106 230 } 107 for(i=mi+1; i<=ni; i++) 231 232 233 /************************************************************************* 234 Linear approximation using least squares method 235 236 The subroutine calculates coefficients of the line approximating given 237 function. 238 239 Input parameters: 240 X - array[0..N-1], it contains a set of abscissas. 241 Y - array[0..N-1], function values. 242 N - number of points, N>=1 243 244 Output parameters: 245 a, b- coefficients of linear approximation a+b*t 246 247 -- ALGLIB -- 248 Copyright by Bochkanov Sergey 249 *************************************************************************/ 250 public static void buildlinearleastsquares(ref double[] x, 251 ref double[] y, 252 int n, 253 ref double a, 254 ref double b) 108 255 { 109 b[i] = 0; 256 double pp = 0; 257 double qq = 0; 258 double pq = 0; 259 double b1 = 0; 260 double b2 = 0; 261 double d1 = 0; 262 double d2 = 0; 263 double t1 = 0; 264 double t2 = 0; 265 double phi = 0; 266 double c = 0; 267 double s = 0; 268 double m = 0; 269 int i = 0; 270 271 pp = n; 272 qq = 0; 273 pq = 0; 274 b1 = 0; 275 b2 = 0; 276 for(i=0; i<=n-1; i++) 277 { 278 pq = pq+x[i]; 279 qq = qq+AP.Math.Sqr(x[i]); 280 b1 = b1+y[i]; 281 b2 = b2+x[i]*y[i]; 282 } 283 phi = Math.Atan2(2*pq, qq-pp)/2; 284 c = Math.Cos(phi); 285 s = Math.Sin(phi); 286 d1 = AP.Math.Sqr(c)*pp+AP.Math.Sqr(s)*qq-2*s*c*pq; 287 d2 = AP.Math.Sqr(s)*pp+AP.Math.Sqr(c)*qq+2*s*c*pq; 288 if( Math.Abs(d1)>Math.Abs(d2) ) 289 { 290 m = Math.Abs(d1); 291 } 292 else 293 { 294 m = Math.Abs(d2); 295 } 296 t1 = c*b1-s*b2; 297 t2 = s*b1+c*b2; 298 if( Math.Abs(d1)>m*AP.Math.MachineEpsilon*1000 ) 299 { 300 t1 = t1/d1; 301 } 302 else 303 { 304 t1 = 0; 305 } 306 if( Math.Abs(d2)>m*AP.Math.MachineEpsilon*1000 ) 307 { 308 t2 = t2/d2; 309 } 310 else 311 { 312 t2 = 0; 313 } 314 a = c*t1+s*t2; 315 b = -(s*t1)+c*t2; 110 316 } 111 for(j=1; j<=ni; j++) 317 318 319 /************************************************************************* 320 Weighted cubic spline approximation using linear least squares 321 322 Input parameters: 323 X - array[0..N-1], abscissas 324 Y - array[0..N-1], function values 325 W - array[0..N-1], weights. 326 A, B- interval to build splines in. 327 N - number of points used. N>=1. 328 M - number of basic splines, M>=2. 329 330 Output parameters: 331 CTbl- coefficients table to be used by SplineInterpolation function. 332 -- ALGLIB -- 333 Copyright by Bochkanov Sergey 334 *************************************************************************/ 335 public static void buildsplineleastsquares(ref double[] x, 336 ref double[] y, 337 ref double[] w, 338 double a, 339 double b, 340 int n, 341 int m, 342 ref double[] ctbl) 112 343 { 113 i1_ = (0) - (1); 114 for(i_=1; i_<=mi;i_++) 115 { 116 a[j,i_] = fmatrix[i_+i1_,j-1]; 117 } 344 int i = 0; 345 int j = 0; 346 double[,] ma = new double[0,0]; 347 double[,] q = new double[0,0]; 348 double[,] vt = new double[0,0]; 349 double[] mb = new double[0]; 350 double[] tau = new double[0]; 351 double[,] b2 = new double[0,0]; 352 double[] tauq = new double[0]; 353 double[] taup = new double[0]; 354 double[] d = new double[0]; 355 double[] e = new double[0]; 356 bool isuppera = new bool(); 357 int mi = 0; 358 int ni = 0; 359 double v = 0; 360 double[] sx = new double[0]; 361 double[] sy = new double[0]; 362 int i_ = 0; 363 364 System.Diagnostics.Debug.Assert(m>=2, "BuildSplineLeastSquares: M is too small!"); 365 mi = n; 366 ni = m; 367 sx = new double[ni-1+1]; 368 sy = new double[ni-1+1]; 369 370 // 371 // Initializing design matrix 372 // Here we are making MI>=NI 373 // 374 ma = new double[ni+1, Math.Max(mi, ni)+1]; 375 mb = new double[Math.Max(mi, ni)+1]; 376 for(j=0; j<=ni-1; j++) 377 { 378 sx[j] = a+(b-a)*j/(ni-1); 379 } 380 for(j=0; j<=ni-1; j++) 381 { 382 for(i=0; i<=ni-1; i++) 383 { 384 sy[i] = 0; 385 } 386 sy[j] = 1; 387 spline3.buildcubicspline(sx, sy, ni, 0, 0.0, 0, 0.0, ref ctbl); 388 for(i=0; i<=mi-1; i++) 389 { 390 ma[j+1,i+1] = w[i]*spline3.splineinterpolation(ref ctbl, x[i]); 391 } 392 } 393 for(j=1; j<=ni; j++) 394 { 395 for(i=mi+1; i<=ni; i++) 396 { 397 ma[j,i] = 0; 398 } 399 } 400 401 // 402 // Initializing right part 403 // 404 for(i=0; i<=mi-1; i++) 405 { 406 mb[i+1] = w[i]*y[i]; 407 } 408 for(i=mi+1; i<=ni; i++) 409 { 410 mb[i] = 0; 411 } 412 mi = Math.Max(mi, ni); 413 414 // 415 // LQ-decomposition of A' 416 // B2 := Q*B 417 // 418 lq.lqdecomposition(ref ma, ni, mi, ref tau); 419 lq.unpackqfromlq(ref ma, ni, mi, ref tau, ni, ref q); 420 b2 = new double[1+1, ni+1]; 421 for(j=1; j<=ni; j++) 422 { 423 b2[1,j] = 0; 424 } 425 for(i=1; i<=ni; i++) 426 { 427 v = 0.0; 428 for(i_=1; i_<=mi;i_++) 429 { 430 v += mb[i_]*q[i,i_]; 431 } 432 b2[1,i] = v; 433 } 434 435 // 436 // Back from A' to A 437 // Making cols(A)=rows(A) 438 // 439 for(i=1; i<=ni-1; i++) 440 { 441 for(i_=i+1; i_<=ni;i_++) 442 { 443 ma[i,i_] = ma[i_,i]; 444 } 445 } 446 for(i=2; i<=ni; i++) 447 { 448 for(j=1; j<=i-1; j++) 449 { 450 ma[i,j] = 0; 451 } 452 } 453 454 // 455 // Bidiagonal decomposition of A 456 // A = Q * d2 * P' 457 // B2 := (Q'*B2')' 458 // 459 bidiagonal.tobidiagonal(ref ma, ni, ni, ref tauq, ref taup); 460 bidiagonal.multiplybyqfrombidiagonal(ref ma, ni, ni, ref tauq, ref b2, 1, ni, true, false); 461 bidiagonal.unpackptfrombidiagonal(ref ma, ni, ni, ref taup, ni, ref vt); 462 bidiagonal.unpackdiagonalsfrombidiagonal(ref ma, ni, ni, ref isuppera, ref d, ref e); 463 464 // 465 // Singular value decomposition of A 466 // A = U * d * V' 467 // B2 := (U'*B2')' 468 // 469 if( !bdsvd.bidiagonalsvddecomposition(ref d, e, ni, isuppera, false, ref b2, 1, ref q, 0, ref vt, ni) ) 470 { 471 for(i=1; i<=ni; i++) 472 { 473 d[i] = 0; 474 b2[1,i] = 0; 475 for(j=1; j<=ni; j++) 476 { 477 if( i==j ) 478 { 479 vt[i,j] = 1; 480 } 481 else 482 { 483 vt[i,j] = 0; 484 } 485 } 486 } 487 b2[1,1] = 1; 488 } 489 490 // 491 // B2 := (d^(-1) * B2')' 492 // 493 for(i=1; i<=ni; i++) 494 { 495 if( d[i]>AP.Math.MachineEpsilon*10*Math.Sqrt(ni)*d[1] ) 496 { 497 b2[1,i] = b2[1,i]/d[i]; 498 } 499 else 500 { 501 b2[1,i] = 0; 502 } 503 } 504 505 // 506 // B := (V * B2')' 507 // 508 for(i=1; i<=ni; i++) 509 { 510 mb[i] = 0; 511 } 512 for(i=1; i<=ni; i++) 513 { 514 v = b2[1,i]; 515 for(i_=1; i_<=ni;i_++) 516 { 517 mb[i_] = mb[i_] + v*vt[i,i_]; 518 } 519 } 520 521 // 522 // Forming result spline 523 // 524 for(i=0; i<=ni-1; i++) 525 { 526 sy[i] = mb[i+1]; 527 } 528 spline3.buildcubicspline(sx, sy, ni, 0, 0.0, 0, 0.0, ref ctbl); 118 529 } 119 for(j=1; j<=ni; j++) 530 531 532 /************************************************************************* 533 Polynomial approximation using least squares method 534 535 The subroutine calculates coefficients of the polynomial approximating 536 given function. It is recommended to use this function only if you need to 537 obtain coefficients of approximation polynomial. If you have to build and 538 calculate polynomial approximation (NOT coefficients), it's better to use 539 BuildChebyshevLeastSquares subroutine in combination with 540 CalculateChebyshevLeastSquares subroutine. The result of Chebyshev 541 polynomial approximation is equivalent to the result obtained using powers 542 of X, but has higher accuracy due to better numerical properties of 543 Chebyshev polynomials. 544 545 Input parameters: 546 X - array[0..N-1], abscissas 547 Y - array[0..N-1], function values 548 N - number of points, N>=1 549 M - order of polynomial required, M>=0 550 551 Output parameters: 552 C - approximating polynomial coefficients, array[0..M], 553 C[i] - coefficient at X^i. 554 555 -- ALGLIB -- 556 Copyright by Bochkanov Sergey 557 *************************************************************************/ 558 public static void buildpolynomialleastsquares(ref double[] x, 559 ref double[] y, 560 int n, 561 int m, 562 ref double[] c) 120 563 { 564 double[] ctbl = new double[0]; 565 double[] w = new double[0]; 566 double[] c1 = new double[0]; 567 double maxx = 0; 568 double minx = 0; 569 int i = 0; 570 int j = 0; 571 int k = 0; 572 double e = 0; 573 double d = 0; 574 double l1 = 0; 575 double l2 = 0; 576 double[] z2 = new double[0]; 577 double[] z1 = new double[0]; 578 579 580 // 581 // Initialize 582 // 583 maxx = x[0]; 584 minx = x[0]; 585 for(i=1; i<=n-1; i++) 586 { 587 if( x[i]>maxx ) 588 { 589 maxx = x[i]; 590 } 591 if( x[i]<minx ) 592 { 593 minx = x[i]; 594 } 595 } 596 if( minx==maxx ) 597 { 598 minx = minx-0.5; 599 maxx = maxx+0.5; 600 } 601 w = new double[n-1+1]; 602 for(i=0; i<=n-1; i++) 603 { 604 w[i] = 1; 605 } 606 607 // 608 // Build Chebyshev approximation 609 // 610 buildchebyshevleastsquares(ref x, ref y, ref w, minx, maxx, n, m, ref ctbl); 611 612 // 613 // From Chebyshev to powers of X 614 // 615 c1 = new double[m+1]; 616 for(i=0; i<=m; i++) 617 { 618 c1[i] = 0; 619 } 620 d = 0; 621 for(i=0; i<=m; i++) 622 { 623 for(k=i; k<=m; k++) 624 { 625 e = c1[k]; 626 c1[k] = 0; 627 if( i<=1 & k==i ) 628 { 629 c1[k] = 1; 630 } 631 else 632 { 633 if( i!=0 ) 634 { 635 c1[k] = 2*d; 636 } 637 if( k>i+1 ) 638 { 639 c1[k] = c1[k]-c1[k-2]; 640 } 641 } 642 d = e; 643 } 644 d = c1[i]; 645 e = 0; 646 k = i; 647 while( k<=m ) 648 { 649 e = e+c1[k]*ctbl[k]; 650 k = k+2; 651 } 652 c1[i] = e; 653 } 654 655 // 656 // Linear translation 657 // 658 l1 = 2/(ctbl[m+2]-ctbl[m+1]); 659 l2 = -(2*ctbl[m+1]/(ctbl[m+2]-ctbl[m+1]))-1; 660 c = new double[m+1]; 661 z2 = new double[m+1]; 662 z1 = new double[m+1]; 663 c[0] = c1[0]; 664 z1[0] = 1; 665 z2[0] = 1; 666 for(i=1; i<=m; i++) 667 { 668 z2[i] = 1; 669 z1[i] = l2*z1[i-1]; 670 c[0] = c[0]+c1[i]*z1[i]; 671 } 672 for(j=1; j<=m; j++) 673 { 674 z2[0] = l1*z2[0]; 675 c[j] = c1[j]*z2[0]; 676 for(i=j+1; i<=m; i++) 677 { 678 k = i-j; 679 z2[k] = l1*z2[k]+z2[k-1]; 680 c[j] = c[j]+c1[i]*z2[k]*z1[k]; 681 } 682 } 683 } 684 685 686 /************************************************************************* 687 Chebyshev polynomial approximation using least squares method. 688 689 The algorithm reduces interval [A, B] to the interval [-1,1], then builds 690 least squares approximation using Chebyshev polynomials. 691 692 Input parameters: 693 X - array[0..N-1], abscissas 694 Y - array[0..N-1], function values 695 W - array[0..N-1], weights 696 A, B- interval to build approximating polynomials in. 697 N - number of points used. N>=1. 698 M - order of polynomial, M>=0. This parameter is passed into 699 CalculateChebyshevLeastSquares function. 700 701 Output parameters: 702 CTbl - coefficient table. This parameter is passed into 703 CalculateChebyshevLeastSquares function. 704 -- ALGLIB -- 705 Copyright by Bochkanov Sergey 706 *************************************************************************/ 707 public static void buildchebyshevleastsquares(ref double[] x, 708 ref double[] y, 709 ref double[] w, 710 double a, 711 double b, 712 int n, 713 int m, 714 ref double[] ctbl) 715 { 716 int i = 0; 717 int j = 0; 718 double[,] ma = new double[0,0]; 719 double[,] q = new double[0,0]; 720 double[,] vt = new double[0,0]; 721 double[] mb = new double[0]; 722 double[] tau = new double[0]; 723 double[,] b2 = new double[0,0]; 724 double[] tauq = new double[0]; 725 double[] taup = new double[0]; 726 double[] d = new double[0]; 727 double[] e = new double[0]; 728 bool isuppera = new bool(); 729 int mi = 0; 730 int ni = 0; 731 double v = 0; 732 int i_ = 0; 733 734 mi = n; 735 ni = m+1; 736 737 // 738 // Initializing design matrix 739 // Here we are making MI>=NI 740 // 741 ma = new double[ni+1, Math.Max(mi, ni)+1]; 742 mb = new double[Math.Max(mi, ni)+1]; 743 for(j=1; j<=ni; j++) 744 { 745 for(i=1; i<=mi; i++) 746 { 747 v = 2*(x[i-1]-a)/(b-a)-1; 748 if( j==1 ) 749 { 750 ma[j,i] = 1.0; 751 } 752 if( j==2 ) 753 { 754 ma[j,i] = v; 755 } 756 if( j>2 ) 757 { 758 ma[j,i] = 2.0*v*ma[j-1,i]-ma[j-2,i]; 759 } 760 } 761 } 762 for(j=1; j<=ni; j++) 763 { 764 for(i=1; i<=mi; i++) 765 { 766 ma[j,i] = w[i-1]*ma[j,i]; 767 } 768 } 769 for(j=1; j<=ni; j++) 770 { 771 for(i=mi+1; i<=ni; i++) 772 { 773 ma[j,i] = 0; 774 } 775 } 776 777 // 778 // Initializing right part 779 // 780 for(i=0; i<=mi-1; i++) 781 { 782 mb[i+1] = w[i]*y[i]; 783 } 121 784 for(i=mi+1; i<=ni; i++) 122 785 { 123 a[j,i] = 0; 124 } 786 mb[i] = 0; 787 } 788 mi = Math.Max(mi, ni); 789 790 // 791 // LQ-decomposition of A' 792 // B2 := Q*B 793 // 794 lq.lqdecomposition(ref ma, ni, mi, ref tau); 795 lq.unpackqfromlq(ref ma, ni, mi, ref tau, ni, ref q); 796 b2 = new double[1+1, ni+1]; 797 for(j=1; j<=ni; j++) 798 { 799 b2[1,j] = 0; 800 } 801 for(i=1; i<=ni; i++) 802 { 803 v = 0.0; 804 for(i_=1; i_<=mi;i_++) 805 { 806 v += mb[i_]*q[i,i_]; 807 } 808 b2[1,i] = v; 809 } 810 811 // 812 // Back from A' to A 813 // Making cols(A)=rows(A) 814 // 815 for(i=1; i<=ni-1; i++) 816 { 817 for(i_=i+1; i_<=ni;i_++) 818 { 819 ma[i,i_] = ma[i_,i]; 820 } 821 } 822 for(i=2; i<=ni; i++) 823 { 824 for(j=1; j<=i-1; j++) 825 { 826 ma[i,j] = 0; 827 } 828 } 829 830 // 831 // Bidiagonal decomposition of A 832 // A = Q * d2 * P' 833 // B2 := (Q'*B2')' 834 // 835 bidiagonal.tobidiagonal(ref ma, ni, ni, ref tauq, ref taup); 836 bidiagonal.multiplybyqfrombidiagonal(ref ma, ni, ni, ref tauq, ref b2, 1, ni, true, false); 837 bidiagonal.unpackptfrombidiagonal(ref ma, ni, ni, ref taup, ni, ref vt); 838 bidiagonal.unpackdiagonalsfrombidiagonal(ref ma, ni, ni, ref isuppera, ref d, ref e); 839 840 // 841 // Singular value decomposition of A 842 // A = U * d * V' 843 // B2 := (U'*B2')' 844 // 845 if( !bdsvd.bidiagonalsvddecomposition(ref d, e, ni, isuppera, false, ref b2, 1, ref q, 0, ref vt, ni) ) 846 { 847 for(i=1; i<=ni; i++) 848 { 849 d[i] = 0; 850 b2[1,i] = 0; 851 for(j=1; j<=ni; j++) 852 { 853 if( i==j ) 854 { 855 vt[i,j] = 1; 856 } 857 else 858 { 859 vt[i,j] = 0; 860 } 861 } 862 } 863 b2[1,1] = 1; 864 } 865 866 // 867 // B2 := (d^(-1) * B2')' 868 // 869 for(i=1; i<=ni; i++) 870 { 871 if( d[i]>AP.Math.MachineEpsilon*10*Math.Sqrt(ni)*d[1] ) 872 { 873 b2[1,i] = b2[1,i]/d[i]; 874 } 875 else 876 { 877 b2[1,i] = 0; 878 } 879 } 880 881 // 882 // B := (V * B2')' 883 // 884 for(i=1; i<=ni; i++) 885 { 886 mb[i] = 0; 887 } 888 for(i=1; i<=ni; i++) 889 { 890 v = b2[1,i]; 891 for(i_=1; i_<=ni;i_++) 892 { 893 mb[i_] = mb[i_] + v*vt[i,i_]; 894 } 895 } 896 897 // 898 // Forming result 899 // 900 ctbl = new double[ni+1+1]; 901 for(i=1; i<=ni; i++) 902 { 903 ctbl[i-1] = mb[i]; 904 } 905 ctbl[ni] = a; 906 ctbl[ni+1] = b; 125 907 } 126 for(j=1; j<=ni; j++) 908 909 910 /************************************************************************* 911 Weighted Chebyshev polynomial constrained least squares approximation. 912 913 The algorithm reduces [A,B] to [-1,1] and builds the Chebyshev polynomials 914 series by approximating a given function using the least squares method. 915 916 Input parameters: 917 X - abscissas, array[0..N-1] 918 Y - function values, array[0..N-1] 919 W - weights, array[0..N-1]. Each item in the squared sum of 920 deviations from given values is multiplied by a square of 921 corresponding weight. 922 A, B- interval in which the approximating polynomials are built. 923 N - number of points, N>0. 924 XC, YC, DC- 925 constraints (see description below)., array[0..NC-1] 926 NC - number of constraints. 0 <= NC < M+1. 927 M - degree of polynomial, M>=0. This parameter is passed into the 928 CalculateChebyshevLeastSquares subroutine. 929 930 Output parameters: 931 CTbl- coefficient table. This parameter is passed into the 932 CalculateChebyshevLeastSquares subroutine. 933 934 Result: 935 True, if the algorithm succeeded. 936 False, if the internal singular value decomposition subroutine hasn't 937 converged or the given constraints could not be met simultaneously (e.g. 938 P(0)=0 è P(0)=1). 939 940 Specifying constraints: 941 This subroutine can solve the problem having constrained function 942 values or its derivatives in several points. NC specifies the number of 943 constraints, DC - the type of constraints, XC and YC - constraints as such. 944 Thus, for each i from 0 to NC-1 the following constraint is given: 945 P(xc[i]) = yc[i], if DC[i]=0 946 or 947 d/dx(P(xc[i])) = yc[i], if DC[i]=1 948 (here P(x) is approximating polynomial). 949 This version of the subroutine supports only either polynomial or its 950 derivative value constraints. If DC[i] is not equal to 0 and 1, the 951 subroutine will be aborted. The number of constraints should be less than 952 the number of degrees of freedom of approximating polynomial - M+1 (at 953 that, it could be equal to 0). 954 955 -- ALGLIB -- 956 Copyright by Bochkanov Sergey 957 *************************************************************************/ 958 public static bool buildchebyshevleastsquaresconstrained(ref double[] x, 959 ref double[] y, 960 ref double[] w, 961 double a, 962 double b, 963 int n, 964 ref double[] xc, 965 ref double[] yc, 966 ref int[] dc, 967 int nc, 968 int m, 969 ref double[] ctbl) 127 970 { 128 for(i=1; i<=mi; i++) 129 { 130 a[j,i] = a[j,i]*w[i-1]; 131 } 971 bool result = new bool(); 972 int i = 0; 973 int j = 0; 974 int reducedsize = 0; 975 double[,] designmatrix = new double[0,0]; 976 double[] rightpart = new double[0]; 977 double[,] cmatrix = new double[0,0]; 978 double[,] c = new double[0,0]; 979 double[,] u = new double[0,0]; 980 double[,] vt = new double[0,0]; 981 double[] d = new double[0]; 982 double[] cr = new double[0]; 983 double[] ws = new double[0]; 984 double[] tj = new double[0]; 985 double[] uj = new double[0]; 986 double[] dtj = new double[0]; 987 double[] tmp = new double[0]; 988 double[] tmp2 = new double[0]; 989 double[,] tmpmatrix = new double[0,0]; 990 double v = 0; 991 int i_ = 0; 992 993 System.Diagnostics.Debug.Assert(n>0); 994 System.Diagnostics.Debug.Assert(m>=0); 995 System.Diagnostics.Debug.Assert(nc>=0 & nc<m+1); 996 result = true; 997 998 // 999 // Initialize design matrix and right part. 1000 // Add fictional rows if needed to ensure that N>=M+1. 1001 // 1002 designmatrix = new double[Math.Max(n, m+1)+1, m+1+1]; 1003 rightpart = new double[Math.Max(n, m+1)+1]; 1004 for(i=1; i<=n; i++) 1005 { 1006 for(j=1; j<=m+1; j++) 1007 { 1008 v = 2*(x[i-1]-a)/(b-a)-1; 1009 if( j==1 ) 1010 { 1011 designmatrix[i,j] = 1.0; 1012 } 1013 if( j==2 ) 1014 { 1015 designmatrix[i,j] = v; 1016 } 1017 if( j>2 ) 1018 { 1019 designmatrix[i,j] = 2.0*v*designmatrix[i,j-1]-designmatrix[i,j-2]; 1020 } 1021 } 1022 } 1023 for(i=1; i<=n; i++) 1024 { 1025 for(j=1; j<=m+1; j++) 1026 { 1027 designmatrix[i,j] = w[i-1]*designmatrix[i,j]; 1028 } 1029 } 1030 for(i=n+1; i<=m+1; i++) 1031 { 1032 for(j=1; j<=m+1; j++) 1033 { 1034 designmatrix[i,j] = 0; 1035 } 1036 } 1037 for(i=0; i<=n-1; i++) 1038 { 1039 rightpart[i+1] = w[i]*y[i]; 1040 } 1041 for(i=n+1; i<=m+1; i++) 1042 { 1043 rightpart[i] = 0; 1044 } 1045 n = Math.Max(n, m+1); 1046 1047 // 1048 // Now N>=M+1 and we are ready to the next stage. 1049 // Handle constraints. 1050 // Represent feasible set of coefficients as x = C*t + d 1051 // 1052 c = new double[m+1+1, m+1+1]; 1053 d = new double[m+1+1]; 1054 if( nc==0 ) 1055 { 1056 1057 // 1058 // No constraints 1059 // 1060 for(i=1; i<=m+1; i++) 1061 { 1062 for(j=1; j<=m+1; j++) 1063 { 1064 c[i,j] = 0; 1065 } 1066 d[i] = 0; 1067 } 1068 for(i=1; i<=m+1; i++) 1069 { 1070 c[i,i] = 1; 1071 } 1072 reducedsize = m+1; 1073 } 1074 else 1075 { 1076 1077 // 1078 // Constraints are present. 1079 // Fill constraints matrix CMatrix and solve CMatrix*x = cr. 1080 // 1081 cmatrix = new double[nc+1, m+1+1]; 1082 cr = new double[nc+1]; 1083 tj = new double[m+1]; 1084 uj = new double[m+1]; 1085 dtj = new double[m+1]; 1086 for(i=0; i<=nc-1; i++) 1087 { 1088 v = 2*(xc[i]-a)/(b-a)-1; 1089 for(j=0; j<=m; j++) 1090 { 1091 if( j==0 ) 1092 { 1093 tj[j] = 1; 1094 uj[j] = 1; 1095 dtj[j] = 0; 1096 } 1097 if( j==1 ) 1098 { 1099 tj[j] = v; 1100 uj[j] = 2*v; 1101 dtj[j] = 1; 1102 } 1103 if( j>1 ) 1104 { 1105 tj[j] = 2*v*tj[j-1]-tj[j-2]; 1106 uj[j] = 2*v*uj[j-1]-uj[j-2]; 1107 dtj[j] = j*uj[j-1]; 1108 } 1109 System.Diagnostics.Debug.Assert(dc[i]==0 | dc[i]==1); 1110 if( dc[i]==0 ) 1111 { 1112 cmatrix[i+1,j+1] = tj[j]; 1113 } 1114 if( dc[i]==1 ) 1115 { 1116 cmatrix[i+1,j+1] = dtj[j]; 1117 } 1118 } 1119 cr[i+1] = yc[i]; 1120 } 1121 1122 // 1123 // Solve CMatrix*x = cr. 1124 // Fill C and d: 1125 // 1. SVD: CMatrix = U * WS * V^T 1126 // 2. C := V[1:M+1,NC+1:M+1] 1127 // 3. tmp := WS^-1 * U^T * cr 1128 // 4. d := V[1:M+1,1:NC] * tmp 1129 // 1130 if( !svd.svddecomposition(cmatrix, nc, m+1, 2, 2, 2, ref ws, ref u, ref vt) ) 1131 { 1132 result = false; 1133 return result; 1134 } 1135 if( ws[1]==0 | ws[nc]<=AP.Math.MachineEpsilon*10*Math.Sqrt(nc)*ws[1] ) 1136 { 1137 result = false; 1138 return result; 1139 } 1140 c = new double[m+1+1, m+1-nc+1]; 1141 d = new double[m+1+1]; 1142 for(i=1; i<=m+1-nc; i++) 1143 { 1144 for(i_=1; i_<=m+1;i_++) 1145 { 1146 c[i_,i] = vt[nc+i,i_]; 1147 } 1148 } 1149 tmp = new double[nc+1]; 1150 for(i=1; i<=nc; i++) 1151 { 1152 v = 0.0; 1153 for(i_=1; i_<=nc;i_++) 1154 { 1155 v += u[i_,i]*cr[i_]; 1156 } 1157 tmp[i] = v/ws[i]; 1158 } 1159 for(i=1; i<=m+1; i++) 1160 { 1161 d[i] = 0; 1162 } 1163 for(i=1; i<=nc; i++) 1164 { 1165 v = tmp[i]; 1166 for(i_=1; i_<=m+1;i_++) 1167 { 1168 d[i_] = d[i_] + v*vt[i,i_]; 1169 } 1170 } 1171 1172 // 1173 // Reduce problem: 1174 // 1. RightPart := RightPart - DesignMatrix*d 1175 // 2. DesignMatrix := DesignMatrix*C 1176 // 1177 for(i=1; i<=n; i++) 1178 { 1179 v = 0.0; 1180 for(i_=1; i_<=m+1;i_++) 1181 { 1182 v += designmatrix[i,i_]*d[i_]; 1183 } 1184 rightpart[i] = rightpart[i]-v; 1185 } 1186 reducedsize = m+1-nc; 1187 tmpmatrix = new double[n+1, reducedsize+1]; 1188 tmp = new double[n+1]; 1189 blas.matrixmatrixmultiply(ref designmatrix, 1, n, 1, m+1, false, ref c, 1, m+1, 1, reducedsize, false, 1.0, ref tmpmatrix, 1, n, 1, reducedsize, 0.0, ref tmp); 1190 blas.copymatrix(ref tmpmatrix, 1, n, 1, reducedsize, ref designmatrix, 1, n, 1, reducedsize); 1191 } 1192 1193 // 1194 // Solve reduced problem DesignMatrix*t = RightPart. 1195 // 1196 if( !svd.svddecomposition(designmatrix, n, reducedsize, 1, 1, 2, ref ws, ref u, ref vt) ) 1197 { 1198 result = false; 1199 return result; 1200 } 1201 tmp = new double[reducedsize+1]; 1202 tmp2 = new double[reducedsize+1]; 1203 for(i=1; i<=reducedsize; i++) 1204 { 1205 tmp[i] = 0; 1206 } 1207 for(i=1; i<=n; i++) 1208 { 1209 v = rightpart[i]; 1210 for(i_=1; i_<=reducedsize;i_++) 1211 { 1212 tmp[i_] = tmp[i_] + v*u[i,i_]; 1213 } 1214 } 1215 for(i=1; i<=reducedsize; i++) 1216 { 1217 if( ws[i]!=0 & ws[i]>AP.Math.MachineEpsilon*10*Math.Sqrt(nc)*ws[1] ) 1218 { 1219 tmp[i] = tmp[i]/ws[i]; 1220 } 1221 else 1222 { 1223 tmp[i] = 0; 1224 } 1225 } 1226 for(i=1; i<=reducedsize; i++) 1227 { 1228 tmp2[i] = 0; 1229 } 1230 for(i=1; i<=reducedsize; i++) 1231 { 1232 v = tmp[i]; 1233 for(i_=1; i_<=reducedsize;i_++) 1234 { 1235 tmp2[i_] = tmp2[i_] + v*vt[i,i_]; 1236 } 1237 } 1238 1239 // 1240 // Solution is in the tmp2. 1241 // Transform it from t to x. 1242 // 1243 ctbl = new double[m+2+1]; 1244 for(i=1; i<=m+1; i++) 1245 { 1246 v = 0.0; 1247 for(i_=1; i_<=reducedsize;i_++) 1248 { 1249 v += c[i,i_]*tmp2[i_]; 1250 } 1251 ctbl[i-1] = v+d[i]; 1252 } 1253 ctbl[m+1] = a; 1254 ctbl[m+2] = b; 1255 return result; 132 1256 } 133 mi = Math.Max(mi, ni); 134 135 // 136 // LQ-decomposition of A' 137 // B2 := Q*B 138 // 139 lq.lqdecomposition(ref a, ni, mi, ref tau); 140 lq.unpackqfromlq(ref a, ni, mi, ref tau, ni, ref q); 141 b2 = new double[1+1, ni+1]; 142 for(j=1; j<=ni; j++) 1257 1258 1259 /************************************************************************* 1260 Calculation of a Chebyshev polynomial obtained during least squares 1261 approximaion at the given point. 1262 1263 Input parameters: 1264 M - order of polynomial (parameter of the 1265 BuildChebyshevLeastSquares function). 1266 A - coefficient table. 1267 A[0..M] contains coefficients of the i-th Chebyshev polynomial. 1268 A[M+1] contains left boundary of approximation interval. 1269 A[M+2] contains right boundary of approximation interval. 1270 X - point to perform calculations in. 1271 1272 The result is the value at the given point. 1273 1274 It should be noted that array A contains coefficients of the Chebyshev 1275 polynomials defined on interval [-1,1]. Argument is reduced to this 1276 interval before calculating polynomial value. 1277 -- ALGLIB -- 1278 Copyright by Bochkanov Sergey 1279 *************************************************************************/ 1280 public static double calculatechebyshevleastsquares(int m, 1281 ref double[] a, 1282 double x) 143 1283 { 144 b2[1,j] = 0; 145 } 146 for(i=1; i<=ni; i++) 147 { 148 v = 0.0; 149 for(i_=1; i_<=mi;i_++) 150 { 151 v += b[i_]*q[i,i_]; 152 } 153 b2[1,i] = v; 154 } 155 156 // 157 // Back from A' to A 158 // Making cols(A)=rows(A) 159 // 160 for(i=1; i<=ni-1; i++) 161 { 162 for(i_=i+1; i_<=ni;i_++) 163 { 164 a[i,i_] = a[i_,i]; 165 } 166 } 167 for(i=2; i<=ni; i++) 168 { 169 for(j=1; j<=i-1; j++) 170 { 171 a[i,j] = 0; 172 } 173 } 174 175 // 176 // Bidiagonal decomposition of A 177 // A = Q * d2 * P' 178 // B2 := (Q'*B2')' 179 // 180 bidiagonal.tobidiagonal(ref a, ni, ni, ref tauq, ref taup); 181 bidiagonal.multiplybyqfrombidiagonal(ref a, ni, ni, ref tauq, ref b2, 1, ni, true, false); 182 bidiagonal.unpackptfrombidiagonal(ref a, ni, ni, ref taup, ni, ref vt); 183 bidiagonal.unpackdiagonalsfrombidiagonal(ref a, ni, ni, ref isuppera, ref d, ref e); 184 185 // 186 // Singular value decomposition of A 187 // A = U * d * V' 188 // B2 := (U'*B2')' 189 // 190 if( !bdsvd.bidiagonalsvddecomposition(ref d, e, ni, isuppera, false, ref b2, 1, ref q, 0, ref vt, ni) ) 191 { 192 for(i=0; i<=ni-1; i++) 193 { 194 c[i] = 0; 195 } 196 return; 197 } 198 199 // 200 // B2 := (d^(-1) * B2')' 201 // 202 if( d[1]!=0 ) 203 { 204 for(i=1; i<=ni; i++) 205 { 206 if( d[i]>AP.Math.MachineEpsilon*10*Math.Sqrt(ni)*d[1] ) 207 { 208 b2[1,i] = b2[1,i]/d[i]; 209 } 210 else 211 { 212 b2[1,i] = 0; 213 } 214 } 215 } 216 217 // 218 // B := (V * B2')' 219 // 220 for(i=1; i<=ni; i++) 221 { 222 b[i] = 0; 223 } 224 for(i=1; i<=ni; i++) 225 { 226 v = b2[1,i]; 227 for(i_=1; i_<=ni;i_++) 228 { 229 b[i_] = b[i_] + v*vt[i,i_]; 230 } 231 } 232 233 // 234 // Out 235 // 236 for(i=1; i<=ni; i++) 237 { 238 c[i-1] = b[i]; 1284 double result = 0; 1285 double b1 = 0; 1286 double b2 = 0; 1287 int i = 0; 1288 1289 x = 2*(x-a[m+1])/(a[m+2]-a[m+1])-1; 1290 b1 = 0; 1291 b2 = 0; 1292 i = m; 1293 do 1294 { 1295 result = 2*x*b1-b2+a[i]; 1296 b2 = b1; 1297 b1 = result; 1298 i = i-1; 1299 } 1300 while( i>=0 ); 1301 result = result-x*b2; 1302 return result; 239 1303 } 240 1304 } 241 242 243 /*************************************************************************244 Linear approximation using least squares method245 246 The subroutine calculates coefficients of the line approximating given247 function.248 249 Input parameters:250 X - array[0..N-1], it contains a set of abscissas.251 Y - array[0..N-1], function values.252 N - number of points, N>=1253 254 Output parameters:255 a, b- coefficients of linear approximation a+b*t256 257 -- ALGLIB --258 Copyright by Bochkanov Sergey259 *************************************************************************/260 public static void buildlinearleastsquares(ref double[] x,261 ref double[] y,262 int n,263 ref double a,264 ref double b)265 {266 double pp = 0;267 double qq = 0;268 double pq = 0;269 double b1 = 0;270 double b2 = 0;271 double d1 = 0;272 double d2 = 0;273 double t1 = 0;274 double t2 = 0;275 double phi = 0;276 double c = 0;277 double s = 0;278 double m = 0;279 int i = 0;280 281 pp = n;282 qq = 0;283 pq = 0;284 b1 = 0;285 b2 = 0;286 for(i=0; i<=n-1; i++)287 {288 pq = pq+x[i];289 qq = qq+AP.Math.Sqr(x[i]);290 b1 = b1+y[i];291 b2 = b2+x[i]*y[i];292 }293 phi = Math.Atan2(2*pq, qq-pp)/2;294 c = Math.Cos(phi);295 s = Math.Sin(phi);296 d1 = AP.Math.Sqr(c)*pp+AP.Math.Sqr(s)*qq-2*s*c*pq;297 d2 = AP.Math.Sqr(s)*pp+AP.Math.Sqr(c)*qq+2*s*c*pq;298 if( Math.Abs(d1)>Math.Abs(d2) )299 {300 m = Math.Abs(d1);301 }302 else303 {304 m = Math.Abs(d2);305 }306 t1 = c*b1-s*b2;307 t2 = s*b1+c*b2;308 if( Math.Abs(d1)>m*AP.Math.MachineEpsilon*1000 )309 {310 t1 = t1/d1;311 }312 else313 {314 t1 = 0;315 }316 if( Math.Abs(d2)>m*AP.Math.MachineEpsilon*1000 )317 {318 t2 = t2/d2;319 }320 else321 {322 t2 = 0;323 }324 a = c*t1+s*t2;325 b = -(s*t1)+c*t2;326 }327 328 329 /*************************************************************************330 Weighted cubic spline approximation using linear least squares331 332 Input parameters:333 X - array[0..N-1], abscissas334 Y - array[0..N-1], function values335 W - array[0..N-1], weights.336 A, B- interval to build splines in.337 N - number of points used. N>=1.338 M - number of basic splines, M>=2.339 340 Output parameters:341 CTbl- coefficients table to be used by SplineInterpolation function.342 -- ALGLIB --343 Copyright by Bochkanov Sergey344 *************************************************************************/345 public static void buildsplineleastsquares(ref double[] x,346 ref double[] y,347 ref double[] w,348 double a,349 double b,350 int n,351 int m,352 ref double[] ctbl)353 {354 int i = 0;355 int j = 0;356 double[,] ma = new double[0,0];357 double[,] q = new double[0,0];358 double[,] vt = new double[0,0];359 double[] mb = new double[0];360 double[] tau = new double[0];361 double[,] b2 = new double[0,0];362 double[] tauq = new double[0];363 double[] taup = new double[0];364 double[] d = new double[0];365 double[] e = new double[0];366 bool isuppera = new bool();367 int mi = 0;368 int ni = 0;369 double v = 0;370 double[] sx = new double[0];371 double[] sy = new double[0];372 int i_ = 0;373 374 System.Diagnostics.Debug.Assert(m>=2, "BuildSplineLeastSquares: M is too small!");375 mi = n;376 ni = m;377 sx = new double[ni-1+1];378 sy = new double[ni-1+1];379 380 //381 // Initializing design matrix382 // Here we are making MI>=NI383 //384 ma = new double[ni+1, Math.Max(mi, ni)+1];385 mb = new double[Math.Max(mi, ni)+1];386 for(j=0; j<=ni-1; j++)387 {388 sx[j] = a+(b-a)*j/(ni-1);389 }390 for(j=0; j<=ni-1; j++)391 {392 for(i=0; i<=ni-1; i++)393 {394 sy[i] = 0;395 }396 sy[j] = 1;397 spline3.buildcubicspline(sx, sy, ni, 0, 0.0, 0, 0.0, ref ctbl);398 for(i=0; i<=mi-1; i++)399 {400 ma[j+1,i+1] = w[i]*spline3.splineinterpolation(ref ctbl, x[i]);401 }402 }403 for(j=1; j<=ni; j++)404 {405 for(i=mi+1; i<=ni; i++)406 {407 ma[j,i] = 0;408 }409 }410 411 //412 // Initializing right part413 //414 for(i=0; i<=mi-1; i++)415 {416 mb[i+1] = w[i]*y[i];417 }418 for(i=mi+1; i<=ni; i++)419 {420 mb[i] = 0;421 }422 mi = Math.Max(mi, ni);423 424 //425 // LQ-decomposition of A'426 // B2 := Q*B427 //428 lq.lqdecomposition(ref ma, ni, mi, ref tau);429 lq.unpackqfromlq(ref ma, ni, mi, ref tau, ni, ref q);430 b2 = new double[1+1, ni+1];431 for(j=1; j<=ni; j++)432 {433 b2[1,j] = 0;434 }435 for(i=1; i<=ni; i++)436 {437 v = 0.0;438 for(i_=1; i_<=mi;i_++)439 {440 v += mb[i_]*q[i,i_];441 }442 b2[1,i] = v;443 }444 445 //446 // Back from A' to A447 // Making cols(A)=rows(A)448 //449 for(i=1; i<=ni-1; i++)450 {451 for(i_=i+1; i_<=ni;i_++)452 {453 ma[i,i_] = ma[i_,i];454 }455 }456 for(i=2; i<=ni; i++)457 {458 for(j=1; j<=i-1; j++)459 {460 ma[i,j] = 0;461 }462 }463 464 //465 // Bidiagonal decomposition of A466 // A = Q * d2 * P'467 // B2 := (Q'*B2')'468 //469 bidiagonal.tobidiagonal(ref ma, ni, ni, ref tauq, ref taup);470 bidiagonal.multiplybyqfrombidiagonal(ref ma, ni, ni, ref tauq, ref b2, 1, ni, true, false);471 bidiagonal.unpackptfrombidiagonal(ref ma, ni, ni, ref taup, ni, ref vt);472 bidiagonal.unpackdiagonalsfrombidiagonal(ref ma, ni, ni, ref isuppera, ref d, ref e);473 474 //475 // Singular value decomposition of A476 // A = U * d * V'477 // B2 := (U'*B2')'478 //479 if( !bdsvd.bidiagonalsvddecomposition(ref d, e, ni, isuppera, false, ref b2, 1, ref q, 0, ref vt, ni) )480 {481 for(i=1; i<=ni; i++)482 {483 d[i] = 0;484 b2[1,i] = 0;485 for(j=1; j<=ni; j++)486 {487 if( i==j )488 {489 vt[i,j] = 1;490 }491 else492 {493 vt[i,j] = 0;494 }495 }496 }497 b2[1,1] = 1;498 }499 500 //501 // B2 := (d^(-1) * B2')'502 //503 for(i=1; i<=ni; i++)504 {505 if( d[i]>AP.Math.MachineEpsilon*10*Math.Sqrt(ni)*d[1] )506 {507 b2[1,i] = b2[1,i]/d[i];508 }509 else510 {511 b2[1,i] = 0;512 }513 }514 515 //516 // B := (V * B2')'517 //518 for(i=1; i<=ni; i++)519 {520 mb[i] = 0;521 }522 for(i=1; i<=ni; i++)523 {524 v = b2[1,i];525 for(i_=1; i_<=ni;i_++)526 {527 mb[i_] = mb[i_] + v*vt[i,i_];528 }529 }530 531 //532 // Forming result spline533 //534 for(i=0; i<=ni-1; i++)535 {536 sy[i] = mb[i+1];537 }538 spline3.buildcubicspline(sx, sy, ni, 0, 0.0, 0, 0.0, ref ctbl);539 }540 541 542 /*************************************************************************543 Polynomial approximation using least squares method544 545 The subroutine calculates coefficients of the polynomial approximating546 given function. It is recommended to use this function only if you need to547 obtain coefficients of approximation polynomial. If you have to build and548 calculate polynomial approximation (NOT coefficients), it's better to use549 BuildChebyshevLeastSquares subroutine in combination with550 CalculateChebyshevLeastSquares subroutine. The result of Chebyshev551 polynomial approximation is equivalent to the result obtained using powers552 of X, but has higher accuracy due to better numerical properties of553 Chebyshev polynomials.554 555 Input parameters:556 X - array[0..N-1], abscissas557 Y - array[0..N-1], function values558 N - number of points, N>=1559 M - order of polynomial required, M>=0560 561 Output parameters:562 C - approximating polynomial coefficients, array[0..M],563 C[i] - coefficient at X^i.564 565 -- ALGLIB --566 Copyright by Bochkanov Sergey567 *************************************************************************/568 public static void buildpolynomialleastsquares(ref double[] x,569 ref double[] y,570 int n,571 int m,572 ref double[] c)573 {574 double[] ctbl = new double[0];575 double[] w = new double[0];576 double[] c1 = new double[0];577 double maxx = 0;578 double minx = 0;579 int i = 0;580 int j = 0;581 int k = 0;582 double e = 0;583 double d = 0;584 double l1 = 0;585 double l2 = 0;586 double[] z2 = new double[0];587 double[] z1 = new double[0];588 589 590 //591 // Initialize592 //593 maxx = x[0];594 minx = x[0];595 for(i=1; i<=n-1; i++)596 {597 if( x[i]>maxx )598 {599 maxx = x[i];600 }601 if( x[i]<minx )602 {603 minx = x[i];604 }605 }606 if( minx==maxx )607 {608 minx = minx-0.5;609 maxx = maxx+0.5;610 }611 w = new double[n-1+1];612 for(i=0; i<=n-1; i++)613 {614 w[i] = 1;615 }616 617 //618 // Build Chebyshev approximation619 //620 buildchebyshevleastsquares(ref x, ref y, ref w, minx, maxx, n, m, ref ctbl);621 622 //623 // From Chebyshev to powers of X624 //625 c1 = new double[m+1];626 for(i=0; i<=m; i++)627 {628 c1[i] = 0;629 }630 d = 0;631 for(i=0; i<=m; i++)632 {633 for(k=i; k<=m; k++)634 {635 e = c1[k];636 c1[k] = 0;637 if( i<=1 & k==i )638 {639 c1[k] = 1;640 }641 else642 {643 if( i!=0 )644 {645 c1[k] = 2*d;646 }647 if( k>i+1 )648 {649 c1[k] = c1[k]-c1[k-2];650 }651 }652 d = e;653 }654 d = c1[i];655 e = 0;656 k = i;657 while( k<=m )658 {659 e = e+c1[k]*ctbl[k];660 k = k+2;661 }662 c1[i] = e;663 }664 665 //666 // Linear translation667 //668 l1 = 2/(ctbl[m+2]-ctbl[m+1]);669 l2 = -(2*ctbl[m+1]/(ctbl[m+2]-ctbl[m+1]))-1;670 c = new double[m+1];671 z2 = new double[m+1];672 z1 = new double[m+1];673 c[0] = c1[0];674 z1[0] = 1;675 z2[0] = 1;676 for(i=1; i<=m; i++)677 {678 z2[i] = 1;679 z1[i] = l2*z1[i-1];680 c[0] = c[0]+c1[i]*z1[i];681 }682 for(j=1; j<=m; j++)683 {684 z2[0] = l1*z2[0];685 c[j] = c1[j]*z2[0];686 for(i=j+1; i<=m; i++)687 {688 k = i-j;689 z2[k] = l1*z2[k]+z2[k-1];690 c[j] = c[j]+c1[i]*z2[k]*z1[k];691 }692 }693 }694 695 696 /*************************************************************************697 Chebyshev polynomial approximation using least squares method.698 699 The algorithm reduces interval [A, B] to the interval [-1,1], then builds700 least squares approximation using Chebyshev polynomials.701 702 Input parameters:703 X - array[0..N-1], abscissas704 Y - array[0..N-1], function values705 W - array[0..N-1], weights706 A, B- interval to build approximating polynomials in.707 N - number of points used. N>=1.708 M - order of polynomial, M>=0. This parameter is passed into709 CalculateChebyshevLeastSquares function.710 711 Output parameters:712 CTbl - coefficient table. This parameter is passed into713 CalculateChebyshevLeastSquares function.714 -- ALGLIB --715 Copyright by Bochkanov Sergey716 *************************************************************************/717 public static void buildchebyshevleastsquares(ref double[] x,718 ref double[] y,719 ref double[] w,720 double a,721 double b,722 int n,723 int m,724 ref double[] ctbl)725 {726 int i = 0;727 int j = 0;728 double[,] ma = new double[0,0];729 double[,] q = new double[0,0];730 double[,] vt = new double[0,0];731 double[] mb = new double[0];732 double[] tau = new double[0];733 double[,] b2 = new double[0,0];734 double[] tauq = new double[0];735 double[] taup = new double[0];736 double[] d = new double[0];737 double[] e = new double[0];738 bool isuppera = new bool();739 int mi = 0;740 int ni = 0;741 double v = 0;742 int i_ = 0;743 744 mi = n;745 ni = m+1;746 747 //748 // Initializing design matrix749 // Here we are making MI>=NI750 //751 ma = new double[ni+1, Math.Max(mi, ni)+1];752 mb = new double[Math.Max(mi, ni)+1];753 for(j=1; j<=ni; j++)754 {755 for(i=1; i<=mi; i++)756 {757 v = 2*(x[i-1]-a)/(b-a)-1;758 if( j==1 )759 {760 ma[j,i] = 1.0;761 }762 if( j==2 )763 {764 ma[j,i] = v;765 }766 if( j>2 )767 {768 ma[j,i] = 2.0*v*ma[j-1,i]-ma[j-2,i];769 }770 }771 }772 for(j=1; j<=ni; j++)773 {774 for(i=1; i<=mi; i++)775 {776 ma[j,i] = w[i-1]*ma[j,i];777 }778 }779 for(j=1; j<=ni; j++)780 {781 for(i=mi+1; i<=ni; i++)782 {783 ma[j,i] = 0;784 }785 }786 787 //788 // Initializing right part789 //790 for(i=0; i<=mi-1; i++)791 {792 mb[i+1] = w[i]*y[i];793 }794 for(i=mi+1; i<=ni; i++)795 {796 mb[i] = 0;797 }798 mi = Math.Max(mi, ni);799 800 //801 // LQ-decomposition of A'802 // B2 := Q*B803 //804 lq.lqdecomposition(ref ma, ni, mi, ref tau);805 lq.unpackqfromlq(ref ma, ni, mi, ref tau, ni, ref q);806 b2 = new double[1+1, ni+1];807 for(j=1; j<=ni; j++)808 {809 b2[1,j] = 0;810 }811 for(i=1; i<=ni; i++)812 {813 v = 0.0;814 for(i_=1; i_<=mi;i_++)815 {816 v += mb[i_]*q[i,i_];817 }818 b2[1,i] = v;819 }820 821 //822 // Back from A' to A823 // Making cols(A)=rows(A)824 //825 for(i=1; i<=ni-1; i++)826 {827 for(i_=i+1; i_<=ni;i_++)828 {829 ma[i,i_] = ma[i_,i];830 }831 }832 for(i=2; i<=ni; i++)833 {834 for(j=1; j<=i-1; j++)835 {836 ma[i,j] = 0;837 }838 }839 840 //841 // Bidiagonal decomposition of A842 // A = Q * d2 * P'843 // B2 := (Q'*B2')'844 //845 bidiagonal.tobidiagonal(ref ma, ni, ni, ref tauq, ref taup);846 bidiagonal.multiplybyqfrombidiagonal(ref ma, ni, ni, ref tauq, ref b2, 1, ni, true, false);847 bidiagonal.unpackptfrombidiagonal(ref ma, ni, ni, ref taup, ni, ref vt);848 bidiagonal.unpackdiagonalsfrombidiagonal(ref ma, ni, ni, ref isuppera, ref d, ref e);849 850 //851 // Singular value decomposition of A852 // A = U * d * V'853 // B2 := (U'*B2')'854 //855 if( !bdsvd.bidiagonalsvddecomposition(ref d, e, ni, isuppera, false, ref b2, 1, ref q, 0, ref vt, ni) )856 {857 for(i=1; i<=ni; i++)858 {859 d[i] = 0;860 b2[1,i] = 0;861 for(j=1; j<=ni; j++)862 {863 if( i==j )864 {865 vt[i,j] = 1;866 }867 else868 {869 vt[i,j] = 0;870 }871 }872 }873 b2[1,1] = 1;874 }875 876 //877 // B2 := (d^(-1) * B2')'878 //879 for(i=1; i<=ni; i++)880 {881 if( d[i]>AP.Math.MachineEpsilon*10*Math.Sqrt(ni)*d[1] )882 {883 b2[1,i] = b2[1,i]/d[i];884 }885 else886 {887 b2[1,i] = 0;888 }889 }890 891 //892 // B := (V * B2')'893 //894 for(i=1; i<=ni; i++)895 {896 mb[i] = 0;897 }898 for(i=1; i<=ni; i++)899 {900 v = b2[1,i];901 for(i_=1; i_<=ni;i_++)902 {903 mb[i_] = mb[i_] + v*vt[i,i_];904 }905 }906 907 //908 // Forming result909 //910 ctbl = new double[ni+1+1];911 for(i=1; i<=ni; i++)912 {913 ctbl[i-1] = mb[i];914 }915 ctbl[ni] = a;916 ctbl[ni+1] = b;917 }918 919 920 /*************************************************************************921 Weighted Chebyshev polynomial constrained least squares approximation.922 923 The algorithm reduces [A,B] to [-1,1] and builds the Chebyshev polynomials924 series by approximating a given function using the least squares method.925 926 Input parameters:927 X - abscissas, array[0..N-1]928 Y - function values, array[0..N-1]929 W - weights, array[0..N-1]. Each item in the squared sum of930 deviations from given values is multiplied by a square of931 corresponding weight.932 A, B- interval in which the approximating polynomials are built.933 N - number of points, N>0.934 XC, YC, DC-935 constraints (see description below)., array[0..NC-1]936 NC - number of constraints. 0 <= NC < M+1.937 M - degree of polynomial, M>=0. This parameter is passed into the938 CalculateChebyshevLeastSquares subroutine.939 940 Output parameters:941 CTbl- coefficient table. This parameter is passed into the942 CalculateChebyshevLeastSquares subroutine.943 944 Result:945 True, if the algorithm succeeded.946 False, if the internal singular value decomposition subroutine hasn't947 converged or the given constraints could not be met simultaneously (e.g.948 P(0)=0 è P(0)=1).949 950 Specifying constraints:951 This subroutine can solve the problem having constrained function952 values or its derivatives in several points. NC specifies the number of953 constraints, DC - the type of constraints, XC and YC - constraints as such.954 Thus, for each i from 0 to NC-1 the following constraint is given:955 P(xc[i]) = yc[i], if DC[i]=0956 or957 d/dx(P(xc[i])) = yc[i], if DC[i]=1958 (here P(x) is approximating polynomial).959 This version of the subroutine supports only either polynomial or its960 derivative value constraints. If DC[i] is not equal to 0 and 1, the961 subroutine will be aborted. The number of constraints should be less than962 the number of degrees of freedom of approximating polynomial - M+1 (at963 that, it could be equal to 0).964 965 -- ALGLIB --966 Copyright by Bochkanov Sergey967 *************************************************************************/968 public static bool buildchebyshevleastsquaresconstrained(ref double[] x,969 ref double[] y,970 ref double[] w,971 double a,972 double b,973 int n,974 ref double[] xc,975 ref double[] yc,976 ref int[] dc,977 int nc,978 int m,979 ref double[] ctbl)980 {981 bool result = new bool();982 int i = 0;983 int j = 0;984 int reducedsize = 0;985 double[,] designmatrix = new double[0,0];986 double[] rightpart = new double[0];987 double[,] cmatrix = new double[0,0];988 double[,] c = new double[0,0];989 double[,] u = new double[0,0];990 double[,] vt = new double[0,0];991 double[] d = new double[0];992 double[] cr = new double[0];993 double[] ws = new double[0];994 double[] tj = new double[0];995 double[] uj = new double[0];996 double[] dtj = new double[0];997 double[] tmp = new double[0];998 double[] tmp2 = new double[0];999 double[,] tmpmatrix = new double[0,0];1000 double v = 0;1001 int i_ = 0;1002 1003 System.Diagnostics.Debug.Assert(n>0);1004 System.Diagnostics.Debug.Assert(m>=0);1005 System.Diagnostics.Debug.Assert(nc>=0 & nc<m+1);1006 result = true;1007 1008 //1009 // Initialize design matrix and right part.1010 // Add fictional rows if needed to ensure that N>=M+1.1011 //1012 designmatrix = new double[Math.Max(n, m+1)+1, m+1+1];1013 rightpart = new double[Math.Max(n, m+1)+1];1014 for(i=1; i<=n; i++)1015 {1016 for(j=1; j<=m+1; j++)1017 {1018 v = 2*(x[i-1]-a)/(b-a)-1;1019 if( j==1 )1020 {1021 designmatrix[i,j] = 1.0;1022 }1023 if( j==2 )1024 {1025 designmatrix[i,j] = v;1026 }1027 if( j>2 )1028 {1029 designmatrix[i,j] = 2.0*v*designmatrix[i,j-1]-designmatrix[i,j-2];1030 }1031 }1032 }1033 for(i=1; i<=n; i++)1034 {1035 for(j=1; j<=m+1; j++)1036 {1037 designmatrix[i,j] = w[i-1]*designmatrix[i,j];1038 }1039 }1040 for(i=n+1; i<=m+1; i++)1041 {1042 for(j=1; j<=m+1; j++)1043 {1044 designmatrix[i,j] = 0;1045 }1046 }1047 for(i=0; i<=n-1; i++)1048 {1049 rightpart[i+1] = w[i]*y[i];1050 }1051 for(i=n+1; i<=m+1; i++)1052 {1053 rightpart[i] = 0;1054 }1055 n = Math.Max(n, m+1);1056 1057 //1058 // Now N>=M+1 and we are ready to the next stage.1059 // Handle constraints.1060 // Represent feasible set of coefficients as x = C*t + d1061 //1062 c = new double[m+1+1, m+1+1];1063 d = new double[m+1+1];1064 if( nc==0 )1065 {1066 1067 //1068 // No constraints1069 //1070 for(i=1; i<=m+1; i++)1071 {1072 for(j=1; j<=m+1; j++)1073 {1074 c[i,j] = 0;1075 }1076 d[i] = 0;1077 }1078 for(i=1; i<=m+1; i++)1079 {1080 c[i,i] = 1;1081 }1082 reducedsize = m+1;1083 }1084 else1085 {1086 1087 //1088 // Constraints are present.1089 // Fill constraints matrix CMatrix and solve CMatrix*x = cr.1090 //1091 cmatrix = new double[nc+1, m+1+1];1092 cr = new double[nc+1];1093 tj = new double[m+1];1094 uj = new double[m+1];1095 dtj = new double[m+1];1096 for(i=0; i<=nc-1; i++)1097 {1098 v = 2*(xc[i]-a)/(b-a)-1;1099 for(j=0; j<=m; j++)1100 {1101 if( j==0 )1102 {1103 tj[j] = 1;1104 uj[j] = 1;1105 dtj[j] = 0;1106 }1107 if( j==1 )1108 {1109 tj[j] = v;1110 uj[j] = 2*v;1111 dtj[j] = 1;1112 }1113 if( j>1 )1114 {1115 tj[j] = 2*v*tj[j-1]-tj[j-2];1116 uj[j] = 2*v*uj[j-1]-uj[j-2];1117 dtj[j] = j*uj[j-1];1118 }1119 System.Diagnostics.Debug.Assert(dc[i]==0 | dc[i]==1);1120 if( dc[i]==0 )1121 {1122 cmatrix[i+1,j+1] = tj[j];1123 }1124 if( dc[i]==1 )1125 {1126 cmatrix[i+1,j+1] = dtj[j];1127 }1128 }1129 cr[i+1] = yc[i];1130 }1131 1132 //1133 // Solve CMatrix*x = cr.1134 // Fill C and d:1135 // 1. SVD: CMatrix = U * WS * V^T1136 // 2. C := V[1:M+1,NC+1:M+1]1137 // 3. tmp := WS^-1 * U^T * cr1138 // 4. d := V[1:M+1,1:NC] * tmp1139 //1140 if( !svd.svddecomposition(cmatrix, nc, m+1, 2, 2, 2, ref ws, ref u, ref vt) )1141 {1142 result = false;1143 return result;1144 }1145 if( ws[1]==0 | ws[nc]<=AP.Math.MachineEpsilon*10*Math.Sqrt(nc)*ws[1] )1146 {1147 result = false;1148 return result;1149 }1150 c = new double[m+1+1, m+1-nc+1];1151 d = new double[m+1+1];1152 for(i=1; i<=m+1-nc; i++)1153 {1154 for(i_=1; i_<=m+1;i_++)1155 {1156 c[i_,i] = vt[nc+i,i_];1157 }1158 }1159 tmp = new double[nc+1];1160 for(i=1; i<=nc; i++)1161 {1162 v = 0.0;1163 for(i_=1; i_<=nc;i_++)1164 {1165 v += u[i_,i]*cr[i_];1166 }1167 tmp[i] = v/ws[i];1168 }1169 for(i=1; i<=m+1; i++)1170 {1171 d[i] = 0;1172 }1173 for(i=1; i<=nc; i++)1174 {1175 v = tmp[i];1176 for(i_=1; i_<=m+1;i_++)1177 {1178 d[i_] = d[i_] + v*vt[i,i_];1179 }1180 }1181 1182 //1183 // Reduce problem:1184 // 1. RightPart := RightPart - DesignMatrix*d1185 // 2. DesignMatrix := DesignMatrix*C1186 //1187 for(i=1; i<=n; i++)1188 {1189 v = 0.0;1190 for(i_=1; i_<=m+1;i_++)1191 {1192 v += designmatrix[i,i_]*d[i_];1193 }1194 rightpart[i] = rightpart[i]-v;1195 }1196 reducedsize = m+1-nc;1197 tmpmatrix = new double[n+1, reducedsize+1];1198 tmp = new double[n+1];1199 blas.matrixmatrixmultiply(ref designmatrix, 1, n, 1, m+1, false, ref c, 1, m+1, 1, reducedsize, false, 1.0, ref tmpmatrix, 1, n, 1, reducedsize, 0.0, ref tmp);1200 blas.copymatrix(ref tmpmatrix, 1, n, 1, reducedsize, ref designmatrix, 1, n, 1, reducedsize);1201 }1202 1203 //1204 // Solve reduced problem DesignMatrix*t = RightPart.1205 //1206 if( !svd.svddecomposition(designmatrix, n, reducedsize, 1, 1, 2, ref ws, ref u, ref vt) )1207 {1208 result = false;1209 return result;1210 }1211 tmp = new double[reducedsize+1];1212 tmp2 = new double[reducedsize+1];1213 for(i=1; i<=reducedsize; i++)1214 {1215 tmp[i] = 0;1216 }1217 for(i=1; i<=n; i++)1218 {1219 v = rightpart[i];1220 for(i_=1; i_<=reducedsize;i_++)1221 {1222 tmp[i_] = tmp[i_] + v*u[i,i_];1223 }1224 }1225 for(i=1; i<=reducedsize; i++)1226 {1227 if( ws[i]!=0 & ws[i]>AP.Math.MachineEpsilon*10*Math.Sqrt(nc)*ws[1] )1228 {1229 tmp[i] = tmp[i]/ws[i];1230 }1231 else1232 {1233 tmp[i] = 0;1234 }1235 }1236 for(i=1; i<=reducedsize; i++)1237 {1238 tmp2[i] = 0;1239 }1240 for(i=1; i<=reducedsize; i++)1241 {1242 v = tmp[i];1243 for(i_=1; i_<=reducedsize;i_++)1244 {1245 tmp2[i_] = tmp2[i_] + v*vt[i,i_];1246 }1247 }1248 1249 //1250 // Solution is in the tmp2.1251 // Transform it from t to x.1252 //1253 ctbl = new double[m+2+1];1254 for(i=1; i<=m+1; i++)1255 {1256 v = 0.0;1257 for(i_=1; i_<=reducedsize;i_++)1258 {1259 v += c[i,i_]*tmp2[i_];1260 }1261 ctbl[i-1] = v+d[i];1262 }1263 ctbl[m+1] = a;1264 ctbl[m+2] = b;1265 return result;1266 }1267 1268 1269 /*************************************************************************1270 Calculation of a Chebyshev polynomial obtained during least squares1271 approximaion at the given point.1272 1273 Input parameters:1274 M - order of polynomial (parameter of the1275 BuildChebyshevLeastSquares function).1276 A - coefficient table.1277 A[0..M] contains coefficients of the i-th Chebyshev polynomial.1278 A[M+1] contains left boundary of approximation interval.1279 A[M+2] contains right boundary of approximation interval.1280 X - point to perform calculations in.1281 1282 The result is the value at the given point.1283 1284 It should be noted that array A contains coefficients of the Chebyshev1285 polynomials defined on interval [-1,1]. Argument is reduced to this1286 interval before calculating polynomial value.1287 -- ALGLIB --1288 Copyright by Bochkanov Sergey1289 *************************************************************************/1290 public static double calculatechebyshevleastsquares(int m,1291 ref double[] a,1292 double x)1293 {1294 double result = 0;1295 double b1 = 0;1296 double b2 = 0;1297 int i = 0;1298 1299 x = 2*(x-a[m+1])/(a[m+2]-a[m+1])-1;1300 b1 = 0;1301 b2 = 0;1302 i = m;1303 do1304 {1305 result = 2*x*b1-b2+a[i];1306 b2 = b1;1307 b1 = result;1308 i = i-1;1309 }1310 while( i>=0 );1311 result = result-x*b2;1312 return result;1313 }1314 1305 }
Note: See TracChangeset
for help on using the changeset viewer.