IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jun 22, 2012, 3:33:32 PM (14 years ago)
Author:
eugene
Message:

read the time table in psTime.c; updates to the variance model tests; better trace options for output masks; convert to new PM_SOURCE_PHOTFIT version of weighting

Location:
branches/eam_branches/ipp-20120601/psphot
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20120601/psphot/src/psphot.h

    r34049 r34053  
    6161bool            psphotSetMaskAndVarianceReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
    6262
     63bool            psphotUpdateVariance (pmConfig *config, const pmFPAview *view, const char *filerule);
     64bool            psphotUpdateVarianceReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
     65
    6366bool            psphotModelBackground (pmConfig *config, const pmFPAview *view, const char *filerule);
    6467bool            psphotModelBackgroundReadoutFileIndex (pmConfig *config, const pmFPAview *view, const char *filerule, int index);
    … …  
    99102
    100103bool            psphotFitSourcesLinear (pmConfig *config, const pmFPAview *view, const char *filerule, bool final);
    101 
    102 # if (HAVE_MODEL_VAR)
    103104bool            psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final, pmSourceFitVarMode fitVarMode);
    104 # else
    105 bool            psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final);
    106 # endif
    107105
    108106bool            psphotSourceSize (pmConfig *config, const pmFPAview *view, const char *filerule, bool getPSFsize);

  • branches/eam_branches/ipp-20120601/psphot/src/psphotEfficiency.c

    r34049 r34053  
    420420
    421421    // psphotFitSourcesLinearReadout subtracts the model fits
    422 # if (HAVE_MODEL_VAR)
    423422    if (!psphotFitSourcesLinearReadout(recipe, readout, fakeSourcesAll, psf, true, PM_SOURCE_PHOTFIT_CONST)) {
    424 # else
    425     if (!psphotFitSourcesLinearReadout(recipe, readout, fakeSourcesAll, psf, true)) {
    426 # endif
    427423        psError(PS_ERR_UNKNOWN, false, "Unable to perform linear fit on fake sources.");
    428424        psFree(fakeSources);

  • branches/eam_branches/ipp-20120601/psphot/src/psphotFitSourcesLinear.c

    r34049 r34053  
    1212static bool SetBorderMatrixElements (psSparseBorder *border, pmReadout *readout, psArray *sources, bool constant_weights, int SKY_FIT_ORDER, psImageMaskType markVal);
    1313
    14 # if (HAVE_MODEL_VAR)
    1514bool psphotGenerateModelVariance (pmConfig *config, const pmFPAview *view, pmFPAfile *file, int index, psMetadata *recipe, pmReadout *readout, psArray *sources);
    1615pmSourceFitVarMode psphotGetFitVarMode (psMetadata *recipe);
    1716bool psphotFreeModelVariance (pmReadout *readout, psArray *sources);
    18 # endif
    1917
    2018// for now, let's store the detections on the readout->analysis for each readout
    … …  
    3028    assert (recipe);
    3129
    32 # if (HAVE_MODEL_VAR)
    3330    pmSourceFitVarMode fitVarMode = psphotGetFitVarMode (recipe);
    3431    if (!fitVarMode) {
    … …  
    4037    // do a single pass.
    4138    pmSourceFitVarMode fitVarModePass1 = (fitVarMode == PM_SOURCE_PHOTFIT_MODEL_VAR) ? PM_SOURCE_PHOTFIT_CONST : fitVarMode;
    42 # endif
    4339
    4440    int num = psphotFileruleCount(config, filerule);
    … …  
    6864        psAssert (psf, "missing psf?");
    6965
    70 # if (HAVE_MODEL_VAR)
    71         if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final, fitVarModePass1))
    72 # else
    73         if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final))
    74 # endif
    75         {
     66        if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final, fitVarModePass1)) {
    7667            psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for %s entry %d", filerule, i);
    7768            return false;
    7869        }
    7970
    80 # if (HAVE_MODEL_VAR)
    8171        // the MODEL_VAR weighting scheme requires knowledge of the model fluxes to generate the variance
    8272        // after we have determined the initial set of fits, then we can generate the variance image and
    … …  
    10494            }
    10595        }
    106 # endif
    10796
    10897        psphotVisualShowResidualImage (readout, (num > 0));
    … …  
    113102}
    114103
    115 # if (HAVE_MODEL_VAR)
    116104// look up the fit variance mode from the recipe; older recipes do not have the value
    117105// 'LINEAR_FIT_VARIANCE_MODE'; in those cases, look for 'CONSTANT_PHOTOMETRIC_WEIGHTS' as a boolean and
    … …  
    136124        return PM_SOURCE_PHOTFIT_IMAGE_VAR;
    137125    }
     126    if (!strcasecmp(fitVarModeString, "SKY")   || !strcasecmp(fitVarModeString, "MODEL_SKY")) {
     127        return PM_SOURCE_PHOTFIT_MODEL_SKY;
     128    }
    138129    if (!strcasecmp(fitVarModeString, "MODEL") || !strcasecmp(fitVarModeString, "MODEL_VAR")) {
    139130        return PM_SOURCE_PHOTFIT_MODEL_VAR;
    … …  
    142133    return PM_SOURCE_PHOTFIT_NONE;
    143134}
    144 # endif
    145 
    146 # if (HAVE_MODEL_VAR)
    147 bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final, pmSourceFitVarMode fitVarMode)
    148 # else
    149 bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final)
    150 # endif
    151 {
     135
     136bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final, pmSourceFitVarMode fitVarMode) {
    152137    bool status;
    153138    float x;
    … …  
    185170    if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined");
    186171
    187 # if (!HAVE_MODEL_VAR)
    188     bool CONSTANT_PHOTOMETRIC_WEIGHTS =
    189         psMetadataLookupBool(&status, recipe, "CONSTANT_PHOTOMETRIC_WEIGHTS");
    190     if (!status) {
    191         psAbort("You must provide a value for the BOOL recipe CONSTANT_PHOTOMETRIC_WEIGHTS");
    192     }
    193 # endif
    194172    int SKY_FIT_ORDER = psMetadataLookupS32(&status, recipe, "SKY_FIT_ORDER");
    195173    if (!status) {
    … …  
    319297    psSparseBorder *border = psSparseBorderAlloc (sparse, nBorder);
    320298
    321 # if (HAVE_MODEL_VAR)
    322299    // if fitVarMode is MODEL_VAR, then we need to generate the model image variance
    323300    // XXX we have two possibilities here:
    … …  
    327304
    328305    // 2) do a single pass, and use the model guess to define the model variance (but do I trust the Model Guess?)
    329 # endif
    330306
    331307    // fill out the sparse matrix elements and border elements (B)
    … …  
    336312
    337313        // diagonal elements of the sparse matrix (auto-cross-product)
    338 # if (HAVE_MODEL_VAR)
    339314        f = pmSourceModelDotModel (SRCi, SRCi, fitVarMode, covarFactor, maskVal);
    340 # else
    341         f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
    342 # endif
    343315        psSparseMatrixElement (sparse, i, i, f);
    344316
    345 # if (HAVE_MODEL_VAR)
    346         // if we have used CONSTANT errors, then we need to calculate the value of the parameter error
    347         if (fitVarMode != PM_SOURCE_PHOTFIT_IMAGE_VAR) {
     317        // if we have used CONSTANT errors, then we need to re-calculate the value of the parameter error
     318        if (fitVarMode == PM_SOURCE_PHOTFIT_CONST) {
    348319            float var = pmSourceModelDotModel (SRCi, SRCi, PM_SOURCE_PHOTFIT_IMAGE_VAR, covarFactor, maskVal);
    349320            errors->data.F32[i] = 1.0 / sqrt(var);
    … …  
    351322            errors->data.F32[i] = 1.0 / sqrt(f);
    352323        }
    353 # else
    354         // the formal error depends on the weighting scheme
    355         if (CONSTANT_PHOTOMETRIC_WEIGHTS) {
    356             float var = pmSourceModelDotModel (SRCi, SRCi, false, covarFactor, maskVal);
    357             errors->data.F32[i] = 1.0 / sqrt(var);
    358         } else {
    359             errors->data.F32[i] = 1.0 / sqrt(f);
    360         }
    361 # endif
    362324
    363325        // find the image x model value
    364 # if (HAVE_MODEL_VAR)
    365326        f = pmSourceDataDotModel (SRCi, SRCi, fitVarMode, covarFactor, maskVal);
    366 # else
    367         f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
    368 # endif
    369327        psSparseVectorElement (sparse, i, f);
    370328
    … …  
    372330        switch (SKY_FIT_ORDER) {
    373331          case 1:
    374 # if (HAVE_MODEL_VAR)
    375332            f = pmSourceModelWeight (SRCi, 1, fitVarMode, covarFactor, maskVal);
    376333            psSparseBorderElementB (border, i, 1, f);
    377334            f = pmSourceModelWeight (SRCi, 2, fitVarMode, covarFactor, maskVal);
    378335            psSparseBorderElementB (border, i, 2, f);
    379 # else
    380             f = pmSourceModelWeight (SRCi, 1, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
    381             psSparseBorderElementB (border, i, 1, f);
    382             f = pmSourceModelWeight (SRCi, 2, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
    383             psSparseBorderElementB (border, i, 2, f);
    384 # endif
    385336
    386337          case 0:
    387 # if (HAVE_MODEL_VAR)
    388338            f = pmSourceModelWeight (SRCi, 0, fitVarMode, covarFactor, maskVal);
    389 # else
    390             f = pmSourceModelWeight (SRCi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
    391 # endif
    392339            psSparseBorderElementB (border, i, 0, f);
    393340            break;
    … …  
    409356
    410357            // got an overlap; calculate cross-product and add to output array
    411 # if (HAVE_MODEL_VAR)
    412358            f = pmSourceModelDotModel (SRCi, SRCj, fitVarMode, covarFactor, maskVal);
    413 # else
    414             f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
    415 # endif
    416359            psSparseMatrixElement (sparse, j, i, f);
    417360        }
    … …  
    451394
    452395    // set the sky, sky_x, sky_y components of border matrix
    453 # if (HAVE_MODEL_VAR)
    454396    SetBorderMatrixElements (border, readout, fitSources, (fitVarMode == PM_SOURCE_PHOTFIT_CONST), SKY_FIT_ORDER, markVal);
    455 # else
    456     SetBorderMatrixElements (border, readout, fitSources, CONSTANT_PHOTOMETRIC_WEIGHTS, SKY_FIT_ORDER, markVal);
    457 # endif
    458397
    459398    psSparseConstraint constraint;
    … …  
    614553}
    615554
    616 # if (HAVE_MODEL_VAR)
    617555bool psphotModelBackgroundReadout(psImage *model,  // Model image
    618556                                  psImage *modelStdev, // Model stdev image
    … …  
    650588    psImage *varModelStdev = psImageAlloc(backBinning->nXruff, backBinning->nYruff, PS_TYPE_F32); // Background model
    651589
     590    // generate an image of the mean variance image in DN
    652591    if (!psphotModelBackgroundReadout(varModel, varModelStdev, NULL, readout, backBinning, config, true)) {
    653592        psError(PS_ERR_UNKNOWN, false, "Unable to generate background model");
    … …  
    668607    psFree (varModelStdev);
    669608
    670     // XXX for a test:
    671     psphotSaveImage (NULL, modelVar, "model.bck.fits");
     609    float gain = 1.0;  // accept 1.0 as a default since it is not critical to the analysis
     610    pmCell *cell = readout->parent; // The parent cell
     611    if (cell) {
     612      gain = psMetadataLookupF32(&status, cell->concepts, "CELL.GAIN"); // Cell gain
     613      if (!status) {
     614        gain = 1.0;           // set note above
     615      }
     616    }
     617    if (gain > 2.0) { /* warn? */ }
    672618
    673619    // insert all of the source models
    … …  
    690636
    691637        // add the source model to the model variance image
     638        // XXX note that this should be added with gain applied
     639        // var_DN = flux_DN / gain [e/DN]
     640        // to do this requires an API upgrade...
    692641        pmSourceAdd (source, PM_MODEL_OP_MODELVAR, maskVal);
    693642    }
    694 
    695     // XXX for a test:
    696     psphotSaveImage (NULL, modelVar, "model.var.fits");
    697     psphotSaveImage (NULL, readout->variance, "image.var.fits");
    698643
    699644    // we save the model variance for future reference
    … …  
    724669    return true;
    725670}
    726 # endif

  • branches/eam_branches/ipp-20120601/psphot/src/psphotMaskReadout.c

    r29936 r34053  
    9494
    9595    // test output of files at this stage
    96     if (psTraceGetLevel("psphot") >= 5) {
     96    if (psTraceGetLevel("psphot.imsave") >= 5) {
    9797        psphotSaveImage (NULL, readout->image,  "image.fits");
    9898        psphotSaveImage (NULL, readout->mask,   "mask.fits");
    … …  
    105105    return true;
    106106}
     107
     108// XXX this function and support below was created to test the theory that the faint-end
     109// bias results from the Poisson variation of the background pixels.  This is NOT the
     110// case.  Using the code below maintains the faint-end bias.
     111bool psphotUpdateVariance (pmConfig *config, const pmFPAview *view, const char *filerule) {
     112
     113    bool status = false;
     114
     115    // select the appropriate recipe information
     116    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     117    psAssert (recipe, "missing recipe?");
     118
     119    int num = psphotFileruleCount(config, filerule);
     120
     121    // loop over the available readouts
     122    for (int i = 0; i < num; i++) {
     123
     124        // Generate the mask and weight images, including the user-defined analysis region of interest
     125        if (!psphotUpdateVarianceReadout (config, view, filerule, i, recipe)) {
     126            psError (PSPHOT_ERR_CONFIG, false, "failed to generate mask for %s entry %d", filerule, i);
     127            return false;
     128        }
     129    }
     130    return true;
     131}
     132
     133// determine the mean variance image (equivalent to the background model, but for the variance image)
     134// set the variance image to the MAX(input, mean)
     135bool psphotUpdateVarianceReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) {
     136
     137    bool status;
     138
     139    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest
     140    psAssert (file, "missing file?");
     141
     142    // find the currently selected readout
     143    pmReadout  *readout = pmFPAviewThisReadout (view, file->fpa);
     144    psAssert (readout, "missing readout?");
     145
     146    pmSourceFitVarMode varMode = psphotGetFitVarMode (recipe);
     147    if (varMode == PM_SOURCE_PHOTFIT_NONE) {
     148      psError (PSPHOT_ERR_CONFIG, false, "need valid LINEAR_FIT_VARIANCE_MODE");
     149      return false;
     150    }
     151
     152    // make this an option via the recipe
     153    if (varMode != PM_SOURCE_PHOTFIT_MODEL_SKY) return true;
     154
     155    // create a model variance image (full-scale image to take result of psImageUnbin below)
     156    psImage *modelVar = psImageCopy (NULL, readout->variance, PS_TYPE_F32);
     157
     158    // find the binning information
     159    psImageBinning *backBinning = psphotBackgroundBinning (modelVar, config);
     160    assert (backBinning);
     161   
     162    psImage *varModel = psImageAlloc(backBinning->nXruff, backBinning->nYruff, PS_TYPE_F32); // Background model
     163    psImage *varModelStdev = psImageAlloc(backBinning->nXruff, backBinning->nYruff, PS_TYPE_F32); // Background model
     164
     165    if (!psphotModelBackgroundReadout(varModel, varModelStdev, NULL, readout, backBinning, config, true)) {
     166        psError(PS_ERR_UNKNOWN, false, "Unable to generate background model");
     167        psFree (varModel);
     168        psFree (varModelStdev);
     169        return false;
     170    }
     171
     172    // linear interpolation to full-scale
     173    if (!psImageUnbin (modelVar, varModel, backBinning)) {
     174        psError (PSPHOT_ERR_PROG, true, "inconsistent sizes for unbinning");
     175        psFree (varModel);
     176        psFree (varModelStdev);
     177        return false;
     178    }
     179
     180    // XXX save these?
     181    psFree (varModel);
     182    psFree (varModelStdev);
     183
     184    psImage *im = readout->image;
     185    psImage *wt = readout->variance;
     186    for (int j = 0; j < im->numRows; j++) {
     187      for (int i = 0; i < im->numCols; i++) {
     188        if (!isfinite(im->data.F32[j][i])) continue;
     189        if (!isfinite(wt->data.F32[j][i])) continue;
     190        // XXX for a test, make variance constant wt->data.F32[j][i] = PS_MAX(wt->data.F32[j][i], modelVar->data.F32[j][i]);
     191        wt->data.F32[j][i] = modelVar->data.F32[j][i];
     192      }
     193    }
     194
     195    // test output of files at this stage
     196    if (psTraceGetLevel("psphot.imsave") >= 5) {
     197        psphotSaveImage (NULL, readout->image,  "image.varsky.fits");
     198        psphotSaveImage (NULL, readout->mask,   "mask.varsky.fits");
     199        psphotSaveImage (NULL, readout->variance, "variance.varsky.fits");
     200    }
     201
     202    psFree (modelVar);
     203
     204    return true;
     205}

  • branches/eam_branches/ipp-20120601/psphot/test/tap_psphot_varmodel.pro

    r33963 r34053  
    22# -*-sh-*-
    33
    4 # $KAPA = kapa -noX
    5 
    6 # PSF.CONVOLVE : if true, we insert delta functions (and optionally
    7 #                galaxies) and smooth the image with the psf model
    8 #                (uses a GAUSS regardless of the model). Note that
    9 #                PSF.CONVOLVE = T is faster than F, but (a) only
    10 #                allows Gauss models and (b) only yields quantized
    11 #                locations
    12 
    13 # config for ppImage to generate chip, mask, weight
    14 $ppImageConfig = -recipe PPIMAGE PPIMAGE_N
    15 $ppImageConfig = $ppImageConfig -Db BACKGROUND T
    16 $ppImageConfig = $ppImageConfig -Db CHIP.FITS T
    17 $ppImageConfig = $ppImageConfig -Db CHIP.MASK.FITS T
    18 $ppImageConfig = $ppImageConfig -Db CHIP.VARIANCE.FITS T
    19 $ppImageConfig = $ppImageConfig -Db BASE.FITS F
    20 $ppImageConfig = $ppImageConfig -Db VARIANCE.BUILD T
    21 $ppImageConfig = $ppImageConfig -Db PHOTOM F
    22 
    23 # basic options for the these images (filter, location, obstype)
    24 $BaseOptions = -type OBJECT -filter r -skymags 20.86 -ra 270.70 -dec -23.70 -pa 0.0
    25 $BaseOptions = $BaseOptions -Df PSASTRO:DVO.GETSTAR.MAX.RHO 50000.0
    26 
    27 # create an image with fake sources and insert the resulting cmf file into a dvodb
    28 $RefConfig = -camera SIMTEST
    29 $RefConfig = $RefConfig -recipe PPSIM STACKTEST.MAKE
    30 $RefConfig = $RefConfig -D PSASTRO:PSASTRO.CATDIR catdir.ref
    31 $RefConfig = $RefConfig -Db PSF.CONVOLVE F
    32 
    33 # options for the reference image
    34 $RefOptions = $BaseOptions -exptime 100.0
    35 $RefOptions = $RefOptions -seeing 1.0
    36 $RefOptions = $RefOptions -D PSF.MODEL PS_MODEL_GAUSS
    37 $RefOptions = $RefOptions -Df STARS.DENSITY 10.0
    38 $RefOptions = $RefOptions -Df STARS.SIGMA.LIM 0.5
    39 
    40 # basic config for ppSim with randomly distributed stars and NO galaxies
    41 $RealConfig = -camera SIMTEST
    42 $RealConfig = $RealConfig -recipe PPSIM STACKTEST.RUN
    43 $RealConfig = $RealConfig -D PSASTRO:PSASTRO.CATDIR catdir.ref
    44 $RealConfig = $RealConfig -Db STARS.FAKE F
    45 $RealConfig = $RealConfig -Db STARS.REAL T
    46 $RealConfig = $RealConfig -Db MATCH.DENSITY F
    47 $RealConfig = $RealConfig -Db PSF.CONVOLVE F
    48 $RealConfig = $RealConfig -Df STARS.DENSITY 10.0
    49 $RealConfig = $RealConfig -Df STARS.SIGMA.LIM 2.5
    50 $RealConfig = $RealConfig -Db GALAXY.FAKE F
    51 $RealConfig = $RealConfig -Db GALAXY.GRID F
    52 
    53 # options for the repeated images
    54 $RealOptions = $BaseOptions -exptime 30.0
    55  
    56 $ExtraOptions = -D PSF.MODEL PS_MODEL_GAUSS
    57 
    58 # sample alternate options:
    59 # $ppSimOptions = $FakeOptions -D PSF.MODEL PS_MODEL_PS1_V1
    60 # $ppSimOptions = $FakeOptions -Df PSF.ARATIO 1.2
    61 # $ppSimOptions = $FakeOptions -Df PSF.THETA +30.0
    62 # $ppSimOptions = $FakeOptions -D PSF.MODEL PS_MODEL_GAUSS
    63 
    64 list fwhm
    65  1.0
    66  1.1
    67  1.2
    68  1.5
    69 end
     4# This script includes a set of tests to demonstrate the dependence of the faint-end bias on the weighting scheme
     5# We have 3 weighting schemes :
     6# CONSTANT -- per-pixel weight is fixed (disadvantage: lower S/N, especially for higher sky?)
     7# IMAGE_VAR - per-pixel weight is Poisson from image (disadvantage: faint-end bias)
     8# MODEL_VAR - per-pixel weight is Poisson from model (disadvantage: 2 linear-fit passes)
     9
     10# Functions:
     11#
     12# init  : initialize variables
     13# mkref : generate a fake reference catalog in a DVO database
     14# mkexp : generate a fake exposure from fake catalog and detrend it (saves basename.in.cmf & basename.fits)
     15# runphot : run psphot on an exposure (saves basename.cmf)
     16# ckchip.mags : generate a set of summary plots for the psphot analysis
     17
     18# I would like to produce a grid of tests:
     19# sky (bright, middle, dark)
     20# fwhm (1.0, 1.3, 1.6)
     21# PSF_MODEL input (GAUSS PS1_V1)
     22# PSF_MODEL apply (GAUSS PS1_V1)
     23# variance mode: CONSTANT, IMAGE_VAR, MODEL_VAR
    7024
    7125macro go
    72 
    7326  mkdir test
    7427
    75   $ExtraOptions = -D PSF.MODEL PS_MODEL_GAUSS
    76   mkexp test/image.00 1.0
    77   $ExtraOptions = -D PSF.MODEL PS_MODEL_PGAUSS
    78   mkexp test/image.01 1.0
    79   $ExtraOptions = -D PSF.MODEL PS_MODEL_PS1_V1
    80   mkexp test/image.02 1.0
     28  local sky fwhm psf_in psf_out
     29
     30  # mkref creates refimage.* and catdir.ref
     31  file catdir.ref found
     32  if (not($found))
     33    mkref
     34  end
     35
     36  $KAPA = kapa -noX
     37  if (not($?RefConfig)) init
     38
     39  foreach sky 19.0 20.0 21.0
     40    foreach fwhm 1.0 1.3 1.6
     41      foreach psf_in GAUSS PS1_V1
     42
     43        sprintf name "test/test.%02d.%02d.%s" $sky {10*$fwhm} $psf_in
     44        # mkexp $name $sky $fwhm $psf_in
     45       
     46        foreach psf_out GAUSS PS1_V1
     47          foreach mode CONSTANT IMAGE_VAR MODEL_VAR
     48            # runphot $name $name.$psf_out.$mode "-D LINEAR_FIT_VARIANCE_MODE $mode -D PSF_MODEL PS_MODEL_$psf_out"
     49            ckchip.mags $name.in.cmf $name.$psf_out.$mode.cmf $name.$psf_out.$mode 0.0
     50          end
     51        end
     52      end
     53    end
     54  end
    8155end
    8256
    8357# create a reference database of fake stars to be used by ppSim below
    8458macro mkref
     59  if (not($?RefConfig)) init
     60
    8561  exec rm -rf catdir.ref
    8662  exec rm -f refimage.fits
    … …  
    10076# create a realistic distribution of fake stars, GAUSS PSF
    10177macro mkexp
     78  if ($0 != 5)
     79    echo "USAGE: mkexp basename sky fwhm psf_model"
     80    break
     81  end
     82
     83  local fwhm basename psf_model
     84  $basename = $1
     85  $sky = $2
     86  $fwhm = $3
     87  $psf_model = $4
     88
     89  $ExtraOptions = -D PSF.MODEL PS_MODEL_$psf_model
     90
     91  # create the raw image
     92  echo ppSim -seeing $fwhm -skymags $sky -nx 3000 -ny 3000 $RealConfig $RealOptions $ExtraOptions $basename
     93  exec ppSim -seeing $fwhm -skymags $sky -nx 3000 -ny 3000 $RealConfig $RealOptions $ExtraOptions $basename
     94  exec /bin/mv -f $basename.cmf $basename.in.cmf
     95
     96  # create the chip output
     97  echo ppImage $ppImageConfig -file $basename.fits $basename
     98  exec ppImage $ppImageConfig -file $basename.fits $basename
     99end
     100
     101# create a realistic distribution of fake stars, GAUSS PSF
     102macro mkexp.deep
    102103  if ($0 != 3)
    103104    echo "USAGE: mkexp basename fwhm"
    … …  
    109110  $fwhm = $2
    110111
     112  $RealOptionsDeep = $BaseOptions -exptime 100
     113
    111114  # create the raw image
    112   echo ppSim -seeing $fwhm -nx 3000 -ny 3000 $RealConfig $RealOptions $ExtraOptions $basename
    113   exec ppSim -seeing $fwhm -nx 3000 -ny 3000 $RealConfig $RealOptions $ExtraOptions $basename
     115  echo ppSim -seeing $fwhm -nx 3000 -ny 3000 $RealConfig $RealOptionsDeep $ExtraOptions $basename
     116  exec ppSim -seeing $fwhm -nx 3000 -ny 3000 $RealConfig $RealOptionsDeep $ExtraOptions $basename
     117  exec /bin/mv -f $basename.cmf $basename.in.cmf
     118
     119  # create the chip output
     120  echo ppImage $ppImageConfig -file $basename.fits $basename
     121  exec ppImage $ppImageConfig -file $basename.fits $basename
     122end
     123
     124# create a realistic distribution of fake stars, GAUSS PSF
     125macro mkexp.bright
     126  if ($0 != 3)
     127    echo "USAGE: mkexp basename fwhm"
     128    break
     129  end
     130
     131  local fwhm basename
     132  $basename = $1
     133  $fwhm = $2
     134
     135  # basic options for the these images (filter, location, obstype)
     136  $BaseOptions = -type OBJECT -filter r -skymags 20.0 -ra 270.70 -dec -23.70 -pa 0.0
     137  $BaseOptions = $BaseOptions -Df PSASTRO:DVO.GETSTAR.MAX.RHO 50000.0
     138  $RealOptionsDeep = $BaseOptions -exptime 100
     139
     140  # create the raw image
     141  echo ppSim -seeing $fwhm -nx 3000 -ny 3000 $RealConfig $RealOptionsDeep $ExtraOptions $basename
     142  exec ppSim -seeing $fwhm -nx 3000 -ny 3000 $RealConfig $RealOptionsDeep $ExtraOptions $basename
    114143  exec /bin/mv -f $basename.cmf $basename.in.cmf
    115144
    … …  
    133162  echo psphot -threads 4 -file $basename.ch.fits -mask $basename.ch.mk.fits -variance $basename.ch.wt.fits $outname $options
    134163  exec psphot -threads 4 -file $basename.ch.fits -mask $basename.ch.mk.fits -variance $basename.ch.wt.fits $outname $options
     164end
     165
     166# compare two cmf files with extname Chip.psf
     167# things to compare:
     168# * completeness (which sources in (1) are not detected in (2)
     169# * positions (X_PSF, Y_PSF)
     170# * instrumental psf mags
     171# * position errors (no input errors; use a model?)
     172# * measured FWHM?
     173# * kron mags (fluxes)
     174# * etc, etc
     175macro ckchip.mags
     176  if ($0 != 5)
     177    echo "USAGE: ckchip.mags (raw) (out) (output) (zpt_off)"
     178    break
     179  end
     180
     181list pairs
     182  PSF_INST_MAG_out      M_raw                 PSF_INST_MAG_raw 2 -0.21 0.21 V
     183  AP_MAG_out            M_raw                 PSF_INST_MAG_raw 2 -0.21 0.21 V
     184end
     185
     186  load.cmf $1 Chip.psf raw
     187  load.cmf $2 Chip.psf out
     188
     189  # images generated with convolution will not have the right output positions
     190  set X_raw = int(X_PSF_raw) + 0.5
     191  set Y_raw = int(Y_PSF_raw) + 0.5
     192  set M_raw = PSF_INST_MAG_raw + $4
     193  set K_out = -2.5*log(KRON_FLUX_out)
     194  match2d X_PSF_out Y_PSF_out X_PSF_raw Y_PSF_raw 1.5 -index1 index1 -index2 index2
     195
     196  local i NX NY nx ny N
     197
     198  resize 1000 1000
     199
     200  clear
     201  section a0 0.0 0.0 1.0 0.5
     202  label -fn courier 14;
     203  style -pt 0 -sz 0.4
     204  show.pair 0
     205  delete -q imag_V dmag_V smag_V Dmag_V dmag_Vraw
     206  for imag -11 -6 0.2
     207    subset dmsub = delta if (rv > $imag) && (rv < $imag + 0.2)
     208    vstat -q dmsub
     209    concat $imag imag_V
     210    concat $MEAN dmag_V
     211    concat $MEDIAN Dmag_V
     212    concat $SIGMA smag_V
     213  end
     214
     215  vstat -q dmag_V
     216  $offset = $MEDIAN
     217  set dmag_V = dmag_V - $offset
     218  set Dmag_V = Dmag_V - $offset
     219
     220  section a1 0.0 0.50 1.0 0.25
     221  lim rv -0.025 0.075; label -fn courier 14; box;
     222  plot -c red   -pt 7 -sz 2.0 imag_V dmag_V
     223  plot -c darkgreen -pt 3 -sz 2.0 imag_V Dmag_V
     224  # plot -c gold -pt 4 -sz 2.0 imag_V dmag_Vraw
     225  plot -c blue  -pt 2 -sz 2.0 imag_V smag_V
     226  # label -y "mean (red), median (green), unclipped mean (gold), sigma (blue) of delta"
     227  label -y "mean (red), median (green), sigma (blue) of delta"
     228  sprintf line "OFFSET: %6.3f" $offset
     229  textline -frac 0.1 0.8 -fn courier 24 "$line"
     230
     231  set lChiNorm = log(PSF_CHISQ_out / PSF_NDOF_out)
     232  reindex chi = lChiNorm using index1
     233  section a2 0.0 0.75 1.0 0.25
     234  label -fn courier 14;
     235  lim rv chi; box;
     236  plot -c red -pt 0 -sz 0.5 rv chi
     237
     238  label -y "chisq"
     239  png -name $3.png
     240
     241  # section a1 0.0 0.5 1.0 0.5
     242  # style -pt 0 -sz 0.4
     243  # show.pair 1
     244end
     245
     246macro go.phot
     247 runphot test.00 test.00.varmode.C "-Db SAVE.RESID T -D LINEAR_FIT_VARIANCE_MODE CONSTANT"
     248 runphot test.00 test.00.varmode.I "-Db SAVE.RESID T -D LINEAR_FIT_VARIANCE_MODE IMAGE_VAR"
     249 runphot test.00 test.00.varmode.M "-Db SAVE.RESID T -D LINEAR_FIT_VARIANCE_MODE MODEL_VAR"
     250 runphot test.00 test.00.varmode.S "-Db SAVE.RESID T -D LINEAR_FIT_VARIANCE_MODE MODEL_SKY"
     251end
     252
     253macro go.vars
     254 dev -n varC; ckchip.mags test.00.in.cmf test.00.varmode.C.cmf test.00.varmode.C 0.0
     255 dev -n varI; ckchip.mags test.00.in.cmf test.00.varmode.I.cmf test.00.varmode.I 0.0
     256 dev -n varM; ckchip.mags test.00.in.cmf test.00.varmode.M.cmf test.00.varmode.M 0.0
     257 dev -n varS; ckchip.mags test.00.in.cmf test.00.varmode.S.cmf test.00.varmode.S 0.0
    135258end
    136259
    … …  
    347470    lim rv $word:4 $word:5; box; plot rv delta
    348471  end
    349   label -y '$word:0' -x '$word:2'
     472  $line = '$word:0' - '$word:1'
     473  label -y "$line" -x '$word:2'
    350474end
    351475
    … …  
    637761end
    638762
     763macro init
     764  # config for ppImage to generate chip, mask, weight
     765  $ppImageConfig = -recipe PPIMAGE PPIMAGE_N
     766  $ppImageConfig = $ppImageConfig -Db BACKGROUND T
     767  $ppImageConfig = $ppImageConfig -Db CHIP.FITS T
     768  $ppImageConfig = $ppImageConfig -Db CHIP.MASK.FITS T
     769  $ppImageConfig = $ppImageConfig -Db CHIP.VARIANCE.FITS T
     770  $ppImageConfig = $ppImageConfig -Db BASE.FITS F
     771  $ppImageConfig = $ppImageConfig -Db VARIANCE.BUILD T
     772  $ppImageConfig = $ppImageConfig -Db PHOTOM F
     773 
     774  # basic options for the these images (filter, location, obstype)
     775  $BaseOptions = -type OBJECT -filter r -ra 270.70 -dec -23.70 -pa 0.0
     776  $BaseOptions = $BaseOptions -Df PSASTRO:DVO.GETSTAR.MAX.RHO 50000.0
     777 
     778  # PSF.CONVOLVE : if true, we insert delta functions (and optionally
     779  #                galaxies) and smooth the image with the psf model
     780  #                (uses a GAUSS regardless of the model). Note that
     781  #                PSF.CONVOLVE = T is faster than F, but (a) only
     782  #                allows Gauss models and (b) only yields quantized
     783  #                locations
     784
     785  # create an image with fake sources and insert the resulting cmf file into a dvodb
     786  $RefConfig = -camera SIMTEST
     787  $RefConfig = $RefConfig -recipe PPSIM STACKTEST.MAKE
     788  $RefConfig = $RefConfig -D PSASTRO:PSASTRO.CATDIR catdir.ref
     789  $RefConfig = $RefConfig -Db PSF.CONVOLVE F
     790 
     791  # options for the reference image
     792  $RefOptions = $BaseOptions
     793  $RefOptions = $RefOptions -exptime 100.0
     794  $RefOptions = $RefOptions -seeing 1.0
     795  $RefOptions = $RefOptions -skymags 21.0 
     796  $RefOptions = $RefOptions -D PSF.MODEL PS_MODEL_GAUSS
     797  $RefOptions = $RefOptions -Df STARS.DENSITY 10.0
     798  $RefOptions = $RefOptions -Df STARS.SIGMA.LIM 0.5
     799
     800  # basic config for ppSim with randomly distributed stars and NO galaxies
     801  $RealConfig = -camera SIMTEST
     802  $RealConfig = $RealConfig -recipe PPSIM STACKTEST.RUN
     803  $RealConfig = $RealConfig -D PSASTRO:PSASTRO.CATDIR catdir.ref
     804  $RealConfig = $RealConfig -Db STARS.FAKE F
     805  $RealConfig = $RealConfig -Db STARS.REAL T
     806  $RealConfig = $RealConfig -Db MATCH.DENSITY F
     807  $RealConfig = $RealConfig -Db PSF.CONVOLVE F
     808  $RealConfig = $RealConfig -Df STARS.DENSITY 10.0
     809  $RealConfig = $RealConfig -Df STARS.SIGMA.LIM 2.5
     810  $RealConfig = $RealConfig -Db GALAXY.FAKE F
     811  $RealConfig = $RealConfig -Db GALAXY.GRID F
     812 
     813  # options for the repeated images
     814  $RealOptions = $BaseOptions -exptime 30.0
     815   
     816  # sample alternate options:
     817  # $ppSimOptions = $FakeOptions -D PSF.MODEL PS_MODEL_PS1_V1
     818  # $ppSimOptions = $FakeOptions -Df PSF.ARATIO 1.2
     819  # $ppSimOptions = $FakeOptions -Df PSF.THETA +30.0
     820  # $ppSimOptions = $FakeOptions -D PSF.MODEL PS_MODEL_GAUSS
     821 
     822  list fwhm
     823   1.0
     824   1.1
     825   1.2
     826   1.5
     827  end
     828end
     829
    639830if ($SCRIPT)
    640831  fulltest 4
Note: See TracChangeset for help on using the changeset viewer.