IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Oct 1, 2009, 4:03:58 PM (17 years ago)
Author:
Paul Price
Message:

Merging stack development branch (implement lax and strict PSF model limits) from branches/pap

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk

  • trunk/psModules/src/objects/models/pmModel_SERSIC.c

    r20001 r25738  
    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 = paramsMinLax;
     54static float *paramsMaxUse = paramsMaxLax;
     55static float betaUse[] = { 1000, 3e6, 5, 5, 1.0, 1.0, 0.5, 2.0 };
     56
     57static bool limitsApply = true;         // Apply limits?
    3358
    3459psF32 PM_MODEL_FUNC (psVector *deriv,
    … …  
    91116bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta)
    92117{
    93     float beta_lim = 0, params_min = 0, params_max = 0;
    94     float f1 = 0, f2 = 0, q1 = 0, q2 = 0;
     118    if (!limitsApply) {
     119        return true;
     120    }
     121    psAssert(nParam >= 0 && nParam <= PM_PAR_7, "Parameter index is out of bounds");
    95122
    96123    // we need to calculate the limits for SXY specially
     124    float q2 = NAN;
    97125    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);
     126        float f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
     127        float f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
     128        float q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
    101129        q1 = (q1 < 0.0) ? 0.0 : q1;
    102130        // if q1 < 0.0, f2 ~ f1, we have a very large axis ratio near 45deg..  Saturate at that
    103131        // angle and let f2,f1 fight it out
    104         q2  = 0.5*sqrt (q1);
     132        q2 = 0.5*sqrtf(q1);
    105133    }
    106134
    107135    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:
     136      case PS_MINIMIZE_BETA_LIMIT: {
     137          psAssert(beta, "Require beta to limit beta");
     138          float limit = betaUse[nParam];
     139          if (nParam == PM_PAR_SXY) {
     140              limit *= q2;
     141          }
     142          if (fabs(beta[nParam]) > fabs(limit)) {
     143              beta[nParam] = (beta[nParam] > 0) ? fabs(limit) : -fabs(limit);
     144              psTrace("psModules.objects", 5, "|beta[nParam==%d]| > |beta_lim|; %g v. %g",
     145                      nParam, beta[nParam], limit);
     146              return false;
     147          }
     148          return true;
     149      }
     150      case PS_MINIMIZE_PARAM_MIN: {
     151          psAssert(params, "Require parameters to limit parameters");
     152          psAssert(paramsMinUse, "Require parameter limits to limit parameters");
     153          float limit = paramsMinUse[nParam];
     154          if (nParam == PM_PAR_SXY) {
     155              limit *= q2;
     156          }
     157          if (params[nParam] < limit) {
     158              params[nParam] = limit;
     159              psTrace("psModules.objects", 5, "params[nParam==%d] < params_min; %g v. %g",
     160                      nParam, params[nParam], limit);
     161              return false;
     162          }
     163          return true;
     164      }
     165      case PS_MINIMIZE_PARAM_MAX: {
     166          psAssert(params, "Require parameters to limit parameters");
     167          psAssert(paramsMaxUse, "Require parameter limits to limit parameters");
     168          float limit = paramsMaxUse[nParam];
     169          if (nParam == PM_PAR_SXY) {
     170              limit *= q2;
     171          }
     172          if (params[nParam] > limit) {
     173              params[nParam] = limit;
     174              psTrace("psModules.objects", 5, "params[nParam==%d] > params_max; %g v. %g",
     175                      nParam, params[nParam], limit);
     176              return false;
     177          }
     178          return true;
     179      }
     180      default:
    217181        psAbort("invalid choice for limits");
    218182    }
    … …  
    220184    return false;
    221185}
    222 
    223186
    224187// make an initial guess for parameters
    … …  
    447410
    448411    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]));
     412             dP, PAR[PM_PAR_I0], (dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]));
    450413
    451414    return status;
     415}
     416
     417
     418void PM_MODEL_SET_LIMITS(pmModelLimitsType type)
     419{
     420    switch (type) {
     421      case PM_MODEL_LIMITS_NONE:
     422        paramsMinUse = NULL;
     423        paramsMaxUse = NULL;
     424        limitsApply = true;
     425        break;
     426      case PM_MODEL_LIMITS_IGNORE:
     427        paramsMinUse = NULL;
     428        paramsMaxUse = NULL;
     429        limitsApply = false;
     430        break;
     431      case PM_MODEL_LIMITS_LAX:
     432        paramsMinUse = paramsMinLax;
     433        paramsMaxUse = paramsMaxLax;
     434        limitsApply = true;
     435        break;
     436      case PM_MODEL_LIMITS_STRICT:
     437        paramsMinUse = paramsMinStrict;
     438        paramsMaxUse = paramsMaxStrict;
     439        limitsApply = true;
     440        break;
     441      default:
     442        psAbort("Unrecognised model limits type: %x", type);
     443    }
     444    return;
    452445}
    453446
    … …  
    460453# undef PM_MODEL_PARAMS_FROM_PSF
    461454# undef PM_MODEL_FIT_STATUS
     455# undef PM_MODEL_SET_LIMITS
Note: See TracChangeset for help on using the changeset viewer.