IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 8926


Ignore:
Timestamp:
Sep 24, 2006, 3:06:28 PM (20 years ago)
Author:
magnier
Message:

finished, tested pmShutterCorrection

Location:
trunk/psModules
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/detrend/pmShutterCorrection.c

    r8886 r8926  
    101101    // XXX we could examine the top 2 or 3 values and decide if we
    102102    // extended exptime enough or median clip.
    103     pars->scale = counts->data.F64[N-1];
     103    pars->scale = counts->data.F32[N-1];
    104104
    105105    // fit a line to the lowest three points and extrapolate to 0.0
    106     tmpX = psVectorAlloc (2, PS_TYPE_F64);
    107     tmpY = psVectorAlloc (2, PS_TYPE_F64);
    108 
    109     tmpX->data.F64[0] = exptime->data.F64[0];
    110     tmpX->data.F64[1] = exptime->data.F64[1];
    111     tmpY->data.F64[0] = counts->data.F64[0];
    112     tmpY->data.F64[1] = counts->data.F64[1];
     106    tmpX = psVectorAlloc (2, PS_TYPE_F32);
     107    tmpY = psVectorAlloc (2, PS_TYPE_F32);
     108    tmpX->n = tmpY->n = 2;
     109
     110    tmpX->data.F32[0] = exptime->data.F32[0];
     111    tmpX->data.F32[1] = exptime->data.F32[1];
     112    tmpY->data.F32[0] = counts->data.F32[0];
     113    tmpY->data.F32[1] = counts->data.F32[1];
    113114
    114115    // fit a line and extrapolate the fit to 0.0
    … …  
    125126    float value = pars->scale * (1 + ratio) / 2.0;
    126127
    127     int Nm = psVectorBracket (exptime, value, (ratio < 1.0));
    128     int Np = (Nm == N - 1) ? Nm - 1 : Nm + 1;
    129 
    130     tmpX->data.F64[0] = counts->data.F64[Nm];
    131     tmpX->data.F64[1] = counts->data.F64[Np];
    132     tmpY->data.F64[0] = exptime->data.F64[Nm];
    133     tmpY->data.F64[1] = exptime->data.F64[Np];
     128    int Np;
     129    if (ratio < 1.0) {
     130        Np = psVectorBracket (counts, value, true);
     131    } else {
     132        Np = psVectorBracketDescend (counts, value, true);
     133    }
     134    int Nm = (Np == 0) ? 1 : Np - 1;
     135
     136    tmpX->data.F32[0] = counts->data.F32[Nm];
     137    tmpX->data.F32[1] = counts->data.F32[Np];
     138    tmpY->data.F32[0] = exptime->data.F32[Nm];
     139    tmpY->data.F32[1] = exptime->data.F32[Np];
    134140
    135141    // fit a line and extrapolate the fit to counts = A (1 + dTk/dTo) : exptime = dTo
    … …  
    146152
    147153// linear fit to the counts and exptime, given a value for offref
    148 pmShutterCorrPars *pmShutterCorrectionLinFit (psVector *exptime, psVector *counts, float offref)
     154pmShutterCorrPars *pmShutterCorrectionLinFit (psVector *exptime, psVector *counts, psVector *cntError, float offref)
    149155{
    150156
    … …  
    152158    psVector *x = psVectorAlloc (exptime->n, PS_TYPE_F32);
    153159    psVector *y = psVectorAlloc (exptime->n, PS_TYPE_F32);
     160    x->n = y->n = exptime->n;
    154161
    155162    for (int i = 0; i < exptime->n; i++) {
    … …  
    167174    line->coeff[1][1] = 0;
    168175
    169     psVectorFitPolynomial2D (line, NULL, 0, counts, NULL, x, y);
     176    psVectorFitPolynomial2D (line, NULL, 0, counts, cntError, x, y);
    170177
    171178    pmShutterCorrPars *pars = pmShutterCorrParsAlloc ();
    … …  
    174181    pars->scale  = line->coeff[1][0];
    175182    pars->offset = line->coeff[0][1] / line->coeff[1][0];
     183
     184    psFree (x);
     185    psFree (y);
     186    psFree (line);
    176187
    177188    return (pars);
    … …  
    196207
    197208// non-linear fit to the counts and exptime, given a guess for the three parameters
    198 pmShutterCorrPars *pmShutterCorrectionFullFit (psVector *exptime, psVector *counts, pmShutterCorrPars *guess)
     209pmShutterCorrPars *pmShutterCorrectionFullFit (psVector *exptime, psVector *counts, psVector *cntError, pmShutterCorrPars *guess)
    199210{
    200211
    … …  
    221232    // construct the coordinate and value entries (y is counts)
    222233    psArray *x = psArrayAlloc(exptime->n);
    223     psVector *y = counts;
    224     psVector *yErr = psVectorAlloc(exptime->n, PS_TYPE_F32);
    225     x->n = yErr->n = exptime->n;
     234    x->n = exptime->n;
    226235
    227236    for (int i = 0; i < exptime->n; i++) {
    228237        psVector *coord = psVectorAlloc(1, PS_TYPE_F32);
    229238        coord->n = 1;
    230 
    231239        coord->data.F32[0] = exptime->data.F32[i];
    232 
    233240        x->data[i] = coord;
    234         yErr->data.F32[i] = sqrt(y->data.F32[i]);
    235     }
    236 
    237     psMinimizeLMChi2(myMin, covar, params, constrain, x, y, yErr, pmShutterCorrectionModel);
     241    }
     242
     243    psMinimizeLMChi2(myMin, covar, params, constrain, x, counts, cntError, pmShutterCorrectionModel);
    238244
    239245    pmShutterCorrPars *pars = pmShutterCorrParsAlloc ();
    … …  
    243249    pars->offref = params->data.F32[2];
    244250
    245     return (pars);
    246 }
     251    psFree (myMin);
     252    psFree (params);
     253    psFree (x);
     254
     255    return (pars);
     256}

  • trunk/psModules/src/detrend/pmShutterCorrection.h

    r8884 r8926  
    4848 *  @author Eugene Magnier, IfA
    4949 *
    50  *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
    51  *  @date $Date: 2006-09-22 21:00:53 $
     50 *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
     51 *  @date $Date: 2006-09-25 01:06:28 $
    5252 *
    5353 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    6666pmShutterCorrPars *pmShutterCorrParsAlloc ();
    6767pmShutterCorrPars *pmShutterCorrectionGuess (psVector *exptime, psVector *counts);
    68 pmShutterCorrPars *pmShutterCorrectionLinFit (psVector *exptime, psVector *counts, float offref);
    69 pmShutterCorrPars *pmShutterCorrectionFullFit (psVector *exptime, psVector *counts, pmShutterCorrPars *guess);
     68pmShutterCorrPars *pmShutterCorrectionLinFit (psVector *exptime, psVector *counts, psVector *cntError, float offref);
     69pmShutterCorrPars *pmShutterCorrectionFullFit (psVector *exptime, psVector *counts, psVector *cntError, pmShutterCorrPars *guess);

  • trunk/psModules/test/detrend/tap_pmShutterCorrection.c

    r8919 r8926  
    1313
    1414    // test allocation, free of pmShutterCorrPars
     15    diag("pmShutterCorrParsAlloc tests");
    1516    {
    1617        psMemId id = psMemGetId();
    … …  
    2526    }
    2627
    27     // test parameter guess (linearly spaced exptimes)
     28    // test parameter guess (linearly spaced exptimes, TK/TO < 1)
     29    diag("pmShutterCorrectionGuess tests : coarse linear-spaced exptimes");
    2830    {
    2931        psMemId id = psMemGetId();
    … …  
    3739        exptime->n = counts->n = NPTS;
    3840
    39         for (int i = 0; i < exptime->n; i++)
    40         {
    41             exptime->data.F32[i] = i;
    42             counts->data.F32[i] = AK*(exptime->data.F32[i] + TK) / (exptime->data.F32[i] + TO);
    43         }
    44 
    45         pmShutterCorrPars *pars = pmShutterCorrectionGuess (exptime, counts);
    46 
    47         ok(pars != NULL, "pmShutterCorrPars successfully allocated");
    48         skip_start(pars == NULL, 0, "Skipping tests because pmShutterCorrParsAlloc() failed");
    49         ok(pars->scale == AK, "got the right scale (got %f vs %f)", pars->scale, AK);
    50         ok(pars->offset == TK, "got the right offset (got %f vs %f)", pars->offset, TK);
    51         ok(pars->offref == TO, "got the right offref (got %f vs %f)", pars->offref, TO);
    52         skip_end();
    53 
    54         psFree(pars);
    55         psFree(exptime);
    56         psFree(counts);
     41        for (int i = 0; i < exptime->n; i++) {
     42            exptime->data.F32[i] = i*0.25;
     43            counts->data.F32[i] = AK*(exptime->data.F32[i] + TK) / (exptime->data.F32[i] + TO);
     44        }
     45
     46        pmShutterCorrPars *pars = pmShutterCorrectionGuess (exptime, counts);
     47
     48        ok(pars != NULL, "pmShutterCorrPars successfully allocated");
     49        skip_start(pars == NULL, 0, "Skipping tests because pmShutterCorrParsAlloc() failed");
     50
     51        // with coarse linearly-spaced times large compared to TO and TK,
     52        // we can't expect very accurate guesses.  the exptime guess is fairly good because
     53        // the largest exptime is much longer than TK or TO
     54        ok(fabs(pars->scale  - AK) < 0.5, "scale guess is close enough (got %f vs %f)",  pars->scale, AK);
     55        ok(fabs(pars->offset - TK) < 0.5, "offset guess is close enough (got %f vs %f)", pars->offset, TK);
     56        ok(fabs(pars->offref - TO) < 0.5, "offref guess is close enough (got %f vs %f)", pars->offref, TO);
     57        skip_end();
     58
     59        psFree(pars);
     60        psFree(exptime);
     61        psFree(counts);
     62        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
     63    }
     64
     65    // test parameter guess (linearly spaced exptimes, TK/TO < 1)
     66    diag("pmShutterCorrectionGuess tests : fine linear-spaced exptimes");
     67    {
     68        psMemId id = psMemGetId();
     69
     70        int NPTS = 20;
     71        float AK = 5.0;
     72        float TK = 0.1;
     73        float TO = 0.2;
     74        psVector *exptime = psVectorAlloc (NPTS, PS_TYPE_F32);
     75        psVector *counts  = psVectorAlloc (NPTS, PS_TYPE_F32);
     76        exptime->n = counts->n = NPTS;
     77
     78        for (int i = 0; i < exptime->n; i++) {
     79            exptime->data.F32[i] = i*0.1;
     80            counts->data.F32[i] = AK*(exptime->data.F32[i] + TK) / (exptime->data.F32[i] + TO);
     81        }
     82
     83        pmShutterCorrPars *pars = pmShutterCorrectionGuess (exptime, counts);
     84
     85        ok(pars != NULL, "pmShutterCorrPars successfully allocated");
     86        skip_start(pars == NULL, 0, "Skipping tests because pmShutterCorrParsAlloc() failed");
     87
     88        // with fine linearly-spaced times large compared to TO and TK,
     89        // we get a good guess to TK and TO, but since the largest exptime is not
     90        // many times larger than TO, we don't do very well with the AK
     91        ok(fabs(pars->scale  - AK) < 0.5, "scale guess is close enough (got %f vs %f)",  pars->scale, AK);
     92        ok(fabs(pars->offset - TK) < 0.05, "offset guess is close enough (got %f vs %f)", pars->offset, TK);
     93        ok(fabs(pars->offref - TO) < 0.05, "offref guess is close enough (got %f vs %f)", pars->offref, TO);
     94        skip_end();
     95
     96        psFree(pars);
     97        psFree(exptime);
     98        psFree(counts);
     99        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
     100    }
     101
     102    // test parameter guess (log spaced exptimes, TK/TO < 1)
     103    diag("pmShutterCorrectionGuess tests : log-spaced exptimes");
     104    {
     105        psMemId id = psMemGetId();
     106
     107        int NPTS = 40;
     108        float AK = 5.0;
     109        float TK = 0.1;
     110        float TO = 0.2;
     111        psVector *exptime = psVectorAlloc (NPTS, PS_TYPE_F32);
     112        psVector *counts  = psVectorAlloc (NPTS, PS_TYPE_F32);
     113        exptime->n = counts->n = NPTS;
     114
     115        for (int i = 0; i < exptime->n; i++) {
     116            exptime->data.F32[i] = pow(10.0, -2 + i*0.1);
     117            counts->data.F32[i] = AK*(exptime->data.F32[i] + TK) / (exptime->data.F32[i] + TO);
     118        }
     119
     120        pmShutterCorrPars *pars = pmShutterCorrectionGuess (exptime, counts);
     121
     122        // with fine log-spaced times well-sampling TO and TK,
     123        // we can expect accurate guesses
     124        ok(pars != NULL, "pmShutterCorrPars successfully allocated");
     125        skip_start(pars == NULL, 0, "Skipping tests because pmShutterCorrParsAlloc() failed");
     126        ok(fabs(pars->scale  - AK) < 0.01, "scale guess is close enough (got %f vs %f)",  pars->scale, AK);
     127        ok(fabs(pars->offset - TK) < 0.01, "offset guess is close enough (got %f vs %f)", pars->offset, TK);
     128        ok(fabs(pars->offref - TO) < 0.01, "offref guess is close enough (got %f vs %f)", pars->offref, TO);
     129        skip_end();
     130
     131        psFree(pars);
     132        psFree(exptime);
     133        psFree(counts);
     134        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
     135    }
     136
     137
     138    // test parameter guess (linearly spaced exptimes, TK/TO > 1)
     139    diag("pmShutterCorrectionGuess tests : coarse linear-spaced exptimes, TK/TO > 1");
     140    {
     141        psMemId id = psMemGetId();
     142
     143        int NPTS = 10;
     144        float AK = 5.0;
     145        float TK = 0.2;
     146        float TO = 0.1;
     147        psVector *exptime = psVectorAlloc (NPTS, PS_TYPE_F32);
     148        psVector *counts  = psVectorAlloc (NPTS, PS_TYPE_F32);
     149        exptime->n = counts->n = NPTS;
     150
     151        for (int i = 0; i < exptime->n; i++) {
     152            exptime->data.F32[i] = i*0.25;
     153            counts->data.F32[i] = AK*(exptime->data.F32[i] + TK) / (exptime->data.F32[i] + TO);
     154        }
     155
     156        pmShutterCorrPars *pars = pmShutterCorrectionGuess (exptime, counts);
     157
     158        ok(pars != NULL, "pmShutterCorrPars successfully allocated");
     159        skip_start(pars == NULL, 0, "Skipping tests because pmShutterCorrParsAlloc() failed");
     160
     161        // with coarse linearly-spaced times large compared to TO and TK,
     162        // we can't expect very accurate guesses.  the exptime guess is fairly good because
     163        // the largest exptime is much longer than TK or TO
     164        ok(fabs(pars->scale  - AK) < 0.5, "scale guess is close enough (got %f vs %f)",  pars->scale, AK);
     165        ok(fabs(pars->offset - TK) < 0.5, "offset guess is close enough (got %f vs %f)", pars->offset, TK);
     166        ok(fabs(pars->offref - TO) < 0.5, "offref guess is close enough (got %f vs %f)", pars->offref, TO);
     167        skip_end();
     168
     169        psFree(pars);
     170        psFree(exptime);
     171        psFree(counts);
     172        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
     173    }
     174
     175    // test parameter guess (linearly spaced exptimes, TK/TO > 1)
     176    diag("pmShutterCorrectionGuess tests : fine linear-spaced exptimes, TK/TO > 1");
     177    {
     178        psMemId id = psMemGetId();
     179
     180        int NPTS = 20;
     181        float AK = 5.0;
     182        float TK = 0.2;
     183        float TO = 0.1;
     184        psVector *exptime = psVectorAlloc (NPTS, PS_TYPE_F32);
     185        psVector *counts  = psVectorAlloc (NPTS, PS_TYPE_F32);
     186        exptime->n = counts->n = NPTS;
     187
     188        for (int i = 0; i < exptime->n; i++) {
     189            exptime->data.F32[i] = i*0.1;
     190            counts->data.F32[i] = AK*(exptime->data.F32[i] + TK) / (exptime->data.F32[i] + TO);
     191        }
     192
     193        pmShutterCorrPars *pars = pmShutterCorrectionGuess (exptime, counts);
     194
     195        ok(pars != NULL, "pmShutterCorrPars successfully allocated");
     196        skip_start(pars == NULL, 0, "Skipping tests because pmShutterCorrParsAlloc() failed");
     197
     198        // with fine linearly-spaced times large compared to TO and TK,
     199        // we get a good guess to TK and TO, but since the largest exptime is not
     200        // many times larger than TO, we don't do very well with the AK
     201        ok(fabs(pars->scale  - AK) < 0.5, "scale guess is close enough (got %f vs %f)",  pars->scale, AK);
     202        ok(fabs(pars->offset - TK) < 0.05, "offset guess is close enough (got %f vs %f)", pars->offset, TK);
     203        ok(fabs(pars->offref - TO) < 0.05, "offref guess is close enough (got %f vs %f)", pars->offref, TO);
     204        skip_end();
     205
     206        psFree(pars);
     207        psFree(exptime);
     208        psFree(counts);
     209        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
     210    }
     211
     212    // test parameter guess (log spaced exptimes, TK/TO > 1)
     213    diag("pmShutterCorrectionGuess tests : log-spaced exptimes, TK/TO > 1");
     214    {
     215        psMemId id = psMemGetId();
     216
     217        int NPTS = 40;
     218        float AK = 5.0;
     219        float TK = 0.2;
     220        float TO = 0.1;
     221        psVector *exptime = psVectorAlloc (NPTS, PS_TYPE_F32);
     222        psVector *counts  = psVectorAlloc (NPTS, PS_TYPE_F32);
     223        exptime->n = counts->n = NPTS;
     224
     225        for (int i = 0; i < exptime->n; i++) {
     226            exptime->data.F32[i] = pow(10.0, -2 + i*0.1);
     227            counts->data.F32[i] = AK*(exptime->data.F32[i] + TK) / (exptime->data.F32[i] + TO);
     228        }
     229
     230        pmShutterCorrPars *pars = pmShutterCorrectionGuess (exptime, counts);
     231
     232        // with fine log-spaced times well-sampling TO and TK,
     233        // we can expect accurate guesses
     234        ok(pars != NULL, "pmShutterCorrPars successfully allocated");
     235        skip_start(pars == NULL, 0, "Skipping tests because pmShutterCorrParsAlloc() failed");
     236        ok(fabs(pars->scale  - AK) < 0.01, "scale guess is close enough (got %f vs %f)",  pars->scale, AK);
     237        ok(fabs(pars->offset - TK) < 0.01, "offset guess is close enough (got %f vs %f)", pars->offset, TK);
     238        ok(fabs(pars->offref - TO) < 0.01, "offref guess is close enough (got %f vs %f)", pars->offref, TO);
     239        skip_end();
     240
     241        psFree(pars);
     242        psFree(exptime);
     243        psFree(counts);
     244        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
     245    }
     246
     247    // test non-linear fitting
     248    diag("pmShutterCorrectionFullFit tests : linear-spaced exptimes");
     249    {
     250        psMemId id = psMemGetId();
     251
     252        int NPTS = 20;
     253        float FL = 10000.0;
     254        float AK = 5.0;
     255        float TK = 0.2;
     256        float TO = 0.1;
     257        psVector *exptime = psVectorAlloc (NPTS, PS_TYPE_F32);
     258        psVector *counts  = psVectorAlloc (NPTS, PS_TYPE_F32);
     259        psVector *cntErr  = psVectorAlloc (NPTS, PS_TYPE_F32);
     260        exptime->n = counts->n = cntErr->n = NPTS;
     261
     262        for (int i = 0; i < exptime->n; i++) {
     263            exptime->data.F32[i] = i*0.1;
     264            counts->data.F32[i] = AK*(exptime->data.F32[i] + TK) / (exptime->data.F32[i] + TO);
     265            cntErr->data.F32[i] = AK*sqrt(FL*(exptime->data.F32[i] + TK)) / (exptime->data.F32[i] + TO);
     266        }
     267
     268        pmShutterCorrPars *guess = pmShutterCorrectionGuess (exptime, counts);
     269        skip_start(guess == NULL, 0, "Skipping tests because pmShutterCorrectionGuess() failed");
     270        pmShutterCorrPars *pars = pmShutterCorrectionFullFit (exptime, counts, cntErr, guess);
     271
     272        // with fine log-spaced times well-sampling TO and TK,
     273        // we can expect accurate guesses
     274        ok(pars != NULL, "pmShutterCorrPars successfully allocated by FullFit");
     275        skip_start(pars == NULL, 0, "Skipping tests because pmShutterCorrParsAlloc() failed");
     276        ok(fabs(pars->scale  - AK) < 0.01, "scale fit is close enough (got %f vs %f)",  pars->scale, AK);
     277        ok(fabs(pars->offset - TK) < 0.01, "offset fit is close enough (got %f vs %f)", pars->offset, TK);
     278        ok(fabs(pars->offref - TO) < 0.01, "offref fit is close enough (got %f vs %f)", pars->offref, TO);
     279        skip_end();
     280
     281        psFree(pars);
     282        skip_end();
     283
     284        psFree(guess);
     285        psFree(exptime);
     286        psFree(counts);
     287        psFree(cntErr);
     288        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
     289    }
     290
     291    // test non-linear fitting
     292    diag("pmShutterCorrectionFullFit tests : log-spaced exptimes");
     293    {
     294        psMemId id = psMemGetId();
     295
     296        int NPTS = 40;
     297        float FL = 10000.0;
     298        float AK = 1.0;
     299        float TK = 0.2;
     300        float TO = 0.1;
     301        psVector *exptime = psVectorAlloc (NPTS, PS_TYPE_F32);
     302        psVector *counts  = psVectorAlloc (NPTS, PS_TYPE_F32);
     303        psVector *cntErr  = psVectorAlloc (NPTS, PS_TYPE_F32);
     304        exptime->n = counts->n = cntErr->n = NPTS;
     305
     306        for (int i = 0; i < exptime->n; i++) {
     307            exptime->data.F32[i] = pow(10.0, -2 + i*0.1);
     308            counts->data.F32[i] = AK*(exptime->data.F32[i] + TK) / (exptime->data.F32[i] + TO);
     309            cntErr->data.F32[i] = AK*sqrt(FL*(exptime->data.F32[i] + TK)) / (exptime->data.F32[i] + TO);
     310        }
     311
     312        pmShutterCorrPars *guess = pmShutterCorrectionGuess (exptime, counts);
     313        skip_start(guess == NULL, 0, "Skipping tests because pmShutterCorrectionGuess() failed");
     314        pmShutterCorrPars *pars = pmShutterCorrectionFullFit (exptime, counts, cntErr, guess);
     315
     316        // with fine log-spaced times well-sampling TO and TK,
     317        // we can expect accurate guesses
     318        ok(pars != NULL, "pmShutterCorrPars successfully allocated by FullFit");
     319        skip_start(pars == NULL, 0, "Skipping tests because pmShutterCorrParsAlloc() failed");
     320        ok(fabs(pars->scale  - AK) < 0.01, "scale fit is close enough (got %f vs %f)",  pars->scale, AK);
     321        ok(fabs(pars->offset - TK) < 0.01, "offset fit is close enough (got %f vs %f)", pars->offset, TK);
     322        ok(fabs(pars->offref - TO) < 0.01, "offref fit is close enough (got %f vs %f)", pars->offref, TO);
     323        skip_end();
     324
     325        psFree(pars);
     326        skip_end();
     327
     328        psFree(guess);
     329        psFree(exptime);
     330        psFree(counts);
     331        psFree(cntErr);
     332        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
     333    }
     334
     335    // XXX should add tests with the input counts scattered with GaussDev...
     336
     337    // test linear fitting
     338    diag("pmShutterCorrectionLinFit tests : linear-spaced exptimes");
     339    {
     340        psMemId id = psMemGetId();
     341
     342        int NPTS = 20;
     343        float FL = 10000.0;
     344        float AK = 5.0;
     345        float TK = 0.2;
     346        float TO = 0.1;
     347        psVector *exptime = psVectorAlloc (NPTS, PS_TYPE_F32);
     348        psVector *counts  = psVectorAlloc (NPTS, PS_TYPE_F32);
     349        psVector *cntErr  = psVectorAlloc (NPTS, PS_TYPE_F32);
     350        exptime->n = counts->n = cntErr->n = NPTS;
     351
     352        for (int i = 0; i < exptime->n; i++) {
     353            exptime->data.F32[i] = i*0.1;
     354            counts->data.F32[i] = AK*(exptime->data.F32[i] + TK) / (exptime->data.F32[i] + TO);
     355            cntErr->data.F32[i] = AK*sqrt(FL*(exptime->data.F32[i] + TK)) / (exptime->data.F32[i] + TO);
     356        }
     357
     358        pmShutterCorrPars *guess = pmShutterCorrectionGuess (exptime, counts);
     359        skip_start(guess == NULL, 0, "Skipping tests because pmShutterCorrectionGuess() failed");
     360        pmShutterCorrPars *full = pmShutterCorrectionFullFit (exptime, counts, cntErr, guess);
     361        pmShutterCorrPars *pars = pmShutterCorrectionLinFit (exptime, counts, cntErr, full->offref);
     362
     363        // with fine log-spaced times well-sampling TO and TK,
     364        // we can expect accurate guesses
     365        ok(pars != NULL, "pmShutterCorrPars successfully allocated by FullFit");
     366        skip_start(pars == NULL, 0, "Skipping tests because pmShutterCorrParsAlloc() failed");
     367        ok(fabs(pars->scale  - AK) < 0.01, "scale fit is close enough (got %f vs %f)",  pars->scale, AK);
     368        ok(fabs(pars->offset - TK) < 0.01, "offset fit is close enough (got %f vs %f)", pars->offset, TK);
     369        ok(fabs(pars->offref - TO) < 0.01, "offref fit is close enough (got %f vs %f)", pars->offref, TO);
     370        skip_end();
     371
     372        psFree(pars);
     373        psFree(full);
     374        skip_end();
     375
     376        psFree(guess);
     377        psFree(exptime);
     378        psFree(counts);
     379        psFree(cntErr);
     380        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
     381    }
     382
     383    // test linear fitting
     384    diag("pmShutterCorrectionLinFit tests : log-spaced exptimes");
     385    {
     386        psMemId id = psMemGetId();
     387
     388        int NPTS = 40;
     389        float FL = 10000.0;
     390        float AK = 1.0;
     391        float TK = 0.2;
     392        float TO = 0.1;
     393        psVector *exptime = psVectorAlloc (NPTS, PS_TYPE_F32);
     394        psVector *counts  = psVectorAlloc (NPTS, PS_TYPE_F32);
     395        psVector *cntErr  = psVectorAlloc (NPTS, PS_TYPE_F32);
     396        exptime->n = counts->n = cntErr->n = NPTS;
     397
     398        for (int i = 0; i < exptime->n; i++) {
     399            exptime->data.F32[i] = pow(10.0, -2 + i*0.1);
     400            counts->data.F32[i] = AK*(exptime->data.F32[i] + TK) / (exptime->data.F32[i] + TO);
     401            cntErr->data.F32[i] = AK*sqrt(FL*(exptime->data.F32[i] + TK)) / (exptime->data.F32[i] + TO);
     402        }
     403
     404        pmShutterCorrPars *guess = pmShutterCorrectionGuess (exptime, counts);
     405        skip_start(guess == NULL, 0, "Skipping tests because pmShutterCorrectionGuess() failed");
     406        pmShutterCorrPars *full = pmShutterCorrectionFullFit (exptime, counts, cntErr, guess);
     407        pmShutterCorrPars *pars = pmShutterCorrectionLinFit (exptime, counts, cntErr, full->offref);
     408
     409        // with fine log-spaced times well-sampling TO and TK,
     410        // we can expect accurate guesses
     411        ok(pars != NULL, "pmShutterCorrPars successfully allocated by FullFit");
     412        skip_start(pars == NULL, 0, "Skipping tests because pmShutterCorrParsAlloc() failed");
     413        ok(fabs(pars->scale  - AK) < 0.01, "scale fit is close enough (got %f vs %f)",  pars->scale, AK);
     414        ok(fabs(pars->offset - TK) < 0.01, "offset fit is close enough (got %f vs %f)", pars->offset, TK);
     415        ok(fabs(pars->offref - TO) < 0.01, "offref fit is close enough (got %f vs %f)", pars->offref, TO);
     416        skip_end();
     417
     418        psFree(pars);
     419        psFree(full);
     420        skip_end();
     421
     422        psFree(guess);
     423        psFree(exptime);
     424        psFree(counts);
     425        psFree(cntErr);
    57426        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
    58427    }
Note: See TracChangeset for help on using the changeset viewer.