IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 3788


Ignore:
Timestamp:
Apr 28, 2005, 11:17:02 PM (21 years ago)
Author:
magnier
Message:

checking in the versions of the files below, incorporating the changes
I made and have already submitted to bugzilla. these files should be
replaced with the versions contributed on the main trunk by MHPCC.
I will tag this state with eam-psphot-1

src/astronomy/psMetadata.c src/astronomy/psMetadata.h
src/collections/psVector.c src/collections/psVector.h
src/dataManip/psMinimize.c src/dataManip/psMinimize.h
src/image/psImage.c src/image/psImage.h src/sysUtils/psTrace.c
src/sysUtils/psTrace.h

Location:
branches/eam-psphot-branch/psLib/src
Files:
18 edited

Legend:

Unmodified
Added
Removed

  • branches/eam-psphot-branch/psLib/src/astronomy/psMetadata.c

    r3541 r3788  
    1212*  @author Ross Harman, MHPCC
    1313*
    14 *  @version $Revision: 1.57 $ $Name: not supported by cvs2svn $
    15 *  @date $Date: 2005-03-29 21:31:54 $
     14*  @version $Revision: 1.57.4.1 $ $Name: not supported by cvs2svn $
     15*  @date $Date: 2005-04-29 09:17:02 $
    1616*
    1717*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    6868/*****************************************************************************/
    6969
    70 static psMetadataItem* makeMetaMulti(psHash* table, const char* key, psMetadataItem* existing)
     70# if (0)
     71    static psMetadataItem* makeMetaMulti(psHash* table, const char* key, psMetadataItem* existing)
    7172{
    7273
    … …  
    9596    return item;
    9697}
     98# endif
    9799
    98100static void metadataItemFree(psMetadataItem* metadataItem)
    … …  
    174176{
    175177    psMetadataItem* metadataItem = NULL;
    176 
     178    char tmp;
     179    int Nbyte;
    177180
    178181    PS_PTR_CHECK_NULL(name,NULL);
    … …  
    185188    p_psMemSetDeallocator(metadataItem, (psFreeFcn) metadataItemFree);
    186189
    187     // Allocate and set metadata item comment
    188     metadataItem->comment = (char *)psAlloc(sizeof(char) * MAX_STRING_LENGTH);
     190    // set metadata item comment
    189191    if (comment == NULL) {
    190192        // Per SDRS, null isn't allowed, must use "" instead
    191         strncpy(metadataItem->comment, "", MAX_STRING_LENGTH);
     193        metadataItem->comment = psStringCopy ("");
    192194    } else {
    193         strncpy(metadataItem->comment, comment, MAX_STRING_LENGTH);
     195        metadataItem->comment = psStringCopy (comment);
    194196    }
    195197
    … …  
    201203
    202204    // Allocate and set metadata item name
    203     metadataItem->name = (char *)psAlloc(sizeof(char) * MAX_STRING_LENGTH);
     205    Nbyte = vsnprintf(&tmp, 0, name, argPtr) + 1;
     206    metadataItem->name = (char *)psAlloc(sizeof(char) * Nbyte);
    204207    vsprintf(metadataItem->name, name, argPtr);
    205208
    … …  
    220223    case PS_META_STR:
    221224        // Perform copy of input strings
    222         metadataItem->data.V = psStringNCopy(va_arg(argPtr, char *), MAX_STRING_LENGTH);
     225        metadataItem->data.V = psStringCopy(va_arg(argPtr, char *));
    223226        break;
    224227    case PS_META_LIST:
    … …  
    230233    case PS_META_ASTROM:
    231234    case PS_META_UNKNOWN:
    232     case PS_META_MULTI:
    233235        // Copy of input data not performed due to variability of data types
    234236        metadataItem->data.V = psMemIncrRefCounter(va_arg(argPtr, psPtr));
     237        break;
     238    case PS_META_MULTI:
     239        // MULTI needs to create a psList entry, value must be NULL
     240        metadataItem->data.list = psListAlloc(NULL);
    235241        break;
    236242    default:
    … …  
    239245        metadataItem = NULL;
    240246    }
    241 
    242247    return metadataItem;
    243248}
    … …  
    262267
    263268    return metadata;
     269}
     270
     271psBool MetadataAddNewItem (psMetadata *md, psMetadataItem *metadataItem, psS32 location)
     272{
     273
     274    char * key = NULL;
     275    psHash *mdTable = NULL;
     276    psList *mdList = NULL;
     277
     278    mdTable = md->table;
     279    mdList = md->list;
     280    key = metadataItem->name;
     281
     282    // OK, this is a new item.
     283    // Add new metadata item to metadata collection's hash
     284    if(!psHashAdd(mdTable, key, metadataItem)) {
     285        psError(PS_ERR_UNKNOWN,false,PS_ERRORTEXT_psMetadata_ADD_TABLE_FAILED,key);
     286        return false;
     287    }
     288
     289    // Add new metadata item to metadata collection's list
     290    if(!psListAdd(mdList, location, metadataItem)) {
     291        psError(PS_ERR_UNKNOWN,false,PS_ERRORTEXT_psMetadata_ADD_COLLECTION_FAILED,key);
     292        return false;
     293    }
     294    // Success
     295    return true;
    264296}
    265297
    … …  
    282314
    283315    // See if key is already in table
    284     existingEntry = (psMetadataItem*)psHashLookup(mdTable, key);
    285 
     316    existingEntry = psMetadataLookup(md, key);
     317
     318    if (existingEntry == NULL) {
     319        MetadataAddNewItem (md, metadataItem, location);
     320        return true;
     321    }
     322
     323    // if existing entry is MULTI, add the new entry to that list
     324    // XXX does this behave correctly for adding a MULTI node to an exiting MULTI node?
     325    if (existingEntry->type == PS_META_MULTI) {
     326        PS_PTR_CHECK_NULL(existingEntry->data.list,NULL);
     327
     328        // Add new metadata item to multi entry's list (XXX is TAIL correct, or location?)
     329        if(!psListAdd(existingEntry->data.list, PS_LIST_TAIL, metadataItem)) {
     330            psError(PS_ERR_UNKNOWN,false,PS_ERRORTEXT_psMetadata_ADD_COLLECTION_FAILED,key);
     331            return false;
     332        }
     333
     334        // Add new metadata item to metadata collection's list
     335        if(!psListAdd(mdList, location, metadataItem)) {
     336            psError(PS_ERR_UNKNOWN,false,PS_ERRORTEXT_psMetadata_ADD_COLLECTION_FAILED,key);
     337            return false;
     338        }
     339        return true;
     340    }
     341
     342    // if existing entry is NOT MULTI, we have several choices:
     343    //   - new entry is MULTI: move existing entry to new MULTI list
     344    //   - DUPLICATE_OK: create a new MULTI and add both old and new to it
     345    //   - REPLACE : replace the existing entry with the new entry
     346    //   - raise an error
     347
     348    // incoming entry is PS_META_MULTI  (DUPLICATE_OK is redundant, REPLACE is non-sequitor)
    286349    if (metadataItem->type == PS_META_MULTI) {
    287         // the incoming entry is PS_META_MULTI
    288 
    289         // force the hash entry to be PS_META_MULTI
    290         existingEntry = makeMetaMulti(mdTable,key,existingEntry);
    291 
    292         // add all the items in the incoming entry to metadata
    293         psList* list = metadataItem->data.list;
    294         if (list != NULL) {
    295             psListIterator* iter = psListIteratorAlloc(list,PS_LIST_HEAD,true);
    296             psMetadataItem* listItem;
    297             while ((listItem=(psMetadataItem*)psListGetAndIncrement(iter)) != NULL) {
    298                 psMetadataAddItem(md,listItem,location,flags);
    299             }
    300             psFree(iter);
    301         }
     350
     351        // replace hash entry with points to existingEntry with new entry
     352        psHashAdd (mdTable, key, metadataItem);
     353
     354        // move the existing hash entry to a node of the new entry
     355        psListAdd (metadataItem->data.list, PS_LIST_TAIL, existingEntry);
    302356
    303357        return true; // all done.
    304358    }
    305359
    306     // how the item is added to the hash depends on prior existence, flags, etc.
    307     if(existingEntry != NULL) { // prior existence
    308         if (existingEntry->type == PS_META_MULTI || (flags & PS_META_DUPLICATE_OK) != 0) {
    309             // duplicate entries allowed - add another entry.
    310 
    311             // make sure the existing entry is PS_META_MULTI
    312             existingEntry = makeMetaMulti(mdTable,key,existingEntry);
    313 
    314             // add to the hash's list of duplicate entries
    315             if (! psListAdd(existingEntry->data.list, PS_LIST_TAIL, metadataItem) ) {
    316                 psError(PS_ERR_UNKNOWN,false,PS_ERRORTEXT_psMetadata_ADD_COLLECTION_FAILED,key);
    317                 return false;
    318             }
    319         } else if ((flags & PS_META_REPLACE) != 0) {
    320             // replace entry instead of creating a duplicate entry.
    321 
    322             // remove the existing entry from metadata
    323             psMetadataRemove(md,0,key);
    324 
    325             // treat as if new (added to list below)
    326             if(!psHashAdd(mdTable, key, metadataItem)) {
    327                 psError(PS_ERR_UNKNOWN,false,PS_ERRORTEXT_psMetadata_ADD_TABLE_FAILED,key);
    328                 return false;
    329             }
    330         } else {
    331             // default is to error on duplicate entry.
    332             psError(PS_ERR_BAD_PARAMETER_VALUE, true,
    333                     PS_ERRORTEXT_psMetadata_DUPLICATE_NOT_ALLOWED);
    334             return false;
    335         }
    336     } else {
    337         // OK, this is a new item.
    338 
    339         // Node doesn't exist - Add new metadata item to metadata collection's hash
    340         if(!psHashAdd(mdTable, key, metadataItem)) {
    341             psError(PS_ERR_UNKNOWN,false,PS_ERRORTEXT_psMetadata_ADD_TABLE_FAILED,key);
    342             return false;
    343         }
    344     }
    345 
    346     if(!psListAdd(mdList, location, metadataItem)) {
    347         psError(PS_ERR_UNKNOWN,false,PS_ERRORTEXT_psMetadata_ADD_COLLECTION_FAILED,key);
    348         return false;
    349     }
    350 
    351     return true;
     360    // duplicate entries allowed - create a MULTI with name and add both to it
     361    if (flags & PS_META_DUPLICATE_OK) {
     362
     363        psMetadataItem *multi;
     364
     365        multi = psMetadataItemAlloc (key, PS_META_MULTI, "multi container", NULL);
     366
     367        // replace hash entry with points to existingEntry with new multi
     368        psHashAdd (mdTable, key, multi);
     369
     370        // move the existing hash entry to a node of the new entry
     371        psListAdd (multi->data.list, PS_LIST_TAIL, existingEntry);
     372
     373        // move the existing hash entry to a node of the new entry
     374        psListAdd (multi->data.list, PS_LIST_TAIL, metadataItem);
     375
     376        // move the existing hash entry to a node of the new entry
     377        psListAdd (mdList, PS_LIST_TAIL, metadataItem);
     378
     379        return true;
     380    }
     381
     382    if (flags & PS_META_REPLACE) {
     383        // replace entry instead of creating a duplicate entry.
     384
     385        // remove the existing entry from metadata
     386        psMetadataRemove(md,0,key);
     387
     388        // treat as if new (added to list below)
     389        MetadataAddNewItem (md, metadataItem, location);
     390        return true;
     391    }
     392
     393    // default is to error on duplicate entry.
     394    psError(PS_ERR_BAD_PARAMETER_VALUE, true,
     395            PS_ERRORTEXT_psMetadata_DUPLICATE_NOT_ALLOWED);
     396    return false;
    352397}
    353398

  • branches/eam-psphot-branch/psLib/src/astronomy/psMetadata.h

    r3544 r3788  
    1111*  @author Ross Harman, MHPCC
    1212*
    13 *  @version $Revision: 1.42 $ $Name: not supported by cvs2svn $
    14 *  @date $Date: 2005-03-29 22:13:53 $
     13*  @version $Revision: 1.42.4.1 $ $Name: not supported by cvs2svn $
     14*  @date $Date: 2005-04-29 09:17:02 $
    1515*
    1616*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii

  • branches/eam-psphot-branch/psLib/src/collections/psVector.c

    r3476 r3788  
    99*  @author Robert DeSonia, MHPCC
    1010*
    11 *  @version $Revision: 1.36 $ $Name: not supported by cvs2svn $
    12 *  @date $Date: 2005-03-22 21:52:49 $
     11*  @version $Revision: 1.36.4.1 $ $Name: not supported by cvs2svn $
     12*  @date $Date: 2005-04-29 09:17:02 $
    1313*
    1414*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    540540    }
    541541}
     542
     543// XXX EAM a utility function
     544bool p_psVectorPrint (FILE *f, psVector *a, char *name)
     545{
     546
     547    fprintf (f, "vector: %s\n", name);
     548
     549    for (int i = 0; i < a[0].n; i++) {
     550        fprintf (f, "%f\n", p_psVectorGetElementF64(a, i));
     551    }
     552    fprintf (f, "\n");
     553    return (true);
     554}
     555

  • branches/eam-psphot-branch/psLib/src/collections/psVector.h

    r3476 r3788  
    1111 *  @author Ross Harman, MHPCC
    1212 *
    13  *  @version $Revision: 1.30 $ $Name: not supported by cvs2svn $
    14  *  @date $Date: 2005-03-22 21:52:49 $
     13 *  @version $Revision: 1.30.4.1 $ $Name: not supported by cvs2svn $
     14 *  @date $Date: 2005-04-29 09:17:02 $
    1515 *
    1616 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    160160);
    161161
     162// XXX EAM a utility function
     163bool p_psVectorPrint (FILE *f, psVector *a, char *name);
    162164
    163165/// @}

  • branches/eam-psphot-branch/psLib/src/dataManip/psMinimize.c

    r3578 r3788  
    99 *  @author GLG, MHPCC
    1010 *
    11  *  @version $Revision: 1.110 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2005-03-31 01:02:15 $
     11 *  @version $Revision: 1.110.6.1 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2005-04-29 09:17:02 $
    1313 *
    1414 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    557557}
    558558
    559 /******************************************************************************
    560 psMinimizeLMChi2():  This routine will take an procedure which calculates
    561 an arbitrary function and it's derivative and minimize the chi-squared match
    562 between that function at the specified coords and the specified value at
    563 those coords.
    564  
    565 XXX: Do this:
    566  After checking that all entries in the paramMask are 1 or 0, when
    567  forming the A matrix from alpha, try this:
    568  
    569      A[i][i] = (1 + lambda*paramask[i]) * alpha[i][i];
    570  
    571 XXX: This is very different from what is specified in the SDR.  Must
    572 coordinate with IfA on new SDR.
    573  
    574 XXX: Do vector/image recycles.
    575  
    576 XXX: probably yErr will be part of the SDR.
    577  
    578 XXX: This must work for both F32 and F64.  F32 is currently implemented.
    579      Note: since the LUD routines are only implemented in F64, then we
    580      will have to convert all F32 input vectors to F64 regardless.  So,
    581      the F64 port might be.
    582  
    583 XXX: Must update the covar matrix.
    584  *****************************************************************************/
     559// XXX EAM this is my re-implementation of MinLM
    585560psBool psMinimizeLMChi2(psMinimization *min,
    586561                        psImage *covar,
    … …  
    601576    PS_PTR_CHECK_NULL(func, NULL);
    602577
     578    // this function has test and current values for several things
     579    // the current best value is in lower case
     580    // the next guess value is in upper case
     581
     582    // allocate internal arrays (current vs Guess)
     583    psImage *alpha   = psImageAlloc (params->n, params->n, PS_TYPE_F64);
     584    psImage *Alpha   = psImageAlloc (params->n, params->n, PS_TYPE_F64);
     585    psVector *beta   = psVectorAlloc (params->n, PS_TYPE_F64);
     586    psVector *Beta   = psVectorAlloc (params->n, PS_TYPE_F64);
     587    psVector *Params = psVectorAlloc (params->n, PS_TYPE_F64);
     588    psVector *dy     = NULL;
     589    psF64 chisq = 0.0;
     590    psF64 Chisq = 0.0;
     591    psF64 lambda = 0.001;
     592
     593    // the initial guess on params is provided by the user
     594    Params = psVectorCopy (Params, params, PS_TYPE_F32);
     595
     596    // the user provides the error or NULL.  we need to convert
     597    // to appropriate weights
     598    dy = psVectorAlloc (y->n, PS_TYPE_F32);
     599    if (yErr != NULL) {
     600        for (int i = 0; i < dy->n; i++) {
     601            dy->data.F32[i] = 1.0 / PS_SQR (yErr->data.F32[i]);
     602        }
     603    } else {
     604        for (int i = 0; i < dy->n; i++) {
     605            dy->data.F32[i] = 1.0;
     606        }
     607    }
     608
     609    // calculate initial alpha and beta, set chisq (min->value)
     610    min->value = p_psMinLM_SetABX (alpha, beta, params, x, y, dy, func);
     611    # ifndef PS_NO_TRACE
     612    // dump some useful info if trace is defined
     613    if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") > 4) {
     614        p_psImagePrint  (psTraceGetDestination(), alpha, "alpha guess");
     615        p_psVectorPrint (psTraceGetDestination(), beta, "beta guess");
     616        p_psVectorPrint (psTraceGetDestination(), params, "params guess");
     617    }
     618    # endif /* PS_NO_TRACE */
     619
     620
     621    // iterate until the tolerance is reached, or give up
     622    while ((min->lastDelta > min->tol) && (min->iter < min->maxIter)) {
     623
     624        // set a new guess for Alpha, Beta, Params
     625        p_psMinLM_GuessABP (Alpha, Beta, Params, alpha, beta, params, lambda);
     626
     627        # ifndef PS_NO_TRACE
     628        // dump some useful info if trace is defined
     629        if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") > 4) {
     630            p_psImagePrint  (psTraceGetDestination(), Alpha, "alpha guess");
     631            p_psVectorPrint (psTraceGetDestination(), Beta, "beta guess");
     632            p_psVectorPrint (psTraceGetDestination(), Params, "params guess");
     633        }
     634        # endif /* PS_NO_TRACE */
     635
     636        // calculate Chisq for new guess, update Alpha & Beta
     637        Chisq = p_psMinLM_SetABX (Alpha, Beta, Params, x, y, dy, func);
     638        psTrace (".psLib.dataManip.psMinimizeLMChi2", 3, "chisq: %f, Chisq %f, delta: %f\n", chisq, Chisq, min->lastDelta);
     639
     640        // accept new guess (if improvement), or increase lambda
     641        if (Chisq < min->value) {
     642            min->lastDelta = (min->value - Chisq) / (dy->n - params->n);
     643            min->value = Chisq;
     644            alpha  = psImageCopy (alpha, Alpha, PS_TYPE_F64);
     645            beta   = psVectorCopy (beta, Beta, PS_TYPE_F64);
     646            params = psVectorCopy (params, Params, PS_TYPE_F32);
     647            lambda *= 0.1;
     648        } else {
     649            lambda *= 10.0;
     650        }
     651        min->iter ++;
     652    }
     653    psTrace (".psLib.dataManip.psMinimizeLMChi2", 3, "chisq: %f, Chisq %f, delta: %f\n", chisq, Chisq, min->lastDelta);
     654
     655    // free the internal temporary data
     656    psFree (alpha);
     657    psFree (Alpha);
     658    psFree (beta);
     659    psFree (Beta);
     660    psFree (Params);
     661    psFree (dy);
     662    return (true);
     663}
     664
     665// XXX EAM: this needs to respect the mask on params
     666// XXX EAM: check not NULL on alpha, beta, params
     667// alpha, beta, params are already allocated
     668psF64 p_psMinLM_SetABX (psImage  *alpha,
     669                        psVector *beta,
     670                        psVector *params,
     671                        const psArray  *x,
     672                        const psVector *y,
     673                        const psVector *dy,
     674                        psMinimizeLMChi2Func func)
     675{
     676
     677    psF64 chisq;
     678    psF64 delta;
     679    psF64 weight;
     680    psF64 ymodel;
     681    psVector *deriv = psVectorAlloc (params->n, PS_TYPE_F32);
     682
     683    // zero alpha and beta for summing below
     684    for (int j = 0; j < params->n; j++) {
     685        for (int k = 0; k < params->n; k++) {
     686            alpha->data.F64[j][k] = 0;
     687        }
     688        beta->data.F64[j] = 0;
     689    }
     690    chisq = 0.0;
     691
     692    // calculate chisq, alpha, beta
     693    for (int i = 0; i < y->n; i++) {
     694        ymodel = func (deriv, params, (psVector *) x->data[i]);
     695
     696        delta = ymodel - y->data.F32[i];
     697        chisq += PS_SQR (delta) * dy->data.F32[i];
     698
     699        for (int j = 0; j < params->n; j++) {
     700            weight = deriv->data.F32[j] * dy->data.F32[i];
     701            for (int k = 0; k <= j; k++) {
     702                alpha->data.F64[j][k] += weight * deriv->data.F32[k];
     703            }
     704            beta->data.F64[j] += weight * delta;
     705        }
     706    }
     707
     708    // calculate lower-left half of alpha
     709    for (int j = 1; j < params->n; j++) {
     710        for (int k = 0; k < j; k++) {
     711            alpha->data.F64[k][j] = alpha->data.F64[j][k];
     712        }
     713    }
     714    psFree (deriv);
     715    return (chisq);
     716}
     717
     718// XXX EAM : can we use static copies of LUv, LUm, A?
     719psBool p_psMinLM_GuessABP (psImage  *Alpha,
     720                           psVector *Beta,
     721                           psVector *Params,
     722                           psImage  *alpha,
     723                           psVector *beta,
     724                           psVector *params,
     725                           psF64 lambda)
     726{
     727
     728    # define USE_LU_DECOMP 1
     729    # if (USE_LU_DECOMP)
     730        psVector *LUv = NULL;
     731    psImage  *LUm = NULL;
     732    psImage  *A   = NULL;
     733    psF32    det;
     734
     735    // LU decomposition version
     736    psTrace (".pslib.dataManip.psMinLM_GuessABP", 3, "using LUD version");
     737
     738    // set new guess values (creates matrix A)
     739    A = psImageCopy (NULL, alpha, PS_TYPE_F64);
     740    for (int j = 0; j < params->n; j++) {
     741        A->data.F64[j][j] = alpha->data.F64[j][j] * (1.0 + lambda);
     742    }
     743
     744    // solve A*beta = Beta (Alpha = 1/A)
     745    // these operations do not modify the input values (creates LUm, LUv)
     746    LUm   = psMatrixLUD (NULL, &LUv, A);
     747    Beta  = psMatrixLUSolve (Beta, LUm, beta, LUv);
     748    Alpha = psMatrixInvert (Alpha, A, &det);
     749
     750    # else
     751        // gauss-jordan version
     752        psTrace (".pslib.dataManip.psMinLM_GuessABP", 3, "using Gauss-J version");
     753
     754    // set new guess values (creates matrix A)
     755    Beta = psVectorCopy (Beta, beta, PS_TYPE_F64);
     756    Alpha = psImageCopy (Alpha, alpha, PS_TYPE_F64);
     757    for (int j = 0; j < params->n; j++) {
     758        Alpha->data.F64[j][j] = alpha->data.F64[j][j] * (1.0 + lambda);
     759    }
     760
     761    psGaussJordan (Alpha, Beta);
     762    # endif
     763
     764    // apply beta to get new params values
     765    for (int j = 0; j < params->n; j++) {
     766        Params->data.F32[j] = params->data.F32[j] - Beta->data.F64[j];
     767    }
     768
     769    # if (USE_LU_DECOMP)
     770        psFree (A);
     771    psFree (LUm);
     772    psFree (LUv);
     773    # endif
     774
     775    return true;
     776}
     777
     778# define SWAP(X,Y) {double tmp=(X); (X) = (Y); (Y) = tmp;}
     779
     780// XXX EAM : temporary gauss-jordan solver based on gene's
     781// version based on the Numerical Recipes version
     782bool psGaussJordan (psImage *a, psVector *b)
     783{
     784
     785    int *indxc,*indxr,*ipiv;
     786    int Nx, icol, irow;
     787    int i, j, k, l, ll;
     788    float big, dum, pivinv;
     789    psF64 *vector;
     790    psF64 **matrix;
     791
     792    Nx = a->numCols;
     793    matrix = a->data.F64;
     794    vector = b->data.F64;
     795
     796    indxc = psAlloc (Nx*sizeof(int));
     797    indxr = psAlloc (Nx*sizeof(int));
     798    ipiv  = psAlloc (Nx*sizeof(int));
     799    for (j = 0; j < Nx; j++)
     800        ipiv[j] = 0;
     801
     802    irow = icol = 0;
     803    big = fabs(matrix[0][0]);
     804
     805    for (i = 0; i < Nx; i++) {
     806        big = 0.0;
     807        for (j = 0; j < Nx; j++) {
     808            if (!finite(matrix[i][j])) {
     809                // XXX EAM: this should use the psError stack
     810                fprintf (stderr, "GAUSSJ: NaN\n");
     811                goto fescape;
     812            }
     813            if (ipiv[j] != 1) {
     814                for (k = 0; k < Nx; k++) {
     815                    if (ipiv[k] == 0) {
     816                        if (fabs (matrix[j][k]) >= big) {
     817                            big  = fabs (matrix[j][k]);
     818                            irow = j;
     819                            icol = k;
     820                        }
     821                    } else {
     822                        if (ipiv[k] > 1) {
     823                            // XXX EAM: this should use the psError stack
     824                            fprintf (stderr, "GAUSSJ: Singular Matrix! (1)\n");
     825                            goto fescape;
     826                        }
     827                    }
     828                }
     829            }
     830        }
     831        ipiv[icol]++;
     832        if (irow != icol) {
     833            for (l = 0; l < Nx; l++) {
     834                SWAP (matrix[irow][l], matrix[icol][l]);
     835            }
     836            SWAP (vector[irow], vector[icol]);
     837        }
     838        indxr[i] = irow;
     839        indxc[i] = icol;
     840        if (matrix[icol][icol] == 0.0) {
     841            // XXX EAM: this should use the psError stack
     842            fprintf (stderr, "GAUSSJ: Singular Matrix! (2)\n");
     843            goto fescape;
     844        }
     845        pivinv = 1.0 / matrix[icol][icol];
     846        matrix[icol][icol] = 1.0;
     847        for (l = 0; l < Nx; l++) {
     848            matrix[icol][l] *= pivinv;
     849        }
     850        vector[icol] *= pivinv;
     851
     852        for (ll = 0; ll < Nx; ll++) {
     853            if (ll != icol) {
     854                dum = matrix[ll][icol];
     855                matrix[ll][icol] = 0.0;
     856                for (l = 0; l < Nx; l++)
     857                    matrix[ll][l] -= matrix[icol][l]*dum;
     858                vector[ll] -= vector[icol]*dum;
     859            }
     860        }
     861    }
     862
     863    for (l = Nx - 1; l >= 0; l--) {
     864        if (indxr[l] != indxc[l])
     865            for (k = 0; k < Nx; k++)
     866                SWAP (matrix[k][indxr[l]], matrix[k][indxc[l]]);
     867    }
     868    psFree (ipiv);
     869    psFree (indxr);
     870    psFree (indxc);
     871    return (true);
     872
     873fescape:
     874    psFree (ipiv);
     875    psFree (indxr);
     876    psFree (indxc);
     877    return (false);
     878}
     879
     880/******************************************************************************
     881psMinimizeLMChi2():  This routine will take an procedure which calculates
     882an arbitrary function and it's derivative and minimize the chi-squared match
     883between that function at the specified coords and the specified value at
     884those coords.
     885 
     886XXX: Do this:
     887 After checking that all entries in the paramMask are 1 or 0, when
     888 forming the A matrix from alpha, try this:
     889 
     890     A[i][i] = (1 + lambda*paramask[i]) * alpha[i][i];
     891 
     892XXX: This is very different from what is specified in the SDR.  Must
     893coordinate with IfA on new SDR.
     894 
     895XXX: Do vector/image recycles.
     896 
     897XXX: probably yErr will be part of the SDR.
     898 
     899XXX: This must work for both F32 and F64.  F32 is currently implemented.
     900     Note: since the LUD routines are only implemented in F64, then we
     901     will have to convert all F32 input vectors to F64 regardless.  So,
     902     the F64 port might be.
     903 
     904XXX: Must update the covar matrix.
     905 *****************************************************************************/
     906psBool psMinimizeLMChi2Old(psMinimization *min,
     907                           psImage *covar,
     908                           psVector *params,
     909                           const psVector *paramMask,
     910                           const psArray *x,
     911                           const psVector *y,
     912                           const psVector *yErr,
     913                           psMinimizeLMChi2Func func)
     914{
     915    PS_PTR_CHECK_NULL(min, NULL);
     916    PS_VECTOR_CHECK_NULL(params, NULL);
     917    PS_VECTOR_CHECK_EMPTY(params, NULL);
     918    PS_PTR_CHECK_NULL(x, NULL);
     919    PS_VECTOR_CHECK_NULL(y, NULL);
     920    PS_VECTOR_CHECK_EMPTY(y, NULL);
     921    PS_VECTOR_CHECK_SIZE_EQUAL(x, y, NULL);
     922    PS_PTR_CHECK_NULL(func, NULL);
     923
    603924    if (paramMask != NULL) {
    604925        PS_VECTOR_CHECK_SIZE_EQUAL(params, paramMask, NULL);
    … …  
    636957    psF32 currChi2 = 0.0;
    637958    psF32 newChi2 = 0.0;
    638     psF32 lamda = 0.00005;
    639     lamda = 0.05;
     959    psF32 lamda = 0.00005;  // XXX EAM : this starting value is VERY small (lamda is mis-spelt)
     960    lamda = 0.05;  // XXX EAM : this starting value is quite large (lamda is mis-spelt)
    640961
    641962    psTrace(".psLib.dataManip.psMinimize", 6,
    … …  
    661982        //
    662983        currChi2 = 0.0;
    663         currValueVec = func(deriv, params, x);
     984        // currValueVec = func(deriv, params, x);
     985
     986        // XXX EAM: use BinaryOp ?
     987        // t1 = BinaryOp (NULL, currValueVec, "-", y);
     988        // t1 = BinaryOp (t1, t1, "*", t1);
     989
     990        // XXX EAM: this ignores yErr
    664991        for (n=0;n<numData;n++) {
    665992            currChi2+= (currValueVec->data.F32[n] - y->data.F32[n]) *
    … …  
    670997        }
    671998
     999        // XXX EAM: this is just for tracing
    6721000        for (p=0;p<numParams;p++) {
    6731001            psTrace(".psLib.dataManip.psMinimize", 6,
    … …  
    6791007        //
    6801008        // Mask elements of the derivative for each data point.
    681         //
     1009        // XXX EAM : is this necessary?  probably not...
    6821010        for (p=0;p<numParams;p++) {
    6831011            if ((paramMask != NULL) && (paramMask->data.U8[p] != 0)) {
    … …  
    6901018        //
    6911019        // Calculate the BETA vector.
    692         //
     1020        // XXX EAM: I think this is wrong
    6931021        for (p=0;p<numParams;p++) {
     1022            if ((paramMask != NULL) && (paramMask->data.U8[p] != 0)) {
     1023                continue;
     1024            }
    6941025            beta->data.F64[p] = 0.0;
    6951026            for (n=0;n<numData;n++) {
    … …  
    7071038        //
    7081039        // Calculate the ALPHA matrix.
    709         //
     1040        // XXX EAM: also wrong? (missing yErr)
    7101041        for (k=0;k<numParams;k++) {
    7111042            for (l=0;l<numParams;l++) {
    … …  
    7731104        //
    7741105        newChi2 = 0.0;
    775         newValueVec = func(deriv, newParams, x);
     1106        // newValueVec = func(deriv, newParams, x);
    7761107        for (n=0;n<numData;n++) {
    7771108            newChi2+= (newValueVec->data.F32[n] - y->data.F32[n]) *
    … …  
    10981429    min->value = 0.0;
    10991430    min->iter = 0;
    1100     min->lastDelta = 0.0;
     1431    min->lastDelta = tol + 1;
    11011432
    11021433    return(min);

  • branches/eam-psphot-branch/psLib/src/dataManip/psMinimize.h

    r3264 r3788  
    88 *  @author GLG, MHPCC
    99 *
    10  *  @version $Revision: 1.39 $ $Name: not supported by cvs2svn $
    11  *  @date $Date: 2005-02-17 19:26:23 $
     10 *  @version $Revision: 1.39.6.1 $ $Name: not supported by cvs2svn $
     11 *  @date $Date: 2005-04-29 09:17:02 $
    1212 *
    1313 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    4242    psF32 tol;                         ///< Error Tolerance
    4343    psF32 value;                       ///< Value of function at minimum
    44     psS32 iter;                          ///< Number of iterations required
     44    psS32 iter;                          ///< Number of iterations to date
    4545    psF32 lastDelta;                   ///< The last difference for the fit
    4646}
    … …  
    7777
    7878typedef
    79 psVector* (*psMinimizeLMChi2Func)(psImage *deriv,
    80                                   const psVector *params,
    81                                   const psArray *x);
     79psF64 (*psMinimizeLMChi2Func)(psVector *deriv,
     80                              psVector *params,
     81                              psVector *x);
    8282
    8383psBool psMinimizeLMChi2(psMinimization *min,
    … …  
    8888                        const psVector *y,
    8989                        const psVector *yErr,
     90                        psMinimizeLMChi2Func func);
     91
     92psBool p_psMinLM_GuessABP (psImage  *Alpha,
     93                           psVector *Beta,
     94                           psVector *Params,
     95                           psImage  *alpha,
     96                           psVector *beta,
     97                           psVector *params,
     98                           psF64 lambda);
     99
     100psF64 p_psMinLM_SetABX (psImage  *alpha,
     101                        psVector *beta,
     102                        psVector *params,
     103                        const psArray  *x,
     104                        const psVector *y,
     105                        const psVector *dy,
    90106                        psMinimizeLMChi2Func func);
    91107
    … …  
    119135                            psMinimizeChi2PowellFunc func);
    120136
    121 
     137// XXX EAM : psGaussJordan provided as an alternate to LU Decomp for psMinimizeLMChi2
     138bool psGaussJordan (psImage *a, psVector *b);
    122139
    123140/* \} */// End of MathGroup Functions

  • branches/eam-psphot-branch/psLib/src/image/psImage.c

    r3446 r3788  
    99 *  @author Ross Harman, MHPCC
    1010 *
    11  *  @version $Revision: 1.61 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2005-03-18 02:35:14 $
     11 *  @version $Revision: 1.61.4.1 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2005-04-29 09:17:02 $
    1313 *
    1414 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    412412
    413413    return numFreed;
     414}
     415
     416// XXX EAM some utility functions
     417bool p_psImagePrint (FILE *f, psImage *a, char *name)
     418{
     419
     420    fprintf (f, "matrix: %s\n", name);
     421
     422    for (int j = 0; j < a[0].numRows; j++) {
     423        for (int i = 0; i < a[0].numCols; i++) {
     424            fprintf (f, "%f  ", p_psImageGetElementF64(a, i, j));
     425        }
     426        fprintf (f, "\n");
     427    }
     428    fprintf (f, "\n");
     429    return (true);
     430}
     431
     432psF64 p_psImageGetElementF64(psImage* image,
     433                             int col,
     434                             int row)
     435{
     436    if (image == NULL) {
     437        return NAN;
     438    }
     439    if (col < 0 || col >= image->numCols) {
     440        return NAN;
     441    }
     442    if (row < 0 || row >= image->numRows) {
     443        return NAN;
     444    }
     445
     446    switch (image->type.type) {
     447    case PS_TYPE_U8:
     448        return image->data.U8[row][col];
     449        break;
     450    case PS_TYPE_U16:
     451        return image->data.U16[row][col];
     452        break;
     453    case PS_TYPE_U32:
     454        return image->data.U32[row][col];
     455        break;
     456    case PS_TYPE_U64:
     457        return image->data.U64[row][col];
     458        break;
     459    case PS_TYPE_S8:
     460        return image->data.S8[row][col];
     461        break;
     462    case PS_TYPE_S16:
     463        return image->data.S16[row][col];
     464        break;
     465    case PS_TYPE_S32:
     466        return image->data.S32[row][col];
     467        break;
     468    case PS_TYPE_S64:
     469        return image->data.S64[row][col];
     470        break;
     471    case PS_TYPE_F32:
     472        return image->data.F32[row][col];
     473        break;
     474    case PS_TYPE_F64:
     475        return image->data.F64[row][col];
     476    default:
     477        return NAN;
     478    }
    414479}
    415480

  • branches/eam-psphot-branch/psLib/src/image/psImage.h

    r3446 r3788  
    1111 *  @author Ross Harman, MHPCC
    1212 *
    13  *  @version $Revision: 1.49 $ $Name: not supported by cvs2svn $
    14  *  @date $Date: 2005-03-18 02:35:14 $
     13 *  @version $Revision: 1.49.4.1 $ $Name: not supported by cvs2svn $
     14 *  @date $Date: 2005-04-29 09:17:02 $
    1515 *
    1616 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    2020
    2121#include <complex.h>
     22#include <stdio.h>
    2223
    2324#include "psType.h"
    … …  
    191192);
    192193
     194
     195// XXX EAM some utility functions
     196
     197psF64 p_psImageGetElementF64(psImage *a, int i, int j);
     198bool p_psImagePrint (FILE *f, psImage *a, char *name);
     199
    193200#define PIXEL_INTERPOLATE_FCN_PROTOTYPE(SUFFIX, RETURNTYPE) \
    194201inline RETURNTYPE p_psImagePixelInterpolate##SUFFIX( \

  • branches/eam-psphot-branch/psLib/src/math/psMinimize.c

    r3578 r3788  
    99 *  @author GLG, MHPCC
    1010 *
    11  *  @version $Revision: 1.110 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2005-03-31 01:02:15 $
     11 *  @version $Revision: 1.110.6.1 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2005-04-29 09:17:02 $
    1313 *
    1414 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    557557}
    558558
    559 /******************************************************************************
    560 psMinimizeLMChi2():  This routine will take an procedure which calculates
    561 an arbitrary function and it's derivative and minimize the chi-squared match
    562 between that function at the specified coords and the specified value at
    563 those coords.
    564  
    565 XXX: Do this:
    566  After checking that all entries in the paramMask are 1 or 0, when
    567  forming the A matrix from alpha, try this:
    568  
    569      A[i][i] = (1 + lambda*paramask[i]) * alpha[i][i];
    570  
    571 XXX: This is very different from what is specified in the SDR.  Must
    572 coordinate with IfA on new SDR.
    573  
    574 XXX: Do vector/image recycles.
    575  
    576 XXX: probably yErr will be part of the SDR.
    577  
    578 XXX: This must work for both F32 and F64.  F32 is currently implemented.
    579      Note: since the LUD routines are only implemented in F64, then we
    580      will have to convert all F32 input vectors to F64 regardless.  So,
    581      the F64 port might be.
    582  
    583 XXX: Must update the covar matrix.
    584  *****************************************************************************/
     559// XXX EAM this is my re-implementation of MinLM
    585560psBool psMinimizeLMChi2(psMinimization *min,
    586561                        psImage *covar,
    … …  
    601576    PS_PTR_CHECK_NULL(func, NULL);
    602577
     578    // this function has test and current values for several things
     579    // the current best value is in lower case
     580    // the next guess value is in upper case
     581
     582    // allocate internal arrays (current vs Guess)
     583    psImage *alpha   = psImageAlloc (params->n, params->n, PS_TYPE_F64);
     584    psImage *Alpha   = psImageAlloc (params->n, params->n, PS_TYPE_F64);
     585    psVector *beta   = psVectorAlloc (params->n, PS_TYPE_F64);
     586    psVector *Beta   = psVectorAlloc (params->n, PS_TYPE_F64);
     587    psVector *Params = psVectorAlloc (params->n, PS_TYPE_F64);
     588    psVector *dy     = NULL;
     589    psF64 chisq = 0.0;
     590    psF64 Chisq = 0.0;
     591    psF64 lambda = 0.001;
     592
     593    // the initial guess on params is provided by the user
     594    Params = psVectorCopy (Params, params, PS_TYPE_F32);
     595
     596    // the user provides the error or NULL.  we need to convert
     597    // to appropriate weights
     598    dy = psVectorAlloc (y->n, PS_TYPE_F32);
     599    if (yErr != NULL) {
     600        for (int i = 0; i < dy->n; i++) {
     601            dy->data.F32[i] = 1.0 / PS_SQR (yErr->data.F32[i]);
     602        }
     603    } else {
     604        for (int i = 0; i < dy->n; i++) {
     605            dy->data.F32[i] = 1.0;
     606        }
     607    }
     608
     609    // calculate initial alpha and beta, set chisq (min->value)
     610    min->value = p_psMinLM_SetABX (alpha, beta, params, x, y, dy, func);
     611    # ifndef PS_NO_TRACE
     612    // dump some useful info if trace is defined
     613    if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") > 4) {
     614        p_psImagePrint  (psTraceGetDestination(), alpha, "alpha guess");
     615        p_psVectorPrint (psTraceGetDestination(), beta, "beta guess");
     616        p_psVectorPrint (psTraceGetDestination(), params, "params guess");
     617    }
     618    # endif /* PS_NO_TRACE */
     619
     620
     621    // iterate until the tolerance is reached, or give up
     622    while ((min->lastDelta > min->tol) && (min->iter < min->maxIter)) {
     623
     624        // set a new guess for Alpha, Beta, Params
     625        p_psMinLM_GuessABP (Alpha, Beta, Params, alpha, beta, params, lambda);
     626
     627        # ifndef PS_NO_TRACE
     628        // dump some useful info if trace is defined
     629        if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") > 4) {
     630            p_psImagePrint  (psTraceGetDestination(), Alpha, "alpha guess");
     631            p_psVectorPrint (psTraceGetDestination(), Beta, "beta guess");
     632            p_psVectorPrint (psTraceGetDestination(), Params, "params guess");
     633        }
     634        # endif /* PS_NO_TRACE */
     635
     636        // calculate Chisq for new guess, update Alpha & Beta
     637        Chisq = p_psMinLM_SetABX (Alpha, Beta, Params, x, y, dy, func);
     638        psTrace (".psLib.dataManip.psMinimizeLMChi2", 3, "chisq: %f, Chisq %f, delta: %f\n", chisq, Chisq, min->lastDelta);
     639
     640        // accept new guess (if improvement), or increase lambda
     641        if (Chisq < min->value) {
     642            min->lastDelta = (min->value - Chisq) / (dy->n - params->n);
     643            min->value = Chisq;
     644            alpha  = psImageCopy (alpha, Alpha, PS_TYPE_F64);
     645            beta   = psVectorCopy (beta, Beta, PS_TYPE_F64);
     646            params = psVectorCopy (params, Params, PS_TYPE_F32);
     647            lambda *= 0.1;
     648        } else {
     649            lambda *= 10.0;
     650        }
     651        min->iter ++;
     652    }
     653    psTrace (".psLib.dataManip.psMinimizeLMChi2", 3, "chisq: %f, Chisq %f, delta: %f\n", chisq, Chisq, min->lastDelta);
     654
     655    // free the internal temporary data
     656    psFree (alpha);
     657    psFree (Alpha);
     658    psFree (beta);
     659    psFree (Beta);
     660    psFree (Params);
     661    psFree (dy);
     662    return (true);
     663}
     664
     665// XXX EAM: this needs to respect the mask on params
     666// XXX EAM: check not NULL on alpha, beta, params
     667// alpha, beta, params are already allocated
     668psF64 p_psMinLM_SetABX (psImage  *alpha,
     669                        psVector *beta,
     670                        psVector *params,
     671                        const psArray  *x,
     672                        const psVector *y,
     673                        const psVector *dy,
     674                        psMinimizeLMChi2Func func)
     675{
     676
     677    psF64 chisq;
     678    psF64 delta;
     679    psF64 weight;
     680    psF64 ymodel;
     681    psVector *deriv = psVectorAlloc (params->n, PS_TYPE_F32);
     682
     683    // zero alpha and beta for summing below
     684    for (int j = 0; j < params->n; j++) {
     685        for (int k = 0; k < params->n; k++) {
     686            alpha->data.F64[j][k] = 0;
     687        }
     688        beta->data.F64[j] = 0;
     689    }
     690    chisq = 0.0;
     691
     692    // calculate chisq, alpha, beta
     693    for (int i = 0; i < y->n; i++) {
     694        ymodel = func (deriv, params, (psVector *) x->data[i]);
     695
     696        delta = ymodel - y->data.F32[i];
     697        chisq += PS_SQR (delta) * dy->data.F32[i];
     698
     699        for (int j = 0; j < params->n; j++) {
     700            weight = deriv->data.F32[j] * dy->data.F32[i];
     701            for (int k = 0; k <= j; k++) {
     702                alpha->data.F64[j][k] += weight * deriv->data.F32[k];
     703            }
     704            beta->data.F64[j] += weight * delta;
     705        }
     706    }
     707
     708    // calculate lower-left half of alpha
     709    for (int j = 1; j < params->n; j++) {
     710        for (int k = 0; k < j; k++) {
     711            alpha->data.F64[k][j] = alpha->data.F64[j][k];
     712        }
     713    }
     714    psFree (deriv);
     715    return (chisq);
     716}
     717
     718// XXX EAM : can we use static copies of LUv, LUm, A?
     719psBool p_psMinLM_GuessABP (psImage  *Alpha,
     720                           psVector *Beta,
     721                           psVector *Params,
     722                           psImage  *alpha,
     723                           psVector *beta,
     724                           psVector *params,
     725                           psF64 lambda)
     726{
     727
     728    # define USE_LU_DECOMP 1
     729    # if (USE_LU_DECOMP)
     730        psVector *LUv = NULL;
     731    psImage  *LUm = NULL;
     732    psImage  *A   = NULL;
     733    psF32    det;
     734
     735    // LU decomposition version
     736    psTrace (".pslib.dataManip.psMinLM_GuessABP", 3, "using LUD version");
     737
     738    // set new guess values (creates matrix A)
     739    A = psImageCopy (NULL, alpha, PS_TYPE_F64);
     740    for (int j = 0; j < params->n; j++) {
     741        A->data.F64[j][j] = alpha->data.F64[j][j] * (1.0 + lambda);
     742    }
     743
     744    // solve A*beta = Beta (Alpha = 1/A)
     745    // these operations do not modify the input values (creates LUm, LUv)
     746    LUm   = psMatrixLUD (NULL, &LUv, A);
     747    Beta  = psMatrixLUSolve (Beta, LUm, beta, LUv);
     748    Alpha = psMatrixInvert (Alpha, A, &det);
     749
     750    # else
     751        // gauss-jordan version
     752        psTrace (".pslib.dataManip.psMinLM_GuessABP", 3, "using Gauss-J version");
     753
     754    // set new guess values (creates matrix A)
     755    Beta = psVectorCopy (Beta, beta, PS_TYPE_F64);
     756    Alpha = psImageCopy (Alpha, alpha, PS_TYPE_F64);
     757    for (int j = 0; j < params->n; j++) {
     758        Alpha->data.F64[j][j] = alpha->data.F64[j][j] * (1.0 + lambda);
     759    }
     760
     761    psGaussJordan (Alpha, Beta);
     762    # endif
     763
     764    // apply beta to get new params values
     765    for (int j = 0; j < params->n; j++) {
     766        Params->data.F32[j] = params->data.F32[j] - Beta->data.F64[j];
     767    }
     768
     769    # if (USE_LU_DECOMP)
     770        psFree (A);
     771    psFree (LUm);
     772    psFree (LUv);
     773    # endif
     774
     775    return true;
     776}
     777
     778# define SWAP(X,Y) {double tmp=(X); (X) = (Y); (Y) = tmp;}
     779
     780// XXX EAM : temporary gauss-jordan solver based on gene's
     781// version based on the Numerical Recipes version
     782bool psGaussJordan (psImage *a, psVector *b)
     783{
     784
     785    int *indxc,*indxr,*ipiv;
     786    int Nx, icol, irow;
     787    int i, j, k, l, ll;
     788    float big, dum, pivinv;
     789    psF64 *vector;
     790    psF64 **matrix;
     791
     792    Nx = a->numCols;
     793    matrix = a->data.F64;
     794    vector = b->data.F64;
     795
     796    indxc = psAlloc (Nx*sizeof(int));
     797    indxr = psAlloc (Nx*sizeof(int));
     798    ipiv  = psAlloc (Nx*sizeof(int));
     799    for (j = 0; j < Nx; j++)
     800        ipiv[j] = 0;
     801
     802    irow = icol = 0;
     803    big = fabs(matrix[0][0]);
     804
     805    for (i = 0; i < Nx; i++) {
     806        big = 0.0;
     807        for (j = 0; j < Nx; j++) {
     808            if (!finite(matrix[i][j])) {
     809                // XXX EAM: this should use the psError stack
     810                fprintf (stderr, "GAUSSJ: NaN\n");
     811                goto fescape;
     812            }
     813            if (ipiv[j] != 1) {
     814                for (k = 0; k < Nx; k++) {
     815                    if (ipiv[k] == 0) {
     816                        if (fabs (matrix[j][k]) >= big) {
     817                            big  = fabs (matrix[j][k]);
     818                            irow = j;
     819                            icol = k;
     820                        }
     821                    } else {
     822                        if (ipiv[k] > 1) {
     823                            // XXX EAM: this should use the psError stack
     824                            fprintf (stderr, "GAUSSJ: Singular Matrix! (1)\n");
     825                            goto fescape;
     826                        }
     827                    }
     828                }
     829            }
     830        }
     831        ipiv[icol]++;
     832        if (irow != icol) {
     833            for (l = 0; l < Nx; l++) {
     834                SWAP (matrix[irow][l], matrix[icol][l]);
     835            }
     836            SWAP (vector[irow], vector[icol]);
     837        }
     838        indxr[i] = irow;
     839        indxc[i] = icol;
     840        if (matrix[icol][icol] == 0.0) {
     841            // XXX EAM: this should use the psError stack
     842            fprintf (stderr, "GAUSSJ: Singular Matrix! (2)\n");
     843            goto fescape;
     844        }
     845        pivinv = 1.0 / matrix[icol][icol];
     846        matrix[icol][icol] = 1.0;
     847        for (l = 0; l < Nx; l++) {
     848            matrix[icol][l] *= pivinv;
     849        }
     850        vector[icol] *= pivinv;
     851
     852        for (ll = 0; ll < Nx; ll++) {
     853            if (ll != icol) {
     854                dum = matrix[ll][icol];
     855                matrix[ll][icol] = 0.0;
     856                for (l = 0; l < Nx; l++)
     857                    matrix[ll][l] -= matrix[icol][l]*dum;
     858                vector[ll] -= vector[icol]*dum;
     859            }
     860        }
     861    }
     862
     863    for (l = Nx - 1; l >= 0; l--) {
     864        if (indxr[l] != indxc[l])
     865            for (k = 0; k < Nx; k++)
     866                SWAP (matrix[k][indxr[l]], matrix[k][indxc[l]]);
     867    }
     868    psFree (ipiv);
     869    psFree (indxr);
     870    psFree (indxc);
     871    return (true);
     872
     873fescape:
     874    psFree (ipiv);
     875    psFree (indxr);
     876    psFree (indxc);
     877    return (false);
     878}
     879
     880/******************************************************************************
     881psMinimizeLMChi2():  This routine will take an procedure which calculates
     882an arbitrary function and it's derivative and minimize the chi-squared match
     883between that function at the specified coords and the specified value at
     884those coords.
     885 
     886XXX: Do this:
     887 After checking that all entries in the paramMask are 1 or 0, when
     888 forming the A matrix from alpha, try this:
     889 
     890     A[i][i] = (1 + lambda*paramask[i]) * alpha[i][i];
     891 
     892XXX: This is very different from what is specified in the SDR.  Must
     893coordinate with IfA on new SDR.
     894 
     895XXX: Do vector/image recycles.
     896 
     897XXX: probably yErr will be part of the SDR.
     898 
     899XXX: This must work for both F32 and F64.  F32 is currently implemented.
     900     Note: since the LUD routines are only implemented in F64, then we
     901     will have to convert all F32 input vectors to F64 regardless.  So,
     902     the F64 port might be.
     903 
     904XXX: Must update the covar matrix.
     905 *****************************************************************************/
     906psBool psMinimizeLMChi2Old(psMinimization *min,
     907                           psImage *covar,
     908                           psVector *params,
     909                           const psVector *paramMask,
     910                           const psArray *x,
     911                           const psVector *y,
     912                           const psVector *yErr,
     913                           psMinimizeLMChi2Func func)
     914{
     915    PS_PTR_CHECK_NULL(min, NULL);
     916    PS_VECTOR_CHECK_NULL(params, NULL);
     917    PS_VECTOR_CHECK_EMPTY(params, NULL);
     918    PS_PTR_CHECK_NULL(x, NULL);
     919    PS_VECTOR_CHECK_NULL(y, NULL);
     920    PS_VECTOR_CHECK_EMPTY(y, NULL);
     921    PS_VECTOR_CHECK_SIZE_EQUAL(x, y, NULL);
     922    PS_PTR_CHECK_NULL(func, NULL);
     923
    603924    if (paramMask != NULL) {
    604925        PS_VECTOR_CHECK_SIZE_EQUAL(params, paramMask, NULL);
    … …  
    636957    psF32 currChi2 = 0.0;
    637958    psF32 newChi2 = 0.0;
    638     psF32 lamda = 0.00005;
    639     lamda = 0.05;
     959    psF32 lamda = 0.00005;  // XXX EAM : this starting value is VERY small (lamda is mis-spelt)
     960    lamda = 0.05;  // XXX EAM : this starting value is quite large (lamda is mis-spelt)
    640961
    641962    psTrace(".psLib.dataManip.psMinimize", 6,
    … …  
    661982        //
    662983        currChi2 = 0.0;
    663         currValueVec = func(deriv, params, x);
     984        // currValueVec = func(deriv, params, x);
     985
     986        // XXX EAM: use BinaryOp ?
     987        // t1 = BinaryOp (NULL, currValueVec, "-", y);
     988        // t1 = BinaryOp (t1, t1, "*", t1);
     989
     990        // XXX EAM: this ignores yErr
    664991        for (n=0;n<numData;n++) {
    665992            currChi2+= (currValueVec->data.F32[n] - y->data.F32[n]) *
    … …  
    670997        }
    671998
     999        // XXX EAM: this is just for tracing
    6721000        for (p=0;p<numParams;p++) {
    6731001            psTrace(".psLib.dataManip.psMinimize", 6,
    … …  
    6791007        //
    6801008        // Mask elements of the derivative for each data point.
    681         //
     1009        // XXX EAM : is this necessary?  probably not...
    6821010        for (p=0;p<numParams;p++) {
    6831011            if ((paramMask != NULL) && (paramMask->data.U8[p] != 0)) {
    … …  
    6901018        //
    6911019        // Calculate the BETA vector.
    692         //
     1020        // XXX EAM: I think this is wrong
    6931021        for (p=0;p<numParams;p++) {
     1022            if ((paramMask != NULL) && (paramMask->data.U8[p] != 0)) {
     1023                continue;
     1024            }
    6941025            beta->data.F64[p] = 0.0;
    6951026            for (n=0;n<numData;n++) {
    … …  
    7071038        //
    7081039        // Calculate the ALPHA matrix.
    709         //
     1040        // XXX EAM: also wrong? (missing yErr)
    7101041        for (k=0;k<numParams;k++) {
    7111042            for (l=0;l<numParams;l++) {
    … …  
    7731104        //
    7741105        newChi2 = 0.0;
    775         newValueVec = func(deriv, newParams, x);
     1106        // newValueVec = func(deriv, newParams, x);
    7761107        for (n=0;n<numData;n++) {
    7771108            newChi2+= (newValueVec->data.F32[n] - y->data.F32[n]) *
    … …  
    10981429    min->value = 0.0;
    10991430    min->iter = 0;
    1100     min->lastDelta = 0.0;
     1431    min->lastDelta = tol + 1;
    11011432
    11021433    return(min);

  • branches/eam-psphot-branch/psLib/src/math/psMinimize.h

    r3264 r3788  
    88 *  @author GLG, MHPCC
    99 *
    10  *  @version $Revision: 1.39 $ $Name: not supported by cvs2svn $
    11  *  @date $Date: 2005-02-17 19:26:23 $
     10 *  @version $Revision: 1.39.6.1 $ $Name: not supported by cvs2svn $
     11 *  @date $Date: 2005-04-29 09:17:02 $
    1212 *
    1313 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    4242    psF32 tol;                         ///< Error Tolerance
    4343    psF32 value;                       ///< Value of function at minimum
    44     psS32 iter;                          ///< Number of iterations required
     44    psS32 iter;                          ///< Number of iterations to date
    4545    psF32 lastDelta;                   ///< The last difference for the fit
    4646}
    … …  
    7777
    7878typedef
    79 psVector* (*psMinimizeLMChi2Func)(psImage *deriv,
    80                                   const psVector *params,
    81                                   const psArray *x);
     79psF64 (*psMinimizeLMChi2Func)(psVector *deriv,
     80                              psVector *params,
     81                              psVector *x);
    8282
    8383psBool psMinimizeLMChi2(psMinimization *min,
    … …  
    8888                        const psVector *y,
    8989                        const psVector *yErr,
     90                        psMinimizeLMChi2Func func);
     91
     92psBool p_psMinLM_GuessABP (psImage  *Alpha,
     93                           psVector *Beta,
     94                           psVector *Params,
     95                           psImage  *alpha,
     96                           psVector *beta,
     97                           psVector *params,
     98                           psF64 lambda);
     99
     100psF64 p_psMinLM_SetABX (psImage  *alpha,
     101                        psVector *beta,
     102                        psVector *params,
     103                        const psArray  *x,
     104                        const psVector *y,
     105                        const psVector *dy,
    90106                        psMinimizeLMChi2Func func);
    91107
    … …  
    119135                            psMinimizeChi2PowellFunc func);
    120136
    121 
     137// XXX EAM : psGaussJordan provided as an alternate to LU Decomp for psMinimizeLMChi2
     138bool psGaussJordan (psImage *a, psVector *b);
    122139
    123140/* \} */// End of MathGroup Functions

  • branches/eam-psphot-branch/psLib/src/mathtypes/psImage.c

    r3446 r3788  
    99 *  @author Ross Harman, MHPCC
    1010 *
    11  *  @version $Revision: 1.61 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2005-03-18 02:35:14 $
     11 *  @version $Revision: 1.61.4.1 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2005-04-29 09:17:02 $
    1313 *
    1414 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    412412
    413413    return numFreed;
     414}
     415
     416// XXX EAM some utility functions
     417bool p_psImagePrint (FILE *f, psImage *a, char *name)
     418{
     419
     420    fprintf (f, "matrix: %s\n", name);
     421
     422    for (int j = 0; j < a[0].numRows; j++) {
     423        for (int i = 0; i < a[0].numCols; i++) {
     424            fprintf (f, "%f  ", p_psImageGetElementF64(a, i, j));
     425        }
     426        fprintf (f, "\n");
     427    }
     428    fprintf (f, "\n");
     429    return (true);
     430}
     431
     432psF64 p_psImageGetElementF64(psImage* image,
     433                             int col,
     434                             int row)
     435{
     436    if (image == NULL) {
     437        return NAN;
     438    }
     439    if (col < 0 || col >= image->numCols) {
     440        return NAN;
     441    }
     442    if (row < 0 || row >= image->numRows) {
     443        return NAN;
     444    }
     445
     446    switch (image->type.type) {
     447    case PS_TYPE_U8:
     448        return image->data.U8[row][col];
     449        break;
     450    case PS_TYPE_U16:
     451        return image->data.U16[row][col];
     452        break;
     453    case PS_TYPE_U32:
     454        return image->data.U32[row][col];
     455        break;
     456    case PS_TYPE_U64:
     457        return image->data.U64[row][col];
     458        break;
     459    case PS_TYPE_S8:
     460        return image->data.S8[row][col];
     461        break;
     462    case PS_TYPE_S16:
     463        return image->data.S16[row][col];
     464        break;
     465    case PS_TYPE_S32:
     466        return image->data.S32[row][col];
     467        break;
     468    case PS_TYPE_S64:
     469        return image->data.S64[row][col];
     470        break;
     471    case PS_TYPE_F32:
     472        return image->data.F32[row][col];
     473        break;
     474    case PS_TYPE_F64:
     475        return image->data.F64[row][col];
     476    default:
     477        return NAN;
     478    }
    414479}
    415480

  • branches/eam-psphot-branch/psLib/src/mathtypes/psImage.h

    r3446 r3788  
    1111 *  @author Ross Harman, MHPCC
    1212 *
    13  *  @version $Revision: 1.49 $ $Name: not supported by cvs2svn $
    14  *  @date $Date: 2005-03-18 02:35:14 $
     13 *  @version $Revision: 1.49.4.1 $ $Name: not supported by cvs2svn $
     14 *  @date $Date: 2005-04-29 09:17:02 $
    1515 *
    1616 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    2020
    2121#include <complex.h>
     22#include <stdio.h>
    2223
    2324#include "psType.h"
    … …  
    191192);
    192193
     194
     195// XXX EAM some utility functions
     196
     197psF64 p_psImageGetElementF64(psImage *a, int i, int j);
     198bool p_psImagePrint (FILE *f, psImage *a, char *name);
     199
    193200#define PIXEL_INTERPOLATE_FCN_PROTOTYPE(SUFFIX, RETURNTYPE) \
    194201inline RETURNTYPE p_psImagePixelInterpolate##SUFFIX( \

  • branches/eam-psphot-branch/psLib/src/mathtypes/psVector.c

    r3476 r3788  
    99*  @author Robert DeSonia, MHPCC
    1010*
    11 *  @version $Revision: 1.36 $ $Name: not supported by cvs2svn $
    12 *  @date $Date: 2005-03-22 21:52:49 $
     11*  @version $Revision: 1.36.4.1 $ $Name: not supported by cvs2svn $
     12*  @date $Date: 2005-04-29 09:17:02 $
    1313*
    1414*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    540540    }
    541541}
     542
     543// XXX EAM a utility function
     544bool p_psVectorPrint (FILE *f, psVector *a, char *name)
     545{
     546
     547    fprintf (f, "vector: %s\n", name);
     548
     549    for (int i = 0; i < a[0].n; i++) {
     550        fprintf (f, "%f\n", p_psVectorGetElementF64(a, i));
     551    }
     552    fprintf (f, "\n");
     553    return (true);
     554}
     555

  • branches/eam-psphot-branch/psLib/src/mathtypes/psVector.h

    r3476 r3788  
    1111 *  @author Ross Harman, MHPCC
    1212 *
    13  *  @version $Revision: 1.30 $ $Name: not supported by cvs2svn $
    14  *  @date $Date: 2005-03-22 21:52:49 $
     13 *  @version $Revision: 1.30.4.1 $ $Name: not supported by cvs2svn $
     14 *  @date $Date: 2005-04-29 09:17:02 $
    1515 *
    1616 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    160160);
    161161
     162// XXX EAM a utility function
     163bool p_psVectorPrint (FILE *f, psVector *a, char *name);
    162164
    163165/// @}

  • branches/eam-psphot-branch/psLib/src/sys/psTrace.c

    r3541 r3788  
    99 *  @author GLG, MHPCC
    1010 *
    11  *  @version $Revision: 1.44 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2005-03-29 21:31:54 $
     11 *  @version $Revision: 1.44.4.1 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2005-04-29 09:17:02 $
    1313 *
    1414 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    252252    if (cRoot == NULL) {
    253253        initTrace();
     254    }
     255
     256    if (traceFP == NULL) {
     257        traceFP = stderr;
    254258    }
    255259
    … …  
    425429    psS32 i = 0;
    426430
     431    if (traceFP == NULL) {
     432        traceFP = stderr;
     433    }
     434
    427435    if (comp->name[0] == '\0') {
    428436        return;
    429437    } else {
    430438        if (comp->level == PS_DEFAULT_TRACE_LEVEL) {
    431             if (traceFP == NULL) {
    432                 printf("%*s%-*s %d\n", depth, "", 20 - depth, comp->name,
    433                        defLevel);
    434             } else {
    435                 fprintf(traceFP,"%*s%-*s %d\n", depth, "", 20 - depth, comp->name,
    436                         defLevel);
    437             }
     439            fprintf(traceFP,"%*s%-*s %d\n", depth, "", 20 - depth, comp->name, defLevel);
    438440        } else {
    439             if (traceFP == NULL) {
    440                 printf("%*s%-*s %d\n", depth, "", 20 - depth, comp->name,
    441                        comp->level);
    442             } else {
    443                 fprintf(traceFP, "%*s%-*s %d\n", depth, "", 20 - depth, comp->name,
    444                         comp->level);
    445             }
     441            fprintf(traceFP, "%*s%-*s %d\n", depth, "", 20 - depth, comp->name, comp->level);
    446442        }
    447443    }
    … …  
    498494    psS32 i = 0;
    499495
     496    if (traceFP == NULL) {
     497        traceFP = stderr;
     498    }
     499
    500500    if (NULL == comp) {
    501501        psError(PS_ERR_BAD_PARAMETER_NULL, true,
    … …  
    514514        fmt = va_arg(ap, char *);
    515515
    516         if (traceFP != NULL) {
    517             // We indent each message one space for each level of the message.
    518             for (i = 0; i < level; i++) {
    519                 fprintf(traceFP, " ");
    520             }
    521             vfprintf(traceFP, fmt, ap);
    522         } else {
    523             // We indent each message one space for each level of the message.
    524             for (i = 0; i < level; i++) {
    525                 putchar(' ');
    526             }
    527             vprintf(fmt, ap);
    528         }
     516        // We indent each message one space for each level of the message.
     517        for (i = 0; i < level; i++) {
     518            fprintf(traceFP, " ");
     519        }
     520        vfprintf(traceFP, fmt, ap);
    529521        va_end(ap);
    530522    }
    531     // NOTE: should we free *fmt as well? Read the man page.
    532 }
    533 
     523}
     524
     525// XXX EAM : I've added code to close the old traceFP (safely)
    534526void psTraceSetDestination(FILE * fp)
    535527{
     528
     529    bool special;
     530
     531    // XXX EAM perhaps return an error?
     532    if (fp == NULL) {
     533        return;
     534    }
     535
     536    // cannot close traceFP if one of the special FILE ptrs
     537    special  = (traceFP == NULL);
     538    special |= (traceFP == stdin);
     539    special |= (traceFP == stdout);
     540    special |= (traceFP == stderr);
     541
     542    if (!special) {
     543        fclose (traceFP);
     544    }
    536545    traceFP = fp;
    537546}
    538547
     548FILE *psTraceGetDestination()
     549{
     550    if (traceFP == NULL) {
     551        traceFP = stderr;
     552    }
     553    return traceFP;
     554}
     555
    539556#endif

  • branches/eam-psphot-branch/psLib/src/sys/psTrace.h

    r3264 r3788  
    99 *  @author GLG, MHPCC
    1010 *
    11  *  @version $Revision: 1.31 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2005-02-17 19:26:24 $
     11 *  @version $Revision: 1.31.6.1 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2005-04-29 09:17:02 $
    1313 *
    1414 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    2929//#define PS_NO_TRACE 1   ///< to turn off all tracing
    3030
    31 #if defined(PS_NO_TRACE)
    32 #        define psTrace(facil, level, ...) (void)0
    33 /* do nothing */
    34 #        define p_psTrace(facil, level, ...)  (void)0
    35 /* do nothing */
    36 #        define psTraceSetLevel(facil,level) 0
    37 #        define psTraceGetLevel(facil) 0
    38 #        define psTraceReset() (void)0     /* do nothing */
    39 #        define psTraceFree() (void)0      /* do nothing */
    40 #        define psTracePrintLevels() (void)0
    41 /* do nothing */
    42 #        define psTraceSetDestination(fp) (void)0
    43 /* do nothing */
    44 #    else
     31// XXX EAM : the old 'empty' values of (void) 0 are dangerous
     32# if defined(PS_NO_TRACE)
     33    #   define psTrace(facil, level, ...)   /* do nothing */
     34    #   define p_psTrace(facil, level, ...) /* do nothing */
     35    #   define psTraceSetLevel(facil,level) /* do nothing */
     36    #   define psTraceGetLevel(facil)       /* do nothing */
     37    #   define psTraceReset()               /* do nothing */
     38    #   define psTraceFree()                /* do nothing */
     39    #   define psTracePrintLevels()         /* do nothing */
     40    #   define psTraceSetDestination(fp)    /* do nothing */
     41    #   define psTraceSetDestination()      /* do nothing */
     42    #   define PS_TRACE_ON 0
    4543
    46     /** Basic structure for the component tree.  A component is a string of the
    47         form aaa.bbb.ccc, and may itself contain further subcomponents.  The
    48         Component structure doesn't in fact contain it's full name, but only the
    49         last part. */
    50     typedef struct p_psComponent
    51     {
    52         const char *name;           // last part of name of component
    53         psS32 level;                  // trace level for this component
    54         bool p_psSpecified;
    55         psS32 n;                      // number of subcomponents
    56         struct p_psComponent* *subcomp;     // next level of subcomponents
    57     }
     44    # else /* PS_NO_TRACE */
     45        #   define PS_TRACE_ON 1
     46
     47        /** Basic structure for the component tree.  A component is a string of the
     48            form aaa.bbb.ccc, and may itself contain further subcomponents.  The
     49            Component structure doesn't in fact contain it's full name, but only the
     50            last part. */
     51        typedef struct p_psComponent
     52        {
     53            const char *name;   // last part of name of component
     54            psS32 level;   // trace level for this component
     55            bool p_psSpecified;
     56            psS32 n;    // number of subcomponents
     57            struct p_psComponent* *subcomp;     // next level of subcomponents
     58        }
    5859p_psComponent;
    5960
    … …  
    7273#ifndef SWIG
    7374#define psTrace(facil, level, ...) p_psTrace(facil, level, __VA_ARGS__)
    74 #endif
     75#endif /* SWIG */
    7576
    76 #endif
     77#endif /* DOXYGEN */
    7778
    7879/// Set trace level
    … …  
    8283
    8384/// Get the trace level
    84 psS32 psTraceGetLevel(const char *facil)     ///< facilty of interest
     85psS32 psTraceGetLevel(const char *facil) ///< facilty of interest
    8586;
    8687
    … …  
    9495void psTraceSetDestination(FILE * fp);
    9596
     97/// Get the current destination for trace messages.
     98FILE *psTraceGetDestination(void);
     99
    96100/* \} */// End of SystemGroup Functions
    97101
    98 #endif
     102#endif /* PS_NO_TRACE */
    99103
    100 #endif
     104#endif /* PS_TRACE_H */
     105

  • branches/eam-psphot-branch/psLib/src/sysUtils/psTrace.c

    r3541 r3788  
    99 *  @author GLG, MHPCC
    1010 *
    11  *  @version $Revision: 1.44 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2005-03-29 21:31:54 $
     11 *  @version $Revision: 1.44.4.1 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2005-04-29 09:17:02 $
    1313 *
    1414 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    252252    if (cRoot == NULL) {
    253253        initTrace();
     254    }
     255
     256    if (traceFP == NULL) {
     257        traceFP = stderr;
    254258    }
    255259
    … …  
    425429    psS32 i = 0;
    426430
     431    if (traceFP == NULL) {
     432        traceFP = stderr;
     433    }
     434
    427435    if (comp->name[0] == '\0') {
    428436        return;
    429437    } else {
    430438        if (comp->level == PS_DEFAULT_TRACE_LEVEL) {
    431             if (traceFP == NULL) {
    432                 printf("%*s%-*s %d\n", depth, "", 20 - depth, comp->name,
    433                        defLevel);
    434             } else {
    435                 fprintf(traceFP,"%*s%-*s %d\n", depth, "", 20 - depth, comp->name,
    436                         defLevel);
    437             }
     439            fprintf(traceFP,"%*s%-*s %d\n", depth, "", 20 - depth, comp->name, defLevel);
    438440        } else {
    439             if (traceFP == NULL) {
    440                 printf("%*s%-*s %d\n", depth, "", 20 - depth, comp->name,
    441                        comp->level);
    442             } else {
    443                 fprintf(traceFP, "%*s%-*s %d\n", depth, "", 20 - depth, comp->name,
    444                         comp->level);
    445             }
     441            fprintf(traceFP, "%*s%-*s %d\n", depth, "", 20 - depth, comp->name, comp->level);
    446442        }
    447443    }
    … …  
    498494    psS32 i = 0;
    499495
     496    if (traceFP == NULL) {
     497        traceFP = stderr;
     498    }
     499
    500500    if (NULL == comp) {
    501501        psError(PS_ERR_BAD_PARAMETER_NULL, true,
    … …  
    514514        fmt = va_arg(ap, char *);
    515515
    516         if (traceFP != NULL) {
    517             // We indent each message one space for each level of the message.
    518             for (i = 0; i < level; i++) {
    519                 fprintf(traceFP, " ");
    520             }
    521             vfprintf(traceFP, fmt, ap);
    522         } else {
    523             // We indent each message one space for each level of the message.
    524             for (i = 0; i < level; i++) {
    525                 putchar(' ');
    526             }
    527             vprintf(fmt, ap);
    528         }
     516        // We indent each message one space for each level of the message.
     517        for (i = 0; i < level; i++) {
     518            fprintf(traceFP, " ");
     519        }
     520        vfprintf(traceFP, fmt, ap);
    529521        va_end(ap);
    530522    }
    531     // NOTE: should we free *fmt as well? Read the man page.
    532 }
    533 
     523}
     524
     525// XXX EAM : I've added code to close the old traceFP (safely)
    534526void psTraceSetDestination(FILE * fp)
    535527{
     528
     529    bool special;
     530
     531    // XXX EAM perhaps return an error?
     532    if (fp == NULL) {
     533        return;
     534    }
     535
     536    // cannot close traceFP if one of the special FILE ptrs
     537    special  = (traceFP == NULL);
     538    special |= (traceFP == stdin);
     539    special |= (traceFP == stdout);
     540    special |= (traceFP == stderr);
     541
     542    if (!special) {
     543        fclose (traceFP);
     544    }
    536545    traceFP = fp;
    537546}
    538547
     548FILE *psTraceGetDestination()
     549{
     550    if (traceFP == NULL) {
     551        traceFP = stderr;
     552    }
     553    return traceFP;
     554}
     555
    539556#endif

  • branches/eam-psphot-branch/psLib/src/sysUtils/psTrace.h

    r3264 r3788  
    99 *  @author GLG, MHPCC
    1010 *
    11  *  @version $Revision: 1.31 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2005-02-17 19:26:24 $
     11 *  @version $Revision: 1.31.6.1 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2005-04-29 09:17:02 $
    1313 *
    1414 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    2929//#define PS_NO_TRACE 1   ///< to turn off all tracing
    3030
    31 #if defined(PS_NO_TRACE)
    32 #        define psTrace(facil, level, ...) (void)0
    33 /* do nothing */
    34 #        define p_psTrace(facil, level, ...)  (void)0
    35 /* do nothing */
    36 #        define psTraceSetLevel(facil,level) 0
    37 #        define psTraceGetLevel(facil) 0
    38 #        define psTraceReset() (void)0     /* do nothing */
    39 #        define psTraceFree() (void)0      /* do nothing */
    40 #        define psTracePrintLevels() (void)0
    41 /* do nothing */
    42 #        define psTraceSetDestination(fp) (void)0
    43 /* do nothing */
    44 #    else
     31// XXX EAM : the old 'empty' values of (void) 0 are dangerous
     32# if defined(PS_NO_TRACE)
     33    #   define psTrace(facil, level, ...)   /* do nothing */
     34    #   define p_psTrace(facil, level, ...) /* do nothing */
     35    #   define psTraceSetLevel(facil,level) /* do nothing */
     36    #   define psTraceGetLevel(facil)       /* do nothing */
     37    #   define psTraceReset()               /* do nothing */
     38    #   define psTraceFree()                /* do nothing */
     39    #   define psTracePrintLevels()         /* do nothing */
     40    #   define psTraceSetDestination(fp)    /* do nothing */
     41    #   define psTraceSetDestination()      /* do nothing */
     42    #   define PS_TRACE_ON 0
    4543
    46     /** Basic structure for the component tree.  A component is a string of the
    47         form aaa.bbb.ccc, and may itself contain further subcomponents.  The
    48         Component structure doesn't in fact contain it's full name, but only the
    49         last part. */
    50     typedef struct p_psComponent
    51     {
    52         const char *name;           // last part of name of component
    53         psS32 level;                  // trace level for this component
    54         bool p_psSpecified;
    55         psS32 n;                      // number of subcomponents
    56         struct p_psComponent* *subcomp;     // next level of subcomponents
    57     }
     44    # else /* PS_NO_TRACE */
     45        #   define PS_TRACE_ON 1
     46
     47        /** Basic structure for the component tree.  A component is a string of the
     48            form aaa.bbb.ccc, and may itself contain further subcomponents.  The
     49            Component structure doesn't in fact contain it's full name, but only the
     50            last part. */
     51        typedef struct p_psComponent
     52        {
     53            const char *name;   // last part of name of component
     54            psS32 level;   // trace level for this component
     55            bool p_psSpecified;
     56            psS32 n;    // number of subcomponents
     57            struct p_psComponent* *subcomp;     // next level of subcomponents
     58        }
    5859p_psComponent;
    5960
    … …  
    7273#ifndef SWIG
    7374#define psTrace(facil, level, ...) p_psTrace(facil, level, __VA_ARGS__)
    74 #endif
     75#endif /* SWIG */
    7576
    76 #endif
     77#endif /* DOXYGEN */
    7778
    7879/// Set trace level
    … …  
    8283
    8384/// Get the trace level
    84 psS32 psTraceGetLevel(const char *facil)     ///< facilty of interest
     85psS32 psTraceGetLevel(const char *facil) ///< facilty of interest
    8586;
    8687
    … …  
    9495void psTraceSetDestination(FILE * fp);
    9596
     97/// Get the current destination for trace messages.
     98FILE *psTraceGetDestination(void);
     99
    96100/* \} */// End of SystemGroup Functions
    97101
    98 #endif
     102#endif /* PS_NO_TRACE */
    99103
    100 #endif
     104#endif /* PS_TRACE_H */
     105
Note: See TracChangeset for help on using the changeset viewer.