Index: /branches/eam-psphot-branch/psLib/src/astronomy/psMetadata.c =================================================================== --- /branches/eam-psphot-branch/psLib/src/astronomy/psMetadata.c (revision 3787) +++ /branches/eam-psphot-branch/psLib/src/astronomy/psMetadata.c (revision 3788) @@ -12,6 +12,6 @@ * @author Ross Harman, MHPCC * -* @version $Revision: 1.57 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-03-29 21:31:54 $ +* @version $Revision: 1.57.4.1 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-04-29 09:17:02 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -68,5 +68,6 @@ /*****************************************************************************/ -static psMetadataItem* makeMetaMulti(psHash* table, const char* key, psMetadataItem* existing) +# if (0) + static psMetadataItem* makeMetaMulti(psHash* table, const char* key, psMetadataItem* existing) { @@ -95,4 +96,5 @@ return item; } +# endif static void metadataItemFree(psMetadataItem* metadataItem) @@ -174,5 +176,6 @@ { psMetadataItem* metadataItem = NULL; - + char tmp; + int Nbyte; PS_PTR_CHECK_NULL(name,NULL); @@ -185,11 +188,10 @@ p_psMemSetDeallocator(metadataItem, (psFreeFcn) metadataItemFree); - // Allocate and set metadata item comment - metadataItem->comment = (char *)psAlloc(sizeof(char) * MAX_STRING_LENGTH); + // set metadata item comment if (comment == NULL) { // Per SDRS, null isn't allowed, must use "" instead - strncpy(metadataItem->comment, "", MAX_STRING_LENGTH); + metadataItem->comment = psStringCopy (""); } else { - strncpy(metadataItem->comment, comment, MAX_STRING_LENGTH); + metadataItem->comment = psStringCopy (comment); } @@ -201,5 +203,6 @@ // Allocate and set metadata item name - metadataItem->name = (char *)psAlloc(sizeof(char) * MAX_STRING_LENGTH); + Nbyte = vsnprintf(&tmp, 0, name, argPtr) + 1; + metadataItem->name = (char *)psAlloc(sizeof(char) * Nbyte); vsprintf(metadataItem->name, name, argPtr); @@ -220,5 +223,5 @@ case PS_META_STR: // Perform copy of input strings - metadataItem->data.V = psStringNCopy(va_arg(argPtr, char *), MAX_STRING_LENGTH); + metadataItem->data.V = psStringCopy(va_arg(argPtr, char *)); break; case PS_META_LIST: @@ -230,7 +233,10 @@ case PS_META_ASTROM: case PS_META_UNKNOWN: - case PS_META_MULTI: // Copy of input data not performed due to variability of data types metadataItem->data.V = psMemIncrRefCounter(va_arg(argPtr, psPtr)); + break; + case PS_META_MULTI: + // MULTI needs to create a psList entry, value must be NULL + metadataItem->data.list = psListAlloc(NULL); break; default: @@ -239,5 +245,4 @@ metadataItem = NULL; } - return metadataItem; } @@ -262,4 +267,31 @@ return metadata; +} + +psBool MetadataAddNewItem (psMetadata *md, psMetadataItem *metadataItem, psS32 location) +{ + + char * key = NULL; + psHash *mdTable = NULL; + psList *mdList = NULL; + + mdTable = md->table; + mdList = md->list; + key = metadataItem->name; + + // OK, this is a new item. + // Add new metadata item to metadata collection's hash + if(!psHashAdd(mdTable, key, metadataItem)) { + psError(PS_ERR_UNKNOWN,false,PS_ERRORTEXT_psMetadata_ADD_TABLE_FAILED,key); + return false; + } + + // Add new metadata item to metadata collection's list + if(!psListAdd(mdList, location, metadataItem)) { + psError(PS_ERR_UNKNOWN,false,PS_ERRORTEXT_psMetadata_ADD_COLLECTION_FAILED,key); + return false; + } + // Success + return true; } @@ -282,72 +314,85 @@ // See if key is already in table - existingEntry = (psMetadataItem*)psHashLookup(mdTable, key); - + existingEntry = psMetadataLookup(md, key); + + if (existingEntry == NULL) { + MetadataAddNewItem (md, metadataItem, location); + return true; + } + + // if existing entry is MULTI, add the new entry to that list + // XXX does this behave correctly for adding a MULTI node to an exiting MULTI node? + if (existingEntry->type == PS_META_MULTI) { + PS_PTR_CHECK_NULL(existingEntry->data.list,NULL); + + // Add new metadata item to multi entry's list (XXX is TAIL correct, or location?) + if(!psListAdd(existingEntry->data.list, PS_LIST_TAIL, metadataItem)) { + psError(PS_ERR_UNKNOWN,false,PS_ERRORTEXT_psMetadata_ADD_COLLECTION_FAILED,key); + return false; + } + + // Add new metadata item to metadata collection's list + if(!psListAdd(mdList, location, metadataItem)) { + psError(PS_ERR_UNKNOWN,false,PS_ERRORTEXT_psMetadata_ADD_COLLECTION_FAILED,key); + return false; + } + return true; + } + + // if existing entry is NOT MULTI, we have several choices: + // - new entry is MULTI: move existing entry to new MULTI list + // - DUPLICATE_OK: create a new MULTI and add both old and new to it + // - REPLACE : replace the existing entry with the new entry + // - raise an error + + // incoming entry is PS_META_MULTI (DUPLICATE_OK is redundant, REPLACE is non-sequitor) if (metadataItem->type == PS_META_MULTI) { - // the incoming entry is PS_META_MULTI - - // force the hash entry to be PS_META_MULTI - existingEntry = makeMetaMulti(mdTable,key,existingEntry); - - // add all the items in the incoming entry to metadata - psList* list = metadataItem->data.list; - if (list != NULL) { - psListIterator* iter = psListIteratorAlloc(list,PS_LIST_HEAD,true); - psMetadataItem* listItem; - while ((listItem=(psMetadataItem*)psListGetAndIncrement(iter)) != NULL) { - psMetadataAddItem(md,listItem,location,flags); - } - psFree(iter); - } + + // replace hash entry with points to existingEntry with new entry + psHashAdd (mdTable, key, metadataItem); + + // move the existing hash entry to a node of the new entry + psListAdd (metadataItem->data.list, PS_LIST_TAIL, existingEntry); return true; // all done. } - // how the item is added to the hash depends on prior existence, flags, etc. - if(existingEntry != NULL) { // prior existence - if (existingEntry->type == PS_META_MULTI || (flags & PS_META_DUPLICATE_OK) != 0) { - // duplicate entries allowed - add another entry. - - // make sure the existing entry is PS_META_MULTI - existingEntry = makeMetaMulti(mdTable,key,existingEntry); - - // add to the hash's list of duplicate entries - if (! psListAdd(existingEntry->data.list, PS_LIST_TAIL, metadataItem) ) { - psError(PS_ERR_UNKNOWN,false,PS_ERRORTEXT_psMetadata_ADD_COLLECTION_FAILED,key); - return false; - } - } else if ((flags & PS_META_REPLACE) != 0) { - // replace entry instead of creating a duplicate entry. - - // remove the existing entry from metadata - psMetadataRemove(md,0,key); - - // treat as if new (added to list below) - if(!psHashAdd(mdTable, key, metadataItem)) { - psError(PS_ERR_UNKNOWN,false,PS_ERRORTEXT_psMetadata_ADD_TABLE_FAILED,key); - return false; - } - } else { - // default is to error on duplicate entry. - psError(PS_ERR_BAD_PARAMETER_VALUE, true, - PS_ERRORTEXT_psMetadata_DUPLICATE_NOT_ALLOWED); - return false; - } - } else { - // OK, this is a new item. - - // Node doesn't exist - Add new metadata item to metadata collection's hash - if(!psHashAdd(mdTable, key, metadataItem)) { - psError(PS_ERR_UNKNOWN,false,PS_ERRORTEXT_psMetadata_ADD_TABLE_FAILED,key); - return false; - } - } - - if(!psListAdd(mdList, location, metadataItem)) { - psError(PS_ERR_UNKNOWN,false,PS_ERRORTEXT_psMetadata_ADD_COLLECTION_FAILED,key); - return false; - } - - return true; + // duplicate entries allowed - create a MULTI with name and add both to it + if (flags & PS_META_DUPLICATE_OK) { + + psMetadataItem *multi; + + multi = psMetadataItemAlloc (key, PS_META_MULTI, "multi container", NULL); + + // replace hash entry with points to existingEntry with new multi + psHashAdd (mdTable, key, multi); + + // move the existing hash entry to a node of the new entry + psListAdd (multi->data.list, PS_LIST_TAIL, existingEntry); + + // move the existing hash entry to a node of the new entry + psListAdd (multi->data.list, PS_LIST_TAIL, metadataItem); + + // move the existing hash entry to a node of the new entry + psListAdd (mdList, PS_LIST_TAIL, metadataItem); + + return true; + } + + if (flags & PS_META_REPLACE) { + // replace entry instead of creating a duplicate entry. + + // remove the existing entry from metadata + psMetadataRemove(md,0,key); + + // treat as if new (added to list below) + MetadataAddNewItem (md, metadataItem, location); + return true; + } + + // default is to error on duplicate entry. + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psMetadata_DUPLICATE_NOT_ALLOWED); + return false; } Index: /branches/eam-psphot-branch/psLib/src/astronomy/psMetadata.h =================================================================== --- /branches/eam-psphot-branch/psLib/src/astronomy/psMetadata.h (revision 3787) +++ /branches/eam-psphot-branch/psLib/src/astronomy/psMetadata.h (revision 3788) @@ -11,6 +11,6 @@ * @author Ross Harman, MHPCC * -* @version $Revision: 1.42 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-03-29 22:13:53 $ +* @version $Revision: 1.42.4.1 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-04-29 09:17:02 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii Index: /branches/eam-psphot-branch/psLib/src/collections/psVector.c =================================================================== --- /branches/eam-psphot-branch/psLib/src/collections/psVector.c (revision 3787) +++ /branches/eam-psphot-branch/psLib/src/collections/psVector.c (revision 3788) @@ -9,6 +9,6 @@ * @author Robert DeSonia, MHPCC * -* @version $Revision: 1.36 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-03-22 21:52:49 $ +* @version $Revision: 1.36.4.1 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-04-29 09:17:02 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -540,2 +540,16 @@ } } + +// XXX EAM a utility function +bool p_psVectorPrint (FILE *f, psVector *a, char *name) +{ + + fprintf (f, "vector: %s\n", name); + + for (int i = 0; i < a[0].n; i++) { + fprintf (f, "%f\n", p_psVectorGetElementF64(a, i)); + } + fprintf (f, "\n"); + return (true); +} + Index: /branches/eam-psphot-branch/psLib/src/collections/psVector.h =================================================================== --- /branches/eam-psphot-branch/psLib/src/collections/psVector.h (revision 3787) +++ /branches/eam-psphot-branch/psLib/src/collections/psVector.h (revision 3788) @@ -11,6 +11,6 @@ * @author Ross Harman, MHPCC * - * @version $Revision: 1.30 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-03-22 21:52:49 $ + * @version $Revision: 1.30.4.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-04-29 09:17:02 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -160,4 +160,6 @@ ); +// XXX EAM a utility function +bool p_psVectorPrint (FILE *f, psVector *a, char *name); /// @} Index: /branches/eam-psphot-branch/psLib/src/dataManip/psMinimize.c =================================================================== --- /branches/eam-psphot-branch/psLib/src/dataManip/psMinimize.c (revision 3787) +++ /branches/eam-psphot-branch/psLib/src/dataManip/psMinimize.c (revision 3788) @@ -9,6 +9,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.110 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-03-31 01:02:15 $ + * @version $Revision: 1.110.6.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-04-29 09:17:02 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -557,30 +557,5 @@ } -/****************************************************************************** -psMinimizeLMChi2(): This routine will take an procedure which calculates -an arbitrary function and it's derivative and minimize the chi-squared match -between that function at the specified coords and the specified value at -those coords. - -XXX: Do this: - After checking that all entries in the paramMask are 1 or 0, when - forming the A matrix from alpha, try this: - - A[i][i] = (1 + lambda*paramask[i]) * alpha[i][i]; - -XXX: This is very different from what is specified in the SDR. Must -coordinate with IfA on new SDR. - -XXX: Do vector/image recycles. - -XXX: probably yErr will be part of the SDR. - -XXX: This must work for both F32 and F64. F32 is currently implemented. - Note: since the LUD routines are only implemented in F64, then we - will have to convert all F32 input vectors to F64 regardless. So, - the F64 port might be. - -XXX: Must update the covar matrix. - *****************************************************************************/ +// XXX EAM this is my re-implementation of MinLM psBool psMinimizeLMChi2(psMinimization *min, psImage *covar, @@ -601,4 +576,350 @@ PS_PTR_CHECK_NULL(func, NULL); + // this function has test and current values for several things + // the current best value is in lower case + // the next guess value is in upper case + + // allocate internal arrays (current vs Guess) + psImage *alpha = psImageAlloc (params->n, params->n, PS_TYPE_F64); + psImage *Alpha = psImageAlloc (params->n, params->n, PS_TYPE_F64); + psVector *beta = psVectorAlloc (params->n, PS_TYPE_F64); + psVector *Beta = psVectorAlloc (params->n, PS_TYPE_F64); + psVector *Params = psVectorAlloc (params->n, PS_TYPE_F64); + psVector *dy = NULL; + psF64 chisq = 0.0; + psF64 Chisq = 0.0; + psF64 lambda = 0.001; + + // the initial guess on params is provided by the user + Params = psVectorCopy (Params, params, PS_TYPE_F32); + + // the user provides the error or NULL. we need to convert + // to appropriate weights + dy = psVectorAlloc (y->n, PS_TYPE_F32); + if (yErr != NULL) { + for (int i = 0; i < dy->n; i++) { + dy->data.F32[i] = 1.0 / PS_SQR (yErr->data.F32[i]); + } + } else { + for (int i = 0; i < dy->n; i++) { + dy->data.F32[i] = 1.0; + } + } + + // calculate initial alpha and beta, set chisq (min->value) + min->value = p_psMinLM_SetABX (alpha, beta, params, x, y, dy, func); + # ifndef PS_NO_TRACE + // dump some useful info if trace is defined + if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") > 4) { + p_psImagePrint (psTraceGetDestination(), alpha, "alpha guess"); + p_psVectorPrint (psTraceGetDestination(), beta, "beta guess"); + p_psVectorPrint (psTraceGetDestination(), params, "params guess"); + } + # endif /* PS_NO_TRACE */ + + + // iterate until the tolerance is reached, or give up + while ((min->lastDelta > min->tol) && (min->iter < min->maxIter)) { + + // set a new guess for Alpha, Beta, Params + p_psMinLM_GuessABP (Alpha, Beta, Params, alpha, beta, params, lambda); + + # ifndef PS_NO_TRACE + // dump some useful info if trace is defined + if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") > 4) { + p_psImagePrint (psTraceGetDestination(), Alpha, "alpha guess"); + p_psVectorPrint (psTraceGetDestination(), Beta, "beta guess"); + p_psVectorPrint (psTraceGetDestination(), Params, "params guess"); + } + # endif /* PS_NO_TRACE */ + + // calculate Chisq for new guess, update Alpha & Beta + Chisq = p_psMinLM_SetABX (Alpha, Beta, Params, x, y, dy, func); + psTrace (".psLib.dataManip.psMinimizeLMChi2", 3, "chisq: %f, Chisq %f, delta: %f\n", chisq, Chisq, min->lastDelta); + + // accept new guess (if improvement), or increase lambda + if (Chisq < min->value) { + min->lastDelta = (min->value - Chisq) / (dy->n - params->n); + min->value = Chisq; + alpha = psImageCopy (alpha, Alpha, PS_TYPE_F64); + beta = psVectorCopy (beta, Beta, PS_TYPE_F64); + params = psVectorCopy (params, Params, PS_TYPE_F32); + lambda *= 0.1; + } else { + lambda *= 10.0; + } + min->iter ++; + } + psTrace (".psLib.dataManip.psMinimizeLMChi2", 3, "chisq: %f, Chisq %f, delta: %f\n", chisq, Chisq, min->lastDelta); + + // free the internal temporary data + psFree (alpha); + psFree (Alpha); + psFree (beta); + psFree (Beta); + psFree (Params); + psFree (dy); + return (true); +} + +// XXX EAM: this needs to respect the mask on params +// XXX EAM: check not NULL on alpha, beta, params +// alpha, beta, params are already allocated +psF64 p_psMinLM_SetABX (psImage *alpha, + psVector *beta, + psVector *params, + const psArray *x, + const psVector *y, + const psVector *dy, + psMinimizeLMChi2Func func) +{ + + psF64 chisq; + psF64 delta; + psF64 weight; + psF64 ymodel; + psVector *deriv = psVectorAlloc (params->n, PS_TYPE_F32); + + // zero alpha and beta for summing below + for (int j = 0; j < params->n; j++) { + for (int k = 0; k < params->n; k++) { + alpha->data.F64[j][k] = 0; + } + beta->data.F64[j] = 0; + } + chisq = 0.0; + + // calculate chisq, alpha, beta + for (int i = 0; i < y->n; i++) { + ymodel = func (deriv, params, (psVector *) x->data[i]); + + delta = ymodel - y->data.F32[i]; + chisq += PS_SQR (delta) * dy->data.F32[i]; + + for (int j = 0; j < params->n; j++) { + weight = deriv->data.F32[j] * dy->data.F32[i]; + for (int k = 0; k <= j; k++) { + alpha->data.F64[j][k] += weight * deriv->data.F32[k]; + } + beta->data.F64[j] += weight * delta; + } + } + + // calculate lower-left half of alpha + for (int j = 1; j < params->n; j++) { + for (int k = 0; k < j; k++) { + alpha->data.F64[k][j] = alpha->data.F64[j][k]; + } + } + psFree (deriv); + return (chisq); +} + +// XXX EAM : can we use static copies of LUv, LUm, A? +psBool p_psMinLM_GuessABP (psImage *Alpha, + psVector *Beta, + psVector *Params, + psImage *alpha, + psVector *beta, + psVector *params, + psF64 lambda) +{ + + # define USE_LU_DECOMP 1 + # if (USE_LU_DECOMP) + psVector *LUv = NULL; + psImage *LUm = NULL; + psImage *A = NULL; + psF32 det; + + // LU decomposition version + psTrace (".pslib.dataManip.psMinLM_GuessABP", 3, "using LUD version"); + + // set new guess values (creates matrix A) + A = psImageCopy (NULL, alpha, PS_TYPE_F64); + for (int j = 0; j < params->n; j++) { + A->data.F64[j][j] = alpha->data.F64[j][j] * (1.0 + lambda); + } + + // solve A*beta = Beta (Alpha = 1/A) + // these operations do not modify the input values (creates LUm, LUv) + LUm = psMatrixLUD (NULL, &LUv, A); + Beta = psMatrixLUSolve (Beta, LUm, beta, LUv); + Alpha = psMatrixInvert (Alpha, A, &det); + + # else + // gauss-jordan version + psTrace (".pslib.dataManip.psMinLM_GuessABP", 3, "using Gauss-J version"); + + // set new guess values (creates matrix A) + Beta = psVectorCopy (Beta, beta, PS_TYPE_F64); + Alpha = psImageCopy (Alpha, alpha, PS_TYPE_F64); + for (int j = 0; j < params->n; j++) { + Alpha->data.F64[j][j] = alpha->data.F64[j][j] * (1.0 + lambda); + } + + psGaussJordan (Alpha, Beta); + # endif + + // apply beta to get new params values + for (int j = 0; j < params->n; j++) { + Params->data.F32[j] = params->data.F32[j] - Beta->data.F64[j]; + } + + # if (USE_LU_DECOMP) + psFree (A); + psFree (LUm); + psFree (LUv); + # endif + + return true; +} + +# define SWAP(X,Y) {double tmp=(X); (X) = (Y); (Y) = tmp;} + +// XXX EAM : temporary gauss-jordan solver based on gene's +// version based on the Numerical Recipes version +bool psGaussJordan (psImage *a, psVector *b) +{ + + int *indxc,*indxr,*ipiv; + int Nx, icol, irow; + int i, j, k, l, ll; + float big, dum, pivinv; + psF64 *vector; + psF64 **matrix; + + Nx = a->numCols; + matrix = a->data.F64; + vector = b->data.F64; + + indxc = psAlloc (Nx*sizeof(int)); + indxr = psAlloc (Nx*sizeof(int)); + ipiv = psAlloc (Nx*sizeof(int)); + for (j = 0; j < Nx; j++) + ipiv[j] = 0; + + irow = icol = 0; + big = fabs(matrix[0][0]); + + for (i = 0; i < Nx; i++) { + big = 0.0; + for (j = 0; j < Nx; j++) { + if (!finite(matrix[i][j])) { + // XXX EAM: this should use the psError stack + fprintf (stderr, "GAUSSJ: NaN\n"); + goto fescape; + } + if (ipiv[j] != 1) { + for (k = 0; k < Nx; k++) { + if (ipiv[k] == 0) { + if (fabs (matrix[j][k]) >= big) { + big = fabs (matrix[j][k]); + irow = j; + icol = k; + } + } else { + if (ipiv[k] > 1) { + // XXX EAM: this should use the psError stack + fprintf (stderr, "GAUSSJ: Singular Matrix! (1)\n"); + goto fescape; + } + } + } + } + } + ipiv[icol]++; + if (irow != icol) { + for (l = 0; l < Nx; l++) { + SWAP (matrix[irow][l], matrix[icol][l]); + } + SWAP (vector[irow], vector[icol]); + } + indxr[i] = irow; + indxc[i] = icol; + if (matrix[icol][icol] == 0.0) { + // XXX EAM: this should use the psError stack + fprintf (stderr, "GAUSSJ: Singular Matrix! (2)\n"); + goto fescape; + } + pivinv = 1.0 / matrix[icol][icol]; + matrix[icol][icol] = 1.0; + for (l = 0; l < Nx; l++) { + matrix[icol][l] *= pivinv; + } + vector[icol] *= pivinv; + + for (ll = 0; ll < Nx; ll++) { + if (ll != icol) { + dum = matrix[ll][icol]; + matrix[ll][icol] = 0.0; + for (l = 0; l < Nx; l++) + matrix[ll][l] -= matrix[icol][l]*dum; + vector[ll] -= vector[icol]*dum; + } + } + } + + for (l = Nx - 1; l >= 0; l--) { + if (indxr[l] != indxc[l]) + for (k = 0; k < Nx; k++) + SWAP (matrix[k][indxr[l]], matrix[k][indxc[l]]); + } + psFree (ipiv); + psFree (indxr); + psFree (indxc); + return (true); + +fescape: + psFree (ipiv); + psFree (indxr); + psFree (indxc); + return (false); +} + +/****************************************************************************** +psMinimizeLMChi2(): This routine will take an procedure which calculates +an arbitrary function and it's derivative and minimize the chi-squared match +between that function at the specified coords and the specified value at +those coords. + +XXX: Do this: + After checking that all entries in the paramMask are 1 or 0, when + forming the A matrix from alpha, try this: + + A[i][i] = (1 + lambda*paramask[i]) * alpha[i][i]; + +XXX: This is very different from what is specified in the SDR. Must +coordinate with IfA on new SDR. + +XXX: Do vector/image recycles. + +XXX: probably yErr will be part of the SDR. + +XXX: This must work for both F32 and F64. F32 is currently implemented. + Note: since the LUD routines are only implemented in F64, then we + will have to convert all F32 input vectors to F64 regardless. So, + the F64 port might be. + +XXX: Must update the covar matrix. + *****************************************************************************/ +psBool psMinimizeLMChi2Old(psMinimization *min, + psImage *covar, + psVector *params, + const psVector *paramMask, + const psArray *x, + const psVector *y, + const psVector *yErr, + psMinimizeLMChi2Func func) +{ + PS_PTR_CHECK_NULL(min, NULL); + PS_VECTOR_CHECK_NULL(params, NULL); + PS_VECTOR_CHECK_EMPTY(params, NULL); + PS_PTR_CHECK_NULL(x, NULL); + PS_VECTOR_CHECK_NULL(y, NULL); + PS_VECTOR_CHECK_EMPTY(y, NULL); + PS_VECTOR_CHECK_SIZE_EQUAL(x, y, NULL); + PS_PTR_CHECK_NULL(func, NULL); + if (paramMask != NULL) { PS_VECTOR_CHECK_SIZE_EQUAL(params, paramMask, NULL); @@ -636,6 +957,6 @@ psF32 currChi2 = 0.0; psF32 newChi2 = 0.0; - psF32 lamda = 0.00005; - lamda = 0.05; + psF32 lamda = 0.00005; // XXX EAM : this starting value is VERY small (lamda is mis-spelt) + lamda = 0.05; // XXX EAM : this starting value is quite large (lamda is mis-spelt) psTrace(".psLib.dataManip.psMinimize", 6, @@ -661,5 +982,11 @@ // currChi2 = 0.0; - currValueVec = func(deriv, params, x); + // currValueVec = func(deriv, params, x); + + // XXX EAM: use BinaryOp ? + // t1 = BinaryOp (NULL, currValueVec, "-", y); + // t1 = BinaryOp (t1, t1, "*", t1); + + // XXX EAM: this ignores yErr for (n=0;ndata.F32[n] - y->data.F32[n]) * @@ -670,4 +997,5 @@ } + // XXX EAM: this is just for tracing for (p=0;pdata.U8[p] != 0)) { @@ -690,6 +1018,9 @@ // // Calculate the BETA vector. - // + // XXX EAM: I think this is wrong for (p=0;pdata.U8[p] != 0)) { + continue; + } beta->data.F64[p] = 0.0; for (n=0;ndata.F32[n] - y->data.F32[n]) * @@ -1098,5 +1429,5 @@ min->value = 0.0; min->iter = 0; - min->lastDelta = 0.0; + min->lastDelta = tol + 1; return(min); Index: /branches/eam-psphot-branch/psLib/src/dataManip/psMinimize.h =================================================================== --- /branches/eam-psphot-branch/psLib/src/dataManip/psMinimize.h (revision 3787) +++ /branches/eam-psphot-branch/psLib/src/dataManip/psMinimize.h (revision 3788) @@ -8,6 +8,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.39 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-02-17 19:26:23 $ + * @version $Revision: 1.39.6.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-04-29 09:17:02 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -42,5 +42,5 @@ psF32 tol; ///< Error Tolerance psF32 value; ///< Value of function at minimum - psS32 iter; ///< Number of iterations required + psS32 iter; ///< Number of iterations to date psF32 lastDelta; ///< The last difference for the fit } @@ -77,7 +77,7 @@ typedef -psVector* (*psMinimizeLMChi2Func)(psImage *deriv, - const psVector *params, - const psArray *x); +psF64 (*psMinimizeLMChi2Func)(psVector *deriv, + psVector *params, + psVector *x); psBool psMinimizeLMChi2(psMinimization *min, @@ -88,4 +88,20 @@ const psVector *y, const psVector *yErr, + psMinimizeLMChi2Func func); + +psBool p_psMinLM_GuessABP (psImage *Alpha, + psVector *Beta, + psVector *Params, + psImage *alpha, + psVector *beta, + psVector *params, + psF64 lambda); + +psF64 p_psMinLM_SetABX (psImage *alpha, + psVector *beta, + psVector *params, + const psArray *x, + const psVector *y, + const psVector *dy, psMinimizeLMChi2Func func); @@ -119,5 +135,6 @@ psMinimizeChi2PowellFunc func); - +// XXX EAM : psGaussJordan provided as an alternate to LU Decomp for psMinimizeLMChi2 +bool psGaussJordan (psImage *a, psVector *b); /* \} */// End of MathGroup Functions Index: /branches/eam-psphot-branch/psLib/src/image/psImage.c =================================================================== --- /branches/eam-psphot-branch/psLib/src/image/psImage.c (revision 3787) +++ /branches/eam-psphot-branch/psLib/src/image/psImage.c (revision 3788) @@ -9,6 +9,6 @@ * @author Ross Harman, MHPCC * - * @version $Revision: 1.61 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-03-18 02:35:14 $ + * @version $Revision: 1.61.4.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-04-29 09:17:02 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -412,4 +412,69 @@ return numFreed; +} + +// XXX EAM some utility functions +bool p_psImagePrint (FILE *f, psImage *a, char *name) +{ + + fprintf (f, "matrix: %s\n", name); + + for (int j = 0; j < a[0].numRows; j++) { + for (int i = 0; i < a[0].numCols; i++) { + fprintf (f, "%f ", p_psImageGetElementF64(a, i, j)); + } + fprintf (f, "\n"); + } + fprintf (f, "\n"); + return (true); +} + +psF64 p_psImageGetElementF64(psImage* image, + int col, + int row) +{ + if (image == NULL) { + return NAN; + } + if (col < 0 || col >= image->numCols) { + return NAN; + } + if (row < 0 || row >= image->numRows) { + return NAN; + } + + switch (image->type.type) { + case PS_TYPE_U8: + return image->data.U8[row][col]; + break; + case PS_TYPE_U16: + return image->data.U16[row][col]; + break; + case PS_TYPE_U32: + return image->data.U32[row][col]; + break; + case PS_TYPE_U64: + return image->data.U64[row][col]; + break; + case PS_TYPE_S8: + return image->data.S8[row][col]; + break; + case PS_TYPE_S16: + return image->data.S16[row][col]; + break; + case PS_TYPE_S32: + return image->data.S32[row][col]; + break; + case PS_TYPE_S64: + return image->data.S64[row][col]; + break; + case PS_TYPE_F32: + return image->data.F32[row][col]; + break; + case PS_TYPE_F64: + return image->data.F64[row][col]; + default: + return NAN; + } } Index: /branches/eam-psphot-branch/psLib/src/image/psImage.h =================================================================== --- /branches/eam-psphot-branch/psLib/src/image/psImage.h (revision 3787) +++ /branches/eam-psphot-branch/psLib/src/image/psImage.h (revision 3788) @@ -11,6 +11,6 @@ * @author Ross Harman, MHPCC * - * @version $Revision: 1.49 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-03-18 02:35:14 $ + * @version $Revision: 1.49.4.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-04-29 09:17:02 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -20,4 +20,5 @@ #include +#include #include "psType.h" @@ -191,4 +192,10 @@ ); + +// XXX EAM some utility functions + +psF64 p_psImageGetElementF64(psImage *a, int i, int j); +bool p_psImagePrint (FILE *f, psImage *a, char *name); + #define PIXEL_INTERPOLATE_FCN_PROTOTYPE(SUFFIX, RETURNTYPE) \ inline RETURNTYPE p_psImagePixelInterpolate##SUFFIX( \ Index: /branches/eam-psphot-branch/psLib/src/math/psMinimize.c =================================================================== --- /branches/eam-psphot-branch/psLib/src/math/psMinimize.c (revision 3787) +++ /branches/eam-psphot-branch/psLib/src/math/psMinimize.c (revision 3788) @@ -9,6 +9,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.110 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-03-31 01:02:15 $ + * @version $Revision: 1.110.6.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-04-29 09:17:02 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -557,30 +557,5 @@ } -/****************************************************************************** -psMinimizeLMChi2(): This routine will take an procedure which calculates -an arbitrary function and it's derivative and minimize the chi-squared match -between that function at the specified coords and the specified value at -those coords. - -XXX: Do this: - After checking that all entries in the paramMask are 1 or 0, when - forming the A matrix from alpha, try this: - - A[i][i] = (1 + lambda*paramask[i]) * alpha[i][i]; - -XXX: This is very different from what is specified in the SDR. Must -coordinate with IfA on new SDR. - -XXX: Do vector/image recycles. - -XXX: probably yErr will be part of the SDR. - -XXX: This must work for both F32 and F64. F32 is currently implemented. - Note: since the LUD routines are only implemented in F64, then we - will have to convert all F32 input vectors to F64 regardless. So, - the F64 port might be. - -XXX: Must update the covar matrix. - *****************************************************************************/ +// XXX EAM this is my re-implementation of MinLM psBool psMinimizeLMChi2(psMinimization *min, psImage *covar, @@ -601,4 +576,350 @@ PS_PTR_CHECK_NULL(func, NULL); + // this function has test and current values for several things + // the current best value is in lower case + // the next guess value is in upper case + + // allocate internal arrays (current vs Guess) + psImage *alpha = psImageAlloc (params->n, params->n, PS_TYPE_F64); + psImage *Alpha = psImageAlloc (params->n, params->n, PS_TYPE_F64); + psVector *beta = psVectorAlloc (params->n, PS_TYPE_F64); + psVector *Beta = psVectorAlloc (params->n, PS_TYPE_F64); + psVector *Params = psVectorAlloc (params->n, PS_TYPE_F64); + psVector *dy = NULL; + psF64 chisq = 0.0; + psF64 Chisq = 0.0; + psF64 lambda = 0.001; + + // the initial guess on params is provided by the user + Params = psVectorCopy (Params, params, PS_TYPE_F32); + + // the user provides the error or NULL. we need to convert + // to appropriate weights + dy = psVectorAlloc (y->n, PS_TYPE_F32); + if (yErr != NULL) { + for (int i = 0; i < dy->n; i++) { + dy->data.F32[i] = 1.0 / PS_SQR (yErr->data.F32[i]); + } + } else { + for (int i = 0; i < dy->n; i++) { + dy->data.F32[i] = 1.0; + } + } + + // calculate initial alpha and beta, set chisq (min->value) + min->value = p_psMinLM_SetABX (alpha, beta, params, x, y, dy, func); + # ifndef PS_NO_TRACE + // dump some useful info if trace is defined + if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") > 4) { + p_psImagePrint (psTraceGetDestination(), alpha, "alpha guess"); + p_psVectorPrint (psTraceGetDestination(), beta, "beta guess"); + p_psVectorPrint (psTraceGetDestination(), params, "params guess"); + } + # endif /* PS_NO_TRACE */ + + + // iterate until the tolerance is reached, or give up + while ((min->lastDelta > min->tol) && (min->iter < min->maxIter)) { + + // set a new guess for Alpha, Beta, Params + p_psMinLM_GuessABP (Alpha, Beta, Params, alpha, beta, params, lambda); + + # ifndef PS_NO_TRACE + // dump some useful info if trace is defined + if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") > 4) { + p_psImagePrint (psTraceGetDestination(), Alpha, "alpha guess"); + p_psVectorPrint (psTraceGetDestination(), Beta, "beta guess"); + p_psVectorPrint (psTraceGetDestination(), Params, "params guess"); + } + # endif /* PS_NO_TRACE */ + + // calculate Chisq for new guess, update Alpha & Beta + Chisq = p_psMinLM_SetABX (Alpha, Beta, Params, x, y, dy, func); + psTrace (".psLib.dataManip.psMinimizeLMChi2", 3, "chisq: %f, Chisq %f, delta: %f\n", chisq, Chisq, min->lastDelta); + + // accept new guess (if improvement), or increase lambda + if (Chisq < min->value) { + min->lastDelta = (min->value - Chisq) / (dy->n - params->n); + min->value = Chisq; + alpha = psImageCopy (alpha, Alpha, PS_TYPE_F64); + beta = psVectorCopy (beta, Beta, PS_TYPE_F64); + params = psVectorCopy (params, Params, PS_TYPE_F32); + lambda *= 0.1; + } else { + lambda *= 10.0; + } + min->iter ++; + } + psTrace (".psLib.dataManip.psMinimizeLMChi2", 3, "chisq: %f, Chisq %f, delta: %f\n", chisq, Chisq, min->lastDelta); + + // free the internal temporary data + psFree (alpha); + psFree (Alpha); + psFree (beta); + psFree (Beta); + psFree (Params); + psFree (dy); + return (true); +} + +// XXX EAM: this needs to respect the mask on params +// XXX EAM: check not NULL on alpha, beta, params +// alpha, beta, params are already allocated +psF64 p_psMinLM_SetABX (psImage *alpha, + psVector *beta, + psVector *params, + const psArray *x, + const psVector *y, + const psVector *dy, + psMinimizeLMChi2Func func) +{ + + psF64 chisq; + psF64 delta; + psF64 weight; + psF64 ymodel; + psVector *deriv = psVectorAlloc (params->n, PS_TYPE_F32); + + // zero alpha and beta for summing below + for (int j = 0; j < params->n; j++) { + for (int k = 0; k < params->n; k++) { + alpha->data.F64[j][k] = 0; + } + beta->data.F64[j] = 0; + } + chisq = 0.0; + + // calculate chisq, alpha, beta + for (int i = 0; i < y->n; i++) { + ymodel = func (deriv, params, (psVector *) x->data[i]); + + delta = ymodel - y->data.F32[i]; + chisq += PS_SQR (delta) * dy->data.F32[i]; + + for (int j = 0; j < params->n; j++) { + weight = deriv->data.F32[j] * dy->data.F32[i]; + for (int k = 0; k <= j; k++) { + alpha->data.F64[j][k] += weight * deriv->data.F32[k]; + } + beta->data.F64[j] += weight * delta; + } + } + + // calculate lower-left half of alpha + for (int j = 1; j < params->n; j++) { + for (int k = 0; k < j; k++) { + alpha->data.F64[k][j] = alpha->data.F64[j][k]; + } + } + psFree (deriv); + return (chisq); +} + +// XXX EAM : can we use static copies of LUv, LUm, A? +psBool p_psMinLM_GuessABP (psImage *Alpha, + psVector *Beta, + psVector *Params, + psImage *alpha, + psVector *beta, + psVector *params, + psF64 lambda) +{ + + # define USE_LU_DECOMP 1 + # if (USE_LU_DECOMP) + psVector *LUv = NULL; + psImage *LUm = NULL; + psImage *A = NULL; + psF32 det; + + // LU decomposition version + psTrace (".pslib.dataManip.psMinLM_GuessABP", 3, "using LUD version"); + + // set new guess values (creates matrix A) + A = psImageCopy (NULL, alpha, PS_TYPE_F64); + for (int j = 0; j < params->n; j++) { + A->data.F64[j][j] = alpha->data.F64[j][j] * (1.0 + lambda); + } + + // solve A*beta = Beta (Alpha = 1/A) + // these operations do not modify the input values (creates LUm, LUv) + LUm = psMatrixLUD (NULL, &LUv, A); + Beta = psMatrixLUSolve (Beta, LUm, beta, LUv); + Alpha = psMatrixInvert (Alpha, A, &det); + + # else + // gauss-jordan version + psTrace (".pslib.dataManip.psMinLM_GuessABP", 3, "using Gauss-J version"); + + // set new guess values (creates matrix A) + Beta = psVectorCopy (Beta, beta, PS_TYPE_F64); + Alpha = psImageCopy (Alpha, alpha, PS_TYPE_F64); + for (int j = 0; j < params->n; j++) { + Alpha->data.F64[j][j] = alpha->data.F64[j][j] * (1.0 + lambda); + } + + psGaussJordan (Alpha, Beta); + # endif + + // apply beta to get new params values + for (int j = 0; j < params->n; j++) { + Params->data.F32[j] = params->data.F32[j] - Beta->data.F64[j]; + } + + # if (USE_LU_DECOMP) + psFree (A); + psFree (LUm); + psFree (LUv); + # endif + + return true; +} + +# define SWAP(X,Y) {double tmp=(X); (X) = (Y); (Y) = tmp;} + +// XXX EAM : temporary gauss-jordan solver based on gene's +// version based on the Numerical Recipes version +bool psGaussJordan (psImage *a, psVector *b) +{ + + int *indxc,*indxr,*ipiv; + int Nx, icol, irow; + int i, j, k, l, ll; + float big, dum, pivinv; + psF64 *vector; + psF64 **matrix; + + Nx = a->numCols; + matrix = a->data.F64; + vector = b->data.F64; + + indxc = psAlloc (Nx*sizeof(int)); + indxr = psAlloc (Nx*sizeof(int)); + ipiv = psAlloc (Nx*sizeof(int)); + for (j = 0; j < Nx; j++) + ipiv[j] = 0; + + irow = icol = 0; + big = fabs(matrix[0][0]); + + for (i = 0; i < Nx; i++) { + big = 0.0; + for (j = 0; j < Nx; j++) { + if (!finite(matrix[i][j])) { + // XXX EAM: this should use the psError stack + fprintf (stderr, "GAUSSJ: NaN\n"); + goto fescape; + } + if (ipiv[j] != 1) { + for (k = 0; k < Nx; k++) { + if (ipiv[k] == 0) { + if (fabs (matrix[j][k]) >= big) { + big = fabs (matrix[j][k]); + irow = j; + icol = k; + } + } else { + if (ipiv[k] > 1) { + // XXX EAM: this should use the psError stack + fprintf (stderr, "GAUSSJ: Singular Matrix! (1)\n"); + goto fescape; + } + } + } + } + } + ipiv[icol]++; + if (irow != icol) { + for (l = 0; l < Nx; l++) { + SWAP (matrix[irow][l], matrix[icol][l]); + } + SWAP (vector[irow], vector[icol]); + } + indxr[i] = irow; + indxc[i] = icol; + if (matrix[icol][icol] == 0.0) { + // XXX EAM: this should use the psError stack + fprintf (stderr, "GAUSSJ: Singular Matrix! (2)\n"); + goto fescape; + } + pivinv = 1.0 / matrix[icol][icol]; + matrix[icol][icol] = 1.0; + for (l = 0; l < Nx; l++) { + matrix[icol][l] *= pivinv; + } + vector[icol] *= pivinv; + + for (ll = 0; ll < Nx; ll++) { + if (ll != icol) { + dum = matrix[ll][icol]; + matrix[ll][icol] = 0.0; + for (l = 0; l < Nx; l++) + matrix[ll][l] -= matrix[icol][l]*dum; + vector[ll] -= vector[icol]*dum; + } + } + } + + for (l = Nx - 1; l >= 0; l--) { + if (indxr[l] != indxc[l]) + for (k = 0; k < Nx; k++) + SWAP (matrix[k][indxr[l]], matrix[k][indxc[l]]); + } + psFree (ipiv); + psFree (indxr); + psFree (indxc); + return (true); + +fescape: + psFree (ipiv); + psFree (indxr); + psFree (indxc); + return (false); +} + +/****************************************************************************** +psMinimizeLMChi2(): This routine will take an procedure which calculates +an arbitrary function and it's derivative and minimize the chi-squared match +between that function at the specified coords and the specified value at +those coords. + +XXX: Do this: + After checking that all entries in the paramMask are 1 or 0, when + forming the A matrix from alpha, try this: + + A[i][i] = (1 + lambda*paramask[i]) * alpha[i][i]; + +XXX: This is very different from what is specified in the SDR. Must +coordinate with IfA on new SDR. + +XXX: Do vector/image recycles. + +XXX: probably yErr will be part of the SDR. + +XXX: This must work for both F32 and F64. F32 is currently implemented. + Note: since the LUD routines are only implemented in F64, then we + will have to convert all F32 input vectors to F64 regardless. So, + the F64 port might be. + +XXX: Must update the covar matrix. + *****************************************************************************/ +psBool psMinimizeLMChi2Old(psMinimization *min, + psImage *covar, + psVector *params, + const psVector *paramMask, + const psArray *x, + const psVector *y, + const psVector *yErr, + psMinimizeLMChi2Func func) +{ + PS_PTR_CHECK_NULL(min, NULL); + PS_VECTOR_CHECK_NULL(params, NULL); + PS_VECTOR_CHECK_EMPTY(params, NULL); + PS_PTR_CHECK_NULL(x, NULL); + PS_VECTOR_CHECK_NULL(y, NULL); + PS_VECTOR_CHECK_EMPTY(y, NULL); + PS_VECTOR_CHECK_SIZE_EQUAL(x, y, NULL); + PS_PTR_CHECK_NULL(func, NULL); + if (paramMask != NULL) { PS_VECTOR_CHECK_SIZE_EQUAL(params, paramMask, NULL); @@ -636,6 +957,6 @@ psF32 currChi2 = 0.0; psF32 newChi2 = 0.0; - psF32 lamda = 0.00005; - lamda = 0.05; + psF32 lamda = 0.00005; // XXX EAM : this starting value is VERY small (lamda is mis-spelt) + lamda = 0.05; // XXX EAM : this starting value is quite large (lamda is mis-spelt) psTrace(".psLib.dataManip.psMinimize", 6, @@ -661,5 +982,11 @@ // currChi2 = 0.0; - currValueVec = func(deriv, params, x); + // currValueVec = func(deriv, params, x); + + // XXX EAM: use BinaryOp ? + // t1 = BinaryOp (NULL, currValueVec, "-", y); + // t1 = BinaryOp (t1, t1, "*", t1); + + // XXX EAM: this ignores yErr for (n=0;ndata.F32[n] - y->data.F32[n]) * @@ -670,4 +997,5 @@ } + // XXX EAM: this is just for tracing for (p=0;pdata.U8[p] != 0)) { @@ -690,6 +1018,9 @@ // // Calculate the BETA vector. - // + // XXX EAM: I think this is wrong for (p=0;pdata.U8[p] != 0)) { + continue; + } beta->data.F64[p] = 0.0; for (n=0;ndata.F32[n] - y->data.F32[n]) * @@ -1098,5 +1429,5 @@ min->value = 0.0; min->iter = 0; - min->lastDelta = 0.0; + min->lastDelta = tol + 1; return(min); Index: /branches/eam-psphot-branch/psLib/src/math/psMinimize.h =================================================================== --- /branches/eam-psphot-branch/psLib/src/math/psMinimize.h (revision 3787) +++ /branches/eam-psphot-branch/psLib/src/math/psMinimize.h (revision 3788) @@ -8,6 +8,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.39 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-02-17 19:26:23 $ + * @version $Revision: 1.39.6.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-04-29 09:17:02 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -42,5 +42,5 @@ psF32 tol; ///< Error Tolerance psF32 value; ///< Value of function at minimum - psS32 iter; ///< Number of iterations required + psS32 iter; ///< Number of iterations to date psF32 lastDelta; ///< The last difference for the fit } @@ -77,7 +77,7 @@ typedef -psVector* (*psMinimizeLMChi2Func)(psImage *deriv, - const psVector *params, - const psArray *x); +psF64 (*psMinimizeLMChi2Func)(psVector *deriv, + psVector *params, + psVector *x); psBool psMinimizeLMChi2(psMinimization *min, @@ -88,4 +88,20 @@ const psVector *y, const psVector *yErr, + psMinimizeLMChi2Func func); + +psBool p_psMinLM_GuessABP (psImage *Alpha, + psVector *Beta, + psVector *Params, + psImage *alpha, + psVector *beta, + psVector *params, + psF64 lambda); + +psF64 p_psMinLM_SetABX (psImage *alpha, + psVector *beta, + psVector *params, + const psArray *x, + const psVector *y, + const psVector *dy, psMinimizeLMChi2Func func); @@ -119,5 +135,6 @@ psMinimizeChi2PowellFunc func); - +// XXX EAM : psGaussJordan provided as an alternate to LU Decomp for psMinimizeLMChi2 +bool psGaussJordan (psImage *a, psVector *b); /* \} */// End of MathGroup Functions Index: /branches/eam-psphot-branch/psLib/src/mathtypes/psImage.c =================================================================== --- /branches/eam-psphot-branch/psLib/src/mathtypes/psImage.c (revision 3787) +++ /branches/eam-psphot-branch/psLib/src/mathtypes/psImage.c (revision 3788) @@ -9,6 +9,6 @@ * @author Ross Harman, MHPCC * - * @version $Revision: 1.61 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-03-18 02:35:14 $ + * @version $Revision: 1.61.4.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-04-29 09:17:02 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -412,4 +412,69 @@ return numFreed; +} + +// XXX EAM some utility functions +bool p_psImagePrint (FILE *f, psImage *a, char *name) +{ + + fprintf (f, "matrix: %s\n", name); + + for (int j = 0; j < a[0].numRows; j++) { + for (int i = 0; i < a[0].numCols; i++) { + fprintf (f, "%f ", p_psImageGetElementF64(a, i, j)); + } + fprintf (f, "\n"); + } + fprintf (f, "\n"); + return (true); +} + +psF64 p_psImageGetElementF64(psImage* image, + int col, + int row) +{ + if (image == NULL) { + return NAN; + } + if (col < 0 || col >= image->numCols) { + return NAN; + } + if (row < 0 || row >= image->numRows) { + return NAN; + } + + switch (image->type.type) { + case PS_TYPE_U8: + return image->data.U8[row][col]; + break; + case PS_TYPE_U16: + return image->data.U16[row][col]; + break; + case PS_TYPE_U32: + return image->data.U32[row][col]; + break; + case PS_TYPE_U64: + return image->data.U64[row][col]; + break; + case PS_TYPE_S8: + return image->data.S8[row][col]; + break; + case PS_TYPE_S16: + return image->data.S16[row][col]; + break; + case PS_TYPE_S32: + return image->data.S32[row][col]; + break; + case PS_TYPE_S64: + return image->data.S64[row][col]; + break; + case PS_TYPE_F32: + return image->data.F32[row][col]; + break; + case PS_TYPE_F64: + return image->data.F64[row][col]; + default: + return NAN; + } } Index: /branches/eam-psphot-branch/psLib/src/mathtypes/psImage.h =================================================================== --- /branches/eam-psphot-branch/psLib/src/mathtypes/psImage.h (revision 3787) +++ /branches/eam-psphot-branch/psLib/src/mathtypes/psImage.h (revision 3788) @@ -11,6 +11,6 @@ * @author Ross Harman, MHPCC * - * @version $Revision: 1.49 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-03-18 02:35:14 $ + * @version $Revision: 1.49.4.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-04-29 09:17:02 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -20,4 +20,5 @@ #include +#include #include "psType.h" @@ -191,4 +192,10 @@ ); + +// XXX EAM some utility functions + +psF64 p_psImageGetElementF64(psImage *a, int i, int j); +bool p_psImagePrint (FILE *f, psImage *a, char *name); + #define PIXEL_INTERPOLATE_FCN_PROTOTYPE(SUFFIX, RETURNTYPE) \ inline RETURNTYPE p_psImagePixelInterpolate##SUFFIX( \ Index: /branches/eam-psphot-branch/psLib/src/mathtypes/psVector.c =================================================================== --- /branches/eam-psphot-branch/psLib/src/mathtypes/psVector.c (revision 3787) +++ /branches/eam-psphot-branch/psLib/src/mathtypes/psVector.c (revision 3788) @@ -9,6 +9,6 @@ * @author Robert DeSonia, MHPCC * -* @version $Revision: 1.36 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-03-22 21:52:49 $ +* @version $Revision: 1.36.4.1 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-04-29 09:17:02 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -540,2 +540,16 @@ } } + +// XXX EAM a utility function +bool p_psVectorPrint (FILE *f, psVector *a, char *name) +{ + + fprintf (f, "vector: %s\n", name); + + for (int i = 0; i < a[0].n; i++) { + fprintf (f, "%f\n", p_psVectorGetElementF64(a, i)); + } + fprintf (f, "\n"); + return (true); +} + Index: /branches/eam-psphot-branch/psLib/src/mathtypes/psVector.h =================================================================== --- /branches/eam-psphot-branch/psLib/src/mathtypes/psVector.h (revision 3787) +++ /branches/eam-psphot-branch/psLib/src/mathtypes/psVector.h (revision 3788) @@ -11,6 +11,6 @@ * @author Ross Harman, MHPCC * - * @version $Revision: 1.30 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-03-22 21:52:49 $ + * @version $Revision: 1.30.4.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-04-29 09:17:02 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -160,4 +160,6 @@ ); +// XXX EAM a utility function +bool p_psVectorPrint (FILE *f, psVector *a, char *name); /// @} Index: /branches/eam-psphot-branch/psLib/src/sys/psTrace.c =================================================================== --- /branches/eam-psphot-branch/psLib/src/sys/psTrace.c (revision 3787) +++ /branches/eam-psphot-branch/psLib/src/sys/psTrace.c (revision 3788) @@ -9,6 +9,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.44 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-03-29 21:31:54 $ + * @version $Revision: 1.44.4.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-04-29 09:17:02 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -252,4 +252,8 @@ if (cRoot == NULL) { initTrace(); + } + + if (traceFP == NULL) { + traceFP = stderr; } @@ -425,23 +429,15 @@ psS32 i = 0; + if (traceFP == NULL) { + traceFP = stderr; + } + if (comp->name[0] == '\0') { return; } else { if (comp->level == PS_DEFAULT_TRACE_LEVEL) { - if (traceFP == NULL) { - printf("%*s%-*s %d\n", depth, "", 20 - depth, comp->name, - defLevel); - } else { - fprintf(traceFP,"%*s%-*s %d\n", depth, "", 20 - depth, comp->name, - defLevel); - } + fprintf(traceFP,"%*s%-*s %d\n", depth, "", 20 - depth, comp->name, defLevel); } else { - if (traceFP == NULL) { - printf("%*s%-*s %d\n", depth, "", 20 - depth, comp->name, - comp->level); - } else { - fprintf(traceFP, "%*s%-*s %d\n", depth, "", 20 - depth, comp->name, - comp->level); - } + fprintf(traceFP, "%*s%-*s %d\n", depth, "", 20 - depth, comp->name, comp->level); } } @@ -498,4 +494,8 @@ psS32 i = 0; + if (traceFP == NULL) { + traceFP = stderr; + } + if (NULL == comp) { psError(PS_ERR_BAD_PARAMETER_NULL, true, @@ -514,26 +514,43 @@ fmt = va_arg(ap, char *); - if (traceFP != NULL) { - // We indent each message one space for each level of the message. - for (i = 0; i < level; i++) { - fprintf(traceFP, " "); - } - vfprintf(traceFP, fmt, ap); - } else { - // We indent each message one space for each level of the message. - for (i = 0; i < level; i++) { - putchar(' '); - } - vprintf(fmt, ap); - } + // We indent each message one space for each level of the message. + for (i = 0; i < level; i++) { + fprintf(traceFP, " "); + } + vfprintf(traceFP, fmt, ap); va_end(ap); } - // NOTE: should we free *fmt as well? Read the man page. -} - +} + +// XXX EAM : I've added code to close the old traceFP (safely) void psTraceSetDestination(FILE * fp) { + + bool special; + + // XXX EAM perhaps return an error? + if (fp == NULL) { + return; + } + + // cannot close traceFP if one of the special FILE ptrs + special = (traceFP == NULL); + special |= (traceFP == stdin); + special |= (traceFP == stdout); + special |= (traceFP == stderr); + + if (!special) { + fclose (traceFP); + } traceFP = fp; } +FILE *psTraceGetDestination() +{ + if (traceFP == NULL) { + traceFP = stderr; + } + return traceFP; +} + #endif Index: /branches/eam-psphot-branch/psLib/src/sys/psTrace.h =================================================================== --- /branches/eam-psphot-branch/psLib/src/sys/psTrace.h (revision 3787) +++ /branches/eam-psphot-branch/psLib/src/sys/psTrace.h (revision 3788) @@ -9,6 +9,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.31 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-02-17 19:26:24 $ + * @version $Revision: 1.31.6.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-04-29 09:17:02 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -29,31 +29,32 @@ //#define PS_NO_TRACE 1 ///< to turn off all tracing -#if defined(PS_NO_TRACE) -# define psTrace(facil, level, ...) (void)0 -/* do nothing */ -# define p_psTrace(facil, level, ...) (void)0 -/* do nothing */ -# define psTraceSetLevel(facil,level) 0 -# define psTraceGetLevel(facil) 0 -# define psTraceReset() (void)0 /* do nothing */ -# define psTraceFree() (void)0 /* do nothing */ -# define psTracePrintLevels() (void)0 -/* do nothing */ -# define psTraceSetDestination(fp) (void)0 -/* do nothing */ -# else +// XXX EAM : the old 'empty' values of (void) 0 are dangerous +# if defined(PS_NO_TRACE) + # define psTrace(facil, level, ...) /* do nothing */ + # define p_psTrace(facil, level, ...) /* do nothing */ + # define psTraceSetLevel(facil,level) /* do nothing */ + # define psTraceGetLevel(facil) /* do nothing */ + # define psTraceReset() /* do nothing */ + # define psTraceFree() /* do nothing */ + # define psTracePrintLevels() /* do nothing */ + # define psTraceSetDestination(fp) /* do nothing */ + # define psTraceSetDestination() /* do nothing */ + # define PS_TRACE_ON 0 - /** Basic structure for the component tree. A component is a string of the - form aaa.bbb.ccc, and may itself contain further subcomponents. The - Component structure doesn't in fact contain it's full name, but only the - last part. */ - typedef struct p_psComponent - { - const char *name; // last part of name of component - psS32 level; // trace level for this component - bool p_psSpecified; - psS32 n; // number of subcomponents - struct p_psComponent* *subcomp; // next level of subcomponents - } + # else /* PS_NO_TRACE */ + # define PS_TRACE_ON 1 + + /** Basic structure for the component tree. A component is a string of the + form aaa.bbb.ccc, and may itself contain further subcomponents. The + Component structure doesn't in fact contain it's full name, but only the + last part. */ + typedef struct p_psComponent + { + const char *name; // last part of name of component + psS32 level; // trace level for this component + bool p_psSpecified; + psS32 n; // number of subcomponents + struct p_psComponent* *subcomp; // next level of subcomponents + } p_psComponent; @@ -72,7 +73,7 @@ #ifndef SWIG #define psTrace(facil, level, ...) p_psTrace(facil, level, __VA_ARGS__) -#endif +#endif /* SWIG */ -#endif +#endif /* DOXYGEN */ /// Set trace level @@ -82,5 +83,5 @@ /// Get the trace level -psS32 psTraceGetLevel(const char *facil) ///< facilty of interest +psS32 psTraceGetLevel(const char *facil) ///< facilty of interest ; @@ -94,7 +95,11 @@ void psTraceSetDestination(FILE * fp); +/// Get the current destination for trace messages. +FILE *psTraceGetDestination(void); + /* \} */// End of SystemGroup Functions -#endif +#endif /* PS_NO_TRACE */ -#endif +#endif /* PS_TRACE_H */ + Index: /branches/eam-psphot-branch/psLib/src/sysUtils/psTrace.c =================================================================== --- /branches/eam-psphot-branch/psLib/src/sysUtils/psTrace.c (revision 3787) +++ /branches/eam-psphot-branch/psLib/src/sysUtils/psTrace.c (revision 3788) @@ -9,6 +9,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.44 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-03-29 21:31:54 $ + * @version $Revision: 1.44.4.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-04-29 09:17:02 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -252,4 +252,8 @@ if (cRoot == NULL) { initTrace(); + } + + if (traceFP == NULL) { + traceFP = stderr; } @@ -425,23 +429,15 @@ psS32 i = 0; + if (traceFP == NULL) { + traceFP = stderr; + } + if (comp->name[0] == '\0') { return; } else { if (comp->level == PS_DEFAULT_TRACE_LEVEL) { - if (traceFP == NULL) { - printf("%*s%-*s %d\n", depth, "", 20 - depth, comp->name, - defLevel); - } else { - fprintf(traceFP,"%*s%-*s %d\n", depth, "", 20 - depth, comp->name, - defLevel); - } + fprintf(traceFP,"%*s%-*s %d\n", depth, "", 20 - depth, comp->name, defLevel); } else { - if (traceFP == NULL) { - printf("%*s%-*s %d\n", depth, "", 20 - depth, comp->name, - comp->level); - } else { - fprintf(traceFP, "%*s%-*s %d\n", depth, "", 20 - depth, comp->name, - comp->level); - } + fprintf(traceFP, "%*s%-*s %d\n", depth, "", 20 - depth, comp->name, comp->level); } } @@ -498,4 +494,8 @@ psS32 i = 0; + if (traceFP == NULL) { + traceFP = stderr; + } + if (NULL == comp) { psError(PS_ERR_BAD_PARAMETER_NULL, true, @@ -514,26 +514,43 @@ fmt = va_arg(ap, char *); - if (traceFP != NULL) { - // We indent each message one space for each level of the message. - for (i = 0; i < level; i++) { - fprintf(traceFP, " "); - } - vfprintf(traceFP, fmt, ap); - } else { - // We indent each message one space for each level of the message. - for (i = 0; i < level; i++) { - putchar(' '); - } - vprintf(fmt, ap); - } + // We indent each message one space for each level of the message. + for (i = 0; i < level; i++) { + fprintf(traceFP, " "); + } + vfprintf(traceFP, fmt, ap); va_end(ap); } - // NOTE: should we free *fmt as well? Read the man page. -} - +} + +// XXX EAM : I've added code to close the old traceFP (safely) void psTraceSetDestination(FILE * fp) { + + bool special; + + // XXX EAM perhaps return an error? + if (fp == NULL) { + return; + } + + // cannot close traceFP if one of the special FILE ptrs + special = (traceFP == NULL); + special |= (traceFP == stdin); + special |= (traceFP == stdout); + special |= (traceFP == stderr); + + if (!special) { + fclose (traceFP); + } traceFP = fp; } +FILE *psTraceGetDestination() +{ + if (traceFP == NULL) { + traceFP = stderr; + } + return traceFP; +} + #endif Index: /branches/eam-psphot-branch/psLib/src/sysUtils/psTrace.h =================================================================== --- /branches/eam-psphot-branch/psLib/src/sysUtils/psTrace.h (revision 3787) +++ /branches/eam-psphot-branch/psLib/src/sysUtils/psTrace.h (revision 3788) @@ -9,6 +9,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.31 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-02-17 19:26:24 $ + * @version $Revision: 1.31.6.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-04-29 09:17:02 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -29,31 +29,32 @@ //#define PS_NO_TRACE 1 ///< to turn off all tracing -#if defined(PS_NO_TRACE) -# define psTrace(facil, level, ...) (void)0 -/* do nothing */ -# define p_psTrace(facil, level, ...) (void)0 -/* do nothing */ -# define psTraceSetLevel(facil,level) 0 -# define psTraceGetLevel(facil) 0 -# define psTraceReset() (void)0 /* do nothing */ -# define psTraceFree() (void)0 /* do nothing */ -# define psTracePrintLevels() (void)0 -/* do nothing */ -# define psTraceSetDestination(fp) (void)0 -/* do nothing */ -# else +// XXX EAM : the old 'empty' values of (void) 0 are dangerous +# if defined(PS_NO_TRACE) + # define psTrace(facil, level, ...) /* do nothing */ + # define p_psTrace(facil, level, ...) /* do nothing */ + # define psTraceSetLevel(facil,level) /* do nothing */ + # define psTraceGetLevel(facil) /* do nothing */ + # define psTraceReset() /* do nothing */ + # define psTraceFree() /* do nothing */ + # define psTracePrintLevels() /* do nothing */ + # define psTraceSetDestination(fp) /* do nothing */ + # define psTraceSetDestination() /* do nothing */ + # define PS_TRACE_ON 0 - /** Basic structure for the component tree. A component is a string of the - form aaa.bbb.ccc, and may itself contain further subcomponents. The - Component structure doesn't in fact contain it's full name, but only the - last part. */ - typedef struct p_psComponent - { - const char *name; // last part of name of component - psS32 level; // trace level for this component - bool p_psSpecified; - psS32 n; // number of subcomponents - struct p_psComponent* *subcomp; // next level of subcomponents - } + # else /* PS_NO_TRACE */ + # define PS_TRACE_ON 1 + + /** Basic structure for the component tree. A component is a string of the + form aaa.bbb.ccc, and may itself contain further subcomponents. The + Component structure doesn't in fact contain it's full name, but only the + last part. */ + typedef struct p_psComponent + { + const char *name; // last part of name of component + psS32 level; // trace level for this component + bool p_psSpecified; + psS32 n; // number of subcomponents + struct p_psComponent* *subcomp; // next level of subcomponents + } p_psComponent; @@ -72,7 +73,7 @@ #ifndef SWIG #define psTrace(facil, level, ...) p_psTrace(facil, level, __VA_ARGS__) -#endif +#endif /* SWIG */ -#endif +#endif /* DOXYGEN */ /// Set trace level @@ -82,5 +83,5 @@ /// Get the trace level -psS32 psTraceGetLevel(const char *facil) ///< facilty of interest +psS32 psTraceGetLevel(const char *facil) ///< facilty of interest ; @@ -94,7 +95,11 @@ void psTraceSetDestination(FILE * fp); +/// Get the current destination for trace messages. +FILE *psTraceGetDestination(void); + /* \} */// End of SystemGroup Functions -#endif +#endif /* PS_NO_TRACE */ -#endif +#endif /* PS_TRACE_H */ +