Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.3.0/ALGLIB-2.3.0/polint.cs @ 2806

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

Added plugin for new version of ALGLIB. #875 (Update ALGLIB sources)

File size: 40.8 KB
Line 
1/*************************************************************************
2Copyright (c) 2006-2009, 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 polint
26    {
27        /*************************************************************************
28        Polynomial fitting report:
29            TaskRCond       reciprocal of task's condition number
30            RMSError        RMS error
31            AvgError        average error
32            AvgRelError     average relative error (for non-zero Y[I])
33            MaxError        maximum error
34        *************************************************************************/
35        public struct polynomialfitreport
36        {
37            public double taskrcond;
38            public double rmserror;
39            public double avgerror;
40            public double avgrelerror;
41            public double maxerror;
42        };
43
44
45
46
47        /*************************************************************************
48        Lagrange intepolant: generation of the model on the general grid.
49        This function has O(N^2) complexity.
50
51        INPUT PARAMETERS:
52            X   -   abscissas, array[0..N-1]
53            Y   -   function values, array[0..N-1]
54            N   -   number of points, N>=1
55
56        OIYTPUT PARAMETERS
57            P   -   barycentric model which represents Lagrange interpolant
58                    (see ratint unit info and BarycentricCalc() description for
59                    more information).
60
61          -- ALGLIB --
62             Copyright 02.12.2009 by Bochkanov Sergey
63        *************************************************************************/
64        public static void polynomialbuild(ref double[] x,
65            ref double[] y,
66            int n,
67            ref ratint.barycentricinterpolant p)
68        {
69            int j = 0;
70            int k = 0;
71            double[] w = new double[0];
72            double b = 0;
73            double a = 0;
74            double v = 0;
75            double mx = 0;
76            int i_ = 0;
77
78            System.Diagnostics.Debug.Assert(n>0, "PolIntBuild: N<=0!");
79           
80            //
81            // calculate W[j]
82            // multi-pass algorithm is used to avoid overflow
83            //
84            w = new double[n];
85            a = x[0];
86            b = x[0];
87            for(j=0; j<=n-1; j++)
88            {
89                w[j] = 1;
90                a = Math.Min(a, x[j]);
91                b = Math.Max(b, x[j]);
92            }
93            for(k=0; k<=n-1; k++)
94            {
95               
96                //
97                // W[K] is used instead of 0.0 because
98                // cycle on J does not touch K-th element
99                // and we MUST get maximum from ALL elements
100                //
101                mx = Math.Abs(w[k]);
102                for(j=0; j<=n-1; j++)
103                {
104                    if( j!=k )
105                    {
106                        v = (b-a)/(x[j]-x[k]);
107                        w[j] = w[j]*v;
108                        mx = Math.Max(mx, Math.Abs(w[j]));
109                    }
110                }
111                if( k%5==0 )
112                {
113                   
114                    //
115                    // every 5-th run we renormalize W[]
116                    //
117                    v = 1/mx;
118                    for(i_=0; i_<=n-1;i_++)
119                    {
120                        w[i_] = v*w[i_];
121                    }
122                }
123            }
124            ratint.barycentricbuildxyw(ref x, ref y, ref w, n, ref p);
125        }
126
127
128        /*************************************************************************
129        Lagrange intepolant: generation of the model on equidistant grid.
130        This function has O(N) complexity.
131
132        INPUT PARAMETERS:
133            A   -   left boundary of [A,B]
134            B   -   right boundary of [A,B]
135            Y   -   function values at the nodes, array[0..N-1]
136            N   -   number of points, N>=1
137                    for N=1 a constant model is constructed.
138
139        OIYTPUT PARAMETERS
140            P   -   barycentric model which represents Lagrange interpolant
141                    (see ratint unit info and BarycentricCalc() description for
142                    more information).
143
144          -- ALGLIB --
145             Copyright 03.12.2009 by Bochkanov Sergey
146        *************************************************************************/
147        public static void polynomialbuildeqdist(double a,
148            double b,
149            ref double[] y,
150            int n,
151            ref ratint.barycentricinterpolant p)
152        {
153            int i = 0;
154            double[] w = new double[0];
155            double[] x = new double[0];
156            double v = 0;
157
158            System.Diagnostics.Debug.Assert(n>0, "PolIntBuildEqDist: N<=0!");
159           
160            //
161            // Special case: N=1
162            //
163            if( n==1 )
164            {
165                x = new double[1];
166                w = new double[1];
167                x[0] = 0.5*(b+a);
168                w[0] = 1;
169                ratint.barycentricbuildxyw(ref x, ref y, ref w, 1, ref p);
170                return;
171            }
172           
173            //
174            // general case
175            //
176            x = new double[n];
177            w = new double[n];
178            v = 1;
179            for(i=0; i<=n-1; i++)
180            {
181                w[i] = v;
182                x[i] = a+(b-a)*i/(n-1);
183                v = -(v*(n-1-i));
184                v = v/(i+1);
185            }
186            ratint.barycentricbuildxyw(ref x, ref y, ref w, n, ref p);
187        }
188
189
190        /*************************************************************************
191        Lagrange intepolant on Chebyshev grid (first kind).
192        This function has O(N) complexity.
193
194        INPUT PARAMETERS:
195            A   -   left boundary of [A,B]
196            B   -   right boundary of [A,B]
197            Y   -   function values at the nodes, array[0..N-1],
198                    Y[I] = Y(0.5*(B+A) + 0.5*(B-A)*Cos(PI*(2*i+1)/(2*n)))
199            N   -   number of points, N>=1
200                    for N=1 a constant model is constructed.
201
202        OIYTPUT PARAMETERS
203            P   -   barycentric model which represents Lagrange interpolant
204                    (see ratint unit info and BarycentricCalc() description for
205                    more information).
206
207          -- ALGLIB --
208             Copyright 03.12.2009 by Bochkanov Sergey
209        *************************************************************************/
210        public static void polynomialbuildcheb1(double a,
211            double b,
212            ref double[] y,
213            int n,
214            ref ratint.barycentricinterpolant p)
215        {
216            int i = 0;
217            double[] w = new double[0];
218            double[] x = new double[0];
219            double v = 0;
220            double t = 0;
221
222            System.Diagnostics.Debug.Assert(n>0, "PolIntBuildCheb1: N<=0!");
223           
224            //
225            // Special case: N=1
226            //
227            if( n==1 )
228            {
229                x = new double[1];
230                w = new double[1];
231                x[0] = 0.5*(b+a);
232                w[0] = 1;
233                ratint.barycentricbuildxyw(ref x, ref y, ref w, 1, ref p);
234                return;
235            }
236           
237            //
238            // general case
239            //
240            x = new double[n];
241            w = new double[n];
242            v = 1;
243            for(i=0; i<=n-1; i++)
244            {
245                t = Math.Tan(0.5*Math.PI*(2*i+1)/(2*n));
246                w[i] = 2*v*t/(1+AP.Math.Sqr(t));
247                x[i] = 0.5*(b+a)+0.5*(b-a)*(1-AP.Math.Sqr(t))/(1+AP.Math.Sqr(t));
248                v = -v;
249            }
250            ratint.barycentricbuildxyw(ref x, ref y, ref w, n, ref p);
251        }
252
253
254        /*************************************************************************
255        Lagrange intepolant on Chebyshev grid (second kind).
256        This function has O(N) complexity.
257
258        INPUT PARAMETERS:
259            A   -   left boundary of [A,B]
260            B   -   right boundary of [A,B]
261            Y   -   function values at the nodes, array[0..N-1],
262                    Y[I] = Y(0.5*(B+A) + 0.5*(B-A)*Cos(PI*i/(n-1)))
263            N   -   number of points, N>=1
264                    for N=1 a constant model is constructed.
265
266        OIYTPUT PARAMETERS
267            P   -   barycentric model which represents Lagrange interpolant
268                    (see ratint unit info and BarycentricCalc() description for
269                    more information).
270
271          -- ALGLIB --
272             Copyright 03.12.2009 by Bochkanov Sergey
273        *************************************************************************/
274        public static void polynomialbuildcheb2(double a,
275            double b,
276            ref double[] y,
277            int n,
278            ref ratint.barycentricinterpolant p)
279        {
280            int i = 0;
281            double[] w = new double[0];
282            double[] x = new double[0];
283            double v = 0;
284
285            System.Diagnostics.Debug.Assert(n>0, "PolIntBuildCheb2: N<=0!");
286           
287            //
288            // Special case: N=1
289            //
290            if( n==1 )
291            {
292                x = new double[1];
293                w = new double[1];
294                x[0] = 0.5*(b+a);
295                w[0] = 1;
296                ratint.barycentricbuildxyw(ref x, ref y, ref w, 1, ref p);
297                return;
298            }
299           
300            //
301            // general case
302            //
303            x = new double[n];
304            w = new double[n];
305            v = 1;
306            for(i=0; i<=n-1; i++)
307            {
308                if( i==0 | i==n-1 )
309                {
310                    w[i] = v*0.5;
311                }
312                else
313                {
314                    w[i] = v;
315                }
316                x[i] = 0.5*(b+a)+0.5*(b-a)*Math.Cos(Math.PI*i/(n-1));
317                v = -v;
318            }
319            ratint.barycentricbuildxyw(ref x, ref y, ref w, n, ref p);
320        }
321
322
323        /*************************************************************************
324        Fast equidistant polynomial interpolation function with O(N) complexity
325
326        INPUT PARAMETERS:
327            A   -   left boundary of [A,B]
328            B   -   right boundary of [A,B]
329            F   -   function values, array[0..N-1]
330            N   -   number of points on equidistant grid, N>=1
331                    for N=1 a constant model is constructed.
332            T   -   position where P(x) is calculated
333
334        RESULT
335            value of the Lagrange interpolant at T
336           
337        IMPORTANT
338            this function provides fast interface which is not overflow-safe
339            nor it is very precise.
340            the best option is to use  PolynomialBuildEqDist()/BarycentricCalc()
341            subroutines unless you are pretty sure that your data will not result
342            in overflow.
343
344          -- ALGLIB --
345             Copyright 02.12.2009 by Bochkanov Sergey
346        *************************************************************************/
347        public static double polynomialcalceqdist(double a,
348            double b,
349            ref double[] f,
350            int n,
351            double t)
352        {
353            double result = 0;
354            double s1 = 0;
355            double s2 = 0;
356            double v = 0;
357            double threshold = 0;
358            double s = 0;
359            double h = 0;
360            int i = 0;
361            int j = 0;
362            double w = 0;
363            double x = 0;
364
365            System.Diagnostics.Debug.Assert(n>0, "PolIntEqDist: N<=0!");
366            threshold = Math.Sqrt(AP.Math.MinRealNumber);
367           
368            //
369            // Special case: N=1
370            //
371            if( n==1 )
372            {
373                result = f[0];
374                return result;
375            }
376           
377            //
378            // First, decide: should we use "safe" formula (guarded
379            // against overflow) or fast one?
380            //
381            j = 0;
382            s = t-a;
383            for(i=1; i<=n-1; i++)
384            {
385                x = a+(double)(i)/((double)(n-1))*(b-a);
386                if( (double)(Math.Abs(t-x))<(double)(Math.Abs(s)) )
387                {
388                    s = t-x;
389                    j = i;
390                }
391            }
392            if( (double)(s)==(double)(0) )
393            {
394                result = f[j];
395                return result;
396            }
397            if( (double)(Math.Abs(s))>(double)(threshold) )
398            {
399               
400                //
401                // use fast formula
402                //
403                j = -1;
404                s = 1.0;
405            }
406           
407            //
408            // Calculate using safe or fast barycentric formula
409            //
410            s1 = 0;
411            s2 = 0;
412            w = 1.0;
413            h = (b-a)/(n-1);
414            for(i=0; i<=n-1; i++)
415            {
416                if( i!=j )
417                {
418                    v = s*w/(t-(a+i*h));
419                    s1 = s1+v*f[i];
420                    s2 = s2+v;
421                }
422                else
423                {
424                    v = w;
425                    s1 = s1+v*f[i];
426                    s2 = s2+v;
427                }
428                w = -(w*(n-1-i));
429                w = w/(i+1);
430            }
431            result = s1/s2;
432            return result;
433        }
434
435
436        /*************************************************************************
437        Fast polynomial interpolation function on Chebyshev points (first kind)
438        with O(N) complexity.
439
440        INPUT PARAMETERS:
441            A   -   left boundary of [A,B]
442            B   -   right boundary of [A,B]
443            F   -   function values, array[0..N-1]
444            N   -   number of points on Chebyshev grid (first kind),
445                    X[i] = 0.5*(B+A) + 0.5*(B-A)*Cos(PI*(2*i+1)/(2*n))
446                    for N=1 a constant model is constructed.
447            T   -   position where P(x) is calculated
448
449        RESULT
450            value of the Lagrange interpolant at T
451
452        IMPORTANT
453            this function provides fast interface which is not overflow-safe
454            nor it is very precise.
455            the best option is to use  PolIntBuildCheb1()/BarycentricCalc()
456            subroutines unless you are pretty sure that your data will not result
457            in overflow.
458
459          -- ALGLIB --
460             Copyright 02.12.2009 by Bochkanov Sergey
461        *************************************************************************/
462        public static double polynomialcalccheb1(double a,
463            double b,
464            ref double[] f,
465            int n,
466            double t)
467        {
468            double result = 0;
469            double s1 = 0;
470            double s2 = 0;
471            double v = 0;
472            double threshold = 0;
473            double s = 0;
474            int i = 0;
475            int j = 0;
476            double a0 = 0;
477            double delta = 0;
478            double alpha = 0;
479            double beta = 0;
480            double ca = 0;
481            double sa = 0;
482            double tempc = 0;
483            double temps = 0;
484            double x = 0;
485            double w = 0;
486            double p1 = 0;
487
488            System.Diagnostics.Debug.Assert(n>0, "PolIntCheb1: N<=0!");
489            threshold = Math.Sqrt(AP.Math.MinRealNumber);
490            t = (t-0.5*(a+b))/(0.5*(b-a));
491           
492            //
493            // Fast exit
494            //
495            if( n==1 )
496            {
497                result = f[0];
498                return result;
499            }
500           
501            //
502            // Prepare information for the recurrence formula
503            // used to calculate sin(pi*(2j+1)/(2n+2)) and
504            // cos(pi*(2j+1)/(2n+2)):
505            //
506            // A0    = pi/(2n+2)
507            // Delta = pi/(n+1)
508            // Alpha = 2 sin^2 (Delta/2)
509            // Beta  = sin(Delta)
510            //
511            // so that sin(..) = sin(A0+j*delta) and cos(..) = cos(A0+j*delta).
512            // Then we use
513            //
514            // sin(x+delta) = sin(x) - (alpha*sin(x) - beta*cos(x))
515            // cos(x+delta) = cos(x) - (alpha*cos(x) - beta*sin(x))
516            //
517            // to repeatedly calculate sin(..) and cos(..).
518            //
519            a0 = Math.PI/(2*(n-1)+2);
520            delta = 2*Math.PI/(2*(n-1)+2);
521            alpha = 2*AP.Math.Sqr(Math.Sin(delta/2));
522            beta = Math.Sin(delta);
523           
524            //
525            // First, decide: should we use "safe" formula (guarded
526            // against overflow) or fast one?
527            //
528            ca = Math.Cos(a0);
529            sa = Math.Sin(a0);
530            j = 0;
531            x = ca;
532            s = t-x;
533            for(i=1; i<=n-1; i++)
534            {
535               
536                //
537                // Next X[i]
538                //
539                temps = sa-(alpha*sa-beta*ca);
540                tempc = ca-(alpha*ca+beta*sa);
541                sa = temps;
542                ca = tempc;
543                x = ca;
544               
545                //
546                // Use X[i]
547                //
548                if( (double)(Math.Abs(t-x))<(double)(Math.Abs(s)) )
549                {
550                    s = t-x;
551                    j = i;
552                }
553            }
554            if( (double)(s)==(double)(0) )
555            {
556                result = f[j];
557                return result;
558            }
559            if( (double)(Math.Abs(s))>(double)(threshold) )
560            {
561               
562                //
563                // use fast formula
564                //
565                j = -1;
566                s = 1.0;
567            }
568           
569            //
570            // Calculate using safe or fast barycentric formula
571            //
572            s1 = 0;
573            s2 = 0;
574            ca = Math.Cos(a0);
575            sa = Math.Sin(a0);
576            p1 = 1.0;
577            for(i=0; i<=n-1; i++)
578            {
579               
580                //
581                // Calculate X[i], W[i]
582                //
583                x = ca;
584                w = p1*sa;
585               
586                //
587                // Proceed
588                //
589                if( i!=j )
590                {
591                    v = s*w/(t-x);
592                    s1 = s1+v*f[i];
593                    s2 = s2+v;
594                }
595                else
596                {
597                    v = w;
598                    s1 = s1+v*f[i];
599                    s2 = s2+v;
600                }
601               
602                //
603                // Next CA, SA, P1
604                //
605                temps = sa-(alpha*sa-beta*ca);
606                tempc = ca-(alpha*ca+beta*sa);
607                sa = temps;
608                ca = tempc;
609                p1 = -p1;
610            }
611            result = s1/s2;
612            return result;
613        }
614
615
616        /*************************************************************************
617        Fast polynomial interpolation function on Chebyshev points (second kind)
618        with O(N) complexity.
619
620        INPUT PARAMETERS:
621            A   -   left boundary of [A,B]
622            B   -   right boundary of [A,B]
623            F   -   function values, array[0..N-1]
624            N   -   number of points on Chebyshev grid (second kind),
625                    X[i] = 0.5*(B+A) + 0.5*(B-A)*Cos(PI*i/(n-1))
626                    for N=1 a constant model is constructed.
627            T   -   position where P(x) is calculated
628
629        RESULT
630            value of the Lagrange interpolant at T
631
632        IMPORTANT
633            this function provides fast interface which is not overflow-safe
634            nor it is very precise.
635            the best option is to use PolIntBuildCheb2()/BarycentricCalc()
636            subroutines unless you are pretty sure that your data will not result
637            in overflow.
638
639          -- ALGLIB --
640             Copyright 02.12.2009 by Bochkanov Sergey
641        *************************************************************************/
642        public static double polynomialcalccheb2(double a,
643            double b,
644            ref double[] f,
645            int n,
646            double t)
647        {
648            double result = 0;
649            double s1 = 0;
650            double s2 = 0;
651            double v = 0;
652            double threshold = 0;
653            double s = 0;
654            int i = 0;
655            int j = 0;
656            double a0 = 0;
657            double delta = 0;
658            double alpha = 0;
659            double beta = 0;
660            double ca = 0;
661            double sa = 0;
662            double tempc = 0;
663            double temps = 0;
664            double x = 0;
665            double w = 0;
666            double p1 = 0;
667
668            System.Diagnostics.Debug.Assert(n>0, "PolIntCheb2: N<=0!");
669            threshold = Math.Sqrt(AP.Math.MinRealNumber);
670            t = (t-0.5*(a+b))/(0.5*(b-a));
671           
672            //
673            // Fast exit
674            //
675            if( n==1 )
676            {
677                result = f[0];
678                return result;
679            }
680           
681            //
682            // Prepare information for the recurrence formula
683            // used to calculate sin(pi*i/n) and
684            // cos(pi*i/n):
685            //
686            // A0    = 0
687            // Delta = pi/n
688            // Alpha = 2 sin^2 (Delta/2)
689            // Beta  = sin(Delta)
690            //
691            // so that sin(..) = sin(A0+j*delta) and cos(..) = cos(A0+j*delta).
692            // Then we use
693            //
694            // sin(x+delta) = sin(x) - (alpha*sin(x) - beta*cos(x))
695            // cos(x+delta) = cos(x) - (alpha*cos(x) - beta*sin(x))
696            //
697            // to repeatedly calculate sin(..) and cos(..).
698            //
699            a0 = 0.0;
700            delta = Math.PI/(n-1);
701            alpha = 2*AP.Math.Sqr(Math.Sin(delta/2));
702            beta = Math.Sin(delta);
703           
704            //
705            // First, decide: should we use "safe" formula (guarded
706            // against overflow) or fast one?
707            //
708            ca = Math.Cos(a0);
709            sa = Math.Sin(a0);
710            j = 0;
711            x = ca;
712            s = t-x;
713            for(i=1; i<=n-1; i++)
714            {
715               
716                //
717                // Next X[i]
718                //
719                temps = sa-(alpha*sa-beta*ca);
720                tempc = ca-(alpha*ca+beta*sa);
721                sa = temps;
722                ca = tempc;
723                x = ca;
724               
725                //
726                // Use X[i]
727                //
728                if( (double)(Math.Abs(t-x))<(double)(Math.Abs(s)) )
729                {
730                    s = t-x;
731                    j = i;
732                }
733            }
734            if( (double)(s)==(double)(0) )
735            {
736                result = f[j];
737                return result;
738            }
739            if( (double)(Math.Abs(s))>(double)(threshold) )
740            {
741               
742                //
743                // use fast formula
744                //
745                j = -1;
746                s = 1.0;
747            }
748           
749            //
750            // Calculate using safe or fast barycentric formula
751            //
752            s1 = 0;
753            s2 = 0;
754            ca = Math.Cos(a0);
755            sa = Math.Sin(a0);
756            p1 = 1.0;
757            for(i=0; i<=n-1; i++)
758            {
759               
760                //
761                // Calculate X[i], W[i]
762                //
763                x = ca;
764                if( i==0 | i==n-1 )
765                {
766                    w = 0.5*p1;
767                }
768                else
769                {
770                    w = 1.0*p1;
771                }
772               
773                //
774                // Proceed
775                //
776                if( i!=j )
777                {
778                    v = s*w/(t-x);
779                    s1 = s1+v*f[i];
780                    s2 = s2+v;
781                }
782                else
783                {
784                    v = w;
785                    s1 = s1+v*f[i];
786                    s2 = s2+v;
787                }
788               
789                //
790                // Next CA, SA, P1
791                //
792                temps = sa-(alpha*sa-beta*ca);
793                tempc = ca-(alpha*ca+beta*sa);
794                sa = temps;
795                ca = tempc;
796                p1 = -p1;
797            }
798            result = s1/s2;
799            return result;
800        }
801
802
803        /*************************************************************************
804        Least squares fitting by polynomial.
805
806        This subroutine is "lightweight" alternative for more complex and feature-
807        rich PolynomialFitWC().  See  PolynomialFitWC() for more information about
808        subroutine parameters (we don't duplicate it here because of length)
809
810          -- ALGLIB PROJECT --
811             Copyright 12.10.2009 by Bochkanov Sergey
812        *************************************************************************/
813        public static void polynomialfit(ref double[] x,
814            ref double[] y,
815            int n,
816            int m,
817            ref int info,
818            ref ratint.barycentricinterpolant p,
819            ref polynomialfitreport rep)
820        {
821            int i = 0;
822            double[] w = new double[0];
823            double[] xc = new double[0];
824            double[] yc = new double[0];
825            int[] dc = new int[0];
826
827            if( n>0 )
828            {
829                w = new double[n];
830                for(i=0; i<=n-1; i++)
831                {
832                    w[i] = 1;
833                }
834            }
835            polynomialfitwc(x, y, ref w, n, xc, yc, ref dc, 0, m, ref info, ref p, ref rep);
836        }
837
838
839        /*************************************************************************
840        Weighted  fitting  by  Chebyshev  polynomial  in  barycentric  form,  with
841        constraints on function values or first derivatives.
842
843        Small regularizing term is used when solving constrained tasks (to improve
844        stability).
845
846        Task is linear, so linear least squares solver is used. Complexity of this
847        computational scheme is O(N*M^2), mostly dominated by least squares solver
848
849        SEE ALSO:
850            PolynomialFit()
851
852        INPUT PARAMETERS:
853            X   -   points, array[0..N-1].
854            Y   -   function values, array[0..N-1].
855            W   -   weights, array[0..N-1]
856                    Each summand in square  sum  of  approximation deviations from
857                    given  values  is  multiplied  by  the square of corresponding
858                    weight. Fill it by 1's if you don't  want  to  solve  weighted
859                    task.
860            N   -   number of points, N>0.
861            XC  -   points where polynomial values/derivatives are constrained,
862                    array[0..K-1].
863            YC  -   values of constraints, array[0..K-1]
864            DC  -   array[0..K-1], types of constraints:
865                    * DC[i]=0   means that P(XC[i])=YC[i]
866                    * DC[i]=1   means that P'(XC[i])=YC[i]
867                    SEE BELOW FOR IMPORTANT INFORMATION ON CONSTRAINTS
868            K   -   number of constraints, 0<=K<M.
869                    K=0 means no constraints (XC/YC/DC are not used in such cases)
870            M   -   number of basis functions (= polynomial_degree + 1), M>=1
871
872        OUTPUT PARAMETERS:
873            Info-   same format as in LSFitLinearW() subroutine:
874                    * Info>0    task is solved
875                    * Info<=0   an error occured:
876                                -4 means inconvergence of internal SVD
877                                -3 means inconsistent constraints
878                                -1 means another errors in parameters passed
879                                   (N<=0, for example)
880            P   -   interpolant in barycentric form.
881            Rep -   report, same format as in LSFitLinearW() subroutine.
882                    Following fields are set:
883                    * RMSError      rms error on the (X,Y).
884                    * AvgError      average error on the (X,Y).
885                    * AvgRelError   average relative error on the non-zero Y
886                    * MaxError      maximum error
887                                    NON-WEIGHTED ERRORS ARE CALCULATED
888
889        IMPORTANT:
890            this subroitine doesn't calculate task's condition number for K<>0.
891
892        SETTING CONSTRAINTS - DANGERS AND OPPORTUNITIES:
893
894        Setting constraints can lead  to undesired  results,  like ill-conditioned
895        behavior, or inconsistency being detected. From the other side,  it allows
896        us to improve quality of the fit. Here we summarize  our  experience  with
897        constrained regression splines:
898        * even simple constraints can be inconsistent, see  Wikipedia  article  on
899          this subject: http://en.wikipedia.org/wiki/Birkhoff_interpolation
900        * the  greater  is  M (given  fixed  constraints),  the  more chances that
901          constraints will be consistent
902        * in the general case, consistency of constraints is NOT GUARANTEED.
903        * in the one special cases, however, we can  guarantee  consistency.  This
904          case  is:  M>1  and constraints on the function values (NOT DERIVATIVES)
905
906        Our final recommendation is to use constraints  WHEN  AND  ONLY  when  you
907        can't solve your task without them. Anything beyond  special  cases  given
908        above is not guaranteed and may result in inconsistency.
909
910          -- ALGLIB PROJECT --
911             Copyright 10.12.2009 by Bochkanov Sergey
912        *************************************************************************/
913        public static void polynomialfitwc(double[] x,
914            double[] y,
915            ref double[] w,
916            int n,
917            double[] xc,
918            double[] yc,
919            ref int[] dc,
920            int k,
921            int m,
922            ref int info,
923            ref ratint.barycentricinterpolant p,
924            ref polynomialfitreport rep)
925        {
926            double xa = 0;
927            double xb = 0;
928            double sa = 0;
929            double sb = 0;
930            double[] xoriginal = new double[0];
931            double[] yoriginal = new double[0];
932            double[] y2 = new double[0];
933            double[] w2 = new double[0];
934            double[] tmp = new double[0];
935            double[] tmp2 = new double[0];
936            double[] tmpdiff = new double[0];
937            double[] bx = new double[0];
938            double[] by = new double[0];
939            double[] bw = new double[0];
940            double[,] fmatrix = new double[0,0];
941            double[,] cmatrix = new double[0,0];
942            int i = 0;
943            int j = 0;
944            double mx = 0;
945            double decay = 0;
946            double u = 0;
947            double v = 0;
948            double s = 0;
949            int relcnt = 0;
950            lsfit.lsfitreport lrep = new lsfit.lsfitreport();
951            int i_ = 0;
952
953            x = (double[])x.Clone();
954            y = (double[])y.Clone();
955            xc = (double[])xc.Clone();
956            yc = (double[])yc.Clone();
957
958            if( m<1 | n<1 | k<0 | k>=m )
959            {
960                info = -1;
961                return;
962            }
963            for(i=0; i<=k-1; i++)
964            {
965                info = 0;
966                if( dc[i]<0 )
967                {
968                    info = -1;
969                }
970                if( dc[i]>1 )
971                {
972                    info = -1;
973                }
974                if( info<0 )
975                {
976                    return;
977                }
978            }
979           
980            //
981            // weight decay for correct handling of task which becomes
982            // degenerate after constraints are applied
983            //
984            decay = 10000*AP.Math.MachineEpsilon;
985           
986            //
987            // Scale X, Y, XC, YC
988            //
989            lsfit.lsfitscalexy(ref x, ref y, n, ref xc, ref yc, ref dc, k, ref xa, ref xb, ref sa, ref sb, ref xoriginal, ref yoriginal);
990           
991            //
992            // allocate space, initialize/fill:
993            // * FMatrix-   values of basis functions at X[]
994            // * CMatrix-   values (derivatives) of basis functions at XC[]
995            // * fill constraints matrix
996            // * fill first N rows of design matrix with values
997            // * fill next M rows of design matrix with regularizing term
998            // * append M zeros to Y
999            // * append M elements, mean(abs(W)) each, to W
1000            //
1001            y2 = new double[n+m];
1002            w2 = new double[n+m];
1003            tmp = new double[m];
1004            tmpdiff = new double[m];
1005            fmatrix = new double[n+m, m];
1006            if( k>0 )
1007            {
1008                cmatrix = new double[k, m+1];
1009            }
1010           
1011            //
1012            // Fill design matrix, Y2, W2:
1013            // * first N rows with basis functions for original points
1014            // * next M rows with decay terms
1015            //
1016            for(i=0; i<=n-1; i++)
1017            {
1018               
1019                //
1020                // prepare Ith row
1021                // use Tmp for calculations to avoid multidimensional arrays overhead
1022                //
1023                for(j=0; j<=m-1; j++)
1024                {
1025                    if( j==0 )
1026                    {
1027                        tmp[j] = 1;
1028                    }
1029                    else
1030                    {
1031                        if( j==1 )
1032                        {
1033                            tmp[j] = x[i];
1034                        }
1035                        else
1036                        {
1037                            tmp[j] = 2*x[i]*tmp[j-1]-tmp[j-2];
1038                        }
1039                    }
1040                }
1041                for(i_=0; i_<=m-1;i_++)
1042                {
1043                    fmatrix[i,i_] = tmp[i_];
1044                }
1045            }
1046            for(i=0; i<=m-1; i++)
1047            {
1048                for(j=0; j<=m-1; j++)
1049                {
1050                    if( i==j )
1051                    {
1052                        fmatrix[n+i,j] = decay;
1053                    }
1054                    else
1055                    {
1056                        fmatrix[n+i,j] = 0;
1057                    }
1058                }
1059            }
1060            for(i_=0; i_<=n-1;i_++)
1061            {
1062                y2[i_] = y[i_];
1063            }
1064            for(i_=0; i_<=n-1;i_++)
1065            {
1066                w2[i_] = w[i_];
1067            }
1068            mx = 0;
1069            for(i=0; i<=n-1; i++)
1070            {
1071                mx = mx+Math.Abs(w[i]);
1072            }
1073            mx = mx/n;
1074            for(i=0; i<=m-1; i++)
1075            {
1076                y2[n+i] = 0;
1077                w2[n+i] = mx;
1078            }
1079           
1080            //
1081            // fill constraints matrix
1082            //
1083            for(i=0; i<=k-1; i++)
1084            {
1085               
1086                //
1087                // prepare Ith row
1088                // use Tmp for basis function values,
1089                // TmpDiff for basos function derivatives
1090                //
1091                for(j=0; j<=m-1; j++)
1092                {
1093                    if( j==0 )
1094                    {
1095                        tmp[j] = 1;
1096                        tmpdiff[j] = 0;
1097                    }
1098                    else
1099                    {
1100                        if( j==1 )
1101                        {
1102                            tmp[j] = xc[i];
1103                            tmpdiff[j] = 1;
1104                        }
1105                        else
1106                        {
1107                            tmp[j] = 2*xc[i]*tmp[j-1]-tmp[j-2];
1108                            tmpdiff[j] = 2*(tmp[j-1]+xc[i]*tmpdiff[j-1])-tmpdiff[j-2];
1109                        }
1110                    }
1111                }
1112                if( dc[i]==0 )
1113                {
1114                    for(i_=0; i_<=m-1;i_++)
1115                    {
1116                        cmatrix[i,i_] = tmp[i_];
1117                    }
1118                }
1119                if( dc[i]==1 )
1120                {
1121                    for(i_=0; i_<=m-1;i_++)
1122                    {
1123                        cmatrix[i,i_] = tmpdiff[i_];
1124                    }
1125                }
1126                cmatrix[i,m] = yc[i];
1127            }
1128           
1129            //
1130            // Solve constrained task
1131            //
1132            if( k>0 )
1133            {
1134               
1135                //
1136                // solve using regularization
1137                //
1138                lsfit.lsfitlinearwc(y2, ref w2, ref fmatrix, cmatrix, n+m, m, k, ref info, ref tmp, ref lrep);
1139            }
1140            else
1141            {
1142               
1143                //
1144                // no constraints, no regularization needed
1145                //
1146                lsfit.lsfitlinearwc(y, ref w, ref fmatrix, cmatrix, n, m, 0, ref info, ref tmp, ref lrep);
1147            }
1148            if( info<0 )
1149            {
1150                return;
1151            }
1152           
1153            //
1154            // Generate barycentric model and scale it
1155            // * BX, BY store barycentric model nodes
1156            // * FMatrix is reused (remember - it is at least MxM, what we need)
1157            //
1158            // Model intialization is done in O(M^2). In principle, it can be
1159            // done in O(M*log(M)), but before it we solved task with O(N*M^2)
1160            // complexity, so it is only a small amount of total time spent.
1161            //
1162            bx = new double[m];
1163            by = new double[m];
1164            bw = new double[m];
1165            tmp2 = new double[m];
1166            s = 1;
1167            for(i=0; i<=m-1; i++)
1168            {
1169                if( m!=1 )
1170                {
1171                    u = Math.Cos(Math.PI*i/(m-1));
1172                }
1173                else
1174                {
1175                    u = 0;
1176                }
1177                v = 0;
1178                for(j=0; j<=m-1; j++)
1179                {
1180                    if( j==0 )
1181                    {
1182                        tmp2[j] = 1;
1183                    }
1184                    else
1185                    {
1186                        if( j==1 )
1187                        {
1188                            tmp2[j] = u;
1189                        }
1190                        else
1191                        {
1192                            tmp2[j] = 2*u*tmp2[j-1]-tmp2[j-2];
1193                        }
1194                    }
1195                    v = v+tmp[j]*tmp2[j];
1196                }
1197                bx[i] = u;
1198                by[i] = v;
1199                bw[i] = s;
1200                if( i==0 | i==m-1 )
1201                {
1202                    bw[i] = 0.5*bw[i];
1203                }
1204                s = -s;
1205            }
1206            ratint.barycentricbuildxyw(ref bx, ref by, ref bw, m, ref p);
1207            ratint.barycentriclintransx(ref p, 2/(xb-xa), -((xa+xb)/(xb-xa)));
1208            ratint.barycentriclintransy(ref p, sb-sa, sa);
1209           
1210            //
1211            // Scale absolute errors obtained from LSFitLinearW.
1212            // Relative error should be calculated separately
1213            // (because of shifting/scaling of the task)
1214            //
1215            rep.taskrcond = lrep.taskrcond;
1216            rep.rmserror = lrep.rmserror*(sb-sa);
1217            rep.avgerror = lrep.avgerror*(sb-sa);
1218            rep.maxerror = lrep.maxerror*(sb-sa);
1219            rep.avgrelerror = 0;
1220            relcnt = 0;
1221            for(i=0; i<=n-1; i++)
1222            {
1223                if( (double)(yoriginal[i])!=(double)(0) )
1224                {
1225                    rep.avgrelerror = rep.avgrelerror+Math.Abs(ratint.barycentriccalc(ref p, xoriginal[i])-yoriginal[i])/Math.Abs(yoriginal[i]);
1226                    relcnt = relcnt+1;
1227                }
1228            }
1229            if( relcnt!=0 )
1230            {
1231                rep.avgrelerror = rep.avgrelerror/relcnt;
1232            }
1233        }
1234    }
1235}
Note: See TracBrowser for help on using the repository browser.