[2154] | 1 | /*************************************************************************
|
---|
| 2 | Copyright (c) 2006-2007, Sergey Bochkanov (ALGLIB project).
|
---|
| 3 |
|
---|
[2430] | 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.
|
---|
[2154] | 9 |
|
---|
[2430] | 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.
|
---|
[2154] | 14 |
|
---|
[2430] | 15 | A copy of the GNU General Public License is available at
|
---|
| 16 | http://www.fsf.org/licensing/licenses
|
---|
[2154] | 17 |
|
---|
[2430] | 18 | >>> END OF LICENSE >>>
|
---|
[2154] | 19 | *************************************************************************/
|
---|
| 20 |
|
---|
| 21 | using System;
|
---|
| 22 |
|
---|
[2430] | 23 | namespace alglib
|
---|
[2154] | 24 | {
|
---|
[2430] | 25 | public class leastsquares
|
---|
| 26 | {
|
---|
| 27 | /*************************************************************************
|
---|
| 28 | Weighted approximation by arbitrary function basis in a space of arbitrary
|
---|
| 29 | dimension using linear least squares method.
|
---|
[2154] | 30 |
|
---|
[2430] | 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.
|
---|
[2154] | 37 |
|
---|
[2430] | 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.
|
---|
[2154] | 42 |
|
---|
[2430] | 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.
|
---|
[2154] | 45 |
|
---|
[2430] | 46 | N - number of points used. N>=1.
|
---|
| 47 | M - number of basis functions, M>=1.
|
---|
[2154] | 48 |
|
---|
[2430] | 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.
|
---|
[2154] | 53 |
|
---|
[2430] | 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)
|
---|
| 63 | {
|
---|
| 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;
|
---|
[2154] | 82 |
|
---|
[2430] | 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++)
|
---|
[2154] | 94 | {
|
---|
[2430] | 95 | b[i] = w[i-1]*y[i-1];
|
---|
[2154] | 96 | }
|
---|
| 97 | for(i=mi+1; i<=ni; i++)
|
---|
| 98 | {
|
---|
[2430] | 99 | b[i] = 0;
|
---|
[2154] | 100 | }
|
---|
[2430] | 101 | for(j=1; j<=ni; j++)
|
---|
[2154] | 102 | {
|
---|
[2430] | 103 | i1_ = (0) - (1);
|
---|
| 104 | for(i_=1; i_<=mi;i_++)
|
---|
| 105 | {
|
---|
| 106 | a[j,i_] = fmatrix[i_+i1_,j-1];
|
---|
| 107 | }
|
---|
[2154] | 108 | }
|
---|
[2430] | 109 | for(j=1; j<=ni; j++)
|
---|
[2154] | 110 | {
|
---|
[2430] | 111 | for(i=mi+1; i<=ni; i++)
|
---|
| 112 | {
|
---|
| 113 | a[j,i] = 0;
|
---|
| 114 | }
|
---|
[2154] | 115 | }
|
---|
[2430] | 116 | for(j=1; j<=ni; j++)
|
---|
[2154] | 117 | {
|
---|
[2430] | 118 | for(i=1; i<=mi; i++)
|
---|
| 119 | {
|
---|
| 120 | a[j,i] = a[j,i]*w[i-1];
|
---|
| 121 | }
|
---|
[2154] | 122 | }
|
---|
[2430] | 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++)
|
---|
[2154] | 133 | {
|
---|
[2430] | 134 | b2[1,j] = 0;
|
---|
[2154] | 135 | }
|
---|
[2430] | 136 | for(i=1; i<=ni; i++)
|
---|
[2154] | 137 | {
|
---|
[2430] | 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;
|
---|
[2154] | 144 | }
|
---|
[2430] | 145 |
|
---|
| 146 | //
|
---|
| 147 | // Back from A' to A
|
---|
| 148 | // Making cols(A)=rows(A)
|
---|
| 149 | //
|
---|
| 150 | for(i=1; i<=ni-1; i++)
|
---|
[2154] | 151 | {
|
---|
[2430] | 152 | for(i_=i+1; i_<=ni;i_++)
|
---|
[2154] | 153 | {
|
---|
[2430] | 154 | a[i,i_] = a[i_,i];
|
---|
[2154] | 155 | }
|
---|
[2430] | 156 | }
|
---|
| 157 | for(i=2; i<=ni; i++)
|
---|
| 158 | {
|
---|
| 159 | for(j=1; j<=i-1; j++)
|
---|
[2154] | 160 | {
|
---|
[2430] | 161 | a[i,j] = 0;
|
---|
[2154] | 162 | }
|
---|
| 163 | }
|
---|
[2430] | 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) )
|
---|
[2154] | 181 | {
|
---|
[2430] | 182 | for(i=0; i<=ni-1; i++)
|
---|
| 183 | {
|
---|
| 184 | c[i] = 0;
|
---|
| 185 | }
|
---|
| 186 | return;
|
---|
[2154] | 187 | }
|
---|
[2430] | 188 |
|
---|
| 189 | //
|
---|
| 190 | // B2 := (d^(-1) * B2')'
|
---|
| 191 | //
|
---|
[2563] | 192 | if( (double)(d[1])!=(double)(0) )
|
---|
[2430] | 193 | {
|
---|
| 194 | for(i=1; i<=ni; i++)
|
---|
| 195 | {
|
---|
[2563] | 196 | if( (double)(d[i])>(double)(AP.Math.MachineEpsilon*10*Math.Sqrt(ni)*d[1]) )
|
---|
[2430] | 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 | }
|
---|
[2154] | 230 | }
|
---|
| 231 |
|
---|
| 232 |
|
---|
[2430] | 233 | /*************************************************************************
|
---|
| 234 | Linear approximation using least squares method
|
---|
[2154] | 235 |
|
---|
[2430] | 236 | The subroutine calculates coefficients of the line approximating given
|
---|
| 237 | function.
|
---|
[2154] | 238 |
|
---|
[2430] | 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
|
---|
[2154] | 243 |
|
---|
[2430] | 244 | Output parameters:
|
---|
| 245 | a, b- coefficients of linear approximation a+b*t
|
---|
[2154] | 246 |
|
---|
[2430] | 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)
|
---|
| 255 | {
|
---|
| 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;
|
---|
[2154] | 270 |
|
---|
[2430] | 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;
|
---|
[2563] | 288 | if( (double)(Math.Abs(d1))>(double)(Math.Abs(d2)) )
|
---|
[2430] | 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;
|
---|
[2563] | 298 | if( (double)(Math.Abs(d1))>(double)(m*AP.Math.MachineEpsilon*1000) )
|
---|
[2430] | 299 | {
|
---|
| 300 | t1 = t1/d1;
|
---|
| 301 | }
|
---|
| 302 | else
|
---|
| 303 | {
|
---|
| 304 | t1 = 0;
|
---|
| 305 | }
|
---|
[2563] | 306 | if( (double)(Math.Abs(d2))>(double)(m*AP.Math.MachineEpsilon*1000) )
|
---|
[2430] | 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;
|
---|
[2154] | 316 | }
|
---|
| 317 |
|
---|
| 318 |
|
---|
[2430] | 319 | /*************************************************************************
|
---|
| 320 | Weighted cubic spline approximation using linear least squares
|
---|
[2154] | 321 |
|
---|
[2430] | 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.
|
---|
[2154] | 329 |
|
---|
[2430] | 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)
|
---|
| 343 | {
|
---|
| 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;
|
---|
[2154] | 363 |
|
---|
[2430] | 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++)
|
---|
[2154] | 377 | {
|
---|
[2430] | 378 | sx[j] = a+(b-a)*j/(ni-1);
|
---|
[2154] | 379 | }
|
---|
[2430] | 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 | //
|
---|
[2154] | 404 | for(i=0; i<=mi-1; i++)
|
---|
| 405 | {
|
---|
[2430] | 406 | mb[i+1] = w[i]*y[i];
|
---|
[2154] | 407 | }
|
---|
| 408 | for(i=mi+1; i<=ni; i++)
|
---|
| 409 | {
|
---|
[2430] | 410 | mb[i] = 0;
|
---|
[2154] | 411 | }
|
---|
[2430] | 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++)
|
---|
[2154] | 422 | {
|
---|
[2430] | 423 | b2[1,j] = 0;
|
---|
[2154] | 424 | }
|
---|
[2430] | 425 | for(i=1; i<=ni; i++)
|
---|
[2154] | 426 | {
|
---|
[2430] | 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;
|
---|
[2154] | 433 | }
|
---|
[2430] | 434 |
|
---|
| 435 | //
|
---|
| 436 | // Back from A' to A
|
---|
| 437 | // Making cols(A)=rows(A)
|
---|
| 438 | //
|
---|
| 439 | for(i=1; i<=ni-1; i++)
|
---|
[2154] | 440 | {
|
---|
[2430] | 441 | for(i_=i+1; i_<=ni;i_++)
|
---|
| 442 | {
|
---|
| 443 | ma[i,i_] = ma[i_,i];
|
---|
| 444 | }
|
---|
[2154] | 445 | }
|
---|
[2430] | 446 | for(i=2; i<=ni; i++)
|
---|
[2154] | 447 | {
|
---|
[2430] | 448 | for(j=1; j<=i-1; j++)
|
---|
[2154] | 449 | {
|
---|
[2430] | 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++)
|
---|
[2154] | 476 | {
|
---|
[2430] | 477 | if( i==j )
|
---|
| 478 | {
|
---|
| 479 | vt[i,j] = 1;
|
---|
| 480 | }
|
---|
| 481 | else
|
---|
| 482 | {
|
---|
| 483 | vt[i,j] = 0;
|
---|
| 484 | }
|
---|
[2154] | 485 | }
|
---|
| 486 | }
|
---|
[2430] | 487 | b2[1,1] = 1;
|
---|
[2154] | 488 | }
|
---|
[2430] | 489 |
|
---|
| 490 | //
|
---|
| 491 | // B2 := (d^(-1) * B2')'
|
---|
| 492 | //
|
---|
| 493 | for(i=1; i<=ni; i++)
|
---|
[2154] | 494 | {
|
---|
[2563] | 495 | if( (double)(d[i])>(double)(AP.Math.MachineEpsilon*10*Math.Sqrt(ni)*d[1]) )
|
---|
[2430] | 496 | {
|
---|
| 497 | b2[1,i] = b2[1,i]/d[i];
|
---|
| 498 | }
|
---|
| 499 | else
|
---|
| 500 | {
|
---|
| 501 | b2[1,i] = 0;
|
---|
| 502 | }
|
---|
[2154] | 503 | }
|
---|
[2430] | 504 |
|
---|
| 505 | //
|
---|
| 506 | // B := (V * B2')'
|
---|
| 507 | //
|
---|
| 508 | for(i=1; i<=ni; i++)
|
---|
[2154] | 509 | {
|
---|
[2430] | 510 | mb[i] = 0;
|
---|
[2154] | 511 | }
|
---|
[2430] | 512 | for(i=1; i<=ni; i++)
|
---|
[2154] | 513 | {
|
---|
[2430] | 514 | v = b2[1,i];
|
---|
| 515 | for(i_=1; i_<=ni;i_++)
|
---|
| 516 | {
|
---|
| 517 | mb[i_] = mb[i_] + v*vt[i,i_];
|
---|
| 518 | }
|
---|
[2154] | 519 | }
|
---|
[2430] | 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);
|
---|
[2154] | 529 | }
|
---|
| 530 |
|
---|
| 531 |
|
---|
[2430] | 532 | /*************************************************************************
|
---|
| 533 | Polynomial approximation using least squares method
|
---|
[2154] | 534 |
|
---|
[2430] | 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.
|
---|
[2154] | 544 |
|
---|
[2430] | 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
|
---|
[2154] | 550 |
|
---|
[2430] | 551 | Output parameters:
|
---|
| 552 | C - approximating polynomial coefficients, array[0..M],
|
---|
| 553 | C[i] - coefficient at X^i.
|
---|
[2154] | 554 |
|
---|
[2430] | 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)
|
---|
| 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];
|
---|
[2154] | 578 |
|
---|
[2430] | 579 |
|
---|
| 580 | //
|
---|
| 581 | // Initialize
|
---|
| 582 | //
|
---|
| 583 | maxx = x[0];
|
---|
| 584 | minx = x[0];
|
---|
| 585 | for(i=1; i<=n-1; i++)
|
---|
[2154] | 586 | {
|
---|
[2563] | 587 | if( (double)(x[i])>(double)(maxx) )
|
---|
[2430] | 588 | {
|
---|
| 589 | maxx = x[i];
|
---|
| 590 | }
|
---|
[2563] | 591 | if( (double)(x[i])<(double)(minx) )
|
---|
[2430] | 592 | {
|
---|
| 593 | minx = x[i];
|
---|
| 594 | }
|
---|
[2154] | 595 | }
|
---|
[2563] | 596 | if( (double)(minx)==(double)(maxx) )
|
---|
[2154] | 597 | {
|
---|
[2430] | 598 | minx = minx-0.5;
|
---|
| 599 | maxx = maxx+0.5;
|
---|
[2154] | 600 | }
|
---|
[2430] | 601 | w = new double[n-1+1];
|
---|
| 602 | for(i=0; i<=n-1; i++)
|
---|
[2154] | 603 | {
|
---|
[2430] | 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++)
|
---|
[2154] | 624 | {
|
---|
[2430] | 625 | e = c1[k];
|
---|
| 626 | c1[k] = 0;
|
---|
| 627 | if( i<=1 & k==i )
|
---|
[2154] | 628 | {
|
---|
[2430] | 629 | c1[k] = 1;
|
---|
[2154] | 630 | }
|
---|
[2430] | 631 | else
|
---|
[2154] | 632 | {
|
---|
[2430] | 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 | }
|
---|
[2154] | 641 | }
|
---|
[2430] | 642 | d = e;
|
---|
[2154] | 643 | }
|
---|
[2430] | 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;
|
---|
[2154] | 653 | }
|
---|
[2430] | 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++)
|
---|
[2154] | 667 | {
|
---|
[2430] | 668 | z2[i] = 1;
|
---|
| 669 | z1[i] = l2*z1[i-1];
|
---|
| 670 | c[0] = c[0]+c1[i]*z1[i];
|
---|
[2154] | 671 | }
|
---|
[2430] | 672 | for(j=1; j<=m; j++)
|
---|
[2154] | 673 | {
|
---|
[2430] | 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 | }
|
---|
[2154] | 682 | }
|
---|
| 683 | }
|
---|
| 684 |
|
---|
| 685 |
|
---|
[2430] | 686 | /*************************************************************************
|
---|
| 687 | Chebyshev polynomial approximation using least squares method.
|
---|
[2154] | 688 |
|
---|
[2430] | 689 | The algorithm reduces interval [A, B] to the interval [-1,1], then builds
|
---|
| 690 | least squares approximation using Chebyshev polynomials.
|
---|
[2154] | 691 |
|
---|
[2430] | 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.
|
---|
[2154] | 700 |
|
---|
[2430] | 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;
|
---|
[2154] | 733 |
|
---|
[2430] | 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++)
|
---|
[2154] | 744 | {
|
---|
[2430] | 745 | for(i=1; i<=mi; i++)
|
---|
[2154] | 746 | {
|
---|
[2430] | 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 | }
|
---|
[2154] | 760 | }
|
---|
[2430] | 761 | }
|
---|
| 762 | for(j=1; j<=ni; j++)
|
---|
| 763 | {
|
---|
| 764 | for(i=1; i<=mi; i++)
|
---|
[2154] | 765 | {
|
---|
[2430] | 766 | ma[j,i] = w[i-1]*ma[j,i];
|
---|
[2154] | 767 | }
|
---|
[2430] | 768 | }
|
---|
| 769 | for(j=1; j<=ni; j++)
|
---|
| 770 | {
|
---|
| 771 | for(i=mi+1; i<=ni; i++)
|
---|
[2154] | 772 | {
|
---|
[2430] | 773 | ma[j,i] = 0;
|
---|
[2154] | 774 | }
|
---|
| 775 | }
|
---|
[2430] | 776 |
|
---|
| 777 | //
|
---|
| 778 | // Initializing right part
|
---|
| 779 | //
|
---|
| 780 | for(i=0; i<=mi-1; i++)
|
---|
[2154] | 781 | {
|
---|
[2430] | 782 | mb[i+1] = w[i]*y[i];
|
---|
[2154] | 783 | }
|
---|
| 784 | for(i=mi+1; i<=ni; i++)
|
---|
| 785 | {
|
---|
[2430] | 786 | mb[i] = 0;
|
---|
[2154] | 787 | }
|
---|
[2430] | 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++)
|
---|
[2154] | 798 | {
|
---|
[2430] | 799 | b2[1,j] = 0;
|
---|
[2154] | 800 | }
|
---|
[2430] | 801 | for(i=1; i<=ni; i++)
|
---|
[2154] | 802 | {
|
---|
[2430] | 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;
|
---|
[2154] | 809 | }
|
---|
[2430] | 810 |
|
---|
| 811 | //
|
---|
| 812 | // Back from A' to A
|
---|
| 813 | // Making cols(A)=rows(A)
|
---|
| 814 | //
|
---|
| 815 | for(i=1; i<=ni-1; i++)
|
---|
[2154] | 816 | {
|
---|
[2430] | 817 | for(i_=i+1; i_<=ni;i_++)
|
---|
| 818 | {
|
---|
| 819 | ma[i,i_] = ma[i_,i];
|
---|
| 820 | }
|
---|
[2154] | 821 | }
|
---|
[2430] | 822 | for(i=2; i<=ni; i++)
|
---|
[2154] | 823 | {
|
---|
[2430] | 824 | for(j=1; j<=i-1; j++)
|
---|
[2154] | 825 | {
|
---|
[2430] | 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++)
|
---|
[2154] | 852 | {
|
---|
[2430] | 853 | if( i==j )
|
---|
| 854 | {
|
---|
| 855 | vt[i,j] = 1;
|
---|
| 856 | }
|
---|
| 857 | else
|
---|
| 858 | {
|
---|
| 859 | vt[i,j] = 0;
|
---|
| 860 | }
|
---|
[2154] | 861 | }
|
---|
| 862 | }
|
---|
[2430] | 863 | b2[1,1] = 1;
|
---|
[2154] | 864 | }
|
---|
[2430] | 865 |
|
---|
| 866 | //
|
---|
| 867 | // B2 := (d^(-1) * B2')'
|
---|
| 868 | //
|
---|
| 869 | for(i=1; i<=ni; i++)
|
---|
[2154] | 870 | {
|
---|
[2563] | 871 | if( (double)(d[i])>(double)(AP.Math.MachineEpsilon*10*Math.Sqrt(ni)*d[1]) )
|
---|
[2430] | 872 | {
|
---|
| 873 | b2[1,i] = b2[1,i]/d[i];
|
---|
| 874 | }
|
---|
| 875 | else
|
---|
| 876 | {
|
---|
| 877 | b2[1,i] = 0;
|
---|
| 878 | }
|
---|
[2154] | 879 | }
|
---|
[2430] | 880 |
|
---|
| 881 | //
|
---|
| 882 | // B := (V * B2')'
|
---|
| 883 | //
|
---|
| 884 | for(i=1; i<=ni; i++)
|
---|
[2154] | 885 | {
|
---|
[2430] | 886 | mb[i] = 0;
|
---|
[2154] | 887 | }
|
---|
[2430] | 888 | for(i=1; i<=ni; i++)
|
---|
[2154] | 889 | {
|
---|
[2430] | 890 | v = b2[1,i];
|
---|
| 891 | for(i_=1; i_<=ni;i_++)
|
---|
| 892 | {
|
---|
| 893 | mb[i_] = mb[i_] + v*vt[i,i_];
|
---|
| 894 | }
|
---|
[2154] | 895 | }
|
---|
[2430] | 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;
|
---|
[2154] | 907 | }
|
---|
| 908 |
|
---|
| 909 |
|
---|
[2430] | 910 | /*************************************************************************
|
---|
| 911 | Weighted Chebyshev polynomial constrained least squares approximation.
|
---|
[2154] | 912 |
|
---|
[2430] | 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.
|
---|
[2154] | 915 |
|
---|
[2430] | 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.
|
---|
[2154] | 929 |
|
---|
[2430] | 930 | Output parameters:
|
---|
| 931 | CTbl- coefficient table. This parameter is passed into the
|
---|
| 932 | CalculateChebyshevLeastSquares subroutine.
|
---|
[2154] | 933 |
|
---|
[2430] | 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).
|
---|
[2154] | 939 |
|
---|
[2430] | 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).
|
---|
[2154] | 954 |
|
---|
[2430] | 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)
|
---|
| 970 | {
|
---|
| 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;
|
---|
[2154] | 992 |
|
---|
[2430] | 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++)
|
---|
[2154] | 1005 | {
|
---|
[2430] | 1006 | for(j=1; j<=m+1; j++)
|
---|
[2154] | 1007 | {
|
---|
[2430] | 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 | }
|
---|
[2154] | 1021 | }
|
---|
[2430] | 1022 | }
|
---|
| 1023 | for(i=1; i<=n; i++)
|
---|
| 1024 | {
|
---|
| 1025 | for(j=1; j<=m+1; j++)
|
---|
[2154] | 1026 | {
|
---|
[2430] | 1027 | designmatrix[i,j] = w[i-1]*designmatrix[i,j];
|
---|
[2154] | 1028 | }
|
---|
[2430] | 1029 | }
|
---|
| 1030 | for(i=n+1; i<=m+1; i++)
|
---|
| 1031 | {
|
---|
| 1032 | for(j=1; j<=m+1; j++)
|
---|
[2154] | 1033 | {
|
---|
[2430] | 1034 | designmatrix[i,j] = 0;
|
---|
[2154] | 1035 | }
|
---|
| 1036 | }
|
---|
[2430] | 1037 | for(i=0; i<=n-1; i++)
|
---|
[2154] | 1038 | {
|
---|
[2430] | 1039 | rightpart[i+1] = w[i]*y[i];
|
---|
[2154] | 1040 | }
|
---|
[2430] | 1041 | for(i=n+1; i<=m+1; i++)
|
---|
[2154] | 1042 | {
|
---|
[2430] | 1043 | rightpart[i] = 0;
|
---|
[2154] | 1044 | }
|
---|
[2430] | 1045 | n = Math.Max(n, m+1);
|
---|
[2154] | 1046 |
|
---|
| 1047 | //
|
---|
[2430] | 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
|
---|
[2154] | 1051 | //
|
---|
[2430] | 1052 | c = new double[m+1+1, m+1+1];
|
---|
| 1053 | d = new double[m+1+1];
|
---|
| 1054 | if( nc==0 )
|
---|
[2154] | 1055 | {
|
---|
[2430] | 1056 |
|
---|
| 1057 | //
|
---|
| 1058 | // No constraints
|
---|
| 1059 | //
|
---|
| 1060 | for(i=1; i<=m+1; i++)
|
---|
[2154] | 1061 | {
|
---|
[2430] | 1062 | for(j=1; j<=m+1; j++)
|
---|
| 1063 | {
|
---|
| 1064 | c[i,j] = 0;
|
---|
| 1065 | }
|
---|
| 1066 | d[i] = 0;
|
---|
[2154] | 1067 | }
|
---|
[2430] | 1068 | for(i=1; i<=m+1; i++)
|
---|
| 1069 | {
|
---|
| 1070 | c[i,i] = 1;
|
---|
| 1071 | }
|
---|
| 1072 | reducedsize = m+1;
|
---|
[2154] | 1073 | }
|
---|
[2430] | 1074 | else
|
---|
[2154] | 1075 | {
|
---|
[2430] | 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++)
|
---|
[2154] | 1087 | {
|
---|
[2430] | 1088 | v = 2*(xc[i]-a)/(b-a)-1;
|
---|
| 1089 | for(j=0; j<=m; j++)
|
---|
[2154] | 1090 | {
|
---|
[2430] | 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 | }
|
---|
[2154] | 1118 | }
|
---|
[2430] | 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 | }
|
---|
[2563] | 1135 | if( (double)(ws[1])==(double)(0) | (double)(ws[nc])<=(double)(AP.Math.MachineEpsilon*10*Math.Sqrt(nc)*ws[1]) )
|
---|
[2430] | 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_++)
|
---|
[2154] | 1145 | {
|
---|
[2430] | 1146 | c[i_,i] = vt[nc+i,i_];
|
---|
[2154] | 1147 | }
|
---|
[2430] | 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_++)
|
---|
[2154] | 1154 | {
|
---|
[2430] | 1155 | v += u[i_,i]*cr[i_];
|
---|
[2154] | 1156 | }
|
---|
[2430] | 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_++)
|
---|
[2154] | 1167 | {
|
---|
[2430] | 1168 | d[i_] = d[i_] + v*vt[i,i_];
|
---|
[2154] | 1169 | }
|
---|
[2430] | 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_++)
|
---|
[2154] | 1181 | {
|
---|
[2430] | 1182 | v += designmatrix[i,i_]*d[i_];
|
---|
[2154] | 1183 | }
|
---|
[2430] | 1184 | rightpart[i] = rightpart[i]-v;
|
---|
[2154] | 1185 | }
|
---|
[2430] | 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);
|
---|
[2154] | 1191 | }
|
---|
| 1192 |
|
---|
| 1193 | //
|
---|
[2430] | 1194 | // Solve reduced problem DesignMatrix*t = RightPart.
|
---|
[2154] | 1195 | //
|
---|
[2430] | 1196 | if( !svd.svddecomposition(designmatrix, n, reducedsize, 1, 1, 2, ref ws, ref u, ref vt) )
|
---|
[2154] | 1197 | {
|
---|
| 1198 | result = false;
|
---|
| 1199 | return result;
|
---|
| 1200 | }
|
---|
[2430] | 1201 | tmp = new double[reducedsize+1];
|
---|
| 1202 | tmp2 = new double[reducedsize+1];
|
---|
| 1203 | for(i=1; i<=reducedsize; i++)
|
---|
[2154] | 1204 | {
|
---|
[2430] | 1205 | tmp[i] = 0;
|
---|
[2154] | 1206 | }
|
---|
[2430] | 1207 | for(i=1; i<=n; i++)
|
---|
[2154] | 1208 | {
|
---|
[2430] | 1209 | v = rightpart[i];
|
---|
| 1210 | for(i_=1; i_<=reducedsize;i_++)
|
---|
[2154] | 1211 | {
|
---|
[2430] | 1212 | tmp[i_] = tmp[i_] + v*u[i,i_];
|
---|
[2154] | 1213 | }
|
---|
| 1214 | }
|
---|
[2430] | 1215 | for(i=1; i<=reducedsize; i++)
|
---|
[2154] | 1216 | {
|
---|
[2563] | 1217 | if( (double)(ws[i])!=(double)(0) & (double)(ws[i])>(double)(AP.Math.MachineEpsilon*10*Math.Sqrt(nc)*ws[1]) )
|
---|
[2154] | 1218 | {
|
---|
[2430] | 1219 | tmp[i] = tmp[i]/ws[i];
|
---|
[2154] | 1220 | }
|
---|
[2430] | 1221 | else
|
---|
| 1222 | {
|
---|
| 1223 | tmp[i] = 0;
|
---|
| 1224 | }
|
---|
[2154] | 1225 | }
|
---|
[2430] | 1226 | for(i=1; i<=reducedsize; i++)
|
---|
[2154] | 1227 | {
|
---|
[2430] | 1228 | tmp2[i] = 0;
|
---|
[2154] | 1229 | }
|
---|
[2430] | 1230 | for(i=1; i<=reducedsize; i++)
|
---|
[2154] | 1231 | {
|
---|
| 1232 | v = tmp[i];
|
---|
[2430] | 1233 | for(i_=1; i_<=reducedsize;i_++)
|
---|
[2154] | 1234 | {
|
---|
[2430] | 1235 | tmp2[i_] = tmp2[i_] + v*vt[i,i_];
|
---|
[2154] | 1236 | }
|
---|
| 1237 | }
|
---|
| 1238 |
|
---|
| 1239 | //
|
---|
[2430] | 1240 | // Solution is in the tmp2.
|
---|
| 1241 | // Transform it from t to x.
|
---|
[2154] | 1242 | //
|
---|
[2430] | 1243 | ctbl = new double[m+2+1];
|
---|
| 1244 | for(i=1; i<=m+1; i++)
|
---|
[2154] | 1245 | {
|
---|
| 1246 | v = 0.0;
|
---|
[2430] | 1247 | for(i_=1; i_<=reducedsize;i_++)
|
---|
[2154] | 1248 | {
|
---|
[2430] | 1249 | v += c[i,i_]*tmp2[i_];
|
---|
[2154] | 1250 | }
|
---|
[2430] | 1251 | ctbl[i-1] = v+d[i];
|
---|
[2154] | 1252 | }
|
---|
[2430] | 1253 | ctbl[m+1] = a;
|
---|
| 1254 | ctbl[m+2] = b;
|
---|
[2154] | 1255 | return result;
|
---|
| 1256 | }
|
---|
| 1257 |
|
---|
| 1258 |
|
---|
[2430] | 1259 | /*************************************************************************
|
---|
| 1260 | Calculation of a Chebyshev polynomial obtained during least squares
|
---|
| 1261 | approximaion at the given point.
|
---|
[2154] | 1262 |
|
---|
[2430] | 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.
|
---|
[2154] | 1271 |
|
---|
[2430] | 1272 | The result is the value at the given point.
|
---|
[2154] | 1273 |
|
---|
[2430] | 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)
|
---|
| 1283 | {
|
---|
| 1284 | double result = 0;
|
---|
| 1285 | double b1 = 0;
|
---|
| 1286 | double b2 = 0;
|
---|
| 1287 | int i = 0;
|
---|
[2154] | 1288 |
|
---|
[2430] | 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;
|
---|
[2154] | 1303 | }
|
---|
| 1304 | }
|
---|
| 1305 | }
|
---|