IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Sep 23, 2009, 5:12:22 PM (17 years ago)
Author:
Paul Price
Message:

Added multiple sets of limits for models. Involved reorganisation of how the models are sucked into the pmModelClasses: now using header files instead of including the source for each model. This allows static variables to be properly respected.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/pap/psModules/src/objects/models/pmModel_SERSIC.c

    r20001 r25521  
    2323   *****************************************************************************/
    2424
     25#include <stdio.h>
     26#include <pslib.h>
     27
     28#include "pmMoments.h"
     29#include "pmPeaks.h"
     30#include "pmSource.h"
     31#include "pmModel.h"
     32#include "pmModel_SERSIC.h"
     33
    2534# define PM_MODEL_FUNC            pmModelFunc_SERSIC
    2635# define PM_MODEL_FLUX            pmModelFlux_SERSIC
    … …  
    3140# define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_SERSIC
    3241# define PM_MODEL_FIT_STATUS      pmModelFitStatus_SERSIC
     42# define PM_MODEL_SET_LIMITS      pmModelSetLimits_SERSIC
     43
     44// Lax parameter limits
     45static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.05, 0.05, -1.0, 0.05 };
     46static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 4.0 };
     47
     48// Strict parameter limits
     49static float *paramsMinStrict = paramsMinLax;
     50static float *paramsMaxStrict = paramsMaxLax;
     51
     52// Parameter limits to use
     53static float *paramsMinUse = NULL;
     54static float *paramsMaxUse = NULL;
     55static float betaUse[] = { 1000, 3e6, 5, 5, 1.0, 1.0, 0.5, 2.0 };
    3356
    3457psF32 PM_MODEL_FUNC (psVector *deriv,
    … …  
    91114bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta)
    92115{
    93     float beta_lim = 0, params_min = 0, params_max = 0;
    94     float f1 = 0, f2 = 0, q1 = 0, q2 = 0;
     116    psAssert(nParam >= 0 && nParam <= PM_PAR_7, "Parameter index is out of bounds");
    95117
    96118    // we need to calculate the limits for SXY specially
     119    float q2 = NAN;
    97120    if (nParam == PM_PAR_SXY) {
    98         f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
    99         f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
    100         q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
     121        float f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
     122        float f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
     123        float q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
    101124        q1 = (q1 < 0.0) ? 0.0 : q1;
    102125        // if q1 < 0.0, f2 ~ f1, we have a very large axis ratio near 45deg..  Saturate at that
    103126        // angle and let f2,f1 fight it out
    104         q2  = 0.5*sqrt (q1);
     127        q2 = 0.5*sqrtf(q1);
    105128    }
    106129
    107130    switch (mode) {
    108     case PS_MINIMIZE_BETA_LIMIT:
    109         switch (nParam) {
    110         case PM_PAR_SKY:
    111             beta_lim = 1000;
    112             break;
    113         case PM_PAR_I0:
    114             beta_lim = 3e6;
    115             break;
    116         case PM_PAR_XPOS:
    117             beta_lim = 5;
    118             break;
    119         case PM_PAR_YPOS:
    120             beta_lim = 5;
    121             break;
    122         case PM_PAR_SXX:
    123             beta_lim = 1.0;
    124             break;
    125         case PM_PAR_SYY:
    126             beta_lim = 1.0;
    127             break;
    128         case PM_PAR_SXY:
    129             beta_lim =  0.5*q2;
    130             break;
    131         case PM_PAR_7:
    132             beta_lim = 2.0;
    133             break;
    134         default:
    135             psAbort("invalid parameter %d for beta test", nParam);
    136         }
    137         if (fabs(beta[nParam]) > fabs(beta_lim)) {
    138             beta[nParam] = (beta[nParam] > 0) ? fabs(beta_lim) : -fabs(beta_lim);
    139             psTrace ("psModules.objects", 5, "|beta[nParam==%d]| > |beta_lim|; %g v. %g",
    140                      nParam, beta[nParam], beta_lim);
    141             return false;
    142         }
    143         return true;
    144     case PS_MINIMIZE_PARAM_MIN:
    145         switch (nParam) {
    146         case PM_PAR_SKY:
    147             params_min = -1000;
    148             break;
    149         case PM_PAR_I0:
    150             params_min =     0.01;
    151             break;
    152         case PM_PAR_XPOS:
    153             params_min =  -100;
    154             break;
    155         case PM_PAR_YPOS:
    156             params_min =  -100;
    157             break;
    158         case PM_PAR_SXX:
    159             params_min =   0.05;
    160             break;
    161         case PM_PAR_SYY:
    162             params_min =   0.05;
    163             break;
    164         case PM_PAR_SXY:
    165             params_min =  -q2;
    166             break;
    167         case PM_PAR_7:
    168             params_min =   0.05;
    169             break;
    170         default:
    171             psAbort("invalid parameter %d for param min test", nParam);
    172         }
    173         if (params[nParam] < params_min) {
    174             params[nParam] = params_min;
    175             psTrace ("psModules.objects", 5, "params[nParam==%d] < params_min; %g v. %g",
    176                      nParam, params[nParam], params_min);
    177             return false;
    178         }
    179         return true;
    180     case PS_MINIMIZE_PARAM_MAX:
    181         switch (nParam) {
    182         case PM_PAR_SKY:
    183             params_max =   1e5;
    184             break;
    185         case PM_PAR_I0:
    186             params_max =   1e8;
    187             break;
    188         case PM_PAR_XPOS:
    189             params_max =   1e4;
    190             break;
    191         case PM_PAR_YPOS:
    192             params_max =   1e4;
    193             break;
    194         case PM_PAR_SXX:
    195             params_max =   100;
    196             break;
    197         case PM_PAR_SYY:
    198             params_max =   100;
    199             break;
    200         case PM_PAR_SXY:
    201             params_max =  +q2;
    202             break;
    203         case PM_PAR_7:
    204             params_max =   4.0;
    205             break;
    206         default:
    207             psAbort("invalid parameter %d for param max test", nParam);
    208         }
    209         if (params[nParam] > params_max) {
    210             params[nParam] = params_max;
    211             psTrace ("psModules.objects", 5, "params[nParam==%d] > params_max; %g v. %g",
    212                      nParam, params[nParam], params_max);
    213             return false;
    214         }
    215         return true;
    216     default:
     131      case PS_MINIMIZE_BETA_LIMIT: {
     132          psAssert(beta, "Require beta to limit beta");
     133          float limit = betaUse[nParam];
     134          if (nParam == PM_PAR_SXY) {
     135              limit *= q2;
     136          }
     137          if (fabs(beta[nParam]) > fabs(limit)) {
     138              beta[nParam] = (beta[nParam] > 0) ? fabs(limit) : -fabs(limit);
     139              psTrace("psModules.objects", 5, "|beta[nParam==%d]| > |beta_lim|; %g v. %g",
     140                      nParam, beta[nParam], limit);
     141              return false;
     142          }
     143          return true;
     144      }
     145      case PS_MINIMIZE_PARAM_MIN: {
     146          psAssert(params, "Require parameters to limit parameters");
     147          psAssert(paramsMinUse, "Require parameter limits to limit parameters");
     148          float limit = paramsMinUse[nParam];
     149          if (nParam == PM_PAR_SXY) {
     150              limit *= q2;
     151          }
     152          if (params[nParam] < limit) {
     153              params[nParam] = limit;
     154              psTrace("psModules.objects", 5, "params[nParam==%d] < params_min; %g v. %g",
     155                      nParam, params[nParam], limit);
     156              return false;
     157          }
     158          return true;
     159      }
     160      case PS_MINIMIZE_PARAM_MAX: {
     161          psAssert(params, "Require parameters to limit parameters");
     162          psAssert(paramsMaxUse, "Require parameter limits to limit parameters");
     163          float limit = paramsMaxUse[nParam];
     164          if (nParam == PM_PAR_SXY) {
     165              limit *= q2;
     166          }
     167          if (params[nParam] > limit) {
     168              params[nParam] = limit;
     169              psTrace("psModules.objects", 5, "params[nParam==%d] > params_max; %g v. %g",
     170                      nParam, params[nParam], limit);
     171              return false;
     172          }
     173          return true;
     174      }
     175      default:
    217176        psAbort("invalid choice for limits");
    218177    }
    … …  
    220179    return false;
    221180}
    222 
    223181
    224182// make an initial guess for parameters
    … …  
    447405
    448406    fprintf (stderr, "SERSIC status pars: dP: %f, I0: %f, S/N: %f\n",
    449              dP, PAR[PM_PAR_I0], (dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]));
     407             dP, PAR[PM_PAR_I0], (dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]));
    450408
    451409    return status;
     410}
     411
     412void PM_MODEL_SET_LIMITS(pmModelLimitsType type)
     413{
     414    switch (type) {
     415      case PM_MODEL_LIMITS_NONE:
     416        paramsMinUse = NULL;
     417        paramsMaxUse = NULL;
     418      case PM_MODEL_LIMITS_LAX:
     419        paramsMinUse = paramsMinLax;
     420        paramsMaxUse = paramsMaxLax;
     421        break;
     422      case PM_MODEL_LIMITS_STRICT:
     423        paramsMinUse = paramsMinStrict;
     424        paramsMaxUse = paramsMaxStrict;
     425        break;
     426      default:
     427        psAbort("Unrecognised model limits type: %x", type);
     428    }
     429    return;
    452430}
    453431
    … …  
    460438# undef PM_MODEL_PARAMS_FROM_PSF
    461439# undef PM_MODEL_FIT_STATUS
     440# undef PM_MODEL_SET_LIMITS
Note: See TracChangeset for help on using the changeset viewer.