IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jan 28, 2009, 2:33:51 PM (17 years ago)
Author:
Paul Price
Message:

Changing pmReadout.weight to variance. Adding pmReadout.covariance which will carry around a covariance pseudo-matrix, allowing us to calculate the pixel-to-pixel variance. pmReadout.covariance now exists and is set, but no mechanism yet to read/write, or use it.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/pap_branch_20090128/psModules/src/objects/pmSourceFitSet.c

    r21183 r21211  
    66 *  @author GLG, MHPCC
    77 *
    8  *  @version $Revision: 1.14 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2009-01-27 06:39:38 $
     8 *  @version $Revision: 1.14.2.1 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2009-01-29 00:33:51 $
    1010 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    1111 *
    … …  
    5555
    5656    if (!fitSets) {
    57         fitSets = psArrayAlloc (PS_MAX (1, nThreads));
    58     }
    59 
    60     // the allocated elements should be NULL on psArrayAlloc, 
    61     // and a previously allocated array of fitSets should have been cleared 
     57        fitSets = psArrayAlloc (PS_MAX (1, nThreads));
     58    }
     59
     60    // the allocated elements should be NULL on psArrayAlloc,
     61    // and a previously allocated array of fitSets should have been cleared
    6262    // before pmSourceFitSetInit is called
    6363    for (int i = 0; i < fitSets->n; i++) {
    64         psAssert (fitSets->data[i] == NULL, "failure to init or clear fitSets?");
     64        psAssert (fitSets->data[i] == NULL, "failure to init or clear fitSets?");
    6565    }
    6666    return true;
    … …  
    125125    // we do not need to lock to do this....
    126126    for (int i = 0; i < fitSets->n; i++) {
    127         thisSet = fitSets->data[i];
    128         if (!thisSet) continue;
    129         if (thisSet->thread == id) {
    130             break;
    131         }
    132         thisSet = NULL;
     127        thisSet = fitSets->data[i];
     128        if (!thisSet) continue;
     129        if (thisSet->thread == id) {
     130            break;
     131        }
     132        thisSet = NULL;
    133133    }
    134134    psAssert (thisSet == NULL, "pmSourceFitSetDataSet() called but previous entry not cleared");
    … …  
    137137
    138138    pthread_mutex_lock (&fitSetInitMutex);
    139          
     139
    140140    // find an entry that is NULL and place it there
    141141    for (int i = 0; i < fitSets->n; i++) {
    142         if (fitSets->data[i]) continue;
    143         fitSets->data[i] = thisSet;
    144         pthread_mutex_unlock (&fitSetInitMutex);
    145         return thisSet;
     142        if (fitSets->data[i]) continue;
     143        fitSets->data[i] = thisSet;
     144        pthread_mutex_unlock (&fitSetInitMutex);
     145        return thisSet;
    146146    }
    147147    pthread_mutex_unlock (&fitSetInitMutex);
    … …  
    160160    pmSourceFitSetData *thisSet = NULL;
    161161    for (int i = 0; i < fitSets->n; i++) {
    162         thisSet = fitSets->data[i];
    163         if (!thisSet) continue;
    164         if (thisSet->thread == id) {
    165             break;
    166         }
    167         thisSet = NULL;
     162        thisSet = fitSets->data[i];
     163        if (!thisSet) continue;
     164        if (thisSet->thread == id) {
     165            break;
     166        }
     167        thisSet = NULL;
    168168    }
    169169    psAssert (thisSet != NULL, "pmSourceFitSetDataGet() called, but no entry found");
    … …  
    184184    pmSourceFitSetData *thisSet = NULL;
    185185    for (i = 0; i < fitSets->n; i++) {
    186         thisSet = fitSets->data[i];
    187         if (!thisSet) continue;
    188         if (thisSet->thread == id) {
    189             break;
    190         }
    191         thisSet = NULL;
     186        thisSet = fitSets->data[i];
     187        if (!thisSet) continue;
     188        if (thisSet->thread == id) {
     189            break;
     190        }
     191        thisSet = NULL;
    192192    }
    193193    psAssert (thisSet != NULL, "pmSourceFitSetDataClear() called, but no entry found");
    … …  
    256256        psVector *derivOne = set->derivSet->data[i];
    257257
    258         // one or the other (param or deriv) must be set
    259         assert ((deriv != NULL) || (param != NULL));
    260 
    261         // if we are setting derive, derivOne and paramOne must be same length
     258        // one or the other (param or deriv) must be set
     259        assert ((deriv != NULL) || (param != NULL));
     260
     261        // if we are setting derive, derivOne and paramOne must be same length
    262262        assert ((deriv == NULL) || (paramOne->n == derivOne->n));
    263        
     263
    264264        for (int j = 0; j < paramOne->n; j++, n++) {
    265             if (param) {
    266                 param->data.F32[n] = paramOne->data.F32[j];
    267             }
     265            if (param) {
     266                param->data.F32[n] = paramOne->data.F32[j];
     267            }
    268268            if (deriv) {
    269269                deriv->data.F32[n] = derivOne->data.F32[j];
    … …  
    275275
    276276// distribute parameters from single param and deriv vectors into FitSet models
    277 bool pmSourceFitSetSplit (pmSourceFitSetData *set, const psVector *deriv, const psVector *param) 
     277bool pmSourceFitSetSplit (pmSourceFitSetData *set, const psVector *deriv, const psVector *param)
    278278{
    279279    PS_ASSERT_VECTOR_NON_NULL(param, false);
    … …  
    301301
    302302// set the model parameters for this fit set
    303 bool pmSourceFitSetValues (pmSourceFitSetData *set, const psVector *dparam, 
     303bool pmSourceFitSetValues (pmSourceFitSetData *set, const psVector *dparam,
    304304                           const psVector *param, pmSource *source,
    305305                           psMinimization *myMin, int nPix, bool fitStatus)
    … …  
    420420
    421421    for (int i = 0; i < modelSet->n; i++) {
    422         pmModel *model = modelSet->data[i];
    423         int nParams = pmModelClassParameterCount (model->type);
    424         for (int j = 0; j < nParams; j++) {
    425             fprintf (file, "%d %d  : %f %f\n", i, j, model->params->data.F32[j], model->dparams->data.F32[j]);
    426         }
     422        pmModel *model = modelSet->data[i];
     423        int nParams = pmModelClassParameterCount (model->type);
     424        for (int j = 0; j < nParams; j++) {
     425            fprintf (file, "%d %d  : %f %f\n", i, j, model->params->data.F32[j], model->dparams->data.F32[j]);
     426        }
    427427    }
    428428    return true;
    … …  
    435435        psVector *derivOne = set->derivSet->data[i];
    436436        for (int j = 0; j < paramOne->n; j++) {
    437             fprintf (file, "%d %d  : %f %f\n", i, j, paramOne->data.F32[j], derivOne->data.F32[j]);
     437            fprintf (file, "%d %d  : %f %f\n", i, j, paramOne->data.F32[j], derivOne->data.F32[j]);
    438438        }
    439439    }
    … …  
    450450    PS_ASSERT_PTR_NON_NULL(source->pixels, false);
    451451    PS_ASSERT_PTR_NON_NULL(source->maskObj, false);
    452     PS_ASSERT_PTR_NON_NULL(source->weight, false);
     452    PS_ASSERT_PTR_NON_NULL(source->variance, false);
    453453
    454454    bool fitStatus = true;
    … …  
    471471                continue;
    472472            }
    473             // skip zero-weight points
    474             if (source->weight->data.F32[i][j] == 0) {
     473            // skip zero-variance points
     474            if (source->variance->data.F32[i][j] == 0) {
    475475                continue;
    476476            }
    … …  
    489489
    490490            // psMinimizeLMChi2 takes wt = 1/dY^2.  suggestion from RHL is to use the local sky
    491             // as weight to avoid the bias from systematic errors here we would just use the
     491            // as variance to avoid the bias from systematic errors here we would just use the
    492492            // source sky variance
    493493            if (PM_SOURCE_FIT_MODEL_PIX_WEIGHTS) {
    494                 yErr->data.F32[nPix] = 1.0 / source->weight->data.F32[i][j];
     494                yErr->data.F32[nPix] = 1.0 / source->variance->data.F32[i][j];
    495495            } else {
    496496                yErr->data.F32[nPix] = 1.0 / PM_SOURCE_FIT_MODEL_WEIGHT;
    … …  
    543543        psFree (params);
    544544        psFree(constraint);
    545         pmSourceFitSetDataClear(); // frees thisSet and removes if from the array of fitSets
     545        pmSourceFitSetDataClear(); // frees thisSet and removes if from the array of fitSets
    546546        return(false);
    547547    }
Note: See TracChangeset for help on using the changeset viewer.