Main Page   Data Structures   File List   Data Fields   Globals
spline.c
Go to the documentation of this file.

00001 /*******************************************************************
00002   RTN SPLINE: Fits cubic smoothing spline to time series
00003
00004   Derived from IMSL routines by Edward R Cook, Tree Ring Laboratory,
00005   Lamont-Doherty Earth Observatory, Palisades, New York, USA
00006
00007   Four routines combined into one by
00008   Richard L Holmes, University of Arizona, Tucson, Arizona, USA
00009   Modified copyright (C) 10 AUG 1998
00010
00011 ********************************************************************/
00012
00013 #include <stdio.h>
00014 #include <stdlib.h>
00015 #include <math.h>
00016
00017 #define PI 3.1415926535897935
00018
00019
00020
00021 /* DYNAMICALLY ALLOCATE A 2-D ARRAY */
00022 /* Assumption:  nrows*ncols*element_size, rounded up to a multiple   */
00023 /* of sizeof(long double), must fit in a long type.  If not, then    */
00024 /* the "i += ..." step could overflow.                               */
00025
00026 void **MATRIX(int nrows, int ncols, int first_row, int first_col,int element_size)
00027 {
00028     void **p;
00029     int alignment;
00030     long i;
00031
00032     if(nrows < 1 || ncols < 1) return(NULL);
00033     i = nrows*sizeof(void *);
00034     /* align the addr of the data to be a multiple of sizeof(long double) */
00035     alignment = i % sizeof(long double);
00036     if(alignment != 0) alignment = sizeof(long double) - alignment;
00037     i += nrows*ncols*element_size+alignment;
00038     if((p = (void **)malloc((size_t)i)) != NULL)
00039     {
00040         /* compute the address of matrix[first_row][0] */
00041         p[0] = (char *)(p+nrows)+alignment-first_col*element_size;
00042         for(i = 1; i < nrows; i++)
00043             /* compute the address of matrix[first_row+i][0] */
00044             p[i] = (char *)(p[i-1])+ncols*element_size;
00045         /* compute the address of matrix[0][0] */
00046         p -= first_row;
00047     }
00048     return(p);
00049 }
00050
00051
00052
00053 /* This function is called from SPLINE when :  */
00054 /* 1. Series is too short to compute spline    */
00055 /* 2. Matrix A is not positive definite        */
00056
00057 void errorAction(int N, double *Y, float *ZF)
00058 {
00059     int k;
00060     double ZN;
00061
00062     ZN = 0.0;
00063     for(k = 1; k <= N; k++)
00064         ZN = ZN + Y[k];
00065
00066     if (N > 0)
00067     {
00068         ZN = ZN/(float)N;
00069         for(k = 1; k <= N; k++)
00070             ZF[k-1] = ZN;
00071     }
00072
00073     return;
00074 }
00075
00076
00077
00078 /* Function  SPLINE: Fits cubic smoothing spline to time series               */
00079 /* Arguments:                                                                 */
00080 /*                                                                            */
00081 /* N:   Number of values in time series                                       */
00082 /* Z:   Time series array to be modeled with spline                           */
00083 /* ZF:  Computed cubic spline function array fit to time series               */
00084 /* ZSP: Length (rigidity) of spline to be used to model series                */
00085 /* ZPV: Portion of variance at wavelength ZSP contained in spline (0<ZPV<1)   */
00086 /*                                                                            */
00087 /* Arguments Z, ZF, ZSP and ZPV are single precision;                         */
00088 /* computation is done entirely in double-precision arithmetic                */
00089
00090 void SPLINE(int N, float *Z, float *ZF, float ZSP, float ZPV)
00091 {
00092     int i, j, k, l, m;
00093     int NC, NC1, NCP1, IMNCP1, I1, I2, JM1, IW, KL, N1, K1;
00094     double RSP, VPV, PD, RN, D1, D2, SUM;
00095     double **A, *F, *Y, C1[5], C2[4];
00096
00097     C1[0] =  0.0;
00098     C1[1] =  1.0;
00099     C1[2] = -4.0;
00100     C1[3] =  6.0;
00101     C1[4] = -2.0;
00102
00103     C2[0] =  0.0;
00104     C2[1] =  0.0;
00105     C2[2] =  0.33333333333333;
00106     C2[3] =  1.33333333333333;
00107
00108     /* Allocate arrays to store intermeediate results */
00109     A = (double **)MATRIX(N+1, 5, 0, 0, sizeof(double));
00110     F = (double *)malloc((N+1)*sizeof(double));
00111     Y = (double *)malloc((N+1)*sizeof(double));
00112     if (A == NULL || F == NULL || Y == NULL)
00113     {
00114         printf("\nSPLINE >> Unable to allocate memory\n");
00115         return;
00116     }
00117
00118     /* Check whether series is too short to compute spline */
00119     if (N < 4)
00120     {
00121         errorAction(N, Y, ZF);
00122         return;
00123     }
00124
00125     /* Copy time series into double precision array */
00126     for(j = 1; j <= N; j++)
00127         Y[j] = (double)Z[j-1];
00128
00129     /* Compute Lagrange multiplier, which defines frequency response of spline */
00130     RSP = (double)ZSP;
00131     VPV = (double)ZPV;
00132     PD = ((1.0/(1.0-VPV)-1.0)*6.0* pow((cos(PI*2.0/RSP)-1.0),2.0))/(cos(PI*2.0/RSP)+2.0);
00133     for(i = 1; i <= N-2; i++)
00134         for(j = 1; j <= 3; j++)
00135         {
00136             A[i][j] = C1[j] + PD * C2[j];
00137             A[i][4] = Y[i] + C1[4] * Y[i+1] + Y[i+2];
00138         }
00139
00140     A[1][1] = C2[1];
00141     A[1][2] = C2[1];
00142     A[2][1] = C2[1];
00143     NC = 2;
00144
00145     /* Begin LUDAPB */
00146     RN = 1.0/((double)(N-2) * 16.0);
00147     D1 = 1.0;
00148     D2 = 0.0;
00149     NCP1 = NC + 1;
00150
00151     /* Initialize zero elements */
00152     if (NC != 0)
00153     {
00154         for(i = 1; i <= NC; i++)
00155             for(j = i; j <= NC; j++)
00156             {
00157                 k = NCP1 - j;
00158                 A[i][k] = 0.0;
00159             }
00160     }
00161
00162     /* i: row index of element being computed */
00163     /* j: column index of element being computed */
00164     /* l: row index of previously computed vector being used to compute inner product */
00165     /* m: column index */
00166     for(i = 1; i <= N-2; i++)
00167     {
00168         IMNCP1 = i - NCP1;
00169         I1 = (1 < 1 - IMNCP1? 1 - IMNCP1: 1);
00170         for(j = I1; j <= NCP1; j++)
00171         {
00172             l = IMNCP1 + j;
00173             I2 = NCP1 - j;
00174             SUM = A[i][j];
00175             JM1 = j - 1;
00176
00177             if (JM1 > 0)
00178             {
00179                 for(k = 1; k <= JM1; k++)
00180                 {
00181                     m = I2 + k;
00182                     SUM = SUM - (A[i][k]*A[l][m]);
00183                 }
00184             }
00185
00186             /* Matrix not positive definite */
00187             if (j == NCP1)
00188             {
00189                 if (A[i][j]+SUM*RN <= A[i][j])
00190                 {
00191                     printf("\nSPLINE >> Matrix not positive definite\n");
00192                     errorAction(N, Y, ZF);
00193                     return;
00194                 }
00195
00196                 A[i][j] = 1.0/sqrt(SUM);
00197
00198                 /* Update determinant */
00199                 D1 = D1*SUM;
00200                 while (fabs(D1) > 1.0)
00201                 {
00202                     D1 = D1*0.0625;
00203                     D2 = D2+4.0;
00204                 }
00205
00206                 while (fabs(D1) <= 0.0625)
00207                 {
00208                     D1 = D1*16.0;
00209                     D2 = D2-4.0;
00210                 }
00211                 continue;
00212             }
00213             A[i][j] = SUM*A[l][NCP1];
00214         }
00215     }
00216     /* End LUDAPB */
00217
00218     /* Begin LUELPB */
00219     /* Solution LY = B */
00220     NC1 = NC + 1;
00221     IW = 0;
00222     l = 0;
00223
00224     for(i = 1; i <= N-2; i++)
00225     {
00226         SUM = A[i][4];
00227         if (NC > 0) {
00228             if (IW != 0) {
00229                 l = l + 1;
00230                 if (l > NC)
00231                 {
00232                     l = NC;
00233                 }
00234                 k = NC1 - l;
00235                 KL = i - l;
00236
00237                 for(j = k; j <= NC; j++)
00238                 {
00239                     SUM = SUM - (A[KL][4]*A[i][j]);
00240                     KL = KL + 1;
00241                 }
00242             } else if (SUM != 0.0)
00243             {
00244                 IW = 1;
00245             }
00246         }
00247         A[i][4] = SUM*A[i][NC1];
00248     }
00249
00250     /* Solution UX = Y */
00251     A[N-2][4] = A[N-2][4]*A[N-2][NC1];
00252     N1 = N-2+1;
00253
00254     for(i = 2; i <= N-2; i++)
00255     {
00256         k = N1 - i;
00257         SUM = A[k][4];
00258         if (NC > 0)
00259         {
00260             KL = k +1;
00261             K1 = (N-2 < k+NC? N-2: k+NC);
00262             l = 1;
00263             for(j = KL; j <= K1; j++)
00264             {
00265                 SUM = SUM - A[j][4]*A[j][NC1-l];
00266                 l = l + 1;
00267             }
00268         }
00269         A[k][4] = SUM*A[k][NC1];
00270     }
00271     /* End LUELPB */
00272
00273     /* Calculate Spline Curve */
00274     for(i = 3; i <= N-2; i++)
00275         F[i] = A[i-2][4]+C1[4]*A[i-1][4]+A[i][4];
00276
00277     F[1] = A[1][4];
00278     F[2] = C1[4]*A[1][4]+A[2][4];
00279     F[N-1] = A[N-2-1][4]+C1[4]*A[N-2][4];
00280     F[N] = A[N-2][4];
00281
00282     for(i = 1; i <= N; i++)
00283         ZF[i-1] = Y[i] - F[i];
00284
00285     return;
00286 }
00287
00288
00289
00290
00291

Generated on Wed Apr 9 08:56:14 2003 for TREES by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002