IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Feb 5, 2009, 5:03:33 PM (17 years ago)
Author:
Paul Price
Message:

Merging pap_branch_20090128. Resolved a small number of conflicts. Compiles, but not tested in detail.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/psphotModelWithPSF.c

    r21183 r21366  
    2020        paramMask = constraint->paramMask;
    2121        if (paramMask != NULL) {
    22           PS_ASSERT_VECTOR_TYPE(paramMask, PS_TYPE_VECTOR_MASK, false);
     22          PS_ASSERT_VECTOR_TYPE(paramMask, PS_TYPE_VECTOR_MASK, false);
    2323            PS_ASSERT_VECTORS_SIZE_EQUAL(params, paramMask, false);
    2424        }
    … …  
    4242    // Alpha & Beta only contain elements to represent the unmasked parameters
    4343    if (!psMinLM_AllocAB (&Alpha, &Beta, params, paramMask)) {
    44         psAbort ("programming error: no unmasked parameters to be fit\n");
    45     }
    46    
     44        psAbort ("programming error: no unmasked parameters to be fit\n");
     45    }
     46
    4747    // allocate internal arrays (current vs Guess)
    4848    psImage *alpha   = psImageAlloc(Alpha->numCols, Alpha->numRows, PS_TYPE_F32);
    … …  
    127127            lambda *= 0.25;
    128128
    129             // save the new convolved model image
    130             psFree (source->modelFlux);
    131             source->modelFlux = pmPCMDataSaveImage(pcm);
     129            // save the new convolved model image
     130            psFree (source->modelFlux);
     131            source->modelFlux = pmPCMDataSaveImage(pcm);
    132132        } else {
    133133            lambda *= 10.0;
    … …  
    142142            psTrace ("psphot", 5, "failure to calculate covariance matrix\n");
    143143        }
    144         // set covar values which are not masked
    145         psImageInit (covar, 0.0);
    146         for (int j = 0, J = 0; j < params->n; j++) {
    147             if (paramMask && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[j])) {
    148                 covar->data.F32[j][j] = 1.0;
    149                 continue;
    150             }
    151             for (int k = 0, K = 0; k < params->n; k++) {
    152                 if (paramMask && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[k])) continue;
    153                 covar->data.F32[j][k] = Alpha->data.F32[J][K];
    154                 K++;
    155             }
    156             J++;
    157         }
     144        // set covar values which are not masked
     145        psImageInit (covar, 0.0);
     146        for (int j = 0, J = 0; j < params->n; j++) {
     147            if (paramMask && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[j])) {
     148                covar->data.F32[j][j] = 1.0;
     149                continue;
     150            }
     151            for (int k = 0, K = 0; k < params->n; k++) {
     152                if (paramMask && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[k])) continue;
     153                covar->data.F32[j][k] = Alpha->data.F32[J][K];
     154                K++;
     155            }
     156            J++;
     157        }
    158158    }
    159159
    … …  
    192192    PS_ASSERT_PTR_NON_NULL(source, NAN);
    193193    PS_ASSERT_IMAGE_NON_NULL(source->pixels, NAN);
    194     PS_ASSERT_IMAGE_NON_NULL(source->weight, NAN);
     194    PS_ASSERT_IMAGE_NON_NULL(source->variance, NAN);
    195195    PS_ASSERT_IMAGE_NON_NULL(source->maskObj, NAN);
    196196
    … …  
    210210    psImageInit (pcm->model, 0.0);
    211211    for (int n = 0; n < params->n; n++) {
    212         if (!pcm->dmodels->data[n]) continue;
    213         psImageInit (pcm->dmodels->data[n], 0.0);
     212        if (!pcm->dmodels->data[n]) continue;
     213        psImageInit (pcm->dmodels->data[n], 0.0);
    214214    }
    215215
    … …  
    218218        for (psS32 j = 0; j < source->pixels->numCols; j++) {
    219219
    220             // XXX can we skip some of the data points where the model
    221             // is not going to be fitted??
     220            // XXX can we skip some of the data points where the model
     221            // is not going to be fitted??
    222222
    223223            // skip masked points
    224             // XXX probably should not skipped masked points:
    225             // XXX skip if convolution of unmasked pixels will not see this pixel
     224            // XXX probably should not skipped masked points:
     225            // XXX skip if convolution of unmasked pixels will not see this pixel
    226226            // if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[i][j]) {
    227             // continue;
    228             // }
    229 
    230             // skip zero-weight points
    231             // XXX why is this not masked?
    232             // if (source->weight->data.F32[i][j] == 0) {
    233             // continue;
    234             // }
     227            // continue;
     228            // }
     229
     230            // skip zero-variance points
     231            // XXX why is this not masked?
     232            // if (source->variance->data.F32[i][j] == 0) {
     233            // continue;
     234            // }
    235235            // skip nan value points
    236             // XXX why is this not masked?
     236            // XXX why is this not masked?
    237237            // if (!isfinite(source->pixels->data.F32[i][j])) {
    238             // continue;
    239             // }
     238            // continue;
     239            // }
    240240
    241241            // Convert i/j to image space:
    … …  
    243243            coord->data.F32[1] = (psF32) (i + source->pixels->row0);
    244244
    245             pcm->model->data.F32[i][j] = func (deriv, params, coord);
    246 
    247             for (int n = 0; n < params->n; n++) {
    248                 if ((paramMask != NULL) && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n])) { continue; }
    249                 psImage *dmodel = pcm->dmodels->data[n];
    250                 dmodel->data.F32[i][j] = deriv->data.F32[n];
    251             }
     245            pcm->model->data.F32[i][j] = func (deriv, params, coord);
     246
     247            for (int n = 0; n < params->n; n++) {
     248                if ((paramMask != NULL) && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n])) { continue; }
     249                psImage *dmodel = pcm->dmodels->data[n];
     250                dmodel->data.F32[i][j] = deriv->data.F32[n];
     251            }
    252252        }
    253253    }
    … …  
    258258    psImageConvolveDirect (pcm->modelConv, pcm->model, psf);
    259259    for (int n = 0; n < pcm->dmodels->n; n++) {
    260         if (pcm->dmodels->data[n] == NULL) continue;
    261         psImage *dmodel = pcm->dmodels->data[n];
    262         psImage *dmodelConv = pcm->dmodelsConv->data[n];
    263         psImageConvolveDirect (dmodelConv, dmodel, psf);
     260        if (pcm->dmodels->data[n] == NULL) continue;
     261        psImage *dmodel = pcm->dmodels->data[n];
     262        psImage *dmodelConv = pcm->dmodelsConv->data[n];
     263        psImageConvolveDirect (dmodelConv, dmodel, psf);
    264264    }
    265265
    266266    // XXX TEST : SAVE IMAGES
    267 # if (SAVE_IMAGES) 
     267# if (SAVE_IMAGES)
    268268    psphotSaveImage (NULL, psf->image, "psf.fits");
    269269    psphotSaveImage (NULL, pcm->model, "model.fits");
    … …  
    271271    psphotSaveImage (NULL, source->pixels, "obj.fits");
    272272    psphotSaveImage (NULL, source->maskObj, "mask.fits");
    273     psphotSaveImage (NULL, source->weight, "weight.fits");
     273    psphotSaveImage (NULL, source->variance, "variance.fits");
    274274# endif
    275275
    276     // 2 *** accumulate alpha & beta 
     276    // 2 *** accumulate alpha & beta
    277277
    278278    // zero alpha and beta for summing below
    … …  
    283283    for (psS32 i = 0; i < source->pixels->numRows; i++) {
    284284        for (psS32 j = 0; j < source->pixels->numCols; j++) {
    285             // XXX are we doing the right thing with the mask?
     285            // XXX are we doing the right thing with the mask?
    286286            // skip masked points
    287287            if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[i][j]) {
    288288                continue;
    289289            }
    290             // skip zero-weight points
    291             if (source->weight->data.F32[i][j] == 0) {
     290            // skip zero-variance points
     291            if (source->variance->data.F32[i][j] == 0) {
    292292                continue;
    293293            }
    … …  
    297297            }
    298298
    299             float ymodel  = pcm->modelConv->data.F32[i][j];
    300             float yweight = 1.0 / source->weight->data.F32[i][j];
    301             float delta = ymodel - source->pixels->data.F32[i][j];
    302 
    303             chisq += PS_SQR(delta) * yweight;
    304 
    305             if (isnan(delta)) psAbort("nan in delta");
    306             if (isnan(chisq)) psAbort("nan in chisq");
    307 
    308             // alpha & beta only contain unmasked elements
    309             for (int n1 = 0, N1 = 0; n1 < params->n; n1++) {
    310                 if ((paramMask != NULL) && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n1])) continue;
    311                 psImage *dmodel = pcm->dmodelsConv->data[n1];
    312                 float weight = dmodel->data.F32[i][j] * yweight;
    313                 for (int n2 = 0, N2 = 0; n2 <= n1; n2++) {
    314                     if ((paramMask != NULL) && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n2])) continue;
    315                     dmodel = pcm->dmodelsConv->data[n2];
    316                     alpha->data.F32[N1][N2] += weight * dmodel->data.F32[i][j];
    317                     N2++;
    318                 }
    319                 beta->data.F32[N1] += weight * delta;
    320                 N1++;
    321             }
    322         }
     299            float ymodel  = pcm->modelConv->data.F32[i][j];
     300            float yweight = 1.0 / source->variance->data.F32[i][j];
     301            float delta = ymodel - source->pixels->data.F32[i][j];
     302
     303            chisq += PS_SQR(delta) * yweight;
     304
     305            if (isnan(delta)) psAbort("nan in delta");
     306            if (isnan(chisq)) psAbort("nan in chisq");
     307
     308            // alpha & beta only contain unmasked elements
     309            for (int n1 = 0, N1 = 0; n1 < params->n; n1++) {
     310                if ((paramMask != NULL) && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n1])) continue;
     311                psImage *dmodel = pcm->dmodelsConv->data[n1];
     312                float weight = dmodel->data.F32[i][j] * yweight;
     313                for (int n2 = 0, N2 = 0; n2 <= n1; n2++) {
     314                    if ((paramMask != NULL) && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n2])) continue;
     315                    dmodel = pcm->dmodelsConv->data[n2];
     316                    alpha->data.F32[N1][N2] += weight * dmodel->data.F32[i][j];
     317                    N2++;
     318                }
     319                beta->data.F32[N1] += weight * delta;
     320                N1++;
     321            }
     322        }
    323323    }
    324324
    … …  
    356356    pcm->dmodels = psArrayAlloc (params->n);
    357357    for (psS32 n = 0; n < params->n; n++) {
    358         pcm->dmodels->data[n] = NULL;
    359         if ((paramMask != NULL) && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n])) { continue; }
    360         pcm->dmodels->data[n] = psImageCopy (NULL, source->pixels, PS_TYPE_F32);
     358        pcm->dmodels->data[n] = NULL;
     359        if ((paramMask != NULL) && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n])) { continue; }
     360        pcm->dmodels->data[n] = psImageCopy (NULL, source->pixels, PS_TYPE_F32);
    361361    }
    362362
    … …  
    365365    pcm->dmodelsConv = psArrayAlloc (params->n);
    366366    for (psS32 n = 0; n < params->n; n++) {
    367         pcm->dmodelsConv->data[n] = NULL;
    368         if ((paramMask != NULL) && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n])) { continue; }
    369         pcm->dmodelsConv->data[n] = psImageCopy (NULL, source->pixels, PS_TYPE_F32);
     367        pcm->dmodelsConv->data[n] = NULL;
     368        if ((paramMask != NULL) && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n])) { continue; }
     369        pcm->dmodelsConv->data[n] = psImageCopy (NULL, source->pixels, PS_TYPE_F32);
    370370    }
    371371
    … …  
    376376
    377377    psImage *model = psImageCopy (NULL, pcm->modelConv, PS_TYPE_F32);
    378    
     378
    379379    return model;
    380380}
    … …  
    383383 *
    384384 * we have a function func(param; value)
    385  
     385
    386386 * basic LMM:
    387  
     387
    388388 - fill in the data (x, y)
    389  
     389
    390390 chisq = SetABX (alpha, beta, params, paramMask, x, y, dy, func)
    391  
     391
    392392 while () {
    393393 GuessABP (Alpha, Beta, Params, alpha, beta, params, paramMask, checkLimits, lambda)
    … …  
    396396 convergence tests...
    397397 }
    398  
    399  
     398
     399
    400400
    401401 ** GuessABP:
Note: See TracChangeset for help on using the changeset viewer.