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_RGAUSS.c

    r20001 r25738  
    2020 *****************************************************************************/
    2121
     22#include <stdio.h>
     23#include <pslib.h>
     24
     25#include "pmMoments.h"
     26#include "pmPeaks.h"
     27#include "pmSource.h"
     28#include "pmModel.h"
     29#include "pmModel_RGAUSS.h"
     30
    2231# define PM_MODEL_FUNC            pmModelFunc_RGAUSS
    2332# define PM_MODEL_FLUX            pmModelFlux_RGAUSS
    … …  
    2837# define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_RGAUSS
    2938# define PM_MODEL_FIT_STATUS      pmModelFitStatus_RGAUSS
     39# define PM_MODEL_SET_LIMITS      pmModelSetLimits_RGAUSS
     40
     41// Lax parameter limits
     42static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, 1.25 };
     43static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 4.0 };
     44
     45// Strict parameter limits
     46static float *paramsMinStrict = paramsMinLax;
     47static float *paramsMaxStrict = paramsMaxLax;
     48
     49// Parameter limits to use
     50static float *paramsMinUse = paramsMinLax;
     51static float *paramsMaxUse = paramsMaxLax;
     52static float betaUse[] = { 1000, 3e6, 5, 5, 0.5, 0.5, 0.5, 0.5 };
     53
     54static bool limitsApply = true;         // Apply limits?
    3055
    3156psF32 PM_MODEL_FUNC (psVector *deriv,
    … …  
    7398# define AR_MAX 20.0
    7499# define AR_RATIO 0.99
     100
    75101bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta)
    76102{
    77     float beta_lim = 0, params_min = 0, params_max = 0;
    78     float f1 = 0, f2 = 0, q1 = 0, q2 = 0;
     103    if (!limitsApply) {
     104        return true;
     105    }
     106    psAssert(nParam >= 0 && nParam <= PM_PAR_7, "Parameter index is out of bounds");
    79107
    80108    // we need to calculate the limits for SXY specially
     109    float q2 = NAN;
    81110    if (nParam == PM_PAR_SXY) {
    82         f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
    83         f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
    84         q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
     111        float f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
     112        float f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
     113        float q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
    85114        q1 = (q1 < 0.0) ? 0.0 : q1;
    86115        // if q1 < 0.0, f2 ~ f1, we have a very large axis ratio near 45deg..  Saturate at that
    87116        // angle and let f2,f1 fight it out
    88         q2  = 0.5*sqrt (q1);
     117        q2 = 0.5*sqrtf(q1);
    89118    }
    90119
    91120    switch (mode) {
    92     case PS_MINIMIZE_BETA_LIMIT:
    93         switch (nParam) {
    94         case PM_PAR_SKY:
    95             beta_lim = 1000;
    96             break;
    97         case PM_PAR_I0:
    98             beta_lim = 3e6;
    99             break;
    100         case PM_PAR_XPOS:
    101             beta_lim = 5;
    102             break;
    103         case PM_PAR_YPOS:
    104             beta_lim = 5;
    105             break;
    106         case PM_PAR_SXX:
    107             beta_lim = 0.5;
    108             break;
    109         case PM_PAR_SYY:
    110             beta_lim = 0.5;
    111             break;
    112         case PM_PAR_SXY:
    113             beta_lim =  0.5*q2;
    114             break;
    115         case PM_PAR_7:
    116             beta_lim = 0.5;
    117             break;
    118         default:
    119             psAbort("invalid parameter %d for beta test", nParam);
    120         }
    121         if (fabs(beta[nParam]) > fabs(beta_lim)) {
    122             beta[nParam] = (beta[nParam] > 0) ? fabs(beta_lim) : -fabs(beta_lim);
    123             psTrace ("psModules.objects", 5, "|beta[nParam==%d]| > |beta_lim|; %g v. %g",
    124                      nParam, beta[nParam], beta_lim);
    125             return false;
    126         }
    127         return true;
    128     case PS_MINIMIZE_PARAM_MIN:
    129         switch (nParam) {
    130         case PM_PAR_SKY:
    131             params_min = -1000;
    132             break;
    133         case PM_PAR_I0:
    134             params_min =   0.01;
    135             break;
    136         case PM_PAR_XPOS:
    137             params_min =  -100;
    138             break;
    139         case PM_PAR_YPOS:
    140             params_min =  -100;
    141             break;
    142         case PM_PAR_SXX:
    143             params_min =   0.5;
    144             break;
    145         case PM_PAR_SYY:
    146             params_min =   0.5;
    147             break;
    148         case PM_PAR_SXY:
    149             params_min =  -q2;
    150             break;
    151         case PM_PAR_7:
    152             params_min =   1.25;
    153             break;
    154         default:
    155             psAbort("invalid parameter %d for param min test", nParam);
    156         }
    157         if (params[nParam] < params_min) {
    158             params[nParam] = params_min;
    159             psTrace ("psModules.objects", 5, "params[nParam==%d] < params_min; %g v. %g",
    160                      nParam, params[nParam], params_min);
    161             return false;
    162         }
    163         return true;
    164     case PS_MINIMIZE_PARAM_MAX:
    165         switch (nParam) {
    166         case PM_PAR_SKY:
    167             params_max =   1e5;
    168             break;
    169         case PM_PAR_I0:
    170             params_max =   1e8;
    171             break;
    172         case PM_PAR_XPOS:
    173             params_max =   1e4;
    174             break;
    175         case PM_PAR_YPOS:
    176             params_max =   1e4;
    177             break;
    178         case PM_PAR_SXX:
    179             params_max =   100;
    180             break;
    181         case PM_PAR_SYY:
    182             params_max =   100;
    183             break;
    184         case PM_PAR_SXY:
    185             params_max =  +q2;
    186             break;
    187         case PM_PAR_7:
    188             params_max =  4.0;
    189             break;
    190         default:
    191             psAbort("invalid parameter %d for param max test", nParam);
    192         }
    193         if (params[nParam] > params_max) {
    194             params[nParam] = params_max;
    195             psTrace ("psModules.objects", 5, "params[nParam==%d] > params_max; %g v. %g",
    196                      nParam, params[nParam], params_max);
    197             return false;
    198         }
    199         return true;
    200     default:
     121      case PS_MINIMIZE_BETA_LIMIT: {
     122          psAssert(beta, "Require beta to limit beta");
     123          float limit = betaUse[nParam];
     124          if (nParam == PM_PAR_SXY) {
     125              limit *= q2;
     126          }
     127          if (fabs(beta[nParam]) > fabs(limit)) {
     128              beta[nParam] = (beta[nParam] > 0) ? fabs(limit) : -fabs(limit);
     129              psTrace("psModules.objects", 5, "|beta[nParam==%d]| > |beta_lim|; %g v. %g",
     130                      nParam, beta[nParam], limit);
     131              return false;
     132          }
     133          return true;
     134      }
     135      case PS_MINIMIZE_PARAM_MIN: {
     136          psAssert(params, "Require parameters to limit parameters");
     137          psAssert(paramsMinUse, "Require parameter limits to limit parameters");
     138          float limit = paramsMinUse[nParam];
     139          if (nParam == PM_PAR_SXY) {
     140              limit *= q2;
     141          }
     142          if (params[nParam] < limit) {
     143              params[nParam] = limit;
     144              psTrace("psModules.objects", 5, "params[nParam==%d] < params_min; %g v. %g",
     145                      nParam, params[nParam], limit);
     146              return false;
     147          }
     148          return true;
     149      }
     150      case PS_MINIMIZE_PARAM_MAX: {
     151          psAssert(params, "Require parameters to limit parameters");
     152          psAssert(paramsMaxUse, "Require parameter limits to limit parameters");
     153          float limit = paramsMaxUse[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_max; %g v. %g",
     160                      nParam, params[nParam], limit);
     161              return false;
     162          }
     163          return true;
     164      }
     165      default:
    201166        psAbort("invalid choice for limits");
    202167    }
    … …  
    204169    return false;
    205170}
     171
    206172
    207173// make an initial guess for parameters
    … …  
    456422}
    457423
     424
     425void PM_MODEL_SET_LIMITS(pmModelLimitsType type)
     426{
     427    switch (type) {
     428      case PM_MODEL_LIMITS_NONE:
     429        paramsMinUse = NULL;
     430        paramsMaxUse = NULL;
     431        limitsApply = true;
     432        break;
     433      case PM_MODEL_LIMITS_IGNORE:
     434        paramsMinUse = NULL;
     435        paramsMaxUse = NULL;
     436        limitsApply = false;
     437        break;
     438      case PM_MODEL_LIMITS_LAX:
     439        paramsMinUse = paramsMinLax;
     440        paramsMaxUse = paramsMaxLax;
     441        limitsApply = true;
     442        break;
     443      case PM_MODEL_LIMITS_STRICT:
     444        paramsMinUse = paramsMinStrict;
     445        paramsMaxUse = paramsMaxStrict;
     446        limitsApply = true;
     447        break;
     448      default:
     449        psAbort("Unrecognised model limits type: %x", type);
     450    }
     451    return;
     452}
     453
    458454# undef PM_MODEL_FUNC
    459455# undef PM_MODEL_FLUX
    … …  
    464460# undef PM_MODEL_PARAMS_FROM_PSF
    465461# undef PM_MODEL_FIT_STATUS
     462# undef PM_MODEL_SET_LIMITS
Note: See TracChangeset for help on using the changeset viewer.