Free cookie consent management tool by TermsFeed Policy Generator

source: branches/Persistence Test/ALGLIB/leastsquares.cs @ 2497

Last change on this file since 2497 was 2430, checked in by gkronber, 15 years ago

Moved ALGLIB code into a separate plugin. #783

File size: 42.2 KB
Line 
1/*************************************************************************
2Copyright (c) 2006-2007, Sergey Bochkanov (ALGLIB project).
3
4>>> SOURCE LICENSE >>>
5This program is free software; you can redistribute it and/or modify
6it under the terms of the GNU General Public License as published by
7the Free Software Foundation (www.fsf.org); either version 2 of the
8License, or (at your option) any later version.
9
10This program is distributed in the hope that it will be useful,
11but WITHOUT ANY WARRANTY; without even the implied warranty of
12MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13GNU General Public License for more details.
14
15A copy of the GNU General Public License is available at
16http://www.fsf.org/licensing/licenses
17
18>>> END OF LICENSE >>>
19*************************************************************************/
20
21using System;
22
23namespace alglib
24{
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.
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)
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;
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            }
230        }
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)
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;
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;
316        }
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)
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;
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);
529        }
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)
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            }
784            for(i=mi+1; i<=ni; i++)
785            {
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;
907        }
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)
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;
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;
1256        }
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)
1283        {
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;
1303        }
1304    }
1305}
Note: See TracBrowser for help on using the repository browser.