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

    r23962 r25521  
    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_PS1_V1.h"
     30
    2231# define PM_MODEL_FUNC            pmModelFunc_PS1_V1
    2332# define PM_MODEL_FLUX            pmModelFlux_PS1_V1
    … …  
    2837# define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_PS1_V1
    2938# define PM_MODEL_FIT_STATUS      pmModelFitStatus_PS1_V1
     39# define PM_MODEL_SET_LIMITS      pmModelSetLimits_PS1_V1
    3040
    3141# define ALPHA   1.666
    3242# define ALPHA_M 0.666
     43
     44// Lax parameter limits
     45static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, -1.0 };
     46static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 };
     47
     48// Strict parameter limits
     49// k = PAR_7 < 0 is very undesirable (big divot in the middle)
     50static float paramsMinStrict[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, 0.0 };
     51static float paramsMaxStrict[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 };
     52
     53// Parameter limits to use
     54static float *paramsMinUse = NULL;
     55static float *paramsMaxUse = NULL;
     56static float betaUse[] = { 1000, 3e6, 5, 5, 1.0, 1.0, 0.5, 2.0 };
     57
    3358
    3459psF32 PM_MODEL_FUNC (psVector *deriv,
    … …  
    84109bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta)
    85110{
    86     float beta_lim = 0, params_min = 0, params_max = 0;
    87     float f1 = 0, f2 = 0, q1 = 0, q2 = 0;
     111    psAssert(nParam >= 0 && nParam <= PM_PAR_7, "Parameter index is out of bounds");
    88112
    89113    // we need to calculate the limits for SXY specially
     114    float q2 = NAN;
    90115    if (nParam == PM_PAR_SXY) {
    91         f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
    92         f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
    93         q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
     116        float f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
     117        float f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
     118        float q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
    94119        q1 = (q1 < 0.0) ? 0.0 : q1;
    95120        // if q1 < 0.0, f2 ~ f1, we have a very large axis ratio near 45deg..  Saturate at that
    96121        // angle and let f2,f1 fight it out
    97         q2  = 0.5*sqrt (q1);
     122        q2 = 0.5*sqrtf(q1);
    98123    }
    99124
    100125    switch (mode) {
    101     case PS_MINIMIZE_BETA_LIMIT:
    102         switch (nParam) {
    103         case PM_PAR_SKY:
    104             beta_lim = 1000;
    105             break;
    106         case PM_PAR_I0:
    107             beta_lim = 3e6;
    108             break;
    109         case PM_PAR_XPOS:
    110             beta_lim = 5;
    111             break;
    112         case PM_PAR_YPOS:
    113             beta_lim = 5;
    114             break;
    115         case PM_PAR_SXX:
    116             beta_lim = 1.0;
    117             break;
    118         case PM_PAR_SYY:
    119             beta_lim = 1.0;
    120             break;
    121         case PM_PAR_SXY:
    122             beta_lim =  0.5*q2;
    123             break;
    124         case PM_PAR_7:
    125             beta_lim = 2.0;
    126             break;
    127         default:
    128             psAbort("invalid parameter %d for beta test", nParam);
    129         }
    130         if (fabs(beta[nParam]) > fabs(beta_lim)) {
    131             beta[nParam] = (beta[nParam] > 0) ? fabs(beta_lim) : -fabs(beta_lim);
    132             psTrace ("psModules.objects", 5, "|beta[nParam==%d]| > |beta_lim|; %g v. %g",
    133                      nParam, beta[nParam], beta_lim);
    134             return false;
    135         }
    136         return true;
    137     case PS_MINIMIZE_PARAM_MIN:
    138         switch (nParam) {
    139         case PM_PAR_SKY:
    140             params_min = -1000;
    141             break;
    142         case PM_PAR_I0:
    143             params_min =   0.01;
    144             break;
    145         case PM_PAR_XPOS:
    146             params_min =  -100;
    147             break;
    148         case PM_PAR_YPOS:
    149             params_min =  -100;
    150             break;
    151         case PM_PAR_SXX:
    152             params_min =   0.5;
    153             break;
    154         case PM_PAR_SYY:
    155             params_min =   0.5;
    156             break;
    157         case PM_PAR_SXY:
    158             params_min =  -q2;
    159             break;
    160         case PM_PAR_7:
    161             params_min =  -1.0;
    162             break;
    163         default:
    164             psAbort("invalid parameter %d for param min test", nParam);
    165         }
    166         if (params[nParam] < params_min) {
    167             params[nParam] = params_min;
    168             psTrace ("psModules.objects", 5, "params[nParam==%d] < params_min; %g v. %g",
    169                      nParam, params[nParam], params_min);
    170             return false;
    171         }
    172         return true;
    173     case PS_MINIMIZE_PARAM_MAX:
    174         switch (nParam) {
    175         case PM_PAR_SKY:
    176             params_max =   1e5;
    177             break;
    178         case PM_PAR_I0:
    179             params_max =   1e8;
    180             break;
    181         case PM_PAR_XPOS:
    182             params_max =   1e4;
    183             break;
    184         case PM_PAR_YPOS:
    185             params_max =   1e4;
    186             break;
    187         case PM_PAR_SXX:
    188             params_max =   100;
    189             break;
    190         case PM_PAR_SYY:
    191             params_max =   100;
    192             break;
    193         case PM_PAR_SXY:
    194             params_max =  +q2;
    195             break;
    196         case PM_PAR_7:
    197             params_max =  20.0;
    198             break;
    199         default:
    200             psAbort("invalid parameter %d for param max test", nParam);
    201         }
    202         if (params[nParam] > params_max) {
    203             params[nParam] = params_max;
    204             psTrace ("psModules.objects", 5, "params[nParam==%d] > params_max; %g v. %g",
    205                      nParam, params[nParam], params_max);
    206             return false;
    207         }
    208         return true;
    209     default:
     126      case PS_MINIMIZE_BETA_LIMIT: {
     127          psAssert(beta, "Require beta to limit beta");
     128          float limit = betaUse[nParam];
     129          if (nParam == PM_PAR_SXY) {
     130              limit *= q2;
     131          }
     132          if (fabs(beta[nParam]) > fabs(limit)) {
     133              beta[nParam] = (beta[nParam] > 0) ? fabs(limit) : -fabs(limit);
     134              psTrace("psModules.objects", 5, "|beta[nParam==%d]| > |beta_lim|; %g v. %g",
     135                      nParam, beta[nParam], limit);
     136              return false;
     137          }
     138          return true;
     139      }
     140      case PS_MINIMIZE_PARAM_MIN: {
     141          psAssert(params, "Require parameters to limit parameters");
     142          psAssert(paramsMinUse, "Require parameter limits to limit parameters");
     143          float limit = paramsMinUse[nParam];
     144          if (nParam == PM_PAR_SXY) {
     145              limit *= q2;
     146          }
     147          if (params[nParam] < limit) {
     148              params[nParam] = limit;
     149              psTrace("psModules.objects", 5, "params[nParam==%d] < params_min; %g v. %g",
     150                      nParam, params[nParam], limit);
     151              return false;
     152          }
     153          return true;
     154      }
     155      case PS_MINIMIZE_PARAM_MAX: {
     156          psAssert(params, "Require parameters to limit parameters");
     157          psAssert(paramsMaxUse, "Require parameter limits to limit parameters");
     158          float limit = paramsMaxUse[nParam];
     159          if (nParam == PM_PAR_SXY) {
     160              limit *= q2;
     161          }
     162          if (params[nParam] > limit) {
     163              params[nParam] = limit;
     164              psTrace("psModules.objects", 5, "params[nParam==%d] > params_max; %g v. %g",
     165                      nParam, params[nParam], limit);
     166              return false;
     167          }
     168          return true;
     169      }
     170      default:
    210171        psAbort("invalid choice for limits");
    211172    }
    … …  
    468429}
    469430
     431
     432void PM_MODEL_SET_LIMITS(pmModelLimitsType type)
     433{
     434    switch (type) {
     435      case PM_MODEL_LIMITS_NONE:
     436        paramsMinUse = NULL;
     437        paramsMaxUse = NULL;
     438      case PM_MODEL_LIMITS_LAX:
     439        paramsMinUse = paramsMinLax;
     440        paramsMaxUse = paramsMaxLax;
     441        break;
     442      case PM_MODEL_LIMITS_STRICT:
     443        paramsMinUse = paramsMinStrict;
     444        paramsMaxUse = paramsMaxStrict;
     445        break;
     446      default:
     447        psAbort("Unrecognised model limits type: %x", type);
     448    }
     449    return;
     450}
     451
     452
    470453# undef PM_MODEL_FUNC
    471454# undef PM_MODEL_FLUX
    … …  
    476459# undef PM_MODEL_PARAMS_FROM_PSF
    477460# undef PM_MODEL_FIT_STATUS
     461# undef PM_MODEL_SET_LIMITS
    478462# undef ALPHA
    479463# undef ALPHA_M
Note: See TracChangeset for help on using the changeset viewer.