Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/ALGLIB/spline1d.cs @ 2601

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

Updated ALGLIB to latest version. #751 (Plugin for for data-modeling with ANN (integrated into CEDMA))

File size: 60.9 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 spline1d
26    {
27        /*************************************************************************
28        1-dimensional spline inteprolant
29        *************************************************************************/
30        public struct spline1dinterpolant
31        {
32            public int n;
33            public int k;
34            public double[] x;
35            public double[] c;
36        };
37
38
39        /*************************************************************************
40        Spline fitting report:
41            TaskRCond       reciprocal of task's condition number
42            RMSError        RMS error
43            AvgError        average error
44            AvgRelError     average relative error (for non-zero Y[I])
45            MaxError        maximum error
46        *************************************************************************/
47        public struct spline1dfitreport
48        {
49            public double taskrcond;
50            public double rmserror;
51            public double avgerror;
52            public double avgrelerror;
53            public double maxerror;
54        };
55
56
57
58
59        public const int spline1dvnum = 11;
60
61
62        /*************************************************************************
63        This subroutine builds linear spline interpolant
64
65        INPUT PARAMETERS:
66            X   -   spline nodes, array[0..N-1]
67            Y   -   function values, array[0..N-1]
68            N   -   points count, N>=2
69           
70        OUTPUT PARAMETERS:
71            C   -   spline interpolant
72
73          -- ALGLIB PROJECT --
74             Copyright 24.06.2007 by Bochkanov Sergey
75        *************************************************************************/
76        public static void spline1dbuildlinear(double[] x,
77            double[] y,
78            int n,
79            ref spline1dinterpolant c)
80        {
81            int i = 0;
82
83            x = (double[])x.Clone();
84            y = (double[])y.Clone();
85
86            System.Diagnostics.Debug.Assert(n>1, "Spline1DBuildLinear: N<2!");
87           
88            //
89            // Sort points
90            //
91            heapsortpoints(ref x, ref y, n);
92           
93            //
94            // Build
95            //
96            c.n = n;
97            c.k = 3;
98            c.x = new double[n];
99            c.c = new double[4*(n-1)];
100            for(i=0; i<=n-1; i++)
101            {
102                c.x[i] = x[i];
103            }
104            for(i=0; i<=n-2; i++)
105            {
106                c.c[4*i+0] = y[i];
107                c.c[4*i+1] = (y[i+1]-y[i])/(x[i+1]-x[i]);
108                c.c[4*i+2] = 0;
109                c.c[4*i+3] = 0;
110            }
111        }
112
113
114        /*************************************************************************
115        This subroutine builds cubic spline interpolant.
116
117        INPUT PARAMETERS:
118            X           -   spline nodes, array[0..N-1]
119            Y           -   function values, array[0..N-1]
120            N           -   points count, N>=2
121            BoundLType  -   boundary condition type for the left boundary
122            BoundL      -   left boundary condition (first or second derivative,
123                            depending on the BoundLType)
124            BoundRType  -   boundary condition type for the right boundary
125            BoundR      -   right boundary condition (first or second derivative,
126                            depending on the BoundRType)
127
128        OUTPUT PARAMETERS:
129            C           -   spline interpolant
130                           
131        The BoundLType/BoundRType parameters can have the following values:
132            * 0, which  corresponds  to  the  parabolically   terminated  spline
133                 (BoundL/BoundR are ignored).
134            * 1, which corresponds to the first derivative boundary condition
135            * 2, which corresponds to the second derivative boundary condition
136
137          -- ALGLIB PROJECT --
138             Copyright 23.06.2007 by Bochkanov Sergey
139        *************************************************************************/
140        public static void spline1dbuildcubic(double[] x,
141            double[] y,
142            int n,
143            int boundltype,
144            double boundl,
145            int boundrtype,
146            double boundr,
147            ref spline1dinterpolant c)
148        {
149            double[] a1 = new double[0];
150            double[] a2 = new double[0];
151            double[] a3 = new double[0];
152            double[] b = new double[0];
153            double[] d = new double[0];
154            int i = 0;
155            int tblsize = 0;
156            double delta = 0;
157            double delta2 = 0;
158            double delta3 = 0;
159
160            x = (double[])x.Clone();
161            y = (double[])y.Clone();
162
163            System.Diagnostics.Debug.Assert(n>=2, "BuildCubicSpline: N<2!");
164            System.Diagnostics.Debug.Assert(boundltype==0 | boundltype==1 | boundltype==2, "BuildCubicSpline: incorrect BoundLType!");
165            System.Diagnostics.Debug.Assert(boundrtype==0 | boundrtype==1 | boundrtype==2, "BuildCubicSpline: incorrect BoundRType!");
166            a1 = new double[n];
167            a2 = new double[n];
168            a3 = new double[n];
169            b = new double[n];
170           
171            //
172            // Special case:
173            // * N=2
174            // * parabolic terminated boundary condition on both ends
175            //
176            if( n==2 & boundltype==0 & boundrtype==0 )
177            {
178               
179                //
180                // Change task type
181                //
182                boundltype = 2;
183                boundl = 0;
184                boundrtype = 2;
185                boundr = 0;
186            }
187           
188            //
189            //
190            // Sort points
191            //
192            heapsortpoints(ref x, ref y, n);
193           
194            //
195            // Left boundary conditions
196            //
197            if( boundltype==0 )
198            {
199                a1[0] = 0;
200                a2[0] = 1;
201                a3[0] = 1;
202                b[0] = 2*(y[1]-y[0])/(x[1]-x[0]);
203            }
204            if( boundltype==1 )
205            {
206                a1[0] = 0;
207                a2[0] = 1;
208                a3[0] = 0;
209                b[0] = boundl;
210            }
211            if( boundltype==2 )
212            {
213                a1[0] = 0;
214                a2[0] = 2;
215                a3[0] = 1;
216                b[0] = 3*(y[1]-y[0])/(x[1]-x[0])-0.5*boundl*(x[1]-x[0]);
217            }
218           
219            //
220            // Central conditions
221            //
222            for(i=1; i<=n-2; i++)
223            {
224                a1[i] = x[i+1]-x[i];
225                a2[i] = 2*(x[i+1]-x[i-1]);
226                a3[i] = x[i]-x[i-1];
227                b[i] = 3*(y[i]-y[i-1])/(x[i]-x[i-1])*(x[i+1]-x[i])+3*(y[i+1]-y[i])/(x[i+1]-x[i])*(x[i]-x[i-1]);
228            }
229           
230            //
231            // Right boundary conditions
232            //
233            if( boundrtype==0 )
234            {
235                a1[n-1] = 1;
236                a2[n-1] = 1;
237                a3[n-1] = 0;
238                b[n-1] = 2*(y[n-1]-y[n-2])/(x[n-1]-x[n-2]);
239            }
240            if( boundrtype==1 )
241            {
242                a1[n-1] = 0;
243                a2[n-1] = 1;
244                a3[n-1] = 0;
245                b[n-1] = boundr;
246            }
247            if( boundrtype==2 )
248            {
249                a1[n-1] = 1;
250                a2[n-1] = 2;
251                a3[n-1] = 0;
252                b[n-1] = 3*(y[n-1]-y[n-2])/(x[n-1]-x[n-2])+0.5*boundr*(x[n-1]-x[n-2]);
253            }
254           
255            //
256            // Solve
257            //
258            solvetridiagonal(a1, a2, a3, b, n, ref d);
259           
260            //
261            // Now problem is reduced to the cubic Hermite spline
262            //
263            spline1dbuildhermite(x, y, d, n, ref c);
264        }
265
266
267        /*************************************************************************
268        This subroutine builds Hermite spline interpolant.
269
270        INPUT PARAMETERS:
271            X           -   spline nodes, array[0..N-1]
272            Y           -   function values, array[0..N-1]
273            D           -   derivatives, array[0..N-1]
274            N           -   points count, N>=2
275
276        OUTPUT PARAMETERS:
277            C           -   spline interpolant.
278
279          -- ALGLIB PROJECT --
280             Copyright 23.06.2007 by Bochkanov Sergey
281        *************************************************************************/
282        public static void spline1dbuildhermite(double[] x,
283            double[] y,
284            double[] d,
285            int n,
286            ref spline1dinterpolant c)
287        {
288            int i = 0;
289            int tblsize = 0;
290            double delta = 0;
291            double delta2 = 0;
292            double delta3 = 0;
293
294            x = (double[])x.Clone();
295            y = (double[])y.Clone();
296            d = (double[])d.Clone();
297
298            System.Diagnostics.Debug.Assert(n>=2, "BuildHermiteSpline: N<2!");
299           
300            //
301            // Sort points
302            //
303            heapsortdpoints(ref x, ref y, ref d, n);
304           
305            //
306            // Build
307            //
308            c.x = new double[n];
309            c.c = new double[4*(n-1)];
310            c.k = 3;
311            c.n = n;
312            for(i=0; i<=n-1; i++)
313            {
314                c.x[i] = x[i];
315            }
316            for(i=0; i<=n-2; i++)
317            {
318                delta = x[i+1]-x[i];
319                delta2 = AP.Math.Sqr(delta);
320                delta3 = delta*delta2;
321                c.c[4*i+0] = y[i];
322                c.c[4*i+1] = d[i];
323                c.c[4*i+2] = (3*(y[i+1]-y[i])-2*d[i]*delta-d[i+1]*delta)/delta2;
324                c.c[4*i+3] = (2*(y[i]-y[i+1])+d[i]*delta+d[i+1]*delta)/delta3;
325            }
326        }
327
328
329        /*************************************************************************
330        This subroutine builds Akima spline interpolant
331
332        INPUT PARAMETERS:
333            X           -   spline nodes, array[0..N-1]
334            Y           -   function values, array[0..N-1]
335            N           -   points count, N>=5
336
337        OUTPUT PARAMETERS:
338            C           -   spline interpolant
339
340          -- ALGLIB PROJECT --
341             Copyright 24.06.2007 by Bochkanov Sergey
342        *************************************************************************/
343        public static void spline1dbuildakima(double[] x,
344            double[] y,
345            int n,
346            ref spline1dinterpolant c)
347        {
348            int i = 0;
349            double[] d = new double[0];
350            double[] w = new double[0];
351            double[] diff = new double[0];
352
353            x = (double[])x.Clone();
354            y = (double[])y.Clone();
355
356            System.Diagnostics.Debug.Assert(n>=5, "BuildAkimaSpline: N<5!");
357           
358            //
359            // Sort points
360            //
361            heapsortpoints(ref x, ref y, n);
362           
363            //
364            // Prepare W (weights), Diff (divided differences)
365            //
366            w = new double[n-1];
367            diff = new double[n-1];
368            for(i=0; i<=n-2; i++)
369            {
370                diff[i] = (y[i+1]-y[i])/(x[i+1]-x[i]);
371            }
372            for(i=1; i<=n-2; i++)
373            {
374                w[i] = Math.Abs(diff[i]-diff[i-1]);
375            }
376           
377            //
378            // Prepare Hermite interpolation scheme
379            //
380            d = new double[n];
381            for(i=2; i<=n-3; i++)
382            {
383                if( (double)(Math.Abs(w[i-1])+Math.Abs(w[i+1]))!=(double)(0) )
384                {
385                    d[i] = (w[i+1]*diff[i-1]+w[i-1]*diff[i])/(w[i+1]+w[i-1]);
386                }
387                else
388                {
389                    d[i] = ((x[i+1]-x[i])*diff[i-1]+(x[i]-x[i-1])*diff[i])/(x[i+1]-x[i-1]);
390                }
391            }
392            d[0] = diffthreepoint(x[0], x[0], y[0], x[1], y[1], x[2], y[2]);
393            d[1] = diffthreepoint(x[1], x[0], y[0], x[1], y[1], x[2], y[2]);
394            d[n-2] = diffthreepoint(x[n-2], x[n-3], y[n-3], x[n-2], y[n-2], x[n-1], y[n-1]);
395            d[n-1] = diffthreepoint(x[n-1], x[n-3], y[n-3], x[n-2], y[n-2], x[n-1], y[n-1]);
396           
397            //
398            // Build Akima spline using Hermite interpolation scheme
399            //
400            spline1dbuildhermite(x, y, d, n, ref c);
401        }
402
403
404        /*************************************************************************
405        Weighted fitting by cubic  spline,  with constraints on function values or
406        derivatives.
407
408        Equidistant grid with M-2 nodes on [min(x,xc),max(x,xc)] is  used to build
409        basis functions. Basis functions are cubic splines with continuous  second
410        derivatives  and  non-fixed first  derivatives  at  interval  ends.  Small
411        regularizing term is used  when  solving  constrained  tasks  (to  improve
412        stability).
413
414        Task is linear, so linear least squares solver is used. Complexity of this
415        computational scheme is O(N*M^2), mostly dominated by least squares solver
416
417        SEE ALSO
418            Spline1DFitHermiteWC()  -   fitting by Hermite splines (more flexible,
419                                        less smooth)
420            Spline1DFitCubic()      -   "lightweight" fitting  by  cubic  splines,
421                                        without invididual weights and constraints
422
423        INPUT PARAMETERS:
424            X   -   points, array[0..N-1].
425            Y   -   function values, array[0..N-1].
426            W   -   weights, array[0..N-1]
427                    Each summand in square  sum  of  approximation deviations from
428                    given  values  is  multiplied  by  the square of corresponding
429                    weight. Fill it by 1's if you don't  want  to  solve  weighted
430                    task.
431            N   -   number of points, N>0.
432            XC  -   points where spline values/derivatives are constrained,
433                    array[0..K-1].
434            YC  -   values of constraints, array[0..K-1]
435            DC  -   array[0..K-1], types of constraints:
436                    * DC[i]=0   means that S(XC[i])=YC[i]
437                    * DC[i]=1   means that S'(XC[i])=YC[i]
438                    SEE BELOW FOR IMPORTANT INFORMATION ON CONSTRAINTS
439            K   -   number of constraints, 0<=K<M.
440                    K=0 means no constraints (XC/YC/DC are not used in such cases)
441            M   -   number of basis functions ( = number_of_nodes+2), M>=4.
442
443        OUTPUT PARAMETERS:
444            Info-   same format as in LSFitLinearWC() subroutine.
445                    * Info>0    task is solved
446                    * Info<=0   an error occured:
447                                -4 means inconvergence of internal SVD
448                                -3 means inconsistent constraints
449                                -1 means another errors in parameters passed
450                                   (N<=0, for example)
451            S   -   spline interpolant.
452            Rep -   report, same format as in LSFitLinearWC() subroutine.
453                    Following fields are set:
454                    * RMSError      rms error on the (X,Y).
455                    * AvgError      average error on the (X,Y).
456                    * AvgRelError   average relative error on the non-zero Y
457                    * MaxError      maximum error
458                                    NON-WEIGHTED ERRORS ARE CALCULATED
459
460        IMPORTANT:
461            this subroitine doesn't calculate task's condition number for K<>0.
462
463        SETTING CONSTRAINTS - DANGERS AND OPPORTUNITIES:
464
465        Setting constraints can lead  to undesired  results,  like ill-conditioned
466        behavior, or inconsistency being detected. From the other side,  it allows
467        us to improve quality of the fit. Here we summarize  our  experience  with
468        constrained regression splines:
469        * excessive constraints can be inconsistent. Splines are  piecewise  cubic
470          functions, and it is easy to create an example, where  large  number  of
471          constraints  concentrated  in  small  area will result in inconsistency.
472          Just because spline is not flexible enough to satisfy all of  them.  And
473          same constraints spread across the  [min(x),max(x)]  will  be  perfectly
474          consistent.
475        * the more evenly constraints are spread across [min(x),max(x)],  the more
476          chances that they will be consistent
477        * the  greater  is  M (given  fixed  constraints),  the  more chances that
478          constraints will be consistent
479        * in the general case, consistency of constraints IS NOT GUARANTEED.
480        * in the several special cases, however, we CAN guarantee consistency.
481        * one of this cases is constraints  on  the  function  values  AND/OR  its
482          derivatives at the interval boundaries.
483        * another  special  case  is ONE constraint on the function value (OR, but
484          not AND, derivative) anywhere in the interval
485
486        Our final recommendation is to use constraints  WHEN  AND  ONLY  WHEN  you
487        can't solve your task without them. Anything beyond  special  cases  given
488        above is not guaranteed and may result in inconsistency.
489
490
491          -- ALGLIB PROJECT --
492             Copyright 18.08.2009 by Bochkanov Sergey
493        *************************************************************************/
494        public static void spline1dfitcubicwc(ref double[] x,
495            ref double[] y,
496            ref double[] w,
497            int n,
498            ref double[] xc,
499            ref double[] yc,
500            ref int[] dc,
501            int k,
502            int m,
503            ref int info,
504            ref spline1dinterpolant s,
505            ref spline1dfitreport rep)
506        {
507            spline1dfitinternal(0, x, y, ref w, n, xc, yc, ref dc, k, m, ref info, ref s, ref rep);
508        }
509
510
511        /*************************************************************************
512        Weighted  fitting  by Hermite spline,  with constraints on function values
513        or first derivatives.
514
515        Equidistant grid with M nodes on [min(x,xc),max(x,xc)] is  used  to  build
516        basis functions. Basis functions are Hermite splines.  Small  regularizing
517        term is used when solving constrained tasks (to improve stability).
518
519        Task is linear, so linear least squares solver is used. Complexity of this
520        computational scheme is O(N*M^2), mostly dominated by least squares solver
521
522        SEE ALSO
523            Spline1DFitCubicWC()    -   fitting by Cubic splines (less flexible,
524                                        more smooth)
525            Spline1DFitHermite()    -   "lightweight" Hermite fitting, without
526                                        invididual weights and constraints
527
528        INPUT PARAMETERS:
529            X   -   points, array[0..N-1].
530            Y   -   function values, array[0..N-1].
531            W   -   weights, array[0..N-1]
532                    Each summand in square  sum  of  approximation deviations from
533                    given  values  is  multiplied  by  the square of corresponding
534                    weight. Fill it by 1's if you don't  want  to  solve  weighted
535                    task.
536            N   -   number of points, N>0.
537            XC  -   points where spline values/derivatives are constrained,
538                    array[0..K-1].
539            YC  -   values of constraints, array[0..K-1]
540            DC  -   array[0..K-1], types of constraints:
541                    * DC[i]=0   means that S(XC[i])=YC[i]
542                    * DC[i]=1   means that S'(XC[i])=YC[i]
543                    SEE BELOW FOR IMPORTANT INFORMATION ON CONSTRAINTS
544            K   -   number of constraints, 0<=K<M.
545                    K=0 means no constraints (XC/YC/DC are not used in such cases)
546            M   -   number of basis functions (= 2 * number of nodes),
547                    M>=4,
548                    M IS EVEN!
549
550        OUTPUT PARAMETERS:
551            Info-   same format as in LSFitLinearW() subroutine:
552                    * Info>0    task is solved
553                    * Info<=0   an error occured:
554                                -4 means inconvergence of internal SVD
555                                -3 means inconsistent constraints
556                                -2 means odd M was passed (which is not supported)
557                                -1 means another errors in parameters passed
558                                   (N<=0, for example)
559            S   -   spline interpolant.
560            Rep -   report, same format as in LSFitLinearW() subroutine.
561                    Following fields are set:
562                    * RMSError      rms error on the (X,Y).
563                    * AvgError      average error on the (X,Y).
564                    * AvgRelError   average relative error on the non-zero Y
565                    * MaxError      maximum error
566                                    NON-WEIGHTED ERRORS ARE CALCULATED
567
568        IMPORTANT:
569            this subroitine doesn't calculate task's condition number for K<>0.
570
571        IMPORTANT:
572            this subroitine supports only even M's
573
574        SETTING CONSTRAINTS - DANGERS AND OPPORTUNITIES:
575
576        Setting constraints can lead  to undesired  results,  like ill-conditioned
577        behavior, or inconsistency being detected. From the other side,  it allows
578        us to improve quality of the fit. Here we summarize  our  experience  with
579        constrained regression splines:
580        * excessive constraints can be inconsistent. Splines are  piecewise  cubic
581          functions, and it is easy to create an example, where  large  number  of
582          constraints  concentrated  in  small  area will result in inconsistency.
583          Just because spline is not flexible enough to satisfy all of  them.  And
584          same constraints spread across the  [min(x),max(x)]  will  be  perfectly
585          consistent.
586        * the more evenly constraints are spread across [min(x),max(x)],  the more
587          chances that they will be consistent
588        * the  greater  is  M (given  fixed  constraints),  the  more chances that
589          constraints will be consistent
590        * in the general case, consistency of constraints is NOT GUARANTEED.
591        * in the several special cases, however, we can guarantee consistency.
592        * one of this cases is  M>=4  and   constraints  on   the  function  value
593          (AND/OR its derivative) at the interval boundaries.
594        * another special case is M>=4  and  ONE  constraint on the function value
595          (OR, BUT NOT AND, derivative) anywhere in [min(x),max(x)]
596
597        Our final recommendation is to use constraints  WHEN  AND  ONLY  when  you
598        can't solve your task without them. Anything beyond  special  cases  given
599        above is not guaranteed and may result in inconsistency.
600
601          -- ALGLIB PROJECT --
602             Copyright 18.08.2009 by Bochkanov Sergey
603        *************************************************************************/
604        public static void spline1dfithermitewc(ref double[] x,
605            ref double[] y,
606            ref double[] w,
607            int n,
608            ref double[] xc,
609            ref double[] yc,
610            ref int[] dc,
611            int k,
612            int m,
613            ref int info,
614            ref spline1dinterpolant s,
615            ref spline1dfitreport rep)
616        {
617            spline1dfitinternal(1, x, y, ref w, n, xc, yc, ref dc, k, m, ref info, ref s, ref rep);
618        }
619
620
621        /*************************************************************************
622        Least squares fitting by cubic spline.
623
624        This subroutine is "lightweight" alternative for more complex and feature-
625        rich Spline1DFitCubicWC().  See  Spline1DFitCubicWC() for more information
626        about subroutine parameters (we don't duplicate it here because of length)
627
628          -- ALGLIB PROJECT --
629             Copyright 18.08.2009 by Bochkanov Sergey
630        *************************************************************************/
631        public static void spline1dfitcubic(ref double[] x,
632            ref double[] y,
633            int n,
634            int m,
635            ref int info,
636            ref spline1dinterpolant s,
637            ref spline1dfitreport rep)
638        {
639            int i = 0;
640            double[] w = new double[0];
641            double[] xc = new double[0];
642            double[] yc = new double[0];
643            int[] dc = new int[0];
644
645            if( n>0 )
646            {
647                w = new double[n];
648                for(i=0; i<=n-1; i++)
649                {
650                    w[i] = 1;
651                }
652            }
653            spline1dfitcubicwc(ref x, ref y, ref w, n, ref xc, ref yc, ref dc, 0, m, ref info, ref s, ref rep);
654        }
655
656
657        /*************************************************************************
658        Least squares fitting by Hermite spline.
659
660        This subroutine is "lightweight" alternative for more complex and feature-
661        rich Spline1DFitHermiteWC().  See Spline1DFitHermiteWC()  description  for
662        more information about subroutine parameters (we don't duplicate  it  here
663        because of length).
664
665          -- ALGLIB PROJECT --
666             Copyright 18.08.2009 by Bochkanov Sergey
667        *************************************************************************/
668        public static void spline1dfithermite(ref double[] x,
669            ref double[] y,
670            int n,
671            int m,
672            ref int info,
673            ref spline1dinterpolant s,
674            ref spline1dfitreport rep)
675        {
676            int i = 0;
677            double[] w = new double[0];
678            double[] xc = new double[0];
679            double[] yc = new double[0];
680            int[] dc = new int[0];
681
682            if( n>0 )
683            {
684                w = new double[n];
685                for(i=0; i<=n-1; i++)
686                {
687                    w[i] = 1;
688                }
689            }
690            spline1dfithermitewc(ref x, ref y, ref w, n, ref xc, ref yc, ref dc, 0, m, ref info, ref s, ref rep);
691        }
692
693
694        /*************************************************************************
695        This subroutine calculates the value of the spline at the given point X.
696
697        INPUT PARAMETERS:
698            C   -   spline interpolant
699            X   -   point
700
701        Result:
702            S(x)
703
704          -- ALGLIB PROJECT --
705             Copyright 23.06.2007 by Bochkanov Sergey
706        *************************************************************************/
707        public static double spline1dcalc(ref spline1dinterpolant c,
708            double x)
709        {
710            double result = 0;
711            int l = 0;
712            int r = 0;
713            int m = 0;
714
715            System.Diagnostics.Debug.Assert(c.k==3, "Spline1DCalc: internal error");
716           
717            //
718            // Binary search in the [ x[0], ..., x[n-2] ] (x[n-1] is not included)
719            //
720            l = 0;
721            r = c.n-2+1;
722            while( l!=r-1 )
723            {
724                m = (l+r)/2;
725                if( (double)(c.x[m])>=(double)(x) )
726                {
727                    r = m;
728                }
729                else
730                {
731                    l = m;
732                }
733            }
734           
735            //
736            // Interpolation
737            //
738            x = x-c.x[l];
739            m = 4*l;
740            result = c.c[m]+x*(c.c[m+1]+x*(c.c[m+2]+x*c.c[m+3]));
741            return result;
742        }
743
744
745        /*************************************************************************
746        This subroutine differentiates the spline.
747
748        INPUT PARAMETERS:
749            C   -   spline interpolant.
750            X   -   point
751
752        Result:
753            S   -   S(x)
754            DS  -   S'(x)
755            D2S -   S''(x)
756
757          -- ALGLIB PROJECT --
758             Copyright 24.06.2007 by Bochkanov Sergey
759        *************************************************************************/
760        public static void spline1ddiff(ref spline1dinterpolant c,
761            double x,
762            ref double s,
763            ref double ds,
764            ref double d2s)
765        {
766            int l = 0;
767            int r = 0;
768            int m = 0;
769
770            System.Diagnostics.Debug.Assert(c.k==3, "Spline1DCalc: internal error");
771           
772            //
773            // Binary search
774            //
775            l = 0;
776            r = c.n-2+1;
777            while( l!=r-1 )
778            {
779                m = (l+r)/2;
780                if( (double)(c.x[m])>=(double)(x) )
781                {
782                    r = m;
783                }
784                else
785                {
786                    l = m;
787                }
788            }
789           
790            //
791            // Differentiation
792            //
793            x = x-c.x[l];
794            m = 4*l;
795            s = c.c[m]+x*(c.c[m+1]+x*(c.c[m+2]+x*c.c[m+3]));
796            ds = c.c[m+1]+2*x*c.c[m+2]+3*AP.Math.Sqr(x)*c.c[m+3];
797            d2s = 2*c.c[m+2]+6*x*c.c[m+3];
798        }
799
800
801        /*************************************************************************
802        This subroutine makes the copy of the spline.
803
804        INPUT PARAMETERS:
805            C   -   spline interpolant.
806
807        Result:
808            CC  -   spline copy
809
810          -- ALGLIB PROJECT --
811             Copyright 29.06.2007 by Bochkanov Sergey
812        *************************************************************************/
813        public static void spline1dcopy(ref spline1dinterpolant c,
814            ref spline1dinterpolant cc)
815        {
816            int i_ = 0;
817
818            cc.n = c.n;
819            cc.k = c.k;
820            cc.x = new double[cc.n];
821            for(i_=0; i_<=cc.n-1;i_++)
822            {
823                cc.x[i_] = c.x[i_];
824            }
825            cc.c = new double[(cc.k+1)*(cc.n-1)];
826            for(i_=0; i_<=(cc.k+1)*(cc.n-1)-1;i_++)
827            {
828                cc.c[i_] = c.c[i_];
829            }
830        }
831
832
833        /*************************************************************************
834        Serialization of the spline interpolant
835
836        INPUT PARAMETERS:
837            B   -   spline interpolant
838
839        OUTPUT PARAMETERS:
840            RA      -   array of real numbers which contains interpolant,
841                        array[0..RLen-1]
842            RLen    -   RA lenght
843
844          -- ALGLIB --
845             Copyright 17.08.2009 by Bochkanov Sergey
846        *************************************************************************/
847        public static void spline1dserialize(ref spline1dinterpolant c,
848            ref double[] ra,
849            ref int ralen)
850        {
851            int i_ = 0;
852            int i1_ = 0;
853
854            ralen = 2+2+c.n+(c.k+1)*(c.n-1);
855            ra = new double[ralen];
856            ra[0] = ralen;
857            ra[1] = spline1dvnum;
858            ra[2] = c.n;
859            ra[3] = c.k;
860            i1_ = (0) - (4);
861            for(i_=4; i_<=4+c.n-1;i_++)
862            {
863                ra[i_] = c.x[i_+i1_];
864            }
865            i1_ = (0) - (4+c.n);
866            for(i_=4+c.n; i_<=4+c.n+(c.k+1)*(c.n-1)-1;i_++)
867            {
868                ra[i_] = c.c[i_+i1_];
869            }
870        }
871
872
873        /*************************************************************************
874        Unserialization of the spline interpolant
875
876        INPUT PARAMETERS:
877            RA  -   array of real numbers which contains interpolant,
878
879        OUTPUT PARAMETERS:
880            B   -   spline interpolant
881
882          -- ALGLIB --
883             Copyright 17.08.2009 by Bochkanov Sergey
884        *************************************************************************/
885        public static void spline1dunserialize(ref double[] ra,
886            ref spline1dinterpolant c)
887        {
888            int i_ = 0;
889            int i1_ = 0;
890
891            System.Diagnostics.Debug.Assert((int)Math.Round(ra[1])==spline1dvnum, "Spline1DUnserialize: corrupted array!");
892            c.n = (int)Math.Round(ra[2]);
893            c.k = (int)Math.Round(ra[3]);
894            c.x = new double[c.n];
895            c.c = new double[(c.k+1)*(c.n-1)];
896            i1_ = (4) - (0);
897            for(i_=0; i_<=c.n-1;i_++)
898            {
899                c.x[i_] = ra[i_+i1_];
900            }
901            i1_ = (4+c.n) - (0);
902            for(i_=0; i_<=(c.k+1)*(c.n-1)-1;i_++)
903            {
904                c.c[i_] = ra[i_+i1_];
905            }
906        }
907
908
909        /*************************************************************************
910        This subroutine unpacks the spline into the coefficients table.
911
912        INPUT PARAMETERS:
913            C   -   spline interpolant.
914            X   -   point
915
916        Result:
917            Tbl -   coefficients table, unpacked format, array[0..N-2, 0..5].
918                    For I = 0...N-2:
919                        Tbl[I,0] = X[i]
920                        Tbl[I,1] = X[i+1]
921                        Tbl[I,2] = C0
922                        Tbl[I,3] = C1
923                        Tbl[I,4] = C2
924                        Tbl[I,5] = C3
925                    On [x[i], x[i+1]] spline is equals to:
926                        S(x) = C0 + C1*t + C2*t^2 + C3*t^3
927                        t = x-x[i]
928
929          -- ALGLIB PROJECT --
930             Copyright 29.06.2007 by Bochkanov Sergey
931        *************************************************************************/
932        public static void spline1dunpack(ref spline1dinterpolant c,
933            ref int n,
934            ref double[,] tbl)
935        {
936            int i = 0;
937            int j = 0;
938
939            tbl = new double[c.n-2+1, 2+c.k+1];
940            n = c.n;
941           
942            //
943            // Fill
944            //
945            for(i=0; i<=n-2; i++)
946            {
947                tbl[i,0] = c.x[i];
948                tbl[i,1] = c.x[i+1];
949                for(j=0; j<=c.k; j++)
950                {
951                    tbl[i,2+j] = c.c[(c.k+1)*i+j];
952                }
953            }
954        }
955
956
957        /*************************************************************************
958        This subroutine performs linear transformation of the spline argument.
959
960        INPUT PARAMETERS:
961            C   -   spline interpolant.
962            A, B-   transformation coefficients: x = A*t + B
963        Result:
964            C   -   transformed spline
965
966          -- ALGLIB PROJECT --
967             Copyright 30.06.2007 by Bochkanov Sergey
968        *************************************************************************/
969        public static void spline1dlintransx(ref spline1dinterpolant c,
970            double a,
971            double b)
972        {
973            int i = 0;
974            int j = 0;
975            int n = 0;
976            double v = 0;
977            double dv = 0;
978            double d2v = 0;
979            double[] x = new double[0];
980            double[] y = new double[0];
981            double[] d = new double[0];
982
983            n = c.n;
984           
985            //
986            // Special case: A=0
987            //
988            if( (double)(a)==(double)(0) )
989            {
990                v = spline1dcalc(ref c, b);
991                for(i=0; i<=n-2; i++)
992                {
993                    c.c[(c.k+1)*i] = v;
994                    for(j=1; j<=c.k; j++)
995                    {
996                        c.c[(c.k+1)*i+j] = 0;
997                    }
998                }
999                return;
1000            }
1001           
1002            //
1003            // General case: A<>0.
1004            // Unpack, X, Y, dY/dX.
1005            // Scale and pack again.
1006            //
1007            System.Diagnostics.Debug.Assert(c.k==3, "Spline1DLinTransX: internal error");
1008            x = new double[n-1+1];
1009            y = new double[n-1+1];
1010            d = new double[n-1+1];
1011            for(i=0; i<=n-1; i++)
1012            {
1013                x[i] = c.x[i];
1014                spline1ddiff(ref c, x[i], ref v, ref dv, ref d2v);
1015                x[i] = (x[i]-b)/a;
1016                y[i] = v;
1017                d[i] = a*dv;
1018            }
1019            spline1dbuildhermite(x, y, d, n, ref c);
1020        }
1021
1022
1023        /*************************************************************************
1024        This subroutine performs linear transformation of the spline.
1025
1026        INPUT PARAMETERS:
1027            C   -   spline interpolant.
1028            A, B-   transformation coefficients: S2(x) = A*S(x) + B
1029        Result:
1030            C   -   transformed spline
1031
1032          -- ALGLIB PROJECT --
1033             Copyright 30.06.2007 by Bochkanov Sergey
1034        *************************************************************************/
1035        public static void spline1dlintransy(ref spline1dinterpolant c,
1036            double a,
1037            double b)
1038        {
1039            int i = 0;
1040            int j = 0;
1041            int n = 0;
1042
1043            n = c.n;
1044            for(i=0; i<=n-2; i++)
1045            {
1046                c.c[(c.k+1)*i] = a*c.c[(c.k+1)*i]+b;
1047                for(j=1; j<=c.k; j++)
1048                {
1049                    c.c[(c.k+1)*i+j] = a*c.c[(c.k+1)*i+j];
1050                }
1051            }
1052        }
1053
1054
1055        /*************************************************************************
1056        This subroutine integrates the spline.
1057
1058        INPUT PARAMETERS:
1059            C   -   spline interpolant.
1060            X   -   right bound of the integration interval [a, x]
1061        Result:
1062            integral(S(t)dt,a,x)
1063
1064          -- ALGLIB PROJECT --
1065             Copyright 23.06.2007 by Bochkanov Sergey
1066        *************************************************************************/
1067        public static double spline1dintegrate(ref spline1dinterpolant c,
1068            double x)
1069        {
1070            double result = 0;
1071            int n = 0;
1072            int i = 0;
1073            int j = 0;
1074            int l = 0;
1075            int r = 0;
1076            int m = 0;
1077            double w = 0;
1078            double v = 0;
1079
1080            n = c.n;
1081           
1082            //
1083            // Binary search in the [ x[0], ..., x[n-2] ] (x[n-1] is not included)
1084            //
1085            l = 0;
1086            r = n-2+1;
1087            while( l!=r-1 )
1088            {
1089                m = (l+r)/2;
1090                if( (double)(c.x[m])>=(double)(x) )
1091                {
1092                    r = m;
1093                }
1094                else
1095                {
1096                    l = m;
1097                }
1098            }
1099           
1100            //
1101            // Integration
1102            //
1103            result = 0;
1104            for(i=0; i<=l-1; i++)
1105            {
1106                w = c.x[i+1]-c.x[i];
1107                m = (c.k+1)*i;
1108                result = result+c.c[m]*w;
1109                v = w;
1110                for(j=1; j<=c.k; j++)
1111                {
1112                    v = v*w;
1113                    result = result+c.c[m+j]*v/(j+1);
1114                }
1115            }
1116            w = x-c.x[l];
1117            m = (c.k+1)*l;
1118            v = w;
1119            result = result+c.c[m]*w;
1120            for(j=1; j<=c.k; j++)
1121            {
1122                v = v*w;
1123                result = result+c.c[m+j]*v/(j+1);
1124            }
1125            return result;
1126        }
1127
1128
1129        /*************************************************************************
1130        Internal spline fitting subroutine
1131
1132          -- ALGLIB PROJECT --
1133             Copyright 08.09.2009 by Bochkanov Sergey
1134        *************************************************************************/
1135        private static void spline1dfitinternal(int st,
1136            double[] x,
1137            double[] y,
1138            ref double[] w,
1139            int n,
1140            double[] xc,
1141            double[] yc,
1142            ref int[] dc,
1143            int k,
1144            int m,
1145            ref int info,
1146            ref spline1dinterpolant s,
1147            ref spline1dfitreport rep)
1148        {
1149            double[,] fmatrix = new double[0,0];
1150            double[,] cmatrix = new double[0,0];
1151            double[] y2 = new double[0];
1152            double[] w2 = new double[0];
1153            double[] sx = new double[0];
1154            double[] sy = new double[0];
1155            double[] sd = new double[0];
1156            double[] tmp = new double[0];
1157            double[] xoriginal = new double[0];
1158            double[] yoriginal = new double[0];
1159            lsfit.lsfitreport lrep = new lsfit.lsfitreport();
1160            double v = 0;
1161            double v0 = 0;
1162            double v1 = 0;
1163            double v2 = 0;
1164            double mx = 0;
1165            spline1dinterpolant s2 = new spline1dinterpolant();
1166            int i = 0;
1167            int j = 0;
1168            int relcnt = 0;
1169            double xa = 0;
1170            double xb = 0;
1171            double sa = 0;
1172            double sb = 0;
1173            double bl = 0;
1174            double br = 0;
1175            double decay = 0;
1176            int i_ = 0;
1177
1178            x = (double[])x.Clone();
1179            y = (double[])y.Clone();
1180            xc = (double[])xc.Clone();
1181            yc = (double[])yc.Clone();
1182
1183            System.Diagnostics.Debug.Assert(st==0 | st==1, "Spline1DFit: internal error!");
1184            if( st==0 & m<4 )
1185            {
1186                info = -1;
1187                return;
1188            }
1189            if( st==1 & m<4 )
1190            {
1191                info = -1;
1192                return;
1193            }
1194            if( n<1 | k<0 | k>=m )
1195            {
1196                info = -1;
1197                return;
1198            }
1199            for(i=0; i<=k-1; i++)
1200            {
1201                info = 0;
1202                if( dc[i]<0 )
1203                {
1204                    info = -1;
1205                }
1206                if( dc[i]>1 )
1207                {
1208                    info = -1;
1209                }
1210                if( info<0 )
1211                {
1212                    return;
1213                }
1214            }
1215            if( st==1 & m%2!=0 )
1216            {
1217               
1218                //
1219                // Hermite fitter must have even number of basis functions
1220                //
1221                info = -2;
1222                return;
1223            }
1224           
1225            //
1226            // weight decay for correct handling of task which becomes
1227            // degenerate after constraints are applied
1228            //
1229            decay = 10000*AP.Math.MachineEpsilon;
1230           
1231            //
1232            // Scale X, Y, XC, YC
1233            //
1234            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);
1235           
1236            //
1237            // allocate space, initialize:
1238            // * SX     -   grid for basis functions
1239            // * SY     -   values of basis functions at grid points
1240            // * FMatrix-   values of basis functions at X[]
1241            // * CMatrix-   values (derivatives) of basis functions at XC[]
1242            //
1243            y2 = new double[n+m];
1244            w2 = new double[n+m];
1245            fmatrix = new double[n+m, m];
1246            if( k>0 )
1247            {
1248                cmatrix = new double[k, m+1];
1249            }
1250            if( st==0 )
1251            {
1252               
1253                //
1254                // allocate space for cubic spline
1255                //
1256                sx = new double[m-2];
1257                sy = new double[m-2];
1258                for(j=0; j<=m-2-1; j++)
1259                {
1260                    sx[j] = (double)(2*j)/((double)(m-2-1))-1;
1261                }
1262            }
1263            if( st==1 )
1264            {
1265               
1266                //
1267                // allocate space for Hermite spline
1268                //
1269                sx = new double[m/2];
1270                sy = new double[m/2];
1271                sd = new double[m/2];
1272                for(j=0; j<=m/2-1; j++)
1273                {
1274                    sx[j] = (double)(2*j)/((double)(m/2-1))-1;
1275                }
1276            }
1277           
1278            //
1279            // Prepare design and constraints matrices:
1280            // * fill constraints matrix
1281            // * fill first N rows of design matrix with values
1282            // * fill next M rows of design matrix with regularizing term
1283            // * append M zeros to Y
1284            // * append M elements, mean(abs(W)) each, to W
1285            //
1286            for(j=0; j<=m-1; j++)
1287            {
1288               
1289                //
1290                // prepare Jth basis function
1291                //
1292                if( st==0 )
1293                {
1294                   
1295                    //
1296                    // cubic spline basis
1297                    //
1298                    for(i=0; i<=m-2-1; i++)
1299                    {
1300                        sy[i] = 0;
1301                    }
1302                    bl = 0;
1303                    br = 0;
1304                    if( j<m-2 )
1305                    {
1306                        sy[j] = 1;
1307                    }
1308                    if( j==m-2 )
1309                    {
1310                        bl = 1;
1311                    }
1312                    if( j==m-1 )
1313                    {
1314                        br = 1;
1315                    }
1316                    spline1dbuildcubic(sx, sy, m-2, 1, bl, 1, br, ref s2);
1317                }
1318                if( st==1 )
1319                {
1320                   
1321                    //
1322                    // Hermite basis
1323                    //
1324                    for(i=0; i<=m/2-1; i++)
1325                    {
1326                        sy[i] = 0;
1327                        sd[i] = 0;
1328                    }
1329                    if( j%2==0 )
1330                    {
1331                        sy[j/2] = 1;
1332                    }
1333                    else
1334                    {
1335                        sd[j/2] = 1;
1336                    }
1337                    spline1dbuildhermite(sx, sy, sd, m/2, ref s2);
1338                }
1339               
1340                //
1341                // values at X[], XC[]
1342                //
1343                for(i=0; i<=n-1; i++)
1344                {
1345                    fmatrix[i,j] = spline1dcalc(ref s2, x[i]);
1346                }
1347                for(i=0; i<=k-1; i++)
1348                {
1349                    System.Diagnostics.Debug.Assert(dc[i]>=0 & dc[i]<=2, "Spline1DFit: internal error!");
1350                    spline1ddiff(ref s2, xc[i], ref v0, ref v1, ref v2);
1351                    if( dc[i]==0 )
1352                    {
1353                        cmatrix[i,j] = v0;
1354                    }
1355                    if( dc[i]==1 )
1356                    {
1357                        cmatrix[i,j] = v1;
1358                    }
1359                    if( dc[i]==2 )
1360                    {
1361                        cmatrix[i,j] = v2;
1362                    }
1363                }
1364            }
1365            for(i=0; i<=k-1; i++)
1366            {
1367                cmatrix[i,m] = yc[i];
1368            }
1369            for(i=0; i<=m-1; i++)
1370            {
1371                for(j=0; j<=m-1; j++)
1372                {
1373                    if( i==j )
1374                    {
1375                        fmatrix[n+i,j] = decay;
1376                    }
1377                    else
1378                    {
1379                        fmatrix[n+i,j] = 0;
1380                    }
1381                }
1382            }
1383            y2 = new double[n+m];
1384            w2 = new double[n+m];
1385            for(i_=0; i_<=n-1;i_++)
1386            {
1387                y2[i_] = y[i_];
1388            }
1389            for(i_=0; i_<=n-1;i_++)
1390            {
1391                w2[i_] = w[i_];
1392            }
1393            mx = 0;
1394            for(i=0; i<=n-1; i++)
1395            {
1396                mx = mx+Math.Abs(w[i]);
1397            }
1398            mx = mx/n;
1399            for(i=0; i<=m-1; i++)
1400            {
1401                y2[n+i] = 0;
1402                w2[n+i] = mx;
1403            }
1404           
1405            //
1406            // Solve constrained task
1407            //
1408            if( k>0 )
1409            {
1410               
1411                //
1412                // solve using regularization
1413                //
1414                lsfit.lsfitlinearwc(y2, ref w2, ref fmatrix, cmatrix, n+m, m, k, ref info, ref tmp, ref lrep);
1415            }
1416            else
1417            {
1418               
1419                //
1420                // no constraints, no regularization needed
1421                //
1422                lsfit.lsfitlinearwc(y, ref w, ref fmatrix, cmatrix, n, m, k, ref info, ref tmp, ref lrep);
1423            }
1424            if( info<0 )
1425            {
1426                return;
1427            }
1428           
1429            //
1430            // Generate spline and scale it
1431            //
1432            if( st==0 )
1433            {
1434               
1435                //
1436                // cubic spline basis
1437                //
1438                for(i_=0; i_<=m-2-1;i_++)
1439                {
1440                    sy[i_] = tmp[i_];
1441                }
1442                spline1dbuildcubic(sx, sy, m-2, 1, tmp[m-2], 1, tmp[m-1], ref s);
1443            }
1444            if( st==1 )
1445            {
1446               
1447                //
1448                // Hermite basis
1449                //
1450                for(i=0; i<=m/2-1; i++)
1451                {
1452                    sy[i] = tmp[2*i];
1453                    sd[i] = tmp[2*i+1];
1454                }
1455                spline1dbuildhermite(sx, sy, sd, m/2, ref s);
1456            }
1457            spline1dlintransx(ref s, 2/(xb-xa), -((xa+xb)/(xb-xa)));
1458            spline1dlintransy(ref s, sb-sa, sa);
1459           
1460            //
1461            // Scale absolute errors obtained from LSFitLinearW.
1462            // Relative error should be calculated separately
1463            // (because of shifting/scaling of the task)
1464            //
1465            rep.taskrcond = lrep.taskrcond;
1466            rep.rmserror = lrep.rmserror*(sb-sa);
1467            rep.avgerror = lrep.avgerror*(sb-sa);
1468            rep.maxerror = lrep.maxerror*(sb-sa);
1469            rep.avgrelerror = 0;
1470            relcnt = 0;
1471            for(i=0; i<=n-1; i++)
1472            {
1473                if( (double)(yoriginal[i])!=(double)(0) )
1474                {
1475                    rep.avgrelerror = rep.avgrelerror+Math.Abs(spline1dcalc(ref s, xoriginal[i])-yoriginal[i])/Math.Abs(yoriginal[i]);
1476                    relcnt = relcnt+1;
1477                }
1478            }
1479            if( relcnt!=0 )
1480            {
1481                rep.avgrelerror = rep.avgrelerror/relcnt;
1482            }
1483        }
1484
1485
1486        /*************************************************************************
1487        Internal subroutine. Heap sort.
1488        *************************************************************************/
1489        private static void heapsortpoints(ref double[] x,
1490            ref double[] y,
1491            int n)
1492        {
1493            int i = 0;
1494            int j = 0;
1495            int k = 0;
1496            int t = 0;
1497            double tmp = 0;
1498            bool isascending = new bool();
1499            bool isdescending = new bool();
1500
1501           
1502            //
1503            // Test for already sorted set
1504            //
1505            isascending = true;
1506            isdescending = true;
1507            for(i=1; i<=n-1; i++)
1508            {
1509                isascending = isascending & (double)(x[i])>(double)(x[i-1]);
1510                isdescending = isdescending & (double)(x[i])<(double)(x[i-1]);
1511            }
1512            if( isascending )
1513            {
1514                return;
1515            }
1516            if( isdescending )
1517            {
1518                for(i=0; i<=n-1; i++)
1519                {
1520                    j = n-1-i;
1521                    if( j<=i )
1522                    {
1523                        break;
1524                    }
1525                    tmp = x[i];
1526                    x[i] = x[j];
1527                    x[j] = tmp;
1528                    tmp = y[i];
1529                    y[i] = y[j];
1530                    y[j] = tmp;
1531                }
1532                return;
1533            }
1534           
1535            //
1536            // Special case: N=1
1537            //
1538            if( n==1 )
1539            {
1540                return;
1541            }
1542           
1543            //
1544            // General case
1545            //
1546            i = 2;
1547            do
1548            {
1549                t = i;
1550                while( t!=1 )
1551                {
1552                    k = t/2;
1553                    if( (double)(x[k-1])>=(double)(x[t-1]) )
1554                    {
1555                        t = 1;
1556                    }
1557                    else
1558                    {
1559                        tmp = x[k-1];
1560                        x[k-1] = x[t-1];
1561                        x[t-1] = tmp;
1562                        tmp = y[k-1];
1563                        y[k-1] = y[t-1];
1564                        y[t-1] = tmp;
1565                        t = k;
1566                    }
1567                }
1568                i = i+1;
1569            }
1570            while( i<=n );
1571            i = n-1;
1572            do
1573            {
1574                tmp = x[i];
1575                x[i] = x[0];
1576                x[0] = tmp;
1577                tmp = y[i];
1578                y[i] = y[0];
1579                y[0] = tmp;
1580                t = 1;
1581                while( t!=0 )
1582                {
1583                    k = 2*t;
1584                    if( k>i )
1585                    {
1586                        t = 0;
1587                    }
1588                    else
1589                    {
1590                        if( k<i )
1591                        {
1592                            if( (double)(x[k])>(double)(x[k-1]) )
1593                            {
1594                                k = k+1;
1595                            }
1596                        }
1597                        if( (double)(x[t-1])>=(double)(x[k-1]) )
1598                        {
1599                            t = 0;
1600                        }
1601                        else
1602                        {
1603                            tmp = x[k-1];
1604                            x[k-1] = x[t-1];
1605                            x[t-1] = tmp;
1606                            tmp = y[k-1];
1607                            y[k-1] = y[t-1];
1608                            y[t-1] = tmp;
1609                            t = k;
1610                        }
1611                    }
1612                }
1613                i = i-1;
1614            }
1615            while( i>=1 );
1616        }
1617
1618
1619        /*************************************************************************
1620        Internal subroutine. Heap sort.
1621        *************************************************************************/
1622        private static void heapsortdpoints(ref double[] x,
1623            ref double[] y,
1624            ref double[] d,
1625            int n)
1626        {
1627            int i = 0;
1628            int j = 0;
1629            int k = 0;
1630            int t = 0;
1631            double tmp = 0;
1632            bool isascending = new bool();
1633            bool isdescending = new bool();
1634
1635           
1636            //
1637            // Test for already sorted set
1638            //
1639            isascending = true;
1640            isdescending = true;
1641            for(i=1; i<=n-1; i++)
1642            {
1643                isascending = isascending & (double)(x[i])>(double)(x[i-1]);
1644                isdescending = isdescending & (double)(x[i])<(double)(x[i-1]);
1645            }
1646            if( isascending )
1647            {
1648                return;
1649            }
1650            if( isdescending )
1651            {
1652                for(i=0; i<=n-1; i++)
1653                {
1654                    j = n-1-i;
1655                    if( j<=i )
1656                    {
1657                        break;
1658                    }
1659                    tmp = x[i];
1660                    x[i] = x[j];
1661                    x[j] = tmp;
1662                    tmp = y[i];
1663                    y[i] = y[j];
1664                    y[j] = tmp;
1665                    tmp = d[i];
1666                    d[i] = d[j];
1667                    d[j] = tmp;
1668                }
1669                return;
1670            }
1671           
1672            //
1673            // Special case: N=1
1674            //
1675            if( n==1 )
1676            {
1677                return;
1678            }
1679           
1680            //
1681            // General case
1682            //
1683            i = 2;
1684            do
1685            {
1686                t = i;
1687                while( t!=1 )
1688                {
1689                    k = t/2;
1690                    if( (double)(x[k-1])>=(double)(x[t-1]) )
1691                    {
1692                        t = 1;
1693                    }
1694                    else
1695                    {
1696                        tmp = x[k-1];
1697                        x[k-1] = x[t-1];
1698                        x[t-1] = tmp;
1699                        tmp = y[k-1];
1700                        y[k-1] = y[t-1];
1701                        y[t-1] = tmp;
1702                        tmp = d[k-1];
1703                        d[k-1] = d[t-1];
1704                        d[t-1] = tmp;
1705                        t = k;
1706                    }
1707                }
1708                i = i+1;
1709            }
1710            while( i<=n );
1711            i = n-1;
1712            do
1713            {
1714                tmp = x[i];
1715                x[i] = x[0];
1716                x[0] = tmp;
1717                tmp = y[i];
1718                y[i] = y[0];
1719                y[0] = tmp;
1720                tmp = d[i];
1721                d[i] = d[0];
1722                d[0] = tmp;
1723                t = 1;
1724                while( t!=0 )
1725                {
1726                    k = 2*t;
1727                    if( k>i )
1728                    {
1729                        t = 0;
1730                    }
1731                    else
1732                    {
1733                        if( k<i )
1734                        {
1735                            if( (double)(x[k])>(double)(x[k-1]) )
1736                            {
1737                                k = k+1;
1738                            }
1739                        }
1740                        if( (double)(x[t-1])>=(double)(x[k-1]) )
1741                        {
1742                            t = 0;
1743                        }
1744                        else
1745                        {
1746                            tmp = x[k-1];
1747                            x[k-1] = x[t-1];
1748                            x[t-1] = tmp;
1749                            tmp = y[k-1];
1750                            y[k-1] = y[t-1];
1751                            y[t-1] = tmp;
1752                            tmp = d[k-1];
1753                            d[k-1] = d[t-1];
1754                            d[t-1] = tmp;
1755                            t = k;
1756                        }
1757                    }
1758                }
1759                i = i-1;
1760            }
1761            while( i>=1 );
1762        }
1763
1764
1765        /*************************************************************************
1766        Internal subroutine. Tridiagonal solver.
1767        *************************************************************************/
1768        private static void solvetridiagonal(double[] a,
1769            double[] b,
1770            double[] c,
1771            double[] d,
1772            int n,
1773            ref double[] x)
1774        {
1775            int k = 0;
1776            double t = 0;
1777
1778            a = (double[])a.Clone();
1779            b = (double[])b.Clone();
1780            c = (double[])c.Clone();
1781            d = (double[])d.Clone();
1782
1783            x = new double[n-1+1];
1784            a[0] = 0;
1785            c[n-1] = 0;
1786            for(k=1; k<=n-1; k++)
1787            {
1788                t = a[k]/b[k-1];
1789                b[k] = b[k]-t*c[k-1];
1790                d[k] = d[k]-t*d[k-1];
1791            }
1792            x[n-1] = d[n-1]/b[n-1];
1793            for(k=n-2; k>=0; k--)
1794            {
1795                x[k] = (d[k]-c[k]*x[k+1])/b[k];
1796            }
1797        }
1798
1799
1800        /*************************************************************************
1801        Internal subroutine. Three-point differentiation
1802        *************************************************************************/
1803        private static double diffthreepoint(double t,
1804            double x0,
1805            double f0,
1806            double x1,
1807            double f1,
1808            double x2,
1809            double f2)
1810        {
1811            double result = 0;
1812            double a = 0;
1813            double b = 0;
1814
1815            t = t-x0;
1816            x1 = x1-x0;
1817            x2 = x2-x0;
1818            a = (f2-f0-x2/x1*(f1-f0))/(AP.Math.Sqr(x2)-x1*x2);
1819            b = (f1-f0-a*AP.Math.Sqr(x1))/x1;
1820            result = 2*a*t+b;
1821            return result;
1822        }
1823    }
1824}
Note: See TracBrowser for help on using the repository browser.