IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jul 19, 2019, 5:33:23 PM (7 years ago)
Author:
eugene
Message:

avoid NANs; testing fitalt

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ohana.20190329/src/opihi/mana/deimos_fitalt.c

    r40663 r40838  
    11# include "data.h"
    22# include "deimos.h"
    3 
    4 # if (0)
    5 int deimos_fitalt (int argc, char **argv) {
    6   OHANA_UNUSED_PARAM(argc);
    7   OHANA_UNUSED_PARAM(argv);
    8   return FALSE;
    9 }
    10 
    11 # else
    123
    134/*
    … …  
    189
    1910// internal functions to fitobj
    20 static float      deimos_get_chisq (float *buffer, float *model, int Nx, int Ny, int row, int *Npts);
    21 double            deimos_LMM_update (float *data, float *model_ref, float *model_obj, float *model_sky, float *model_bck, float *dObj, float *dSky, float *dBck, int Nx, int Ny, int row, float lambda);
     11static float      deimos_get_chisq (float *buffer, float *vardata, float *model, int Nx, int Ny, int row, int *Npts);
     12double            deimos_LMM_update (float *data, float *vardata, float *model_ref, float *model_obj, float *model_sky, float *model_bck, float *dObj, float *dSky, float *dBck, int Nx, int Ny, int row, float lambda);
    2213static opihi_flt *deimos_make_test (opihi_flt *guess, float sigma, int row, int Nrow, float *dValue);
    2314opihi_flt        *deimos_subvector (opihi_flt *fullVec, int row, int Nrow);
    … …  
    2617
    2718  // input parameters:
    28   // * buffer      : 2D image of slit (full frame or cutout?)
    29   // * slit_trace  : spline fit of slit central x pos vs y-coord
    30   // * psf_trace   : spline fit of slit central x pos vs y-coord
    31   // * profile     : slit window profile (vector)
    32   // * PSF         : point-spread function vector (flux normalized, x-dir)
    33   // * stilt       : slit tilt response : 2D kernel?
     19  // * buffer         : 2D image of slit (full frame or cutout?)
     20  // * slit_trace_red : spline fit of slit central x pos vs y-coord
     21  // * slit_trace_blu : spline fit of slit central x pos vs y-coord
     22  // * psf_trace      : spline fit of slit central x pos vs y-coord
     23  // * profile        : slit window profile (vector)
     24  // * PSF            : point-spread function vector (flux normalized, x-dir)
     25  // * stilt          : slit tilt response : 2D kernel?
    3426
    3527  // in-out parameters:
    … …  
    4032  int N;
    4133
    42   Spline *slit_trace = NULL;
     34  Spline *slit_trace_red = NULL;
     35  Spline *slit_trace_blu = NULL;
    4336  Spline *psf_trace  = NULL;
    44   Vector *profile = NULL;
    45 
    46   Vector *obj     = NULL;
    47   Vector *sky     = NULL;
    48   Vector *bck     = NULL;
    49 
    50   Vector *psf     = NULL;
    51 
    52   Buffer *buffer  = NULL;
    53 
    54   float stilt = 0.0; // angle of the slit
     37
     38  Vector *profile    = NULL;
     39  Vector *psf        = NULL;
     40
     41  Vector *obj        = NULL;
     42  Vector *sky        = NULL;
     43  Vector *bck        = NULL;
     44
     45  Buffer *buffer     = NULL;
     46
     47  float stilt        = 0.0; // angle of the slit
     48
    5549  ohana_gaussdev_init ();
    5650
    … …  
    6559    remove_argument (N, &argc, argv);
    6660    VERBOSE = TRUE;
     61  }
     62
     63  int SKIPFIT = FALSE;
     64  if ((N = get_argument (argc, argv, "-skipfit"))) {
     65    remove_argument (N, &argc, argv);
     66    SKIPFIT = TRUE;
     67  }
     68
     69  // for a red vs blu spline, we need to specify the split point
     70  // XXX this is REALLY ad-hoc for Deimos.  not sure how to make
     71  // this more generic (need to define the ranges somewhere)
     72  int redlimit = 4096;
     73  if ((N = get_argument (argc, argv, "-redlimit"))) {
     74    remove_argument (N, &argc, argv);
     75    redlimit = atoi (argv[N]);
     76    remove_argument (N, &argc, argv);
    6777  }
    6878
    … …  
    108118
    109119  // Input parameters:
    110   if ((N = get_argument (argc, argv, "-slit-trace"))) {
    111     remove_argument (N, &argc, argv);
    112     if ((slit_trace = FindSpline (argv[N])) == NULL) return (FALSE);
    113     remove_argument (N, &argc, argv);
    114   } else { goto usage; }
     120  if ((N = get_argument (argc, argv, "-slit-trace-red"))) {
     121    remove_argument (N, &argc, argv);
     122    if ((slit_trace_red = FindSpline (argv[N])) == NULL) return (FALSE);
     123    remove_argument (N, &argc, argv);
     124  } else { goto usage; }
     125  if ((N = get_argument (argc, argv, "-slit-trace-blu"))) {
     126    remove_argument (N, &argc, argv);
     127    if ((slit_trace_blu = FindSpline (argv[N])) == NULL) return (FALSE);
     128    remove_argument (N, &argc, argv);
     129  } else { goto usage; }
     130
    115131  if ((N = get_argument (argc, argv, "-psf-trace"))) {
    116132    remove_argument (N, &argc, argv);
    … …  
    150166    remove_argument (N, &argc, argv);
    151167  } else { goto usage; }
     168
     169  // save the model
     170  Buffer *model = NULL;
     171  if ((N = get_argument (argc, argv, "-model"))) {
     172    remove_argument (N, &argc, argv);
     173    if ((model = SelectBuffer (argv[N], ANYBUFFER, TRUE)) == NULL) return (FALSE);
     174    remove_argument (N, &argc, argv);
     175  }
     176
     177  // save the model
     178  Buffer *varim = NULL;
     179  if ((N = get_argument (argc, argv, "-variance"))) {
     180    remove_argument (N, &argc, argv);
     181    if ((varim = SelectBuffer (argv[N], OLDBUFFER, TRUE)) == NULL) return (FALSE);
     182    remove_argument (N, &argc, argv);
     183  }
    152184
    153185  if (argc != 3) goto usage;
    … …  
    171203  float *bufVal = (float *) buffer->matrix.buffer;
    172204
     205  if (varim) {
     206    if (Nx != varim[0].matrix.Naxis[0]) { gprint (GP_ERR, "inconsistent sizes for image and varim\n"); return FALSE; }
     207    if (Ny != varim[0].matrix.Naxis[1]) { gprint (GP_ERR, "inconsistent sizes for image and varim\n"); return FALSE; }
     208  }
     209  float *varVal = varim ? (float *) varim->matrix.buffer : NULL;
     210
    173211  // obj, sky, bck must be consistent with data
    174212  if (Ny != obj->Nelements) {
    … …  
    183221    gprint (GP_ERR, "inconsistent wavelength scales (backgnd)\n");
    184222    return FALSE;
     223  }
     224
     225  if (model) {
     226    if (!CreateBuffer (model, Nx, Ny, -32, 1.0, 0.0)) { gprint (GP_ERR, "error generating output model buffer\n"); return FALSE; }
    185227  }
    186228
    … …  
    230272    // use scaled-noise to define test ranges
    231273    // we are going to try each row using a window Nrow wide to measure chisq:
    232     for (int row = dRow; row < Ny - dRow - 1; row++) {
     274    for (int row = dRow; !SKIPFIT && (row < Ny - dRow - 1); row++) {
    233275     
    234276      // I am using a numerical equivalent of Levenberg-Marquardt, but only for a single
    … …  
    251293      // current vector value:
    252294      // for deimos_make_model, row is the starting point of the subimage
    253       float *model_ref = deimos_make_model (objCurr, skyCurr, bckCurr, psf, profile, slit_trace, psf_trace, Nx, Nrow, row_0);
     295      float *model_ref = deimos_make_model (objCurr, skyCurr, bckCurr, psf, profile, slit_trace_red, slit_trace_blu, redlimit, psf_trace, Nx, Nrow, row_0);
    254296
    255297      // in the test functions below, dRow is the central pixel of the subset, e.g., obj[row] => objCurr[dRow]
    … …  
    259301      // delta obj (save dObj)
    260302      opihi_flt *objtest = deimos_make_test (objCurr, objNoise[row], dRow, Nrow, &dObj);
    261       float *model_obj = deimos_make_model (objtest, skyCurr, bckCurr, psf, profile, slit_trace, psf_trace, Nx, Nrow, row_0);
     303      float *model_obj = deimos_make_model (objtest, skyCurr, bckCurr, psf, profile, slit_trace_red, slit_trace_blu, redlimit, psf_trace, Nx, Nrow, row_0);
    262304
    263305      // delta sky (save dSky)
    264306      opihi_flt *skytest = deimos_make_test (skyCurr, skyNoise[row], dRow, Nrow, &dSky);
    265       float *model_sky = deimos_make_model (objCurr, skytest, bckCurr, psf, profile, slit_trace, psf_trace, Nx, Nrow, row_0);
     307      float *model_sky = deimos_make_model (objCurr, skytest, bckCurr, psf, profile, slit_trace_red, slit_trace_blu, redlimit, psf_trace, Nx, Nrow, row_0);
    266308
    267309      // delta bck (save dBck)
    268310      opihi_flt *bcktest = deimos_make_test (bckCurr, bckNoise[row], dRow, Nrow, &dBck);
    269       float *model_bck = deimos_make_model (objCurr, skyCurr, bcktest, psf, profile, slit_trace, psf_trace, Nx, Nrow, row_0);
     311      float *model_bck = deimos_make_model (objCurr, skyCurr, bcktest, psf, profile, slit_trace_red, slit_trace_blu, redlimit, psf_trace, Nx, Nrow, row_0);
    270312
    271313      // dObj, dSky, dBck carry the delta used derivatives in and the delta applied to the values out
    272       float chisq = deimos_LMM_update (bufVal, model_ref, model_obj, model_sky, model_bck, &dObj, &dSky, &dBck, Nx, Nrow, row_0, lambda);
     314      float chisq = deimos_LMM_update (bufVal, varVal, model_ref, model_obj, model_sky, model_bck, &dObj, &dSky, &dBck, Nx, Nrow, row_0, lambda);
     315
     316      if (FALSE && (row > 10) && (row < 20)) {
     317        fprintf (stderr, "%d %d : %f %f %f\n", iter, row, dObj, dSky, dBck);
     318      }
    273319
    274320      if (1) {
    … …  
    301347    }
    302348
    303     float *model_full = deimos_make_model (objVal, skyVal, bckVal, psf, profile, slit_trace, psf_trace, Nx, Ny, 0);
     349    float *model_full = deimos_make_model (objVal, skyVal, bckVal, psf, profile, slit_trace_red, slit_trace_blu, redlimit, psf_trace, Nx, Ny, 0);
    304350
    305351    int Npts = 0;
    306     float chisq_full  = deimos_get_chisq (bufVal, model_full, Nx, Ny, 0, &Npts);
     352    float chisq_full  = deimos_get_chisq (bufVal, varVal, model_full, Nx, Ny, 0, &Npts);
     353
     354    if (model) {
     355      free (model[0].matrix.buffer);
     356      model[0].matrix.buffer = (char *) model_full;
     357    } else {
     358      free (model_full);
     359    }
    307360
    308361    fprintf (stderr, "** Full chisq: %f (%d) : %f)\n", chisq_full, Npts, chisq_full / (1.0 * Npts));
    … …  
    312365
    313366 usage:
    314   gprint (GP_ERR, "USAGE: deimos fitobj (buffer) (Xref) -object vector -sky vector -backgnd vector -slit-trace spline -psf-trace spline -profile vector -psf vector -stilt (angle)\n");
     367  gprint (GP_ERR, "USAGE: deimos fitobj (buffer) (Xref) -object vector -sky vector -backgnd vector -slit-trace-red (spline) -slit-trace-blu (spline) -psf-trace (spline) -profile vector -psf vector -stilt (angle)\n");
    315368  return FALSE;
    316369}
    … …  
    351404
    352405
    353 double deimos_LMM_update (float *data, float *model_ref,
     406double deimos_LMM_update (float *data, float *vardata, float *model_ref,
    354407                          float *model_obj, float *model_sky, float *model_bck,
    355408                          float *dObj, float *dSky, float *dBck, int Nx, int Ny, int row, float lambda) {
    … …  
    365418  double dF_obj = 0.0, dF_sky = 0.0, dF_bck = 0.0;
    366419  double d2F_obj_obj = 0.0, d2F_obj_sky = 0.0, d2F_obj_bck = 0.0, d2F_sky_sky = 0.0, d2F_sky_bck = 0.0, d2F_bck_bck = 0.0;
    367 
    368   float W = 1.0;
    369420
    370421  for (int iy = 0; iy < Ny; iy++) {
    … …  
    380431      if (!isfinite(model_bck[pix])) continue;
    381432
     433      if (vardata && !isfinite(vardata[pix_data])) continue;
     434      float Var = vardata ? vardata[pix_data] : 1.0;
     435      float W = (fabs(Var) < 0.01) ? 100.0 : 1.0 / Var;
     436
    382437      float df = W * (data[pix_data] - model_ref[pix]);
    383438
    … …  
    431486}
    432487
    433 static float deimos_get_chisq (float *buffer, float *model, int Nx, int Ny, int row, int *Npts) {
     488static float deimos_get_chisq (float *buffer, float *vardata, float *model, int Nx, int Ny, int row, int *Npts) {
    434489
    435490  int npts = 0;
    … …  
    445500      if (!isfinite(buffer[pix_buffer])) continue;
    446501      if (!isfinite(model[pix_model])) continue;
    447       chisq += SQ(buffer[pix_buffer] - model[pix_model]);
     502
     503      if (vardata && !isfinite(vardata[pix_buffer])) continue;
     504      float Var = vardata ? vardata[pix_buffer] : 1.0;
     505      float W = (fabs(Var) < 0.01) ? 100.0 : 1.0 / Var;
     506
     507      float dchisq = SQ(buffer[pix_buffer] - model[pix_model]) * W;
     508      if (!isfinite(dchisq)) continue;
     509      chisq += dchisq;
    448510      npts ++;
    449511    }
    … …  
    452514  return chisq;
    453515}
    454 # endif
Note: See TracChangeset for help on using the changeset viewer.