IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 40838


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

avoid NANs; testing fitalt

Location:
branches/eam_branches/ohana.20190329/src/opihi/mana
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ohana.20190329/src/opihi/mana/Makefile

    r40753 r40838  
    2828$(SRC)/deimos_fitobj.$(ARCH).o \
    2929$(SRC)/deimos_fitalt.$(ARCH).o \
    30 $(SRC)/deimos_mkalt.$(ARCH).o \
    3130$(SRC)/deimos_mkmodel.$(ARCH).o \
    3231$(SRC)/deimos_arclines.$(ARCH).o \
    3332$(SRC)/deimos_fitprofile.$(ARCH).o \
    3433$(SRC)/findrowpeaks.$(ARCH).o
     34
     35# $(SRC)/deimos_mkalt.$(ARCH).o
    3536
    3637cmds = \

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

    r40817 r40838  
    5353  deimos_make_LSF_kernel (LSF, stilt, Nx);
    5454
     55  // FILE *Fout = NULL;
     56
    5557  // loop over the rows
    5658  for (int iy = 0; iy < Ny; iy++) {
    5759    float Fsum = 0;
    5860    float Ksum = 0;
     61
     62    // if (iy == 5802) Fout = fopen ("test.5802.dat", "w");
     63    // if (iy == 5815) Fout = fopen ("test.5815.dat", "w");
     64
    5965    // cross-correlation of buffer values and kernel values
    6066    for (int ix = 0; ix < Nx; ix++) {
    … …  
    6773        Fsum += Fin[ix + (iy + ky)*Nx] * kernel[ix + ko*Nx];
    6874        Ksum += kernel[ix + ko*Nx];
     75        // if ((iy == 5815) || (iy == 5802)) {
     76        //   fprintf (Fout, "%d %d %d : %f %f\n", iy, ix, ky, Fin[ix + (iy + ky)*Nx], kernel[ix + ko*Nx]);
     77        // }
    6978      }
    7079    }
     80
     81    // if ((iy == 5815) || (iy == 5802)) fclose (Fout);
     82
    7183    coord->elements.Flt[iy] = iy;
    7284    flux->elements.Flt[iy] = Fsum;

  • 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

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

    r40741 r40838  
    11# include "data.h"
    22# include "deimos.h"
     3
     4XXX: deprecated
    35
    46// Use this to test the deimos_make_object function used by deimos_fitobj

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

    r40663 r40838  
    8686// generate a subimage of the full model for the range row - row + Ny
    8787// obj,sky,bck are subset vectors corresponding to the range row to row+Ny
    88 float *deimos_make_model (opihi_flt *obj, opihi_flt *sky, opihi_flt *bck, Vector *psf, Vector *profile, Spline *slit_trace, Spline *psf_trace, int Nx, int Ny, int row) {
     88// the slit_trace_* and psf_trace refer to the full image starting at row
     89float *deimos_make_model (opihi_flt *obj, opihi_flt *sky, opihi_flt *bck, Vector *psf, Vector *profile, Spline *slit_trace_red, Spline *slit_trace_blu, float redlimit, Spline *psf_trace, int Nx, int Ny, int row) {
    8990
    9091  float *raw = deimos_make_straight_image (obj, sky, psf, psf_trace, Nx, Ny, row);
    … …  
    9293  deimos_apply_profile (profile, new, Nx, Ny);
    9394  deimos_add_background (bck, new, Nx, Ny, row);
    94   float *out = deimos_apply_trace (slit_trace, new, Nx, Ny, row);
     95  float *out = deimos_apply_trace (slit_trace_red, slit_trace_blu, redlimit, new, Nx, Ny, row);
    9596
    9697  free (raw);
    … …  
    113114
    114115  // Xref is a global supplied by the user (defaults to Nx/2 if < 0)
    115   int Xoff = (int)(Xref - Npsf);
     116  int XoffRef = (int)(Xref - Npsf);
    116117
    117118  // Nx : width of output window
    … …  
    123124   
    124125    // we are generating the image for just the row range row to row+Ny
     126    // note that obj & sky are subset vectors for just this range of pixels
    125127    opihi_flt objVy = obj[iy]; // if we shift to using the row offset: obj[iy + row]
    126128    opihi_flt skyVy = sky[iy]; // if we shift to using the row offset: sky[iy + row]
    … …  
    129131    int dxi = 0;
    130132    float dxf = 0.0, dxr = 1.0;
     133    int Xoff = XoffRef;
    131134   
    132135    if (psf_trace) {
    133       float dx = spline_apply_dbl (psf_trace->xk, psf_trace->yk, psf_trace->y2, psf_trace->Nknots, iy);
     136      // the psf_trace is referenced against the full image, so we need to add row to iy:
     137      float dx = spline_apply_dbl (psf_trace->xk, psf_trace->yk, psf_trace->y2, psf_trace->Nknots, iy + row);
    134138      dxi = floor(dx);
    135139      dxf = dx - dxi;  // -1.7 -> +0.3, -0.5 -> +0.5, +0.5 ->+0.5, +1.7 -> +0.7
    136140      dxr = 1 - dxf;
    137       Xref += dxi;
     141      Xoff += dxi;
     142      // fprintf (stderr, "%d: %f : %d\n", iy + row, dx, Xoff);
    138143    }
    139144
    … …  
    153158        if (doInterp) {
    154159          if ((n > 0) && (n < NpsfFull)) {
    155             value += objVy*psfV[n]*dxr + objVy*psfV[n-1]*dxf;
     160            if (isfinite(psfV[n]) && isfinite(psfV[n-1])) {
     161              value += objVy*psfV[n]*dxr + objVy*psfV[n-1]*dxf;
     162            }
    156163          }
    157164          if (n == 0) {
    158             value += objVy*psfV[n]*dxr;
     165            if (isfinite(psfV[n])) {
     166              value += objVy*psfV[n]*dxr;
     167            }
    159168          }
    160169        } else {
    161           value += objVy * psfV[n];
     170          if (isfinite(psfV[n])) {
     171            value += objVy * psfV[n];
     172          }
    162173        }
    163174      }
    … …  
    243254
    244255// shift pixels in x based on the trace
    245 float *deimos_apply_trace (Spline *slit_trace, float *input, int Nx, int Ny, int row) {
     256float *deimos_apply_trace (Spline *slit_trace_red, Spline *slit_trace_blu, float redlimit, float *input, int Nx, int Ny, int row) {
    246257
    247258  ALLOCATE_PTR (output, float, Nx*Ny);
    … …  
    251262    // evaluate the trace spline at this y-coord to find the x-coord offset of the profile center
    252263    // NOTE: spline is evaluated at full dispersion coordinate
    253     float dx = -spline_apply_dbl (slit_trace->xk, slit_trace->yk, slit_trace->y2, slit_trace->Nknots, iy + row);
     264    int iy_full = iy + row;
     265    Spline *slit_trace = (iy_full < redlimit) ? slit_trace_red : slit_trace_blu;
     266    float dx = -spline_apply_dbl (slit_trace->xk, slit_trace->yk, slit_trace->y2, slit_trace->Nknots, iy_full);
    254267     
    255268    // extract the integer pixel offset and the fractional offset

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

    r40817 r40838  
    223223      dxr = 1 - dxf;
    224224      Npof += dxi;
     225      // fprintf (stderr, "%d: %f : %d\n", iy, dx, Npof);
    225226      if (VERBOSE && (iy % 200 == 0)) { gprint (GP_ERR, "%d : %f : %d : %f %f\n", iy, dx, dxi, dxf, dxr); }
    226227    }
    … …  
    241242          if (doInterp) {
    242243            if ((n > 0) && (n < NpsfFull)) {
    243               value += objVy*psfV[n]*dxr + objVy*psfV[n-1]*dxf;
     244              if (isfinite(psfV[n]) && isfinite(psfV[n-1])) {
     245                value += objVy*psfV[n]*dxr + objVy*psfV[n-1]*dxf;
     246              }
    244247            }
    245248            if (n == 0) {
    246               value += objVy*psfV[n]*dxr;
     249              if (isfinite(psfV[n])) {
     250                value += objVy*psfV[n]*dxr;
     251              }
    247252            }
    248253          } else {
    249             value += objVy*psfV[n];
     254            if (isfinite(psfV[n])) {
     255              value += objVy*psfV[n];
     256            }
    250257          }
    251258        }
Note: See TracChangeset for help on using the changeset viewer.