IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 1860


Ignore:
Timestamp:
Sep 22, 2004, 7:43:25 PM (22 years ago)
Author:
gusciora
Message:

Misc comments and XXX notation.

Location:
trunk/psLib/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

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

    r1859 r1860  
    99 *  @author George Gusciora, MHPCC
    1010 *
    11  *  @version $Revision: 1.45 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2004-09-23 04:56:41 $
     11 *  @version $Revision: 1.46 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2004-09-23 05:43:25 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    1515 *
    1616 *  XXX: must follow coding name standards on local functions.
     17 *
     18 *  XXX: Section 4.5.1.1 (predefined functions for Gauss minimization via
     19 *       LMM) is not addressed here.  We are waiting for subsequent SDRs
     20 *       which will redfeine the LMM functions.
    1721 *
    1822 */
    … …  
    3842#include "psImage.h"
    3943#include "psTrace.h"
     44#include "psLogMsg.h"
    4045#include "psError.h"
    4146#include "psAbort.h"
    … …  
    9499/*****************************************************************************/
    95100static psMinimizeChi2PowellFunc Chi2PowellFunc = NULL;
    96 //static psMinimizePowellFunc PowellFunc = NULL;
    97101static psVector *myValue;
    98102static psVector *myError;
    … …  
    121125    int i = 0;
    122126    double xSum = 0.0;
     127    if (sums == NULL) {
     128        sums = psVectorAlloc(polyOrder, PS_TYPE_F64);
     129    }
    123130
    124131    xSum = 1.0;
    … …  
    139146
    140147    tmp = psVectorAlloc(x, PS_TYPE_F32);
    141 
    142148    for (i = 0; i < x; i++) {
    143149        tmp->data.F32[i] = (((float)2 * i) / ((float)x)) - 1.0;
    … …  
    224230}
    225231
    226 
    227 /*****************************************************************************/
    228 /* FUNCTION IMPLEMENTATION - PUBLIC                                          */
    229 /*****************************************************************************/
    230 
    231 /*****************************************************************************
    232 psVectorFitSpline1D(): given a psSpline1D data structure and a set of x/y
    233 vectors, this routine generates the cublic splines which satisfy those data
    234 points.
    235  
    236 The formula for calculating the spline polynomials is derived from Numerical
    237 Recipes in C.  The basic idea is that the polynomial is
    238  (1)     y = (A * y[0]) +
    239  (2)         (B * y[1]) +
    240  (3)         ((((A*A*A)-A) * mySpline->p_psDeriv2[0]) * H^2)/6.0 +
    241  (4)         ((((B*B*B)-B) * mySpline->p_psDeriv2[1]) * H^2)/6.0
    242 Where:
    243  H = x[1]-x[0]
    244  A = (x[1]-x)/H
    245  B = (x-x[0])/H
    246 The bulk of the code in this routine is the expansion of the above equation
    247 into a polynomial in terms of x, and then saving the coefficients of the
    248 powers of x in the spline polynomials.  This gets pretty complicated.
    249  
    250 XXX: usage of yErr is not specified in IfA documentation.
    251  *****************************************************************************/
    252 
    253 psSpline1D *psVectorFitSpline1D(psSpline1D *mySpline,              ///< The spline which will be generated.
    254                                 const psVector* restrict x,        ///< Ordinates (or NULL to just use the indices)
    255                                 const psVector* restrict y,        ///< Coordinates
    256                                 const psVector* restrict yErr)     ///< Errors in coordinates, or NULL
    257 {
    258     psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
    259             "---- psVectorFitSpline1D() begin ----\n");
    260     int numSplines = (y->n)-1;
    261     float tmp;
    262     float H;
    263     int i;
    264 
    265     if (mySpline == NULL) {
    266         //XXX psErrorMsg()
    267     }
    268     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]);
    272 
    273     for (i=0;i<numSplines;i++) {
    274         H = x->data.F32[i+1] - x->data.F32[i];
    275         psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
    276                 "X data (%f - %f) (%f)\n", x->data.F32[i], x->data.F32[i+1], H);
    277         //
    278         // ******** Calulate 0-order term ********
    279         //
    280         // From (1)
    281         (mySpline->spline[i])->coeff[0] = (y->data.F32[i] * x->data.F32[i+1]/H);
    282         // From (2)
    283         ((mySpline->spline[i])->coeff[0])-= ((y->data.F32[i+1] * x->data.F32[i])/H);
    284         // From (3)
    285         tmp = (x->data.F32[i+1] * x->data.F32[i+1] * x->data.F32[i+1]) / (H * H * H);
    286         tmp-= (x->data.F32[i+1] / H);
    287         tmp*= (mySpline->p_psDeriv2)[i] * H * H / 6.0;
    288         ((mySpline->spline[i])->coeff[0])+= tmp;
    289         // From (4)
    290         tmp = -(x->data.F32[i] * x->data.F32[i] * x->data.F32[i]) / (H * H * H);
    291         tmp+= (x->data.F32[i] / H);
    292         tmp*= (mySpline->p_psDeriv2)[i+1] * H * H / 6.0;
    293         ((mySpline->spline[i])->coeff[0])+= tmp;
    294 
    295         //
    296         // ******** Calulate 1-order term ********
    297         //
    298         // From (1)
    299         (mySpline->spline[i])->coeff[1] = -(y->data.F32[i]) / H;
    300         // From (2)
    301         ((mySpline->spline[i])->coeff[1])+= (y->data.F32[i+1] / H);
    302         // From (3)
    303         tmp = -3.0 * (x->data.F32[i+1] * x->data.F32[i+1]) / (H * H * H);
    304         tmp+= (1.0 / H);
    305         tmp*= ((mySpline->p_psDeriv2)[i]) * H * H / 6.0;
    306         ((mySpline->spline[i])->coeff[1])+= tmp;
    307         // From (4)
    308         tmp = 3.0 * (x->data.F32[i] * x->data.F32[i]) / (H * H * H);
    309         tmp-= (1.0 / H);
    310         tmp*= ((mySpline->p_psDeriv2)[i+1]) * H * H / 6.0;
    311         ((mySpline->spline[i])->coeff[1])+= tmp;
    312 
    313         //
    314         // ******** Calulate 2-order term ********
    315         //
    316         // From (3)
    317         (mySpline->spline[i])->coeff[2] = ((mySpline->p_psDeriv2)[i]) * 3.0 * x->data.F32[i+1] / (6.0 * H);
    318         // From (4)
    319         ((mySpline->spline[i])->coeff[2])-= (((mySpline->p_psDeriv2)[i+1]) * 3.0 * x->data.F32[i] / (6.0 * H));
    320 
    321         //
    322         // ******** Calulate 3-order term ********
    323         //
    324         // From (3)
    325         (mySpline->spline[i])->coeff[3] = -((mySpline->p_psDeriv2)[i]) / (6.0 * H);
    326         // From (4)
    327         ((mySpline->spline[i])->coeff[3])+=  ((mySpline->p_psDeriv2)[i+1]) / (6.0 * H);
    328 
    329         psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
    330                 "(mySpline->spline[%d])->coeff[0] is %f\n", i, (mySpline->spline[i])->coeff[0]);
    331         psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
    332                 "(mySpline->spline[%d])->coeff[1] is %f\n", i, (mySpline->spline[i])->coeff[1]);
    333         psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
    334                 "(mySpline->spline[%d])->coeff[2] is %f\n", i, (mySpline->spline[i])->coeff[2]);
    335         psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
    336                 "(mySpline->spline[%d])->coeff[3] is %f\n", i, (mySpline->spline[i])->coeff[3]);
    337 
    338     }
    339 
    340     psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
    341             "---- psVectorFitSpline1D() end ----\n");
    342     return(mySpline);
    343 }
    344232
    345233/******************************************************************************
    … …  
    387275
    388276
     277/*****************************************************************************/
     278/* FUNCTION IMPLEMENTATION - PUBLIC                                          */
     279/*****************************************************************************/
     280
     281/*****************************************************************************
     282psVectorFitSpline1D(): given a psSpline1D data structure and a set of x/y
     283vectors, this routine generates the linear or cublic splines which satisfy
     284those data points.
     285 
     286The formula for calculating the spline polynomials is derived from Numerical
     287Recipes in C.  The basic idea is that the polynomial is
     288 (1)     y = (A * y[0]) +
     289 (2)         (B * y[1]) +
     290 (3)         ((((A*A*A)-A) * mySpline->p_psDeriv2[0]) * H^2)/6.0 +
     291 (4)         ((((B*B*B)-B) * mySpline->p_psDeriv2[1]) * H^2)/6.0
     292Where:
     293 H = x[1]-x[0]
     294 A = (x[1]-x)/H
     295 B = (x-x[0])/H
     296The bulk of the code in this routine is the expansion of the above equation
     297into a polynomial in terms of x, and then saving the coefficients of the
     298powers of x in the spline polynomials.  This gets pretty complicated.
     299 
     300XXX: usage of yErr is not specified in IfA documentation.
     301XXX: Must do: if yErr==NULL, set all errors equal.
     302 
     303XXX: Is the x argument redundant?  What do we do if the x argument is
     304 supplied, but does not equal the domains specified in mySpline?
     305 
     306XXX: can psSpline be NULL?
     307 
     308XXX: Implemented in F32 only, must add F64.
     309 *****************************************************************************/
     310psSpline1D *psVectorFitSpline1D(psSpline1D *mySpline,              ///< The spline which will be generated.
     311                                const psVector* restrict x,        ///< Ordinates (or NULL to just use the indices)
     312                                const psVector* restrict y,        ///< Coordinates
     313                                const psVector* restrict yErr)     ///< Errors in coordinates, or NULL
     314{
     315    psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
     316            "---- psVectorFitSpline1D() begin ----\n");
     317    int numSplines = (y->n)-1;
     318    float tmp;
     319    float H;
     320    int i;
     321    float slope;
     322
     323    if (mySpline == NULL) {
     324        //XXX psErrorMsg()
     325    }
     326    if (y == NULL) {
     327        //XXX psErrorMsg()
     328    }
     329
     330    if (y->n != (1 + mySpline->n)) {
     331        psLogMsg(__func__, PS_LOG_WARN,
     332                 "data size / spline size mismatch (%d %d)\n",
     333                 y->n, mySpline->n);
     334        // XXX: psErrorMsg()
     335    }
     336
     337    // If these are linear splines, which means their polynomials will have
     338    // two coefficients, then we do the simple calculation.
     339    if (2 == (mySpline->spline[0])->n) {
     340        for (i=0;i<mySpline->n;i++) {
     341            slope = (y->data.F32[i+1] - y->data.F32[i]) /
     342                    (mySpline->domains[i+1] - mySpline->domains[i]);
     343            (mySpline->spline[i])->coeff[0] = y->data.F32[i] -
     344                                              (slope * mySpline->domains[i]);
     345
     346            (mySpline->spline[i])->coeff[1] = slope;
     347            psTrace(".psLib.dataManip.psFunctions.psSpline1DGen", 4,
     348                    "---- mySpline %d coeffs are (%f, %f)\n", i,
     349                    (mySpline->spline[i])->coeff[0],
     350                    (mySpline->spline[i])->coeff[1]);
     351        }
     352        psTrace(".psLib.dataManip.psFunctions.psSpline1DGen", 4,
     353                "---- Exiting psSpline1DGen()\n");
     354        return((psSpline1D *) mySpline);
     355    }
     356
     357    // Check if these are cubic splines (n==4).  If not, psError.
     358    if (4 != (mySpline->spline[0])->n) {
     359        psLogMsg(__func__, PS_LOG_WARN,
     360                 "Don't know how to generate %d-order splines.",
     361                 (mySpline->spline[0])->n-1);
     362        // XXX: psErrorMsg()
     363        return(NULL);
     364    }
     365
     366    // If we get here, then we know these are cubic splines.  We first
     367    // generate the second derivatives at each data point.
     368    mySpline->p_psDeriv2 = CalculateSecondDerivs(x, y);
     369    for (i=0;i<y->n;i++)
     370        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
     371                "Second deriv[%d] is %f\n", i, mySpline->p_psDeriv2[i]);
     372
     373    // We generate the coefficients of the spline polynomials.  I can't
     374    // concisely explain how this code works.  See above function comments
     375    // and Numerical Recipes in C.
     376    for (i=0;i<numSplines;i++) {
     377        H = x->data.F32[i+1] - x->data.F32[i];
     378        psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
     379                "X data (%f - %f) (%f)\n", x->data.F32[i], x->data.F32[i+1], H);
     380        //
     381        // ******** Calculate 0-order term ********
     382        //
     383        // From (1)
     384        (mySpline->spline[i])->coeff[0] = (y->data.F32[i] * x->data.F32[i+1]/H);
     385        // From (2)
     386        ((mySpline->spline[i])->coeff[0])-= ((y->data.F32[i+1] * x->data.F32[i])/H);
     387        // From (3)
     388        tmp = (x->data.F32[i+1] * x->data.F32[i+1] * x->data.F32[i+1]) / (H * H * H);
     389        tmp-= (x->data.F32[i+1] / H);
     390        tmp*= (mySpline->p_psDeriv2)[i] * H * H / 6.0;
     391        ((mySpline->spline[i])->coeff[0])+= tmp;
     392        // From (4)
     393        tmp = -(x->data.F32[i] * x->data.F32[i] * x->data.F32[i]) / (H * H * H);
     394        tmp+= (x->data.F32[i] / H);
     395        tmp*= (mySpline->p_psDeriv2)[i+1] * H * H / 6.0;
     396        ((mySpline->spline[i])->coeff[0])+= tmp;
     397
     398        //
     399        // ******** Calculate 1-order term ********
     400        //
     401        // From (1)
     402        (mySpline->spline[i])->coeff[1] = -(y->data.F32[i]) / H;
     403        // From (2)
     404        ((mySpline->spline[i])->coeff[1])+= (y->data.F32[i+1] / H);
     405        // From (3)
     406        tmp = -3.0 * (x->data.F32[i+1] * x->data.F32[i+1]) / (H * H * H);
     407        tmp+= (1.0 / H);
     408        tmp*= ((mySpline->p_psDeriv2)[i]) * H * H / 6.0;
     409        ((mySpline->spline[i])->coeff[1])+= tmp;
     410        // From (4)
     411        tmp = 3.0 * (x->data.F32[i] * x->data.F32[i]) / (H * H * H);
     412        tmp-= (1.0 / H);
     413        tmp*= ((mySpline->p_psDeriv2)[i+1]) * H * H / 6.0;
     414        ((mySpline->spline[i])->coeff[1])+= tmp;
     415
     416        //
     417        // ******** Calculate 2-order term ********
     418        //
     419        // From (3)
     420        (mySpline->spline[i])->coeff[2] = ((mySpline->p_psDeriv2)[i]) * 3.0 * x->data.F32[i+1] / (6.0 * H);
     421        // From (4)
     422        ((mySpline->spline[i])->coeff[2])-= (((mySpline->p_psDeriv2)[i+1]) * 3.0 * x->data.F32[i] / (6.0 * H));
     423
     424        //
     425        // ******** Calculate 3-order term ********
     426        //
     427        // From (3)
     428        (mySpline->spline[i])->coeff[3] = -((mySpline->p_psDeriv2)[i]) / (6.0 * H);
     429        // From (4)
     430        ((mySpline->spline[i])->coeff[3])+=  ((mySpline->p_psDeriv2)[i+1]) / (6.0 * H);
     431
     432        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
     433                "(mySpline->spline[%d])->coeff[0] is %f\n", i, (mySpline->spline[i])->coeff[0]);
     434        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
     435                "(mySpline->spline[%d])->coeff[1] is %f\n", i, (mySpline->spline[i])->coeff[1]);
     436        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
     437                "(mySpline->spline[%d])->coeff[2] is %f\n", i, (mySpline->spline[i])->coeff[2]);
     438        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
     439                "(mySpline->spline[%d])->coeff[3] is %f\n", i, (mySpline->spline[i])->coeff[3]);
     440
     441    }
     442
     443    psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
     444            "---- psVectorFitSpline1D() end ----\n");
     445    return(mySpline);
     446}
     447
    389448/******************************************************************************
    390449psMinimizeLMChi2():  This routine will take an procedure which calculates an
    … …  
    398457 
    399458     A[i][i] = (1 + lambda*paramask[i]) * alpha[i][i];
     459 
     460XXX: This is very different from what is specified in the SDR.  Must
     461coordinate with IfA on new SDR.
     462 
     463XXX: Do vector/image recycles.
     464 
     465XXX: probably yErr will be part of the SDR.
     466 
     467XXX: This must work for both F32 and F64.  F32 is currently implemented.
     468     Note: since the LUD routines are only implemented in F64, then we
     469     will have to convert all F32 input vectors to F64 regardless.  So,
     470     the F64 port might be.
    400471 *****************************************************************************/
    401472bool psMinimizeLMChi2(psMinimization *min,
    … …  
    450521
    451522    min->lastDelta = HUGE;
    452     min->lastDelta = 12345.0;
    453523    min->iter = 0;
    454524
    … …  
    627697polynomial, as well as the error for each data point (yErr).
    628698 
    629 XXX: NOTE: yErr is currently ignored.
     699XXX: yErr is currently ignored.
     700 
     701XXX: must add type F32 (currently F64 only).
     702 
     703XXX: Must do: if x==NULL, use the index vector (???).
     704 
     705XXX: Must do: if yErr==NULL, set all errors equal.
    630706 *****************************************************************************/
    631707psPolynomial1D* psVectorFitPolynomial1D(psPolynomial1D* myPoly,
    … …  
    726802}
    727803
     804
     805/******************************************************************************
     806 *****************************************************************************/
    728807psMinimization *psMinimizationAlloc(int maxIter,
    729808                                    float tol)
    … …  
    747826                                         f(param + c * line).
    748827 
    749     Algorithm: XXX completely ad hoc: start with the user-supplied starting
    750     parameter and call that b.  Calculate a/c as a fractional amount
    751     smaller/larger than b.  Repeat this process until a local minimum is
    752     found.
    753  
    754  
     828Algorithm: XXX completely ad hoc: start with the user-supplied starting
     829parameter and call that b.  Calculate a/c as a fractional amount
     830smaller/larger than b.  Repeat this process until a local minimum is found.
     831 
     832XXX: new algorithm: start at x=0, expand in one direction until the function
     833decreases.  Then you have two points in the bracket.  Keep going until it
     834increases, or x is too large.  If thst does not work, expand in the other
     835direction.
    755836 *****************************************************************************/
    756837psVector *p_psDetermineBracket(psVector *params,
    … …  
    9261007    This routine must minimize a possibly multi-dimensional function
    9271008    along a vector defined by line.
     1009 
     1010XXX: Use a p_psName().
    9281011 *****************************************************************************/
    9291012float psLineMin(psMinimization *min,
    … …  
    11131196XXX: We do not use Brent's method.
    11141197 
    1115 XXX: We do not use the
     1198XXX: The SDR is silent about data types.  F32 is implemented here.
     1199 
     1200XXX: Define constants for default dummy number of iterations and tolerance.
    11161201 *****************************************************************************/
    11171202bool psMinimizePowell(psMinimization *min,

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

    r1859 r1860  
    99 *  @author George Gusciora, MHPCC
    1010 *
    11  *  @version $Revision: 1.45 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2004-09-23 04:56:41 $
     11 *  @version $Revision: 1.46 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2004-09-23 05:43:25 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    1515 *
    1616 *  XXX: must follow coding name standards on local functions.
     17 *
     18 *  XXX: Section 4.5.1.1 (predefined functions for Gauss minimization via
     19 *       LMM) is not addressed here.  We are waiting for subsequent SDRs
     20 *       which will redfeine the LMM functions.
    1721 *
    1822 */
    … …  
    3842#include "psImage.h"
    3943#include "psTrace.h"
     44#include "psLogMsg.h"
    4045#include "psError.h"
    4146#include "psAbort.h"
    … …  
    9499/*****************************************************************************/
    95100static psMinimizeChi2PowellFunc Chi2PowellFunc = NULL;
    96 //static psMinimizePowellFunc PowellFunc = NULL;
    97101static psVector *myValue;
    98102static psVector *myError;
    … …  
    121125    int i = 0;
    122126    double xSum = 0.0;
     127    if (sums == NULL) {
     128        sums = psVectorAlloc(polyOrder, PS_TYPE_F64);
     129    }
    123130
    124131    xSum = 1.0;
    … …  
    139146
    140147    tmp = psVectorAlloc(x, PS_TYPE_F32);
    141 
    142148    for (i = 0; i < x; i++) {
    143149        tmp->data.F32[i] = (((float)2 * i) / ((float)x)) - 1.0;
    … …  
    224230}
    225231
    226 
    227 /*****************************************************************************/
    228 /* FUNCTION IMPLEMENTATION - PUBLIC                                          */
    229 /*****************************************************************************/
    230 
    231 /*****************************************************************************
    232 psVectorFitSpline1D(): given a psSpline1D data structure and a set of x/y
    233 vectors, this routine generates the cublic splines which satisfy those data
    234 points.
    235  
    236 The formula for calculating the spline polynomials is derived from Numerical
    237 Recipes in C.  The basic idea is that the polynomial is
    238  (1)     y = (A * y[0]) +
    239  (2)         (B * y[1]) +
    240  (3)         ((((A*A*A)-A) * mySpline->p_psDeriv2[0]) * H^2)/6.0 +
    241  (4)         ((((B*B*B)-B) * mySpline->p_psDeriv2[1]) * H^2)/6.0
    242 Where:
    243  H = x[1]-x[0]
    244  A = (x[1]-x)/H
    245  B = (x-x[0])/H
    246 The bulk of the code in this routine is the expansion of the above equation
    247 into a polynomial in terms of x, and then saving the coefficients of the
    248 powers of x in the spline polynomials.  This gets pretty complicated.
    249  
    250 XXX: usage of yErr is not specified in IfA documentation.
    251  *****************************************************************************/
    252 
    253 psSpline1D *psVectorFitSpline1D(psSpline1D *mySpline,              ///< The spline which will be generated.
    254                                 const psVector* restrict x,        ///< Ordinates (or NULL to just use the indices)
    255                                 const psVector* restrict y,        ///< Coordinates
    256                                 const psVector* restrict yErr)     ///< Errors in coordinates, or NULL
    257 {
    258     psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
    259             "---- psVectorFitSpline1D() begin ----\n");
    260     int numSplines = (y->n)-1;
    261     float tmp;
    262     float H;
    263     int i;
    264 
    265     if (mySpline == NULL) {
    266         //XXX psErrorMsg()
    267     }
    268     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]);
    272 
    273     for (i=0;i<numSplines;i++) {
    274         H = x->data.F32[i+1] - x->data.F32[i];
    275         psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
    276                 "X data (%f - %f) (%f)\n", x->data.F32[i], x->data.F32[i+1], H);
    277         //
    278         // ******** Calulate 0-order term ********
    279         //
    280         // From (1)
    281         (mySpline->spline[i])->coeff[0] = (y->data.F32[i] * x->data.F32[i+1]/H);
    282         // From (2)
    283         ((mySpline->spline[i])->coeff[0])-= ((y->data.F32[i+1] * x->data.F32[i])/H);
    284         // From (3)
    285         tmp = (x->data.F32[i+1] * x->data.F32[i+1] * x->data.F32[i+1]) / (H * H * H);
    286         tmp-= (x->data.F32[i+1] / H);
    287         tmp*= (mySpline->p_psDeriv2)[i] * H * H / 6.0;
    288         ((mySpline->spline[i])->coeff[0])+= tmp;
    289         // From (4)
    290         tmp = -(x->data.F32[i] * x->data.F32[i] * x->data.F32[i]) / (H * H * H);
    291         tmp+= (x->data.F32[i] / H);
    292         tmp*= (mySpline->p_psDeriv2)[i+1] * H * H / 6.0;
    293         ((mySpline->spline[i])->coeff[0])+= tmp;
    294 
    295         //
    296         // ******** Calulate 1-order term ********
    297         //
    298         // From (1)
    299         (mySpline->spline[i])->coeff[1] = -(y->data.F32[i]) / H;
    300         // From (2)
    301         ((mySpline->spline[i])->coeff[1])+= (y->data.F32[i+1] / H);
    302         // From (3)
    303         tmp = -3.0 * (x->data.F32[i+1] * x->data.F32[i+1]) / (H * H * H);
    304         tmp+= (1.0 / H);
    305         tmp*= ((mySpline->p_psDeriv2)[i]) * H * H / 6.0;
    306         ((mySpline->spline[i])->coeff[1])+= tmp;
    307         // From (4)
    308         tmp = 3.0 * (x->data.F32[i] * x->data.F32[i]) / (H * H * H);
    309         tmp-= (1.0 / H);
    310         tmp*= ((mySpline->p_psDeriv2)[i+1]) * H * H / 6.0;
    311         ((mySpline->spline[i])->coeff[1])+= tmp;
    312 
    313         //
    314         // ******** Calulate 2-order term ********
    315         //
    316         // From (3)
    317         (mySpline->spline[i])->coeff[2] = ((mySpline->p_psDeriv2)[i]) * 3.0 * x->data.F32[i+1] / (6.0 * H);
    318         // From (4)
    319         ((mySpline->spline[i])->coeff[2])-= (((mySpline->p_psDeriv2)[i+1]) * 3.0 * x->data.F32[i] / (6.0 * H));
    320 
    321         //
    322         // ******** Calulate 3-order term ********
    323         //
    324         // From (3)
    325         (mySpline->spline[i])->coeff[3] = -((mySpline->p_psDeriv2)[i]) / (6.0 * H);
    326         // From (4)
    327         ((mySpline->spline[i])->coeff[3])+=  ((mySpline->p_psDeriv2)[i+1]) / (6.0 * H);
    328 
    329         psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
    330                 "(mySpline->spline[%d])->coeff[0] is %f\n", i, (mySpline->spline[i])->coeff[0]);
    331         psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
    332                 "(mySpline->spline[%d])->coeff[1] is %f\n", i, (mySpline->spline[i])->coeff[1]);
    333         psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
    334                 "(mySpline->spline[%d])->coeff[2] is %f\n", i, (mySpline->spline[i])->coeff[2]);
    335         psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
    336                 "(mySpline->spline[%d])->coeff[3] is %f\n", i, (mySpline->spline[i])->coeff[3]);
    337 
    338     }
    339 
    340     psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
    341             "---- psVectorFitSpline1D() end ----\n");
    342     return(mySpline);
    343 }
    344232
    345233/******************************************************************************
    … …  
    387275
    388276
     277/*****************************************************************************/
     278/* FUNCTION IMPLEMENTATION - PUBLIC                                          */
     279/*****************************************************************************/
     280
     281/*****************************************************************************
     282psVectorFitSpline1D(): given a psSpline1D data structure and a set of x/y
     283vectors, this routine generates the linear or cublic splines which satisfy
     284those data points.
     285 
     286The formula for calculating the spline polynomials is derived from Numerical
     287Recipes in C.  The basic idea is that the polynomial is
     288 (1)     y = (A * y[0]) +
     289 (2)         (B * y[1]) +
     290 (3)         ((((A*A*A)-A) * mySpline->p_psDeriv2[0]) * H^2)/6.0 +
     291 (4)         ((((B*B*B)-B) * mySpline->p_psDeriv2[1]) * H^2)/6.0
     292Where:
     293 H = x[1]-x[0]
     294 A = (x[1]-x)/H
     295 B = (x-x[0])/H
     296The bulk of the code in this routine is the expansion of the above equation
     297into a polynomial in terms of x, and then saving the coefficients of the
     298powers of x in the spline polynomials.  This gets pretty complicated.
     299 
     300XXX: usage of yErr is not specified in IfA documentation.
     301XXX: Must do: if yErr==NULL, set all errors equal.
     302 
     303XXX: Is the x argument redundant?  What do we do if the x argument is
     304 supplied, but does not equal the domains specified in mySpline?
     305 
     306XXX: can psSpline be NULL?
     307 
     308XXX: Implemented in F32 only, must add F64.
     309 *****************************************************************************/
     310psSpline1D *psVectorFitSpline1D(psSpline1D *mySpline,              ///< The spline which will be generated.
     311                                const psVector* restrict x,        ///< Ordinates (or NULL to just use the indices)
     312                                const psVector* restrict y,        ///< Coordinates
     313                                const psVector* restrict yErr)     ///< Errors in coordinates, or NULL
     314{
     315    psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
     316            "---- psVectorFitSpline1D() begin ----\n");
     317    int numSplines = (y->n)-1;
     318    float tmp;
     319    float H;
     320    int i;
     321    float slope;
     322
     323    if (mySpline == NULL) {
     324        //XXX psErrorMsg()
     325    }
     326    if (y == NULL) {
     327        //XXX psErrorMsg()
     328    }
     329
     330    if (y->n != (1 + mySpline->n)) {
     331        psLogMsg(__func__, PS_LOG_WARN,
     332                 "data size / spline size mismatch (%d %d)\n",
     333                 y->n, mySpline->n);
     334        // XXX: psErrorMsg()
     335    }
     336
     337    // If these are linear splines, which means their polynomials will have
     338    // two coefficients, then we do the simple calculation.
     339    if (2 == (mySpline->spline[0])->n) {
     340        for (i=0;i<mySpline->n;i++) {
     341            slope = (y->data.F32[i+1] - y->data.F32[i]) /
     342                    (mySpline->domains[i+1] - mySpline->domains[i]);
     343            (mySpline->spline[i])->coeff[0] = y->data.F32[i] -
     344                                              (slope * mySpline->domains[i]);
     345
     346            (mySpline->spline[i])->coeff[1] = slope;
     347            psTrace(".psLib.dataManip.psFunctions.psSpline1DGen", 4,
     348                    "---- mySpline %d coeffs are (%f, %f)\n", i,
     349                    (mySpline->spline[i])->coeff[0],
     350                    (mySpline->spline[i])->coeff[1]);
     351        }
     352        psTrace(".psLib.dataManip.psFunctions.psSpline1DGen", 4,
     353                "---- Exiting psSpline1DGen()\n");
     354        return((psSpline1D *) mySpline);
     355    }
     356
     357    // Check if these are cubic splines (n==4).  If not, psError.
     358    if (4 != (mySpline->spline[0])->n) {
     359        psLogMsg(__func__, PS_LOG_WARN,
     360                 "Don't know how to generate %d-order splines.",
     361                 (mySpline->spline[0])->n-1);
     362        // XXX: psErrorMsg()
     363        return(NULL);
     364    }
     365
     366    // If we get here, then we know these are cubic splines.  We first
     367    // generate the second derivatives at each data point.
     368    mySpline->p_psDeriv2 = CalculateSecondDerivs(x, y);
     369    for (i=0;i<y->n;i++)
     370        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
     371                "Second deriv[%d] is %f\n", i, mySpline->p_psDeriv2[i]);
     372
     373    // We generate the coefficients of the spline polynomials.  I can't
     374    // concisely explain how this code works.  See above function comments
     375    // and Numerical Recipes in C.
     376    for (i=0;i<numSplines;i++) {
     377        H = x->data.F32[i+1] - x->data.F32[i];
     378        psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
     379                "X data (%f - %f) (%f)\n", x->data.F32[i], x->data.F32[i+1], H);
     380        //
     381        // ******** Calculate 0-order term ********
     382        //
     383        // From (1)
     384        (mySpline->spline[i])->coeff[0] = (y->data.F32[i] * x->data.F32[i+1]/H);
     385        // From (2)
     386        ((mySpline->spline[i])->coeff[0])-= ((y->data.F32[i+1] * x->data.F32[i])/H);
     387        // From (3)
     388        tmp = (x->data.F32[i+1] * x->data.F32[i+1] * x->data.F32[i+1]) / (H * H * H);
     389        tmp-= (x->data.F32[i+1] / H);
     390        tmp*= (mySpline->p_psDeriv2)[i] * H * H / 6.0;
     391        ((mySpline->spline[i])->coeff[0])+= tmp;
     392        // From (4)
     393        tmp = -(x->data.F32[i] * x->data.F32[i] * x->data.F32[i]) / (H * H * H);
     394        tmp+= (x->data.F32[i] / H);
     395        tmp*= (mySpline->p_psDeriv2)[i+1] * H * H / 6.0;
     396        ((mySpline->spline[i])->coeff[0])+= tmp;
     397
     398        //
     399        // ******** Calculate 1-order term ********
     400        //
     401        // From (1)
     402        (mySpline->spline[i])->coeff[1] = -(y->data.F32[i]) / H;
     403        // From (2)
     404        ((mySpline->spline[i])->coeff[1])+= (y->data.F32[i+1] / H);
     405        // From (3)
     406        tmp = -3.0 * (x->data.F32[i+1] * x->data.F32[i+1]) / (H * H * H);
     407        tmp+= (1.0 / H);
     408        tmp*= ((mySpline->p_psDeriv2)[i]) * H * H / 6.0;
     409        ((mySpline->spline[i])->coeff[1])+= tmp;
     410        // From (4)
     411        tmp = 3.0 * (x->data.F32[i] * x->data.F32[i]) / (H * H * H);
     412        tmp-= (1.0 / H);
     413        tmp*= ((mySpline->p_psDeriv2)[i+1]) * H * H / 6.0;
     414        ((mySpline->spline[i])->coeff[1])+= tmp;
     415
     416        //
     417        // ******** Calculate 2-order term ********
     418        //
     419        // From (3)
     420        (mySpline->spline[i])->coeff[2] = ((mySpline->p_psDeriv2)[i]) * 3.0 * x->data.F32[i+1] / (6.0 * H);
     421        // From (4)
     422        ((mySpline->spline[i])->coeff[2])-= (((mySpline->p_psDeriv2)[i+1]) * 3.0 * x->data.F32[i] / (6.0 * H));
     423
     424        //
     425        // ******** Calculate 3-order term ********
     426        //
     427        // From (3)
     428        (mySpline->spline[i])->coeff[3] = -((mySpline->p_psDeriv2)[i]) / (6.0 * H);
     429        // From (4)
     430        ((mySpline->spline[i])->coeff[3])+=  ((mySpline->p_psDeriv2)[i+1]) / (6.0 * H);
     431
     432        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
     433                "(mySpline->spline[%d])->coeff[0] is %f\n", i, (mySpline->spline[i])->coeff[0]);
     434        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
     435                "(mySpline->spline[%d])->coeff[1] is %f\n", i, (mySpline->spline[i])->coeff[1]);
     436        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
     437                "(mySpline->spline[%d])->coeff[2] is %f\n", i, (mySpline->spline[i])->coeff[2]);
     438        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
     439                "(mySpline->spline[%d])->coeff[3] is %f\n", i, (mySpline->spline[i])->coeff[3]);
     440
     441    }
     442
     443    psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
     444            "---- psVectorFitSpline1D() end ----\n");
     445    return(mySpline);
     446}
     447
    389448/******************************************************************************
    390449psMinimizeLMChi2():  This routine will take an procedure which calculates an
    … …  
    398457 
    399458     A[i][i] = (1 + lambda*paramask[i]) * alpha[i][i];
     459 
     460XXX: This is very different from what is specified in the SDR.  Must
     461coordinate with IfA on new SDR.
     462 
     463XXX: Do vector/image recycles.
     464 
     465XXX: probably yErr will be part of the SDR.
     466 
     467XXX: This must work for both F32 and F64.  F32 is currently implemented.
     468     Note: since the LUD routines are only implemented in F64, then we
     469     will have to convert all F32 input vectors to F64 regardless.  So,
     470     the F64 port might be.
    400471 *****************************************************************************/
    401472bool psMinimizeLMChi2(psMinimization *min,
    … …  
    450521
    451522    min->lastDelta = HUGE;
    452     min->lastDelta = 12345.0;
    453523    min->iter = 0;
    454524
    … …  
    627697polynomial, as well as the error for each data point (yErr).
    628698 
    629 XXX: NOTE: yErr is currently ignored.
     699XXX: yErr is currently ignored.
     700 
     701XXX: must add type F32 (currently F64 only).
     702 
     703XXX: Must do: if x==NULL, use the index vector (???).
     704 
     705XXX: Must do: if yErr==NULL, set all errors equal.
    630706 *****************************************************************************/
    631707psPolynomial1D* psVectorFitPolynomial1D(psPolynomial1D* myPoly,
    … …  
    726802}
    727803
     804
     805/******************************************************************************
     806 *****************************************************************************/
    728807psMinimization *psMinimizationAlloc(int maxIter,
    729808                                    float tol)
    … …  
    747826                                         f(param + c * line).
    748827 
    749     Algorithm: XXX completely ad hoc: start with the user-supplied starting
    750     parameter and call that b.  Calculate a/c as a fractional amount
    751     smaller/larger than b.  Repeat this process until a local minimum is
    752     found.
    753  
    754  
     828Algorithm: XXX completely ad hoc: start with the user-supplied starting
     829parameter and call that b.  Calculate a/c as a fractional amount
     830smaller/larger than b.  Repeat this process until a local minimum is found.
     831 
     832XXX: new algorithm: start at x=0, expand in one direction until the function
     833decreases.  Then you have two points in the bracket.  Keep going until it
     834increases, or x is too large.  If thst does not work, expand in the other
     835direction.
    755836 *****************************************************************************/
    756837psVector *p_psDetermineBracket(psVector *params,
    … …  
    9261007    This routine must minimize a possibly multi-dimensional function
    9271008    along a vector defined by line.
     1009 
     1010XXX: Use a p_psName().
    9281011 *****************************************************************************/
    9291012float psLineMin(psMinimization *min,
    … …  
    11131196XXX: We do not use Brent's method.
    11141197 
    1115 XXX: We do not use the
     1198XXX: The SDR is silent about data types.  F32 is implemented here.
     1199 
     1200XXX: Define constants for default dummy number of iterations and tolerance.
    11161201 *****************************************************************************/
    11171202bool psMinimizePowell(psMinimization *min,
Note: See TracChangeset for help on using the changeset viewer.