IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Sep 22, 2004, 6:56:41 PM (22 years ago)
Author:
gusciora
Message:

Implemented the generation spline functions.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/dataManip/psMinimize.c

    r1848 r1859  
    99 *  @author George Gusciora, MHPCC
    1010 *
    11  *  @version $Revision: 1.44 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2004-09-22 02:42:49 $
     11 *  @version $Revision: 1.45 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2004-09-23 04:56:41 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    4343#include "psMinimize.h"
    4444#include "psMatrix.h"
    45 #include "psTrace.h"
     45#include "psConstants.h"
    4646
    4747/*****************************************************************************/
    4848/* DEFINE STATEMENTS                                                         */
    4949/*****************************************************************************/
    50 
    51 #define MAX_LMM_ITERATIONS 100
    52 #define MAX_MINIMIZE_ITERATIONS 100
    5350
    5451/** Preprocessor macro to generate error on a NULL 1DPolynomial */
    … …  
    9693/* GLOBAL VARIABLES                                                          */
    9794/*****************************************************************************/
    98 
    9995static psMinimizeChi2PowellFunc Chi2PowellFunc = NULL;
    10096//static psMinimizePowellFunc PowellFunc = NULL;
    … …  
    116112input parameter "x" between 0 and input parameter polyOrder.  The result is
    117113returned as a psVector sums.
     114 
     115XXX: change name
    118116 *****************************************************************************/
    119 void p_psBuildSums1D(double x, int polyOrder, psVector* sums)
     117void p_psBuildSums1D(double x,
     118                     int polyOrder,
     119                     psVector* sums)
    120120{
    121121    int i = 0;
    … …  
    148148
    149149/*****************************************************************************
    150  
    151150CalculateSecondDerivs(): Given a set of x/y vectors corresponding to a
    152 tabulated function at n points, this routine calcualtes the second
     151tabulated function at n points, this routine calculates the second
    153152derivatives of the interpolating cubic splines at those n points.
    154153 
    155154The first and second derivatives at the endpoints, undefined in the SDR, are
    156 here defined to be 0.0.
    157  
    158 XXX: This algorithm is very similar to that in Numerical Recipes.
     155here defined to be 0.0.  They can be modified via ypo and yp1.
     156 
     157XXX: This algorithm is derived from the Numerical Recipes.
    159158 *****************************************************************************/
    160159float *CalculateSecondDerivs(const psVector* restrict x,        ///< Ordinates (or NULL to just use the indices)
    … …  
    174173    float *X = (float *) & (x->data.F32[0]);
    175174    float *Y = (float *) & (y->data.F32[0]);
     175    float qn;
    176176
    177177    if (x == NULL) {
    … …  
    183183    }
    184184
    185     // XXX: The first and second derivatives at the endpoints, undefined in
    186     // the SDR, are here defined to be 0.0.
    187     u[0]= 0.0;
    188     u[n-1]= 0.0;
    189     derivs2[0] = 0.0;
    190     derivs2[n-1] = 0.0;
     185    // XXX: The second derivatives at the endpoints, undefined in the SDR,
     186    // are set in psConstants.h: LEFT_SPLINE_DERIV, RIGHT_SPLINE_DERIV.
     187    derivs2[0] = -0.5;
     188    u[0]= (3.0/(X[1]-X[0])) * ((Y[1]-Y[0])/(X[1]-X[0]) - LEFT_SPLINE_DERIV);
    191189
    192190    for (i=1;i<=(n-2);i++) {
    … …  
    194192        p = sig * derivs2[i-1] + 2.0;
    195193        derivs2[i] = (sig - 1.0) / p;
    196         u[i] = (Y[i+1] - Y[i])/(X[i+1]-X[i]) - (Y[i]-Y[i-1])/(X[i]-X[i-1]);
    197         u[i] = (6.0 * u[i] / (X[i+1] - X[i-1]) - sig * u[i-1]) / p;
     194        u[i] = ((Y[i+1] - Y[i])/(X[i+1]-X[i])) - ((Y[i]-Y[i-1])/(X[i]-X[i-1]));
     195        u[i] = ((6.0 * u[i] / (X[i+1] - X[i-1])) - (sig * u[i-1])) / p;
    198196
    199197        psTrace(".psLib.dataManip.CalculateSecondDerivs", 6,
    … …  
    205203    }
    206204
     205    qn = 0.5;
     206    u[n-1] = (3.0/(X[n-1]-X[n-2])) * (RIGHT_SPLINE_DERIV - (Y[n-1]-Y[n-2])/(X[n-1]-X[n-2]));
     207    derivs2[n-1] = (u[n-1] - (qn * u[n-2])) / ((qn * derivs2[n-2]) + 1.0);
     208
    207209    for (k=(n-2);k>=0;k--) {
    208210        derivs2[k] = derivs2[k] * derivs2[k+1] + u[k];
     211
    209212        psTrace(".psLib.dataManip.CalculateSecondDerivs", 6,
    210213                "derivs2[%d] is %f\n", k, derivs2[k]);
    211214    }
     215
    212216    if (mustFreeX == true) {
    213217        psFree(X);
    … …  
    224228/* FUNCTION IMPLEMENTATION - PUBLIC                                          */
    225229/*****************************************************************************/
    226 
    227 
    228230
    229231/*****************************************************************************
    … …  
    238240 (3)         ((((A*A*A)-A) * mySpline->p_psDeriv2[0]) * H^2)/6.0 +
    239241 (4)         ((((B*B*B)-B) * mySpline->p_psDeriv2[1]) * H^2)/6.0
    240 Where
    241  H = x[1]-x[0], A = (x[1]-x)/H, B = (x-x[0])/H
     242Where:
     243 H = x[1]-x[0]
     244 A = (x[1]-x)/H
     245 B = (x-x[0])/H
    242246The bulk of the code in this routine is the expansion of the above equation
    243 into a polynomial in terms of x.  This gets pretty complicated.
     247into a polynomial in terms of x, and then saving the coefficients of the
     248powers of x in the spline polynomials.  This gets pretty complicated.
    244249 
    245250XXX: usage of yErr is not specified in IfA documentation.
    … …  
    262267    }
    263268    mySpline->p_psDeriv2 = CalculateSecondDerivs(x, y);
     269    for (i=0;i<y->n;i++)
     270        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
     271                "Second deriv[%d] is %f\n", i, mySpline->p_psDeriv2[i]);
    264272
    265273    for (i=0;i<numSplines;i++) {
    … …  
    271279        //
    272280        // From (1)
    273         (mySpline->spline[i])->coeff[0] = y->data.F32[i] * (x->data.F32[i+1]/H);
     281        (mySpline->spline[i])->coeff[0] = (y->data.F32[i] * x->data.F32[i+1]/H);
    274282        // From (2)
    275         ((mySpline->spline[i])->coeff[0])-= (y->data.F32[i+1] * x->data.F32[i]/H);
     283        ((mySpline->spline[i])->coeff[0])-= ((y->data.F32[i+1] * x->data.F32[i])/H);
    276284        // From (3)
    277285        tmp = (x->data.F32[i+1] * x->data.F32[i+1] * x->data.F32[i+1]) / (H * H * H);
    … …  
    289297        //
    290298        // From (1)
    291         (mySpline->spline[i])->coeff[1] = - ((mySpline->p_psDeriv2)[i]) / H;
     299        (mySpline->spline[i])->coeff[1] = -(y->data.F32[i]) / H;
    292300        // From (2)
    293         (mySpline->spline[i])->coeff[1]-= ((mySpline->p_psDeriv2)[i+1]) / H;
     301        ((mySpline->spline[i])->coeff[1])+= (y->data.F32[i+1] / H);
    294302        // From (3)
    295         tmp = - (x->data.F32[i+1] * x->data.F32[i+1]) / (H * H * H);
     303        tmp = -3.0 * (x->data.F32[i+1] * x->data.F32[i+1]) / (H * H * H);
    296304        tmp+= (1.0 / H);
    297305        tmp*= ((mySpline->p_psDeriv2)[i]) * H * H / 6.0;
    298         (mySpline->spline[i])->coeff[1]+= tmp;
     306        ((mySpline->spline[i])->coeff[1])+= tmp;
    299307        // From (4)
    300         tmp = (x->data.F32[i] + x->data.F32[i]) / (H * H * H);
     308        tmp = 3.0 * (x->data.F32[i] * x->data.F32[i]) / (H * H * H);
    301309        tmp-= (1.0 / H);
    302310        tmp*= ((mySpline->p_psDeriv2)[i+1]) * H * H / 6.0;
    303         (mySpline->spline[i])->coeff[1]+= tmp;
     311        ((mySpline->spline[i])->coeff[1])+= tmp;
    304312
    305313        //
    … …  
    307315        //
    308316        // From (3)
    309         (mySpline->spline[i])->coeff[2] = ((mySpline->p_psDeriv2)[i]) * 3.0 * x->data.F32[i+1] / 6.0;
     317        (mySpline->spline[i])->coeff[2] = ((mySpline->p_psDeriv2)[i]) * 3.0 * x->data.F32[i+1] / (6.0 * H);
    310318        // From (4)
    311         (mySpline->spline[i])->coeff[2]-= ((mySpline->p_psDeriv2)[i+1]) * 3.0 * x->data.F32[i] / 6.0;
     319        ((mySpline->spline[i])->coeff[2])-= (((mySpline->p_psDeriv2)[i+1]) * 3.0 * x->data.F32[i] / (6.0 * H));
    312320
    313321        //
    … …  
    317325        (mySpline->spline[i])->coeff[3] = -((mySpline->p_psDeriv2)[i]) / (6.0 * H);
    318326        // From (4)
    319         (mySpline->spline[i])->coeff[3]+=  ((mySpline->p_psDeriv2)[i+1]) / (6.0 * H);
     327        ((mySpline->spline[i])->coeff[3])+=  ((mySpline->p_psDeriv2)[i+1]) / (6.0 * H);
    320328
    321329        psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
    … …  
    334342    return(mySpline);
    335343}
     344
     345/******************************************************************************
     346p_psNRSpline1DEval(): This routine does NR-style evaluation of cubic splines.
     347It takes advantage of the 2nd derivatives of the cubic splines, which are
     348stored in the psSPline1D data structure, and computes the interpolated value
     349directly, without computing (or using) the interpolating cubic spline
     350polynomial.
     351 
     352This routine is here mostly for a sanity check on the psLib function
     353evalSpline() which computes the interpolated value based on the cubic spline
     354polynomials which are stored in psSpline1D.
     355 *****************************************************************************/
     356float p_psNRSpline1DEval(psSpline1D *spline,
     357                         const psVector* restrict x,
     358                         const psVector* restrict y,
     359                         float X)
     360{
     361    int n;
     362    int klo;
     363    int khi;
     364    float H;
     365    float A;
     366    float B;
     367    float C;
     368    float D;
     369    float Y;
     370
     371    n = spline->n;
     372    klo = p_psVectorBinDisectF32(spline->domains, (spline->n)+1, X);
     373    khi = klo + 1;
     374    H = (spline->domains)[khi] - (spline->domains)[klo];
     375    A = ((spline->domains)[khi] - X) / H;
     376    B = (X - (spline->domains)[klo]) / H;
     377    C = ((A*A*A)-A) * (H*H/6.0);
     378    D = ((B*B*B)-B) * (H*H/6.0);
     379
     380    Y = (A * y->data.F32[klo]) +
     381        (B * y->data.F32[khi]) +
     382        (C * (spline->p_psDeriv2)[klo]) +
     383        (D * (spline->p_psDeriv2)[khi]);
     384
     385    return(Y);
     386}
     387
    336388
    337389/******************************************************************************
    … …  
    687739}
    688740
    689 
    690 #define STEP_SIZE 0.10
    691741/******************************************************************************
    692     This routine takes as input an arbitrary function, and the parameter to
    693     vary, and the line along which it must vary.  This function produces as
    694     output a bracket [a, b, c] such that
     742p_psDetermineBracket():  This routine takes as input an arbitrary function,
     743and the parameter to vary, and the line along which it must vary.  This
     744function produces as output a bracket [a, b, c] such that
     745 
    695746        f(param + b * line) is less than f(param + a * line) and
    696747                                         f(param + c * line).
    … …  
    700751    smaller/larger than b.  Repeat this process until a local minimum is
    701752    found.
     753 
     754 
    702755 *****************************************************************************/
    703756psVector *p_psDetermineBracket(psVector *params,
    … …  
    719772    float new_cDir = 0.0;
    720773    psVector *bracket = psVectorAlloc(3, PS_TYPE_F32);
    721     float stepSize = STEP_SIZE;
     774    float stepSize = DETERMINE_BRACKET_STEP_SIZE;
    722775    psVector *tmp = NULL;
    723776    int i = 0;
    … …  
    869922    return(NULL);
    870923}
    871 
    872 
    873 
    874 
    875 
    876 
    877 
    878 
    879 
    880 
    881 
    882924
    883925/******************************************************************************
    … …  
    10611103
    10621104
    1063 
    1064 
    1065 
    1066 
    1067 
    1068 
    1069 
    1070 
    1071 
    1072 
    1073 
    1074 
    1075 
    10761105/******************************************************************************
    10771106This routine must minimize a possibly multi-dimensional function.  The
Note: See TracChangeset for help on using the changeset viewer.