/************************************************************************* Copyright (c) 2006-2007, Sergey Bochkanov (ALGLIB project). >>> SOURCE LICENSE >>> This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation (www.fsf.org); either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. A copy of the GNU General Public License is available at http://www.fsf.org/licensing/licenses >>> END OF LICENSE >>> *************************************************************************/ using System; namespace alglib { public class leastsquares { /************************************************************************* Weighted approximation by arbitrary function basis in a space of arbitrary dimension using linear least squares method. Input parameters: Y - array[0..N-1] It contains a set of function values in N points. Space dimension and points don't matter. Procedure works with function values in these points and values of basis functions only. W - array[0..N-1] It contains weights corresponding to function values. Each summand in square sum of approximation deviations from given values is multiplied by the square of corresponding weight. FMatrix-a table of basis functions values, array[0..N-1, 0..M-1]. FMatrix[I, J] - value of J-th basis function in I-th point. N - number of points used. N>=1. M - number of basis functions, M>=1. Output parameters: C - decomposition coefficients. Array of real numbers whose index goes from 0 to M-1. C[j] - j-th basis function coefficient. -- ALGLIB -- Copyright by Bochkanov Sergey *************************************************************************/ public static void buildgeneralleastsquares(ref double[] y, ref double[] w, ref double[,] fmatrix, int n, int m, ref double[] c) { int i = 0; int j = 0; double[,] a = new double[0,0]; double[,] q = new double[0,0]; double[,] vt = new double[0,0]; double[] b = new double[0]; double[] tau = new double[0]; double[,] b2 = new double[0,0]; double[] tauq = new double[0]; double[] taup = new double[0]; double[] d = new double[0]; double[] e = new double[0]; bool isuppera = new bool(); int mi = 0; int ni = 0; double v = 0; int i_ = 0; int i1_ = 0; mi = n; ni = m; c = new double[ni-1+1]; // // Initialize design matrix. // Here we are making MI>=NI. // a = new double[ni+1, Math.Max(mi, ni)+1]; b = new double[Math.Max(mi, ni)+1]; for(i=1; i<=mi; i++) { b[i] = w[i-1]*y[i-1]; } for(i=mi+1; i<=ni; i++) { b[i] = 0; } for(j=1; j<=ni; j++) { i1_ = (0) - (1); for(i_=1; i_<=mi;i_++) { a[j,i_] = fmatrix[i_+i1_,j-1]; } } for(j=1; j<=ni; j++) { for(i=mi+1; i<=ni; i++) { a[j,i] = 0; } } for(j=1; j<=ni; j++) { for(i=1; i<=mi; i++) { a[j,i] = a[j,i]*w[i-1]; } } mi = Math.Max(mi, ni); // // LQ-decomposition of A' // B2 := Q*B // lq.lqdecomposition(ref a, ni, mi, ref tau); lq.unpackqfromlq(ref a, ni, mi, ref tau, ni, ref q); b2 = new double[1+1, ni+1]; for(j=1; j<=ni; j++) { b2[1,j] = 0; } for(i=1; i<=ni; i++) { v = 0.0; for(i_=1; i_<=mi;i_++) { v += b[i_]*q[i,i_]; } b2[1,i] = v; } // // Back from A' to A // Making cols(A)=rows(A) // for(i=1; i<=ni-1; i++) { for(i_=i+1; i_<=ni;i_++) { a[i,i_] = a[i_,i]; } } for(i=2; i<=ni; i++) { for(j=1; j<=i-1; j++) { a[i,j] = 0; } } // // Bidiagonal decomposition of A // A = Q * d2 * P' // B2 := (Q'*B2')' // bidiagonal.tobidiagonal(ref a, ni, ni, ref tauq, ref taup); bidiagonal.multiplybyqfrombidiagonal(ref a, ni, ni, ref tauq, ref b2, 1, ni, true, false); bidiagonal.unpackptfrombidiagonal(ref a, ni, ni, ref taup, ni, ref vt); bidiagonal.unpackdiagonalsfrombidiagonal(ref a, ni, ni, ref isuppera, ref d, ref e); // // Singular value decomposition of A // A = U * d * V' // B2 := (U'*B2')' // if( !bdsvd.bidiagonalsvddecomposition(ref d, e, ni, isuppera, false, ref b2, 1, ref q, 0, ref vt, ni) ) { for(i=0; i<=ni-1; i++) { c[i] = 0; } return; } // // B2 := (d^(-1) * B2')' // if( (double)(d[1])!=(double)(0) ) { for(i=1; i<=ni; i++) { if( (double)(d[i])>(double)(AP.Math.MachineEpsilon*10*Math.Sqrt(ni)*d[1]) ) { b2[1,i] = b2[1,i]/d[i]; } else { b2[1,i] = 0; } } } // // B := (V * B2')' // for(i=1; i<=ni; i++) { b[i] = 0; } for(i=1; i<=ni; i++) { v = b2[1,i]; for(i_=1; i_<=ni;i_++) { b[i_] = b[i_] + v*vt[i,i_]; } } // // Out // for(i=1; i<=ni; i++) { c[i-1] = b[i]; } } /************************************************************************* Linear approximation using least squares method The subroutine calculates coefficients of the line approximating given function. Input parameters: X - array[0..N-1], it contains a set of abscissas. Y - array[0..N-1], function values. N - number of points, N>=1 Output parameters: a, b- coefficients of linear approximation a+b*t -- ALGLIB -- Copyright by Bochkanov Sergey *************************************************************************/ public static void buildlinearleastsquares(ref double[] x, ref double[] y, int n, ref double a, ref double b) { double pp = 0; double qq = 0; double pq = 0; double b1 = 0; double b2 = 0; double d1 = 0; double d2 = 0; double t1 = 0; double t2 = 0; double phi = 0; double c = 0; double s = 0; double m = 0; int i = 0; pp = n; qq = 0; pq = 0; b1 = 0; b2 = 0; for(i=0; i<=n-1; i++) { pq = pq+x[i]; qq = qq+AP.Math.Sqr(x[i]); b1 = b1+y[i]; b2 = b2+x[i]*y[i]; } phi = Math.Atan2(2*pq, qq-pp)/2; c = Math.Cos(phi); s = Math.Sin(phi); d1 = AP.Math.Sqr(c)*pp+AP.Math.Sqr(s)*qq-2*s*c*pq; d2 = AP.Math.Sqr(s)*pp+AP.Math.Sqr(c)*qq+2*s*c*pq; if( (double)(Math.Abs(d1))>(double)(Math.Abs(d2)) ) { m = Math.Abs(d1); } else { m = Math.Abs(d2); } t1 = c*b1-s*b2; t2 = s*b1+c*b2; if( (double)(Math.Abs(d1))>(double)(m*AP.Math.MachineEpsilon*1000) ) { t1 = t1/d1; } else { t1 = 0; } if( (double)(Math.Abs(d2))>(double)(m*AP.Math.MachineEpsilon*1000) ) { t2 = t2/d2; } else { t2 = 0; } a = c*t1+s*t2; b = -(s*t1)+c*t2; } /************************************************************************* Weighted cubic spline approximation using linear least squares Input parameters: X - array[0..N-1], abscissas Y - array[0..N-1], function values W - array[0..N-1], weights. A, B- interval to build splines in. N - number of points used. N>=1. M - number of basic splines, M>=2. Output parameters: CTbl- coefficients table to be used by SplineInterpolation function. -- ALGLIB -- Copyright by Bochkanov Sergey *************************************************************************/ public static void buildsplineleastsquares(ref double[] x, ref double[] y, ref double[] w, double a, double b, int n, int m, ref double[] ctbl) { int i = 0; int j = 0; double[,] ma = new double[0,0]; double[,] q = new double[0,0]; double[,] vt = new double[0,0]; double[] mb = new double[0]; double[] tau = new double[0]; double[,] b2 = new double[0,0]; double[] tauq = new double[0]; double[] taup = new double[0]; double[] d = new double[0]; double[] e = new double[0]; bool isuppera = new bool(); int mi = 0; int ni = 0; double v = 0; double[] sx = new double[0]; double[] sy = new double[0]; int i_ = 0; System.Diagnostics.Debug.Assert(m>=2, "BuildSplineLeastSquares: M is too small!"); mi = n; ni = m; sx = new double[ni-1+1]; sy = new double[ni-1+1]; // // Initializing design matrix // Here we are making MI>=NI // ma = new double[ni+1, Math.Max(mi, ni)+1]; mb = new double[Math.Max(mi, ni)+1]; for(j=0; j<=ni-1; j++) { sx[j] = a+(b-a)*j/(ni-1); } for(j=0; j<=ni-1; j++) { for(i=0; i<=ni-1; i++) { sy[i] = 0; } sy[j] = 1; spline3.buildcubicspline(sx, sy, ni, 0, 0.0, 0, 0.0, ref ctbl); for(i=0; i<=mi-1; i++) { ma[j+1,i+1] = w[i]*spline3.splineinterpolation(ref ctbl, x[i]); } } for(j=1; j<=ni; j++) { for(i=mi+1; i<=ni; i++) { ma[j,i] = 0; } } // // Initializing right part // for(i=0; i<=mi-1; i++) { mb[i+1] = w[i]*y[i]; } for(i=mi+1; i<=ni; i++) { mb[i] = 0; } mi = Math.Max(mi, ni); // // LQ-decomposition of A' // B2 := Q*B // lq.lqdecomposition(ref ma, ni, mi, ref tau); lq.unpackqfromlq(ref ma, ni, mi, ref tau, ni, ref q); b2 = new double[1+1, ni+1]; for(j=1; j<=ni; j++) { b2[1,j] = 0; } for(i=1; i<=ni; i++) { v = 0.0; for(i_=1; i_<=mi;i_++) { v += mb[i_]*q[i,i_]; } b2[1,i] = v; } // // Back from A' to A // Making cols(A)=rows(A) // for(i=1; i<=ni-1; i++) { for(i_=i+1; i_<=ni;i_++) { ma[i,i_] = ma[i_,i]; } } for(i=2; i<=ni; i++) { for(j=1; j<=i-1; j++) { ma[i,j] = 0; } } // // Bidiagonal decomposition of A // A = Q * d2 * P' // B2 := (Q'*B2')' // bidiagonal.tobidiagonal(ref ma, ni, ni, ref tauq, ref taup); bidiagonal.multiplybyqfrombidiagonal(ref ma, ni, ni, ref tauq, ref b2, 1, ni, true, false); bidiagonal.unpackptfrombidiagonal(ref ma, ni, ni, ref taup, ni, ref vt); bidiagonal.unpackdiagonalsfrombidiagonal(ref ma, ni, ni, ref isuppera, ref d, ref e); // // Singular value decomposition of A // A = U * d * V' // B2 := (U'*B2')' // if( !bdsvd.bidiagonalsvddecomposition(ref d, e, ni, isuppera, false, ref b2, 1, ref q, 0, ref vt, ni) ) { for(i=1; i<=ni; i++) { d[i] = 0; b2[1,i] = 0; for(j=1; j<=ni; j++) { if( i==j ) { vt[i,j] = 1; } else { vt[i,j] = 0; } } } b2[1,1] = 1; } // // B2 := (d^(-1) * B2')' // for(i=1; i<=ni; i++) { if( (double)(d[i])>(double)(AP.Math.MachineEpsilon*10*Math.Sqrt(ni)*d[1]) ) { b2[1,i] = b2[1,i]/d[i]; } else { b2[1,i] = 0; } } // // B := (V * B2')' // for(i=1; i<=ni; i++) { mb[i] = 0; } for(i=1; i<=ni; i++) { v = b2[1,i]; for(i_=1; i_<=ni;i_++) { mb[i_] = mb[i_] + v*vt[i,i_]; } } // // Forming result spline // for(i=0; i<=ni-1; i++) { sy[i] = mb[i+1]; } spline3.buildcubicspline(sx, sy, ni, 0, 0.0, 0, 0.0, ref ctbl); } /************************************************************************* Polynomial approximation using least squares method The subroutine calculates coefficients of the polynomial approximating given function. It is recommended to use this function only if you need to obtain coefficients of approximation polynomial. If you have to build and calculate polynomial approximation (NOT coefficients), it's better to use BuildChebyshevLeastSquares subroutine in combination with CalculateChebyshevLeastSquares subroutine. The result of Chebyshev polynomial approximation is equivalent to the result obtained using powers of X, but has higher accuracy due to better numerical properties of Chebyshev polynomials. Input parameters: X - array[0..N-1], abscissas Y - array[0..N-1], function values N - number of points, N>=1 M - order of polynomial required, M>=0 Output parameters: C - approximating polynomial coefficients, array[0..M], C[i] - coefficient at X^i. -- ALGLIB -- Copyright by Bochkanov Sergey *************************************************************************/ public static void buildpolynomialleastsquares(ref double[] x, ref double[] y, int n, int m, ref double[] c) { double[] ctbl = new double[0]; double[] w = new double[0]; double[] c1 = new double[0]; double maxx = 0; double minx = 0; int i = 0; int j = 0; int k = 0; double e = 0; double d = 0; double l1 = 0; double l2 = 0; double[] z2 = new double[0]; double[] z1 = new double[0]; // // Initialize // maxx = x[0]; minx = x[0]; for(i=1; i<=n-1; i++) { if( (double)(x[i])>(double)(maxx) ) { maxx = x[i]; } if( (double)(x[i])<(double)(minx) ) { minx = x[i]; } } if( (double)(minx)==(double)(maxx) ) { minx = minx-0.5; maxx = maxx+0.5; } w = new double[n-1+1]; for(i=0; i<=n-1; i++) { w[i] = 1; } // // Build Chebyshev approximation // buildchebyshevleastsquares(ref x, ref y, ref w, minx, maxx, n, m, ref ctbl); // // From Chebyshev to powers of X // c1 = new double[m+1]; for(i=0; i<=m; i++) { c1[i] = 0; } d = 0; for(i=0; i<=m; i++) { for(k=i; k<=m; k++) { e = c1[k]; c1[k] = 0; if( i<=1 & k==i ) { c1[k] = 1; } else { if( i!=0 ) { c1[k] = 2*d; } if( k>i+1 ) { c1[k] = c1[k]-c1[k-2]; } } d = e; } d = c1[i]; e = 0; k = i; while( k<=m ) { e = e+c1[k]*ctbl[k]; k = k+2; } c1[i] = e; } // // Linear translation // l1 = 2/(ctbl[m+2]-ctbl[m+1]); l2 = -(2*ctbl[m+1]/(ctbl[m+2]-ctbl[m+1]))-1; c = new double[m+1]; z2 = new double[m+1]; z1 = new double[m+1]; c[0] = c1[0]; z1[0] = 1; z2[0] = 1; for(i=1; i<=m; i++) { z2[i] = 1; z1[i] = l2*z1[i-1]; c[0] = c[0]+c1[i]*z1[i]; } for(j=1; j<=m; j++) { z2[0] = l1*z2[0]; c[j] = c1[j]*z2[0]; for(i=j+1; i<=m; i++) { k = i-j; z2[k] = l1*z2[k]+z2[k-1]; c[j] = c[j]+c1[i]*z2[k]*z1[k]; } } } /************************************************************************* Chebyshev polynomial approximation using least squares method. The algorithm reduces interval [A, B] to the interval [-1,1], then builds least squares approximation using Chebyshev polynomials. Input parameters: X - array[0..N-1], abscissas Y - array[0..N-1], function values W - array[0..N-1], weights A, B- interval to build approximating polynomials in. N - number of points used. N>=1. M - order of polynomial, M>=0. This parameter is passed into CalculateChebyshevLeastSquares function. Output parameters: CTbl - coefficient table. This parameter is passed into CalculateChebyshevLeastSquares function. -- ALGLIB -- Copyright by Bochkanov Sergey *************************************************************************/ public static void buildchebyshevleastsquares(ref double[] x, ref double[] y, ref double[] w, double a, double b, int n, int m, ref double[] ctbl) { int i = 0; int j = 0; double[,] ma = new double[0,0]; double[,] q = new double[0,0]; double[,] vt = new double[0,0]; double[] mb = new double[0]; double[] tau = new double[0]; double[,] b2 = new double[0,0]; double[] tauq = new double[0]; double[] taup = new double[0]; double[] d = new double[0]; double[] e = new double[0]; bool isuppera = new bool(); int mi = 0; int ni = 0; double v = 0; int i_ = 0; mi = n; ni = m+1; // // Initializing design matrix // Here we are making MI>=NI // ma = new double[ni+1, Math.Max(mi, ni)+1]; mb = new double[Math.Max(mi, ni)+1]; for(j=1; j<=ni; j++) { for(i=1; i<=mi; i++) { v = 2*(x[i-1]-a)/(b-a)-1; if( j==1 ) { ma[j,i] = 1.0; } if( j==2 ) { ma[j,i] = v; } if( j>2 ) { ma[j,i] = 2.0*v*ma[j-1,i]-ma[j-2,i]; } } } for(j=1; j<=ni; j++) { for(i=1; i<=mi; i++) { ma[j,i] = w[i-1]*ma[j,i]; } } for(j=1; j<=ni; j++) { for(i=mi+1; i<=ni; i++) { ma[j,i] = 0; } } // // Initializing right part // for(i=0; i<=mi-1; i++) { mb[i+1] = w[i]*y[i]; } for(i=mi+1; i<=ni; i++) { mb[i] = 0; } mi = Math.Max(mi, ni); // // LQ-decomposition of A' // B2 := Q*B // lq.lqdecomposition(ref ma, ni, mi, ref tau); lq.unpackqfromlq(ref ma, ni, mi, ref tau, ni, ref q); b2 = new double[1+1, ni+1]; for(j=1; j<=ni; j++) { b2[1,j] = 0; } for(i=1; i<=ni; i++) { v = 0.0; for(i_=1; i_<=mi;i_++) { v += mb[i_]*q[i,i_]; } b2[1,i] = v; } // // Back from A' to A // Making cols(A)=rows(A) // for(i=1; i<=ni-1; i++) { for(i_=i+1; i_<=ni;i_++) { ma[i,i_] = ma[i_,i]; } } for(i=2; i<=ni; i++) { for(j=1; j<=i-1; j++) { ma[i,j] = 0; } } // // Bidiagonal decomposition of A // A = Q * d2 * P' // B2 := (Q'*B2')' // bidiagonal.tobidiagonal(ref ma, ni, ni, ref tauq, ref taup); bidiagonal.multiplybyqfrombidiagonal(ref ma, ni, ni, ref tauq, ref b2, 1, ni, true, false); bidiagonal.unpackptfrombidiagonal(ref ma, ni, ni, ref taup, ni, ref vt); bidiagonal.unpackdiagonalsfrombidiagonal(ref ma, ni, ni, ref isuppera, ref d, ref e); // // Singular value decomposition of A // A = U * d * V' // B2 := (U'*B2')' // if( !bdsvd.bidiagonalsvddecomposition(ref d, e, ni, isuppera, false, ref b2, 1, ref q, 0, ref vt, ni) ) { for(i=1; i<=ni; i++) { d[i] = 0; b2[1,i] = 0; for(j=1; j<=ni; j++) { if( i==j ) { vt[i,j] = 1; } else { vt[i,j] = 0; } } } b2[1,1] = 1; } // // B2 := (d^(-1) * B2')' // for(i=1; i<=ni; i++) { if( (double)(d[i])>(double)(AP.Math.MachineEpsilon*10*Math.Sqrt(ni)*d[1]) ) { b2[1,i] = b2[1,i]/d[i]; } else { b2[1,i] = 0; } } // // B := (V * B2')' // for(i=1; i<=ni; i++) { mb[i] = 0; } for(i=1; i<=ni; i++) { v = b2[1,i]; for(i_=1; i_<=ni;i_++) { mb[i_] = mb[i_] + v*vt[i,i_]; } } // // Forming result // ctbl = new double[ni+1+1]; for(i=1; i<=ni; i++) { ctbl[i-1] = mb[i]; } ctbl[ni] = a; ctbl[ni+1] = b; } /************************************************************************* Weighted Chebyshev polynomial constrained least squares approximation. The algorithm reduces [A,B] to [-1,1] and builds the Chebyshev polynomials series by approximating a given function using the least squares method. Input parameters: X - abscissas, array[0..N-1] Y - function values, array[0..N-1] W - weights, array[0..N-1]. Each item in the squared sum of deviations from given values is multiplied by a square of corresponding weight. A, B- interval in which the approximating polynomials are built. N - number of points, N>0. XC, YC, DC- constraints (see description below)., array[0..NC-1] NC - number of constraints. 0 <= NC < M+1. M - degree of polynomial, M>=0. This parameter is passed into the CalculateChebyshevLeastSquares subroutine. Output parameters: CTbl- coefficient table. This parameter is passed into the CalculateChebyshevLeastSquares subroutine. Result: True, if the algorithm succeeded. False, if the internal singular value decomposition subroutine hasn't converged or the given constraints could not be met simultaneously (e.g. P(0)=0 è P(0)=1). Specifying constraints: This subroutine can solve the problem having constrained function values or its derivatives in several points. NC specifies the number of constraints, DC - the type of constraints, XC and YC - constraints as such. Thus, for each i from 0 to NC-1 the following constraint is given: P(xc[i]) = yc[i], if DC[i]=0 or d/dx(P(xc[i])) = yc[i], if DC[i]=1 (here P(x) is approximating polynomial). This version of the subroutine supports only either polynomial or its derivative value constraints. If DC[i] is not equal to 0 and 1, the subroutine will be aborted. The number of constraints should be less than the number of degrees of freedom of approximating polynomial - M+1 (at that, it could be equal to 0). -- ALGLIB -- Copyright by Bochkanov Sergey *************************************************************************/ public static bool buildchebyshevleastsquaresconstrained(ref double[] x, ref double[] y, ref double[] w, double a, double b, int n, ref double[] xc, ref double[] yc, ref int[] dc, int nc, int m, ref double[] ctbl) { bool result = new bool(); int i = 0; int j = 0; int reducedsize = 0; double[,] designmatrix = new double[0,0]; double[] rightpart = new double[0]; double[,] cmatrix = new double[0,0]; double[,] c = new double[0,0]; double[,] u = new double[0,0]; double[,] vt = new double[0,0]; double[] d = new double[0]; double[] cr = new double[0]; double[] ws = new double[0]; double[] tj = new double[0]; double[] uj = new double[0]; double[] dtj = new double[0]; double[] tmp = new double[0]; double[] tmp2 = new double[0]; double[,] tmpmatrix = new double[0,0]; double v = 0; int i_ = 0; System.Diagnostics.Debug.Assert(n>0); System.Diagnostics.Debug.Assert(m>=0); System.Diagnostics.Debug.Assert(nc>=0 & nc=M+1. // designmatrix = new double[Math.Max(n, m+1)+1, m+1+1]; rightpart = new double[Math.Max(n, m+1)+1]; for(i=1; i<=n; i++) { for(j=1; j<=m+1; j++) { v = 2*(x[i-1]-a)/(b-a)-1; if( j==1 ) { designmatrix[i,j] = 1.0; } if( j==2 ) { designmatrix[i,j] = v; } if( j>2 ) { designmatrix[i,j] = 2.0*v*designmatrix[i,j-1]-designmatrix[i,j-2]; } } } for(i=1; i<=n; i++) { for(j=1; j<=m+1; j++) { designmatrix[i,j] = w[i-1]*designmatrix[i,j]; } } for(i=n+1; i<=m+1; i++) { for(j=1; j<=m+1; j++) { designmatrix[i,j] = 0; } } for(i=0; i<=n-1; i++) { rightpart[i+1] = w[i]*y[i]; } for(i=n+1; i<=m+1; i++) { rightpart[i] = 0; } n = Math.Max(n, m+1); // // Now N>=M+1 and we are ready to the next stage. // Handle constraints. // Represent feasible set of coefficients as x = C*t + d // c = new double[m+1+1, m+1+1]; d = new double[m+1+1]; if( nc==0 ) { // // No constraints // for(i=1; i<=m+1; i++) { for(j=1; j<=m+1; j++) { c[i,j] = 0; } d[i] = 0; } for(i=1; i<=m+1; i++) { c[i,i] = 1; } reducedsize = m+1; } else { // // Constraints are present. // Fill constraints matrix CMatrix and solve CMatrix*x = cr. // cmatrix = new double[nc+1, m+1+1]; cr = new double[nc+1]; tj = new double[m+1]; uj = new double[m+1]; dtj = new double[m+1]; for(i=0; i<=nc-1; i++) { v = 2*(xc[i]-a)/(b-a)-1; for(j=0; j<=m; j++) { if( j==0 ) { tj[j] = 1; uj[j] = 1; dtj[j] = 0; } if( j==1 ) { tj[j] = v; uj[j] = 2*v; dtj[j] = 1; } if( j>1 ) { tj[j] = 2*v*tj[j-1]-tj[j-2]; uj[j] = 2*v*uj[j-1]-uj[j-2]; dtj[j] = j*uj[j-1]; } System.Diagnostics.Debug.Assert(dc[i]==0 | dc[i]==1); if( dc[i]==0 ) { cmatrix[i+1,j+1] = tj[j]; } if( dc[i]==1 ) { cmatrix[i+1,j+1] = dtj[j]; } } cr[i+1] = yc[i]; } // // Solve CMatrix*x = cr. // Fill C and d: // 1. SVD: CMatrix = U * WS * V^T // 2. C := V[1:M+1,NC+1:M+1] // 3. tmp := WS^-1 * U^T * cr // 4. d := V[1:M+1,1:NC] * tmp // if( !svd.svddecomposition(cmatrix, nc, m+1, 2, 2, 2, ref ws, ref u, ref vt) ) { result = false; return result; } if( (double)(ws[1])==(double)(0) | (double)(ws[nc])<=(double)(AP.Math.MachineEpsilon*10*Math.Sqrt(nc)*ws[1]) ) { result = false; return result; } c = new double[m+1+1, m+1-nc+1]; d = new double[m+1+1]; for(i=1; i<=m+1-nc; i++) { for(i_=1; i_<=m+1;i_++) { c[i_,i] = vt[nc+i,i_]; } } tmp = new double[nc+1]; for(i=1; i<=nc; i++) { v = 0.0; for(i_=1; i_<=nc;i_++) { v += u[i_,i]*cr[i_]; } tmp[i] = v/ws[i]; } for(i=1; i<=m+1; i++) { d[i] = 0; } for(i=1; i<=nc; i++) { v = tmp[i]; for(i_=1; i_<=m+1;i_++) { d[i_] = d[i_] + v*vt[i,i_]; } } // // Reduce problem: // 1. RightPart := RightPart - DesignMatrix*d // 2. DesignMatrix := DesignMatrix*C // for(i=1; i<=n; i++) { v = 0.0; for(i_=1; i_<=m+1;i_++) { v += designmatrix[i,i_]*d[i_]; } rightpart[i] = rightpart[i]-v; } reducedsize = m+1-nc; tmpmatrix = new double[n+1, reducedsize+1]; tmp = new double[n+1]; 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); blas.copymatrix(ref tmpmatrix, 1, n, 1, reducedsize, ref designmatrix, 1, n, 1, reducedsize); } // // Solve reduced problem DesignMatrix*t = RightPart. // if( !svd.svddecomposition(designmatrix, n, reducedsize, 1, 1, 2, ref ws, ref u, ref vt) ) { result = false; return result; } tmp = new double[reducedsize+1]; tmp2 = new double[reducedsize+1]; for(i=1; i<=reducedsize; i++) { tmp[i] = 0; } for(i=1; i<=n; i++) { v = rightpart[i]; for(i_=1; i_<=reducedsize;i_++) { tmp[i_] = tmp[i_] + v*u[i,i_]; } } for(i=1; i<=reducedsize; i++) { if( (double)(ws[i])!=(double)(0) & (double)(ws[i])>(double)(AP.Math.MachineEpsilon*10*Math.Sqrt(nc)*ws[1]) ) { tmp[i] = tmp[i]/ws[i]; } else { tmp[i] = 0; } } for(i=1; i<=reducedsize; i++) { tmp2[i] = 0; } for(i=1; i<=reducedsize; i++) { v = tmp[i]; for(i_=1; i_<=reducedsize;i_++) { tmp2[i_] = tmp2[i_] + v*vt[i,i_]; } } // // Solution is in the tmp2. // Transform it from t to x. // ctbl = new double[m+2+1]; for(i=1; i<=m+1; i++) { v = 0.0; for(i_=1; i_<=reducedsize;i_++) { v += c[i,i_]*tmp2[i_]; } ctbl[i-1] = v+d[i]; } ctbl[m+1] = a; ctbl[m+2] = b; return result; } /************************************************************************* Calculation of a Chebyshev polynomial obtained during least squares approximaion at the given point. Input parameters: M - order of polynomial (parameter of the BuildChebyshevLeastSquares function). A - coefficient table. A[0..M] contains coefficients of the i-th Chebyshev polynomial. A[M+1] contains left boundary of approximation interval. A[M+2] contains right boundary of approximation interval. X - point to perform calculations in. The result is the value at the given point. It should be noted that array A contains coefficients of the Chebyshev polynomials defined on interval [-1,1]. Argument is reduced to this interval before calculating polynomial value. -- ALGLIB -- Copyright by Bochkanov Sergey *************************************************************************/ public static double calculatechebyshevleastsquares(int m, ref double[] a, double x) { double result = 0; double b1 = 0; double b2 = 0; int i = 0; x = 2*(x-a[m+1])/(a[m+2]-a[m+1])-1; b1 = 0; b2 = 0; i = m; do { result = 2*x*b1-b2+a[i]; b2 = b1; b1 = result; i = i-1; } while( i>=0 ); result = result-x*b2; return result; } } }