Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.1.2.2591/ALGLIB-2.1.2.2591/ratint.cs @ 2645

Last change on this file since 2645 was 2645, checked in by mkommend, 14 years ago

extracted external libraries and adapted dependent plugins (ticket #837)

File size: 46.9 KB
Line 
1/*************************************************************************
2Copyright (c) 2007-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 ratint
26    {
27        /*************************************************************************
28        Barycentric interpolant.
29        *************************************************************************/
30        public struct barycentricinterpolant
31        {
32            public int n;
33            public double sy;
34            public double[] x;
35            public double[] y;
36            public double[] w;
37        };
38
39
40        /*************************************************************************
41        Barycentric fitting report:
42            TaskRCond       reciprocal of task's condition number
43            RMSError        RMS error
44            AvgError        average error
45            AvgRelError     average relative error (for non-zero Y[I])
46            MaxError        maximum error
47        *************************************************************************/
48        public struct barycentricfitreport
49        {
50            public double taskrcond;
51            public double dbest;
52            public double rmserror;
53            public double avgerror;
54            public double avgrelerror;
55            public double maxerror;
56        };
57
58
59
60
61        public const int brcvnum = 10;
62
63
64        /*************************************************************************
65        Rational interpolation using barycentric formula
66
67        F(t) = SUM(i=0,n-1,w[i]*f[i]/(t-x[i])) / SUM(i=0,n-1,w[i]/(t-x[i]))
68
69        Input parameters:
70            B   -   barycentric interpolant built with one of model building
71                    subroutines.
72            T   -   interpolation point
73
74        Result:
75            barycentric interpolant F(t)
76
77          -- ALGLIB --
78             Copyright 17.08.2009 by Bochkanov Sergey
79        *************************************************************************/
80        public static double barycentriccalc(ref barycentricinterpolant b,
81            double t)
82        {
83            double result = 0;
84            double s1 = 0;
85            double s2 = 0;
86            double s = 0;
87            double v = 0;
88            int i = 0;
89
90           
91            //
92            // special case: N=1
93            //
94            if( b.n==1 )
95            {
96                result = b.sy*b.y[0];
97                return result;
98            }
99           
100            //
101            // Here we assume that task is normalized, i.e.:
102            // 1. abs(Y[i])<=1
103            // 2. abs(W[i])<=1
104            // 3. X[] is ordered
105            //
106            s = Math.Abs(t-b.x[0]);
107            for(i=0; i<=b.n-1; i++)
108            {
109                v = b.x[i];
110                if( (double)(v)==(double)(t) )
111                {
112                    result = b.sy*b.y[i];
113                    return result;
114                }
115                v = Math.Abs(t-v);
116                if( (double)(v)<(double)(s) )
117                {
118                    s = v;
119                }
120            }
121            s1 = 0;
122            s2 = 0;
123            for(i=0; i<=b.n-1; i++)
124            {
125                v = s/(t-b.x[i]);
126                v = v*b.w[i];
127                s1 = s1+v*b.y[i];
128                s2 = s2+v;
129            }
130            result = b.sy*s1/s2;
131            return result;
132        }
133
134
135        /*************************************************************************
136        Differentiation of barycentric interpolant: first derivative.
137
138        Algorithm used in this subroutine is very robust and should not fail until
139        provided with values too close to MaxRealNumber  (usually  MaxRealNumber/N
140        or greater will overflow).
141
142        INPUT PARAMETERS:
143            B   -   barycentric interpolant built with one of model building
144                    subroutines.
145            T   -   interpolation point
146
147        OUTPUT PARAMETERS:
148            F   -   barycentric interpolant at T
149            DF  -   first derivative
150           
151        NOTE
152
153
154          -- ALGLIB --
155             Copyright 17.08.2009 by Bochkanov Sergey
156        *************************************************************************/
157        public static void barycentricdiff1(ref barycentricinterpolant b,
158            double t,
159            ref double f,
160            ref double df)
161        {
162            double v = 0;
163            double vv = 0;
164            int i = 0;
165            int k = 0;
166            double n0 = 0;
167            double n1 = 0;
168            double d0 = 0;
169            double d1 = 0;
170            double s0 = 0;
171            double s1 = 0;
172            double xk = 0;
173            double xi = 0;
174            double xmin = 0;
175            double xmax = 0;
176            double xscale1 = 0;
177            double xoffs1 = 0;
178            double xscale2 = 0;
179            double xoffs2 = 0;
180            double xprev = 0;
181
182           
183            //
184            // special case: N=1
185            //
186            if( b.n==1 )
187            {
188                f = b.sy*b.y[0];
189                df = 0;
190                return;
191            }
192            if( (double)(b.sy)==(double)(0) )
193            {
194                f = 0;
195                df = 0;
196                return;
197            }
198            System.Diagnostics.Debug.Assert((double)(b.sy)>(double)(0), "BarycentricDiff1: internal error");
199           
200            //
201            // We assume than N>1 and B.SY>0. Find:
202            // 1. pivot point (X[i] closest to T)
203            // 2. width of interval containing X[i]
204            //
205            v = Math.Abs(b.x[0]-t);
206            k = 0;
207            xmin = b.x[0];
208            xmax = b.x[0];
209            for(i=1; i<=b.n-1; i++)
210            {
211                vv = b.x[i];
212                if( (double)(Math.Abs(vv-t))<(double)(v) )
213                {
214                    v = Math.Abs(vv-t);
215                    k = i;
216                }
217                xmin = Math.Min(xmin, vv);
218                xmax = Math.Max(xmax, vv);
219            }
220           
221            //
222            // pivot point found, calculate dNumerator and dDenominator
223            //
224            xscale1 = 1/(xmax-xmin);
225            xoffs1 = -(xmin/(xmax-xmin))+1;
226            xscale2 = 2;
227            xoffs2 = -3;
228            t = t*xscale1+xoffs1;
229            t = t*xscale2+xoffs2;
230            xk = b.x[k];
231            xk = xk*xscale1+xoffs1;
232            xk = xk*xscale2+xoffs2;
233            v = t-xk;
234            n0 = 0;
235            n1 = 0;
236            d0 = 0;
237            d1 = 0;
238            xprev = -2;
239            for(i=0; i<=b.n-1; i++)
240            {
241                xi = b.x[i];
242                xi = xi*xscale1+xoffs1;
243                xi = xi*xscale2+xoffs2;
244                System.Diagnostics.Debug.Assert((double)(xi)>(double)(xprev), "BarycentricDiff1: points are too close!");
245                xprev = xi;
246                if( i!=k )
247                {
248                    vv = AP.Math.Sqr(t-xi);
249                    s0 = (t-xk)/(t-xi);
250                    s1 = (xk-xi)/vv;
251                }
252                else
253                {
254                    s0 = 1;
255                    s1 = 0;
256                }
257                vv = b.w[i]*b.y[i];
258                n0 = n0+s0*vv;
259                n1 = n1+s1*vv;
260                vv = b.w[i];
261                d0 = d0+s0*vv;
262                d1 = d1+s1*vv;
263            }
264            f = b.sy*n0/d0;
265            df = (n1*d0-n0*d1)/AP.Math.Sqr(d0);
266            if( (double)(df)!=(double)(0) )
267            {
268                df = Math.Sign(df)*Math.Exp(Math.Log(Math.Abs(df))+Math.Log(b.sy)+Math.Log(xscale1)+Math.Log(xscale2));
269            }
270        }
271
272
273        /*************************************************************************
274        Differentiation of barycentric interpolant: first/second derivatives.
275
276        INPUT PARAMETERS:
277            B   -   barycentric interpolant built with one of model building
278                    subroutines.
279            T   -   interpolation point
280
281        OUTPUT PARAMETERS:
282            F   -   barycentric interpolant at T
283            DF  -   first derivative
284            D2F -   second derivative
285
286        NOTE: this algorithm may fail due to overflow/underflor if  used  on  data
287        whose values are close to MaxRealNumber or MinRealNumber.  Use more robust
288        BarycentricDiff1() subroutine in such cases.
289
290
291          -- ALGLIB --
292             Copyright 17.08.2009 by Bochkanov Sergey
293        *************************************************************************/
294        public static void barycentricdiff2(ref barycentricinterpolant b,
295            double t,
296            ref double f,
297            ref double df,
298            ref double d2f)
299        {
300            double v = 0;
301            double vv = 0;
302            int i = 0;
303            int k = 0;
304            double n0 = 0;
305            double n1 = 0;
306            double n2 = 0;
307            double d0 = 0;
308            double d1 = 0;
309            double d2 = 0;
310            double s0 = 0;
311            double s1 = 0;
312            double s2 = 0;
313            double xk = 0;
314            double xi = 0;
315
316            f = 0;
317            df = 0;
318            d2f = 0;
319           
320            //
321            // special case: N=1
322            //
323            if( b.n==1 )
324            {
325                f = b.sy*b.y[0];
326                df = 0;
327                d2f = 0;
328                return;
329            }
330            if( (double)(b.sy)==(double)(0) )
331            {
332                f = 0;
333                df = 0;
334                d2f = 0;
335                return;
336            }
337            System.Diagnostics.Debug.Assert((double)(b.sy)>(double)(0), "BarycentricDiff: internal error");
338           
339            //
340            // We assume than N>1 and B.SY>0. Find:
341            // 1. pivot point (X[i] closest to T)
342            // 2. width of interval containing X[i]
343            //
344            v = Math.Abs(b.x[0]-t);
345            k = 0;
346            for(i=1; i<=b.n-1; i++)
347            {
348                vv = b.x[i];
349                if( (double)(Math.Abs(vv-t))<(double)(v) )
350                {
351                    v = Math.Abs(vv-t);
352                    k = i;
353                }
354            }
355           
356            //
357            // pivot point found, calculate dNumerator and dDenominator
358            //
359            xk = b.x[k];
360            v = t-xk;
361            n0 = 0;
362            n1 = 0;
363            n2 = 0;
364            d0 = 0;
365            d1 = 0;
366            d2 = 0;
367            for(i=0; i<=b.n-1; i++)
368            {
369                if( i!=k )
370                {
371                    xi = b.x[i];
372                    vv = AP.Math.Sqr(t-xi);
373                    s0 = (t-xk)/(t-xi);
374                    s1 = (xk-xi)/vv;
375                    s2 = -(2*(xk-xi)/(vv*(t-xi)));
376                }
377                else
378                {
379                    s0 = 1;
380                    s1 = 0;
381                    s2 = 0;
382                }
383                vv = b.w[i]*b.y[i];
384                n0 = n0+s0*vv;
385                n1 = n1+s1*vv;
386                n2 = n2+s2*vv;
387                vv = b.w[i];
388                d0 = d0+s0*vv;
389                d1 = d1+s1*vv;
390                d2 = d2+s2*vv;
391            }
392            f = b.sy*n0/d0;
393            df = b.sy*(n1*d0-n0*d1)/AP.Math.Sqr(d0);
394            d2f = b.sy*((n2*d0-n0*d2)*AP.Math.Sqr(d0)-(n1*d0-n0*d1)*2*d0*d1)/AP.Math.Sqr(AP.Math.Sqr(d0));
395        }
396
397
398        /*************************************************************************
399        This subroutine performs linear transformation of the argument.
400
401        INPUT PARAMETERS:
402            B       -   rational interpolant in barycentric form
403            CA, CB  -   transformation coefficients: x = CA*t + CB
404
405        OUTPUT PARAMETERS:
406            B       -   transformed interpolant with X replaced by T
407
408          -- ALGLIB PROJECT --
409             Copyright 19.08.2009 by Bochkanov Sergey
410        *************************************************************************/
411        public static void barycentriclintransx(ref barycentricinterpolant b,
412            double ca,
413            double cb)
414        {
415            int i = 0;
416            int j = 0;
417            double v = 0;
418
419           
420            //
421            // special case, replace by constant F(CB)
422            //
423            if( (double)(ca)==(double)(0) )
424            {
425                b.sy = barycentriccalc(ref b, cb);
426                v = 1;
427                for(i=0; i<=b.n-1; i++)
428                {
429                    b.y[i] = 1;
430                    b.w[i] = v;
431                    v = -v;
432                }
433                return;
434            }
435           
436            //
437            // general case: CA<>0
438            //
439            for(i=0; i<=b.n-1; i++)
440            {
441                b.x[i] = (b.x[i]-cb)/ca;
442            }
443            if( (double)(ca)<(double)(0) )
444            {
445                for(i=0; i<=b.n-1; i++)
446                {
447                    if( i<b.n-1-i )
448                    {
449                        j = b.n-1-i;
450                        v = b.x[i];
451                        b.x[i] = b.x[j];
452                        b.x[j] = v;
453                        v = b.y[i];
454                        b.y[i] = b.y[j];
455                        b.y[j] = v;
456                        v = b.w[i];
457                        b.w[i] = b.w[j];
458                        b.w[j] = v;
459                    }
460                    else
461                    {
462                        break;
463                    }
464                }
465            }
466        }
467
468
469        /*************************************************************************
470        This  subroutine   performs   linear  transformation  of  the  barycentric
471        interpolant.
472
473        INPUT PARAMETERS:
474            B       -   rational interpolant in barycentric form
475            CA, CB  -   transformation coefficients: B2(x) = CA*B(x) + CB
476
477        OUTPUT PARAMETERS:
478            B       -   transformed interpolant
479
480          -- ALGLIB PROJECT --
481             Copyright 19.08.2009 by Bochkanov Sergey
482        *************************************************************************/
483        public static void barycentriclintransy(ref barycentricinterpolant b,
484            double ca,
485            double cb)
486        {
487            int i = 0;
488            double v = 0;
489            int i_ = 0;
490
491            for(i=0; i<=b.n-1; i++)
492            {
493                b.y[i] = ca*b.sy*b.y[i]+cb;
494            }
495            b.sy = 0;
496            for(i=0; i<=b.n-1; i++)
497            {
498                b.sy = Math.Max(b.sy, Math.Abs(b.y[i]));
499            }
500            if( (double)(b.sy)>(double)(0) )
501            {
502                v = 1/b.sy;
503                for(i_=0; i_<=b.n-1;i_++)
504                {
505                    b.y[i_] = v*b.y[i_];
506                }
507            }
508        }
509
510
511        /*************************************************************************
512        Extracts X/Y/W arrays from rational interpolant
513
514        INPUT PARAMETERS:
515            B   -   barycentric interpolant
516
517        OUTPUT PARAMETERS:
518            N   -   nodes count, N>0
519            X   -   interpolation nodes, array[0..N-1]
520            F   -   function values, array[0..N-1]
521            W   -   barycentric weights, array[0..N-1]
522
523          -- ALGLIB --
524             Copyright 17.08.2009 by Bochkanov Sergey
525        *************************************************************************/
526        public static void barycentricunpack(ref barycentricinterpolant b,
527            ref int n,
528            ref double[] x,
529            ref double[] y,
530            ref double[] w)
531        {
532            double v = 0;
533            int i_ = 0;
534
535            n = b.n;
536            x = new double[n];
537            y = new double[n];
538            w = new double[n];
539            v = b.sy;
540            for(i_=0; i_<=n-1;i_++)
541            {
542                x[i_] = b.x[i_];
543            }
544            for(i_=0; i_<=n-1;i_++)
545            {
546                y[i_] = v*b.y[i_];
547            }
548            for(i_=0; i_<=n-1;i_++)
549            {
550                w[i_] = b.w[i_];
551            }
552        }
553
554
555        /*************************************************************************
556        Serialization of the barycentric interpolant
557
558        INPUT PARAMETERS:
559            B   -   barycentric interpolant
560
561        OUTPUT PARAMETERS:
562            RA      -   array of real numbers which contains interpolant,
563                        array[0..RLen-1]
564            RLen    -   RA lenght
565
566          -- ALGLIB --
567             Copyright 17.08.2009 by Bochkanov Sergey
568        *************************************************************************/
569        public static void barycentricserialize(ref barycentricinterpolant b,
570            ref double[] ra,
571            ref int ralen)
572        {
573            int i_ = 0;
574            int i1_ = 0;
575
576            ralen = 2+2+3*b.n;
577            ra = new double[ralen];
578            ra[0] = ralen;
579            ra[1] = brcvnum;
580            ra[2] = b.n;
581            ra[3] = b.sy;
582            i1_ = (0) - (4);
583            for(i_=4; i_<=4+b.n-1;i_++)
584            {
585                ra[i_] = b.x[i_+i1_];
586            }
587            i1_ = (0) - (4+b.n);
588            for(i_=4+b.n; i_<=4+2*b.n-1;i_++)
589            {
590                ra[i_] = b.y[i_+i1_];
591            }
592            i1_ = (0) - (4+2*b.n);
593            for(i_=4+2*b.n; i_<=4+3*b.n-1;i_++)
594            {
595                ra[i_] = b.w[i_+i1_];
596            }
597        }
598
599
600        /*************************************************************************
601        Unserialization of the barycentric interpolant
602
603        INPUT PARAMETERS:
604            RA  -   array of real numbers which contains interpolant,
605
606        OUTPUT PARAMETERS:
607            B   -   barycentric interpolant
608
609          -- ALGLIB --
610             Copyright 17.08.2009 by Bochkanov Sergey
611        *************************************************************************/
612        public static void barycentricunserialize(ref double[] ra,
613            ref barycentricinterpolant b)
614        {
615            int i_ = 0;
616            int i1_ = 0;
617
618            System.Diagnostics.Debug.Assert((int)Math.Round(ra[1])==brcvnum, "BarycentricUnserialize: corrupted array!");
619            b.n = (int)Math.Round(ra[2]);
620            b.sy = ra[3];
621            b.x = new double[b.n];
622            b.y = new double[b.n];
623            b.w = new double[b.n];
624            i1_ = (4) - (0);
625            for(i_=0; i_<=b.n-1;i_++)
626            {
627                b.x[i_] = ra[i_+i1_];
628            }
629            i1_ = (4+b.n) - (0);
630            for(i_=0; i_<=b.n-1;i_++)
631            {
632                b.y[i_] = ra[i_+i1_];
633            }
634            i1_ = (4+2*b.n) - (0);
635            for(i_=0; i_<=b.n-1;i_++)
636            {
637                b.w[i_] = ra[i_+i1_];
638            }
639        }
640
641
642        /*************************************************************************
643        Copying of the barycentric interpolant
644
645        INPUT PARAMETERS:
646            B   -   barycentric interpolant
647
648        OUTPUT PARAMETERS:
649            B2  -   copy(B1)
650
651          -- ALGLIB --
652             Copyright 17.08.2009 by Bochkanov Sergey
653        *************************************************************************/
654        public static void barycentriccopy(ref barycentricinterpolant b,
655            ref barycentricinterpolant b2)
656        {
657            int i_ = 0;
658
659            b2.n = b.n;
660            b2.sy = b.sy;
661            b2.x = new double[b2.n];
662            b2.y = new double[b2.n];
663            b2.w = new double[b2.n];
664            for(i_=0; i_<=b2.n-1;i_++)
665            {
666                b2.x[i_] = b.x[i_];
667            }
668            for(i_=0; i_<=b2.n-1;i_++)
669            {
670                b2.y[i_] = b.y[i_];
671            }
672            for(i_=0; i_<=b2.n-1;i_++)
673            {
674                b2.w[i_] = b.w[i_];
675            }
676        }
677
678
679        /*************************************************************************
680        Rational interpolant from X/Y/W arrays
681
682        F(t) = SUM(i=0,n-1,w[i]*f[i]/(t-x[i])) / SUM(i=0,n-1,w[i]/(t-x[i]))
683
684        INPUT PARAMETERS:
685            X   -   interpolation nodes, array[0..N-1]
686            F   -   function values, array[0..N-1]
687            W   -   barycentric weights, array[0..N-1]
688            N   -   nodes count, N>0
689
690        OUTPUT PARAMETERS:
691            B   -   barycentric interpolant built from (X, Y, W)
692
693          -- ALGLIB --
694             Copyright 17.08.2009 by Bochkanov Sergey
695        *************************************************************************/
696        public static void barycentricbuildxyw(ref double[] x,
697            ref double[] y,
698            ref double[] w,
699            int n,
700            ref barycentricinterpolant b)
701        {
702            int i_ = 0;
703
704            System.Diagnostics.Debug.Assert(n>0, "BarycentricBuildXYW: incorrect N!");
705           
706            //
707            // fill X/Y/W
708            //
709            b.x = new double[n];
710            b.y = new double[n];
711            b.w = new double[n];
712            for(i_=0; i_<=n-1;i_++)
713            {
714                b.x[i_] = x[i_];
715            }
716            for(i_=0; i_<=n-1;i_++)
717            {
718                b.y[i_] = y[i_];
719            }
720            for(i_=0; i_<=n-1;i_++)
721            {
722                b.w[i_] = w[i_];
723            }
724            b.n = n;
725           
726            //
727            // Normalize
728            //
729            barycentricnormalize(ref b);
730        }
731
732
733        /*************************************************************************
734        Rational interpolant without poles
735
736        The subroutine constructs the rational interpolating function without real
737        poles  (see  'Barycentric rational interpolation with no  poles  and  high
738        rates of approximation', Michael S. Floater. and  Kai  Hormann,  for  more
739        information on this subject).
740
741        Input parameters:
742            X   -   interpolation nodes, array[0..N-1].
743            Y   -   function values, array[0..N-1].
744            N   -   number of nodes, N>0.
745            D   -   order of the interpolation scheme, 0 <= D <= N-1.
746                    D<0 will cause an error.
747                    D>=N it will be replaced with D=N-1.
748                    if you don't know what D to choose, use small value about 3-5.
749
750        Output parameters:
751            B   -   barycentric interpolant.
752
753        Note:
754            this algorithm always succeeds and calculates the weights  with  close
755            to machine precision.
756
757          -- ALGLIB PROJECT --
758             Copyright 17.06.2007 by Bochkanov Sergey
759        *************************************************************************/
760        public static void barycentricbuildfloaterhormann(ref double[] x,
761            ref double[] y,
762            int n,
763            int d,
764            ref barycentricinterpolant b)
765        {
766            double s0 = 0;
767            double s = 0;
768            double v = 0;
769            int i = 0;
770            int j = 0;
771            int k = 0;
772            int[] perm = new int[0];
773            double[] wtemp = new double[0];
774            int i_ = 0;
775
776            System.Diagnostics.Debug.Assert(n>0, "BarycentricFloaterHormann: N<=0!");
777            System.Diagnostics.Debug.Assert(d>=0, "BarycentricFloaterHormann: incorrect D!");
778           
779            //
780            // Prepare
781            //
782            if( d>n-1 )
783            {
784                d = n-1;
785            }
786            b.n = n;
787           
788            //
789            // special case: N=1
790            //
791            if( n==1 )
792            {
793                b.x = new double[n];
794                b.y = new double[n];
795                b.w = new double[n];
796                b.x[0] = x[0];
797                b.y[0] = y[0];
798                b.w[0] = 1;
799                barycentricnormalize(ref b);
800                return;
801            }
802           
803            //
804            // Fill X/Y
805            //
806            b.x = new double[n];
807            b.y = new double[n];
808            for(i_=0; i_<=n-1;i_++)
809            {
810                b.x[i_] = x[i_];
811            }
812            for(i_=0; i_<=n-1;i_++)
813            {
814                b.y[i_] = y[i_];
815            }
816            tsort.tagsortfastr(ref b.x, ref b.y, n);
817           
818            //
819            // Calculate Wk
820            //
821            b.w = new double[n];
822            s0 = 1;
823            for(k=1; k<=d; k++)
824            {
825                s0 = -s0;
826            }
827            for(k=0; k<=n-1; k++)
828            {
829               
830                //
831                // Wk
832                //
833                s = 0;
834                for(i=Math.Max(k-d, 0); i<=Math.Min(k, n-1-d); i++)
835                {
836                    v = 1;
837                    for(j=i; j<=i+d; j++)
838                    {
839                        if( j!=k )
840                        {
841                            v = v/Math.Abs(b.x[k]-b.x[j]);
842                        }
843                    }
844                    s = s+v;
845                }
846                b.w[k] = s0*s;
847               
848                //
849                // Next S0
850                //
851                s0 = -s0;
852            }
853           
854            //
855            // Normalize
856            //
857            barycentricnormalize(ref b);
858        }
859
860
861        /*************************************************************************
862        Weghted rational least  squares  fitting  using  Floater-Hormann  rational
863        functions  with  optimal  D  chosen  from  [0,9],  with  constraints   and
864        individual weights.
865
866        Equidistant grid with M+1 node on [min(x),max(x)]  is  used to build basis
867        functions. Different values of D are tried, optimal D (least WEIGHTED root
868        mean square error) is chosen.  Task  is  linear,  so  linear least squares
869        solver  is  used.  Complexity  of  this  computational  scheme is O(N*M^2)
870        (mostly dominated by the least squares solver).
871
872        SEE ALSO
873        * BarycentricFitFloaterHormann(), "lightweight" fitting without invididual
874          weights and constraints.
875
876        INPUT PARAMETERS:
877            X   -   points, array[0..N-1].
878            Y   -   function values, array[0..N-1].
879            W   -   weights, array[0..N-1]
880                    Each summand in square  sum  of  approximation deviations from
881                    given  values  is  multiplied  by  the square of corresponding
882                    weight. Fill it by 1's if you don't  want  to  solve  weighted
883                    task.
884            N   -   number of points, N>0.
885            XC  -   points where function values/derivatives are constrained,
886                    array[0..K-1].
887            YC  -   values of constraints, array[0..K-1]
888            DC  -   array[0..K-1], types of constraints:
889                    * DC[i]=0   means that S(XC[i])=YC[i]
890                    * DC[i]=1   means that S'(XC[i])=YC[i]
891                    SEE BELOW FOR IMPORTANT INFORMATION ON CONSTRAINTS
892            K   -   number of constraints, 0<=K<M.
893                    K=0 means no constraints (XC/YC/DC are not used in such cases)
894            M   -   number of basis functions ( = number_of_nodes), M>=2.
895
896        OUTPUT PARAMETERS:
897            Info-   same format as in LSFitLinearWC() subroutine.
898                    * Info>0    task is solved
899                    * Info<=0   an error occured:
900                                -4 means inconvergence of internal SVD
901                                -3 means inconsistent constraints
902                                -1 means another errors in parameters passed
903                                   (N<=0, for example)
904            B   -   barycentric interpolant.
905            Rep -   report, same format as in LSFitLinearWC() subroutine.
906                    Following fields are set:
907                    * DBest         best value of the D parameter
908                    * RMSError      rms error on the (X,Y).
909                    * AvgError      average error on the (X,Y).
910                    * AvgRelError   average relative error on the non-zero Y
911                    * MaxError      maximum error
912                                    NON-WEIGHTED ERRORS ARE CALCULATED
913
914        IMPORTANT:
915            this subroitine doesn't calculate task's condition number for K<>0.
916
917        SETTING CONSTRAINTS - DANGERS AND OPPORTUNITIES:
918
919        Setting constraints can lead  to undesired  results,  like ill-conditioned
920        behavior, or inconsistency being detected. From the other side,  it allows
921        us to improve quality of the fit. Here we summarize  our  experience  with
922        constrained barycentric interpolants:
923        * excessive  constraints  can  be  inconsistent.   Floater-Hormann   basis
924          functions aren't as flexible as splines (although they are very smooth).
925        * the more evenly constraints are spread across [min(x),max(x)],  the more
926          chances that they will be consistent
927        * the  greater  is  M (given  fixed  constraints),  the  more chances that
928          constraints will be consistent
929        * in the general case, consistency of constraints IS NOT GUARANTEED.
930        * in the several special cases, however, we CAN guarantee consistency.
931        * one of this cases is constraints on the function  VALUES at the interval
932          boundaries. Note that consustency of the  constraints  on  the  function
933          DERIVATIVES is NOT guaranteed (you can use in such cases  cubic  splines
934          which are more flexible).
935        * another  special  case  is ONE constraint on the function value (OR, but
936          not AND, derivative) anywhere in the interval
937
938        Our final recommendation is to use constraints  WHEN  AND  ONLY  WHEN  you
939        can't solve your task without them. Anything beyond  special  cases  given
940        above is not guaranteed and may result in inconsistency.
941
942          -- ALGLIB PROJECT --
943             Copyright 18.08.2009 by Bochkanov Sergey
944        *************************************************************************/
945        public static void barycentricfitfloaterhormannwc(ref double[] x,
946            ref double[] y,
947            ref double[] w,
948            int n,
949            ref double[] xc,
950            ref double[] yc,
951            ref int[] dc,
952            int k,
953            int m,
954            ref int info,
955            ref barycentricinterpolant b,
956            ref barycentricfitreport rep)
957        {
958            int d = 0;
959            int i = 0;
960            double wrmscur = 0;
961            double wrmsbest = 0;
962            barycentricinterpolant locb = new barycentricinterpolant();
963            barycentricfitreport locrep = new barycentricfitreport();
964            int locinfo = 0;
965
966            if( n<1 | m<2 | k<0 | k>=m )
967            {
968                info = -1;
969                return;
970            }
971           
972            //
973            // Find optimal D
974            //
975            // Info is -3 by default (degenerate constraints).
976            // If LocInfo will always be equal to -3, Info will remain equal to -3.
977            // If at least once LocInfo will be -4, Info will be -4.
978            //
979            wrmsbest = AP.Math.MaxRealNumber;
980            rep.dbest = -1;
981            info = -3;
982            for(d=0; d<=Math.Min(9, n-1); d++)
983            {
984                barycentricfitwcfixedd(x, y, ref w, n, xc, yc, ref dc, k, m, d, ref locinfo, ref locb, ref locrep);
985                System.Diagnostics.Debug.Assert(locinfo==-4 | locinfo==-3 | locinfo>0, "BarycentricFitFloaterHormannWC: unexpected result from BarycentricFitWCFixedD!");
986                if( locinfo>0 )
987                {
988                   
989                    //
990                    // Calculate weghted RMS
991                    //
992                    wrmscur = 0;
993                    for(i=0; i<=n-1; i++)
994                    {
995                        wrmscur = wrmscur+AP.Math.Sqr(w[i]*(y[i]-barycentriccalc(ref locb, x[i])));
996                    }
997                    wrmscur = Math.Sqrt(wrmscur/n);
998                    if( (double)(wrmscur)<(double)(wrmsbest) | (double)(rep.dbest)<(double)(0) )
999                    {
1000                        barycentriccopy(ref locb, ref b);
1001                        rep.dbest = d;
1002                        info = 1;
1003                        rep.rmserror = locrep.rmserror;
1004                        rep.avgerror = locrep.avgerror;
1005                        rep.avgrelerror = locrep.avgrelerror;
1006                        rep.maxerror = locrep.maxerror;
1007                        rep.taskrcond = locrep.taskrcond;
1008                        wrmsbest = wrmscur;
1009                    }
1010                }
1011                else
1012                {
1013                    if( locinfo!=-3 & info<0 )
1014                    {
1015                        info = locinfo;
1016                    }
1017                }
1018            }
1019        }
1020
1021
1022        /*************************************************************************
1023        Rational least squares fitting, without weights and constraints.
1024
1025        See BarycentricFitFloaterHormannWC() for more information.
1026
1027          -- ALGLIB PROJECT --
1028             Copyright 18.08.2009 by Bochkanov Sergey
1029        *************************************************************************/
1030        public static void barycentricfitfloaterhormann(ref double[] x,
1031            ref double[] y,
1032            int n,
1033            int m,
1034            ref int info,
1035            ref barycentricinterpolant b,
1036            ref barycentricfitreport rep)
1037        {
1038            double[] w = new double[0];
1039            double[] xc = new double[0];
1040            double[] yc = new double[0];
1041            int[] dc = new int[0];
1042            int i = 0;
1043
1044            if( n<1 )
1045            {
1046                info = -1;
1047                return;
1048            }
1049            w = new double[n];
1050            for(i=0; i<=n-1; i++)
1051            {
1052                w[i] = 1;
1053            }
1054            barycentricfitfloaterhormannwc(ref x, ref y, ref w, n, ref xc, ref yc, ref dc, 0, m, ref info, ref b, ref rep);
1055        }
1056
1057
1058        /*************************************************************************
1059        Normalization of barycentric interpolant.
1060        Internal subroutine
1061        *************************************************************************/
1062        private static void barycentricnormalize(ref barycentricinterpolant b)
1063        {
1064            int[] p1 = new int[0];
1065            int[] p2 = new int[0];
1066            int i = 0;
1067            int j = 0;
1068            int j2 = 0;
1069            double v = 0;
1070            int i_ = 0;
1071
1072           
1073            //
1074            // Normalize task: |Y|<=1, |W|<=1, sort X[]
1075            //
1076            b.sy = 0;
1077            for(i=0; i<=b.n-1; i++)
1078            {
1079                b.sy = Math.Max(b.sy, Math.Abs(b.y[i]));
1080            }
1081            if( (double)(b.sy)>(double)(0) & (double)(Math.Abs(b.sy-1))>(double)(10*AP.Math.MachineEpsilon) )
1082            {
1083                v = 1/b.sy;
1084                for(i_=0; i_<=b.n-1;i_++)
1085                {
1086                    b.y[i_] = v*b.y[i_];
1087                }
1088            }
1089            v = 0;
1090            for(i=0; i<=b.n-1; i++)
1091            {
1092                v = Math.Max(v, Math.Abs(b.w[i]));
1093            }
1094            if( (double)(v)>(double)(0) & (double)(Math.Abs(v-1))>(double)(10*AP.Math.MachineEpsilon) )
1095            {
1096                v = 1/v;
1097                for(i_=0; i_<=b.n-1;i_++)
1098                {
1099                    b.w[i_] = v*b.w[i_];
1100                }
1101            }
1102            for(i=0; i<=b.n-2; i++)
1103            {
1104                if( (double)(b.x[i+1])<(double)(b.x[i]) )
1105                {
1106                    tsort.tagsort(ref b.x, b.n, ref p1, ref p2);
1107                    for(j=0; j<=b.n-1; j++)
1108                    {
1109                        j2 = p2[j];
1110                        v = b.y[j];
1111                        b.y[j] = b.y[j2];
1112                        b.y[j2] = v;
1113                        v = b.w[j];
1114                        b.w[j] = b.w[j2];
1115                        b.w[j2] = v;
1116                    }
1117                    break;
1118                }
1119            }
1120        }
1121
1122
1123        /*************************************************************************
1124        Internal subroutine, calculates barycentric basis functions.
1125        Used for efficient simultaneous calculation of N basis functions.
1126
1127          -- ALGLIB --
1128             Copyright 17.08.2009 by Bochkanov Sergey
1129        *************************************************************************/
1130        private static void barycentriccalcbasis(ref barycentricinterpolant b,
1131            double t,
1132            ref double[] y)
1133        {
1134            double s2 = 0;
1135            double s = 0;
1136            double v = 0;
1137            int i = 0;
1138            int j = 0;
1139            int i_ = 0;
1140
1141           
1142            //
1143            // special case: N=1
1144            //
1145            if( b.n==1 )
1146            {
1147                y[0] = 1;
1148                return;
1149            }
1150           
1151            //
1152            // Here we assume that task is normalized, i.e.:
1153            // 1. abs(Y[i])<=1
1154            // 2. abs(W[i])<=1
1155            // 3. X[] is ordered
1156            //
1157            // First, we decide: should we use "safe" formula (guarded
1158            // against overflow) or fast one?
1159            //
1160            s = Math.Abs(t-b.x[0]);
1161            for(i=0; i<=b.n-1; i++)
1162            {
1163                v = b.x[i];
1164                if( (double)(v)==(double)(t) )
1165                {
1166                    for(j=0; j<=b.n-1; j++)
1167                    {
1168                        y[j] = 0;
1169                    }
1170                    y[i] = 1;
1171                    return;
1172                }
1173                v = Math.Abs(t-v);
1174                if( (double)(v)<(double)(s) )
1175                {
1176                    s = v;
1177                }
1178            }
1179            s2 = 0;
1180            for(i=0; i<=b.n-1; i++)
1181            {
1182                v = s/(t-b.x[i]);
1183                v = v*b.w[i];
1184                y[i] = v;
1185                s2 = s2+v;
1186            }
1187            v = 1/s2;
1188            for(i_=0; i_<=b.n-1;i_++)
1189            {
1190                y[i_] = v*y[i_];
1191            }
1192        }
1193
1194
1195        /*************************************************************************
1196        Internal Floater-Hormann fitting subroutine for fixed D
1197        *************************************************************************/
1198        private static void barycentricfitwcfixedd(double[] x,
1199            double[] y,
1200            ref double[] w,
1201            int n,
1202            double[] xc,
1203            double[] yc,
1204            ref int[] dc,
1205            int k,
1206            int m,
1207            int d,
1208            ref int info,
1209            ref barycentricinterpolant b,
1210            ref barycentricfitreport rep)
1211        {
1212            double[,] fmatrix = new double[0,0];
1213            double[,] cmatrix = new double[0,0];
1214            double[] y2 = new double[0];
1215            double[] w2 = new double[0];
1216            double[] sx = new double[0];
1217            double[] sy = new double[0];
1218            double[] sbf = new double[0];
1219            double[] xoriginal = new double[0];
1220            double[] yoriginal = new double[0];
1221            double[] tmp = new double[0];
1222            lsfit.lsfitreport lrep = new lsfit.lsfitreport();
1223            double v = 0;
1224            double v0 = 0;
1225            double v1 = 0;
1226            double v2 = 0;
1227            double mx = 0;
1228            barycentricinterpolant b2 = new barycentricinterpolant();
1229            int i = 0;
1230            int j = 0;
1231            int relcnt = 0;
1232            double xa = 0;
1233            double xb = 0;
1234            double sa = 0;
1235            double sb = 0;
1236            double decay = 0;
1237            int i_ = 0;
1238
1239            x = (double[])x.Clone();
1240            y = (double[])y.Clone();
1241            xc = (double[])xc.Clone();
1242            yc = (double[])yc.Clone();
1243
1244            if( n<1 | m<2 | k<0 | k>=m )
1245            {
1246                info = -1;
1247                return;
1248            }
1249            for(i=0; i<=k-1; i++)
1250            {
1251                info = 0;
1252                if( dc[i]<0 )
1253                {
1254                    info = -1;
1255                }
1256                if( dc[i]>1 )
1257                {
1258                    info = -1;
1259                }
1260                if( info<0 )
1261                {
1262                    return;
1263                }
1264            }
1265           
1266            //
1267            // weight decay for correct handling of task which becomes
1268            // degenerate after constraints are applied
1269            //
1270            decay = 10000*AP.Math.MachineEpsilon;
1271           
1272            //
1273            // Scale X, Y, XC, YC
1274            //
1275            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);
1276           
1277            //
1278            // allocate space, initialize:
1279            // * FMatrix-   values of basis functions at X[]
1280            // * CMatrix-   values (derivatives) of basis functions at XC[]
1281            //
1282            y2 = new double[n+m];
1283            w2 = new double[n+m];
1284            fmatrix = new double[n+m, m];
1285            if( k>0 )
1286            {
1287                cmatrix = new double[k, m+1];
1288            }
1289            y2 = new double[n+m];
1290            w2 = new double[n+m];
1291           
1292            //
1293            // Prepare design and constraints matrices:
1294            // * fill constraints matrix
1295            // * fill first N rows of design matrix with values
1296            // * fill next M rows of design matrix with regularizing term
1297            // * append M zeros to Y
1298            // * append M elements, mean(abs(W)) each, to W
1299            //
1300            sx = new double[m];
1301            sy = new double[m];
1302            sbf = new double[m];
1303            for(j=0; j<=m-1; j++)
1304            {
1305                sx[j] = (double)(2*j)/((double)(m-1))-1;
1306            }
1307            for(i=0; i<=m-1; i++)
1308            {
1309                sy[i] = 1;
1310            }
1311            barycentricbuildfloaterhormann(ref sx, ref sy, m, d, ref b2);
1312            mx = 0;
1313            for(i=0; i<=n-1; i++)
1314            {
1315                barycentriccalcbasis(ref b2, x[i], ref sbf);
1316                for(i_=0; i_<=m-1;i_++)
1317                {
1318                    fmatrix[i,i_] = sbf[i_];
1319                }
1320                y2[i] = y[i];
1321                w2[i] = w[i];
1322                mx = mx+Math.Abs(w[i])/n;
1323            }
1324            for(i=0; i<=m-1; i++)
1325            {
1326                for(j=0; j<=m-1; j++)
1327                {
1328                    if( i==j )
1329                    {
1330                        fmatrix[n+i,j] = decay;
1331                    }
1332                    else
1333                    {
1334                        fmatrix[n+i,j] = 0;
1335                    }
1336                }
1337                y2[n+i] = 0;
1338                w2[n+i] = mx;
1339            }
1340            if( k>0 )
1341            {
1342                for(j=0; j<=m-1; j++)
1343                {
1344                    for(i=0; i<=m-1; i++)
1345                    {
1346                        sy[i] = 0;
1347                    }
1348                    sy[j] = 1;
1349                    barycentricbuildfloaterhormann(ref sx, ref sy, m, d, ref b2);
1350                    for(i=0; i<=k-1; i++)
1351                    {
1352                        System.Diagnostics.Debug.Assert(dc[i]>=0 & dc[i]<=1, "BarycentricFit: internal error!");
1353                        barycentricdiff1(ref b2, xc[i], ref v0, ref v1);
1354                        if( dc[i]==0 )
1355                        {
1356                            cmatrix[i,j] = v0;
1357                        }
1358                        if( dc[i]==1 )
1359                        {
1360                            cmatrix[i,j] = v1;
1361                        }
1362                    }
1363                }
1364                for(i=0; i<=k-1; i++)
1365                {
1366                    cmatrix[i,m] = yc[i];
1367                }
1368            }
1369           
1370            //
1371            // Solve constrained task
1372            //
1373            if( k>0 )
1374            {
1375               
1376                //
1377                // solve using regularization
1378                //
1379                lsfit.lsfitlinearwc(y2, ref w2, ref fmatrix, cmatrix, n+m, m, k, ref info, ref tmp, ref lrep);
1380            }
1381            else
1382            {
1383               
1384                //
1385                // no constraints, no regularization needed
1386                //
1387                lsfit.lsfitlinearwc(y, ref w, ref fmatrix, cmatrix, n, m, k, ref info, ref tmp, ref lrep);
1388            }
1389            if( info<0 )
1390            {
1391                return;
1392            }
1393           
1394            //
1395            // Generate interpolant and scale it
1396            //
1397            for(i_=0; i_<=m-1;i_++)
1398            {
1399                sy[i_] = tmp[i_];
1400            }
1401            barycentricbuildfloaterhormann(ref sx, ref sy, m, d, ref b);
1402            barycentriclintransx(ref b, 2/(xb-xa), -((xa+xb)/(xb-xa)));
1403            barycentriclintransy(ref b, sb-sa, sa);
1404           
1405            //
1406            // Scale absolute errors obtained from LSFitLinearW.
1407            // Relative error should be calculated separately
1408            // (because of shifting/scaling of the task)
1409            //
1410            rep.taskrcond = lrep.taskrcond;
1411            rep.rmserror = lrep.rmserror*(sb-sa);
1412            rep.avgerror = lrep.avgerror*(sb-sa);
1413            rep.maxerror = lrep.maxerror*(sb-sa);
1414            rep.avgrelerror = 0;
1415            relcnt = 0;
1416            for(i=0; i<=n-1; i++)
1417            {
1418                if( (double)(yoriginal[i])!=(double)(0) )
1419                {
1420                    rep.avgrelerror = rep.avgrelerror+Math.Abs(barycentriccalc(ref b, xoriginal[i])-yoriginal[i])/Math.Abs(yoriginal[i]);
1421                    relcnt = relcnt+1;
1422                }
1423            }
1424            if( relcnt!=0 )
1425            {
1426                rep.avgrelerror = rep.avgrelerror/relcnt;
1427            }
1428        }
1429    }
1430}
Note: See TracBrowser for help on using the repository browser.