Index: /branches/ccl_branches/ipponly-20191108/ppStats/Makefile.am =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppStats/Makefile.am (revision 41082) +++ /branches/ccl_branches/ipponly-20191108/ppStats/Makefile.am (revision 41082) @@ -0,0 +1,9 @@ +SUBDIRS = src + +CLEANFILES = *~ core core.* + +pkgconfigdir = $(libdir)/pkgconfig +pkgconfig_DATA= ppStats.pc + +EXTRA_DIST = \ + ppStats.pc.in Index: /branches/ccl_branches/ipponly-20191108/ppStats/autogen.sh =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppStats/autogen.sh (revision 41082) +++ /branches/ccl_branches/ipponly-20191108/ppStats/autogen.sh (revision 41082) @@ -0,0 +1,103 @@ +#!/bin/sh +# Run this to generate all the initial makefiles, etc. + +srcdir=`dirname $0` +test -z "$srcdir" && srcdir=. + +ORIGDIR=`pwd` +cd $srcdir + +PROJECT=ppStats +TEST_TYPE=-f +# change this to be a unique filename in the top level dir +FILE=autogen.sh + +DIE=0 + +LIBTOOLIZE=libtoolize +ACLOCAL="aclocal $ACLOCAL_FLAGS" +AUTOHEADER=autoheader +AUTOMAKE=automake +AUTOCONF=autoconf + +($LIBTOOLIZE --version) < /dev/null > /dev/null 2>&1 || { + echo + echo "You must have $LIBTOOLIZE installed to compile $PROJECT." + echo "Download the appropriate package for your distribution," + echo "or get the source tarball at http://ftp.gnu.org/gnu/libtool/" + DIE=1 +} + +($ACLOCAL --version) < /dev/null > /dev/null 2>&1 || { + echo + echo "You must have $ACLOCAL installed to compile $PROJECT." + echo "Download the appropriate package for your distribution," + echo "or get the source tarball at http://ftp.gnu.org/gnu/automake/" + DIE=1 +} + +($AUTOHEADER --version) < /dev/null > /dev/null 2>&1 || { + echo + echo "You must have $AUTOHEADER installed to compile $PROJECT." + echo "Download the appropriate package for your distribution," + echo "or get the source tarball at http://ftp.gnu.org/gnu/autoconf/" + DIE=1 +} + +($AUTOMAKE --version) < /dev/null > /dev/null 2>&1 || { + echo + echo "You must have $AUTOMAKE installed to compile $PROJECT." + echo "Download the appropriate package for your distribution," + echo "or get the source tarball at http://ftp.gnu.org/gnu/automake/" + DIE=1 +} + +($AUTOCONF --version) < /dev/null > /dev/null 2>&1 || { + echo + echo "You must have $AUTOCONF installed to compile $PROJECT." + echo "Download the appropriate package for your distribution," + echo "or get the source tarball at http://ftp.gnu.org/gnu/autoconf/" + DIE=1 +} + +if test "$DIE" -eq 1; then + exit 1 +fi + +test $TEST_TYPE $FILE || { + echo "You must run this script in the top-level $PROJECT directory" + exit 1 +} + +if test -z "$*"; then + echo "I am going to run ./configure with no arguments - if you wish " + echo "to pass any to it, please specify them on the $0 command line." +fi + +$LIBTOOLIZE --copy --force || echo "$LIBTOOLIZE failed" +$ACLOCAL || echo "$ACLOCAL failed" +$AUTOHEADER || echo "$AUTOHEADER failed" +$AUTOMAKE --add-missing --force-missing --copy || echo "$AUTOMAKE failed" +$AUTOCONF || echo "$AUTOCONF failed" + +cd $ORIGDIR + +run_configure=true +for arg in $*; do + case $arg in + --no-configure) + run_configure=false + ;; + *) + ;; + esac +done + +if $run_configure; then + $srcdir/configure --enable-maintainer-mode "$@" + echo + echo "Now type 'make' to compile $PROJECT." +else + echo + echo "Now run 'configure' and 'make' to compile $PROJECT." +fi Index: /branches/ccl_branches/ipponly-20191108/ppStats/configure.ac =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppStats/configure.ac (revision 41082) +++ /branches/ccl_branches/ipponly-20191108/ppStats/configure.ac (revision 41082) @@ -0,0 +1,38 @@ +AC_PREREQ(2.61) + +AC_INIT([ppStats], [1.1.1], [ipp-support@ifa.hawaii.edu]) +AC_CONFIG_SRCDIR([src]) + +AM_INIT_AUTOMAKE([1.6 foreign dist-bzip2]) +AM_CONFIG_HEADER([src/config.h]) +AM_MAINTAINER_MODE + +IPP_STDLDFLAGS + +AC_LANG(C) +AC_GNU_SOURCE +AC_PROG_CC +AC_PROG_INSTALL +AC_PROG_LIBTOOL +AC_SYS_LARGEFILE + +PKG_CHECK_MODULES([PSLIB], [pslib >= 1.0.0]) +PKG_CHECK_MODULES([PSMODULE], [psmodules >= 1.0.0]) + +IPP_STDOPTS +IPP_STDCFLAGS + +dnl not clear why the below is needed +CFLAGS="${CFLAGS=} -std=c99" + +IPP_VERSION + +AC_SUBST([PPSTATS_CFLAGS]) +AC_SUBST([PPSTATS_LIBS]) + +AC_CONFIG_FILES([ + Makefile + src/Makefile + ppStats.pc +]) +AC_OUTPUT Index: /branches/ccl_branches/ipponly-20191108/ppStats/ppStats.pc.in =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppStats/ppStats.pc.in (revision 41082) +++ /branches/ccl_branches/ipponly-20191108/ppStats/ppStats.pc.in (revision 41082) @@ -0,0 +1,11 @@ +prefix=@prefix@ +exec_prefix=@exec_prefix@ +libdir=@libdir@ +includedir=@includedir@ + +Name: libppStats +Description: Pan-STARRS Image Statistics +Version: @VERSION@ +Requires: pslib psmodules +Libs: -L${libdir} -lppStats +Cflags: -I${includedir} Index: /branches/ccl_branches/ipponly-20191108/ppStats/src/Makefile.am =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppStats/src/Makefile.am (revision 41082) +++ /branches/ccl_branches/ipponly-20191108/ppStats/src/Makefile.am (revision 41082) @@ -0,0 +1,58 @@ +lib_LTLIBRARIES = libppStats.la + +# Force recompilation of ppStatsVersion.c, since it gets the version information +ppStatsVersion.c: ppStatsVersionDefinitions.h +ppStatsVersionDefinitions.h: ppStatsVersionDefinitions.h.in FORCE + pslib-setsvnversion.pl PPSTATS ppStatsVersionDefinitions.h.in ppStatsVersionDefinitions.h +FORCE: ; + +libppStats_la_CFLAGS = $(PPSTATS_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS) +libppStats_la_LDFLAGS = $(PPSTATS_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS) + +bin_PROGRAMS = ppStats ppStatsFromMetadata + +ppStats_CFLAGS = $(PPSTATS_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS) +ppStats_LDFLAGS = $(PPSTATS_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS) +ppStats_LDADD = libppStats.la + +ppStatsFromMetadata_CFLAGS = $(PPSTATS_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS) +ppStatsFromMetadata_LDFLAGS = $(PPSTATS_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS) + +ppStats_SOURCES = \ + ppStats.c + +ppStatsFromMetadata_SOURCES = \ + ppStatsFromMetadata.c \ + ppStatsFromMetadataEntries.c \ + ppStatsFromMetadataParse.c \ + ppStatsFromMetadataStats.c \ + ppStatsFromMetadataPrint.c + +libppStats_la_SOURCES = \ + ppStatsFPA.c \ + ppStatsData.c \ + ppStatsFringe.c \ + ppStatsLoop.c \ + ppStatsChip.c \ + ppStatsCell.c \ + ppStatsReadout.c \ + ppStatsUtils.c \ + ppStatsPixels.c \ + ppStatsMetadata.c \ + ppStatsSetupFromRecipe.c \ + ppStatsSetupFromArgs.c \ + ppStatsVersion.c + +include_HEADERS = \ + ppStats.h + +noinst_HEADERS = \ + ppStatsInternal.h + +CLEANFILES = *~ + +clean-local: + -rm -f TAGS +# Tags for emacs +tags: + etags `find . -name \*.[ch] -print` Index: /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStats.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStats.c (revision 41082) +++ /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStats.c (revision 41082) @@ -0,0 +1,138 @@ +# include "ppStatsInternal.h" + +int main(int argc, char **argv) { + + psExit status = PS_EXIT_SUCCESS; + + psLibInit(NULL); + psTimerStart("PPSTATS"); + + // Parse the configuration and arguments + pmConfig *config = pmConfigRead(&argc, argv, PPSTATS_RECIPE); + if (!config) { + psErrorStackPrint(stderr, "Unable to read configuration.\n"); + exit(PS_EXIT_CONFIG_ERROR); + } + + // Get the options, open the files + ppStatsData *data = ppStatsSetupFromArgs(&argc, argv, config); + if (!data) { + psErrorStackPrint(stderr, "Unable to parse command-line arguments.\n"); + exit(PS_EXIT_CONFIG_ERROR); + } + + // Output filename is optional + const char *outName = NULL; // Output file name + FILE *outFile = stdout; // Output file + if (argc == 2) { + outName = argv[1]; + psString resolved = pmConfigConvertFilename(outName, config, true, true); // Resolved filename + + if (resolved && strlen(resolved) > 0) { + outFile = fopen(resolved, "w"); + if (!outFile) { + psLogMsg("ppStats", PS_LOG_ERROR, "Unable to open output file %s\n", resolved); + psFree(resolved); + // XXX this could be a system or config error, but not a data error + status = PS_EXIT_CONFIG_ERROR; + goto die; + } + } else { + psErrorStackPrint(stderr, "Unable to open output file %s.\n", resolved); + exit(PS_EXIT_CONFIG_ERROR); + } + psFree(resolved); + } + + // Go through the FPA and do the hard work + psMetadata *results = ppStatsLoop(&status, data, config); + if (status != PS_EXIT_SUCCESS) { + psErrorStackPrint(stderr, "Error in stats loop.\n"); + exit (status); + } + + if (data->showFormat) { + char *formatName = psStringCopy (config->formatName); + if (!formatName) { + formatName = psStringCopy ("UNKNOWN"); + } + psMetadataAddStr(results, PS_LIST_HEAD, "FILE.FORMAT", 0, "File format", formatName); + } + + if (data->showCamera) { + char *cameraName = psStringCopy (config->cameraName); + if (!cameraName) { + cameraName = psStringCopy ("UNKNOWN"); + } + psMetadataAddStr(results, PS_LIST_HEAD, "CAMERA", 0, "camera name", cameraName); + } + + // report on the file disposition + if (data->fileLevel) { + pmFPALevel level = pmFPAPHULevel(config->format); + + const char *levelName = pmFPALevelToName(level); // Level for file + psMetadataAddStr(results, PS_LIST_HEAD, "FILE.LEVEL", 0, "File level", levelName); + + char *classID = NULL; + switch (level) { + case PM_FPA_LEVEL_FPA: + assert (data->fpa != NULL); + assert (data->fileView->chip == -1); + classID = pmFPANameFromRule ("{FPA.NAME}", data->fpa, data->fileView); + break; + case PM_FPA_LEVEL_CHIP: + assert (data->fpa != NULL); + assert (data->fileView->chip != -1); + assert (data->fileView->cell == -1); + classID = pmFPANameFromRule ("{CHIP.NAME}", data->fpa, data->fileView); + break; + case PM_FPA_LEVEL_CELL: + assert (data->fpa != NULL); + assert (data->fileView->chip != -1); + assert (data->fileView->cell != -1); + classID = pmFPANameFromRule ("{CELL.NAME}", data->fpa, data->fileView); + break; + default: + psErrorStackPrint(stderr, "Error in file level.\n"); + exit (PS_EXIT_CONFIG_ERROR); + } + psMetadataAddStr(results, PS_LIST_HEAD, "CLASS.ID", 0, "name for element at file level", classID); + psFree (classID); + } + + // did we actually request any data? + if (psListLength(results->list) == 0) { + psErrorStackPrint(stderr, "No output.\n"); + exit (status); + } + + // Format and print the output + psString output = psMetadataConfigFormat(results); + if (!output) { + psErrorStackPrint(stderr, "Unable to generate configuration file with result.\n"); + psFree(results); + exit(PS_EXIT_CONFIG_ERROR); + } + fprintf(outFile, "%s", output); + psFree(output); + + // Clean up + psFree(results); + if (outName) { + fclose(outFile); + } + + // Common code for the death. +die: + if (status) { + psErrorStackPrint (stderr, "failure in %s", __func__); + } + psFree(data); + psFree(config); + pmConceptsDone(); + pmConfigDone(); + psLibFinalize(); + + return status; +} Index: /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStats.config =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStats.config (revision 41082) +++ /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStats.config (revision 41082) @@ -0,0 +1,25 @@ +### Example ppStats recipe + +MULTI HEADERS # May specify multiple HEADERS +HEADER STR HEADER1,HEADER2 # Read these headers +HEADER STR HEADER3 # Read these headers too + +MULTI STATS # May specify multiple STATS +STAT STR SAMPLE_MEAN,SAMPLE_STDEV # Calculate these statistics +STAT STR ROBUST_MEDIAN # Calculate these statistics too + +MULTI ANALYSIS # metadata blocks to search in chip/cell/readout->analysis +ANALYSIS STR PSPHOT.HEADER +ANALYSIS STR PSASTRO.HEADER + +SAMPLE F32 0.1 # Fraction of cell to sample + +MULTI CHIPS # May specify multiple CHIPS +CHIPS STR chip00, chip01 # Look at these chips +CHIPS STR chip22 # Look at this chip too + +MULTI CELLS # May specify multiple CELLS +CELLS STR LeftCell,RightCell # Look at these cells +CELLS STR UpsideDownCell # Look at this cell too + +MASKVAL U8 0xff Index: /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStats.h =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStats.h (revision 41082) +++ /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStats.h (revision 41082) @@ -0,0 +1,161 @@ + +#ifndef PP_STATS_H +#define PP_STATS_H + +#define PPSTATS_RECIPE "PPSTATS" +#define PPSTATS_MD_RECIPE "PPSTATS_METADATA" + +typedef struct { + // Inputs + psFits *fits; // Input file handle + pmFPA *fpa; // FPA to analyse + pmFPAview *view; // View to analyse + // Stuff to output + psStats *stats; // pixel Statistics to calculate + bool doStats; // Do pixel statistics? + bool fileLevel; // Output file level? + bool showFormat; // Output file format? + bool showCamera; // Output camera name? + pmFPAview *fileView; // View to analyse + + psList *headers; // Headers to read + psList *concepts; // Concepts to read + psList *analysis; // Analysis entries to read + psList *summary; // Summary statistics to calculate + // Options for input data + bool doFirstReadout3D; // for 3D data, use the first readout? + float sample; // Fraction of cell to sample for statistics + psImageMaskType maskVal; // Mask value for images + psList *chips; // Chips to look at + psList *cells; // Cells to look at +} ppStatsData; + +typedef struct { + char *keyword; + psDataType type; + char *statistic; + char *flag; + psMetadataItem *value; + psVector *vector; +} ppStatsEntry; + +// Allocator +ppStatsData *ppStatsDataAlloc(void); + +/// Perform the ppStats steps on the given FPA (optionally for specified view) +psMetadata *ppStatsFPA(psMetadata *out, + pmFPA *fpa, // FPA for which to get statistics + pmFPAview *view, // View for analysis + psImageMaskType maskVal, // Value to mask + pmConfig *config // Configuration + ); + +psExit ppStatsChip(psMetadata *fpaResults, // Metadata holding the fpa results + pmChip *chip, // Chip for which to get statistics + psFits *fits, // FITS file handle + pmFPAview *view, // View for analysis + ppStatsData *data,// The data + pmConfig *config // Configuration + ); + +psExit ppStatsCell(psMetadata *chipResults, // Metadata holding the chip results + pmCell *cell, // Cell for which to get statistics + psFits *fits, // FITS file handle + pmFPAview *view, // View for analysis + ppStatsData *data,// The data + pmConfig *config // Configuration + ); + +/// Supplement the statistics with the fringe solution +bool ppStatsFringe(psMetadata *stats, ///< Statistics metadata to supplement + const pmChip *chip, ///< The chip containing the solution + const char *root, ///< Name of output entry + const char *fringeName ///< Name of the solution in the chip->analysis + ); + + +// measure only the pixel-related statistics (stats, summary) +psMetadata *ppStatsPixels(psMetadata *out, + pmFPA *fpa, // FPA for which to get statistics + pmFPAview *view, // View for analysis + psImageMaskType maskVal, // Value to mask + pmConfig *config // Configuration + ); + +// measure only the non-pixel-related statistics (headers, concepts) +psMetadata *ppStatsMetadata(psMetadata *out, + pmFPA *fpa, // FPA for which to get statistics + pmFPAview *view, // View for analysis + psImageMaskType maskVal, // Value to mask + pmConfig *config // Configuration + ); + +/// Loop over the input image and do all the hard work +psMetadata *ppStatsLoop(psExit *result, + ppStatsData *data, // The data + pmConfig *config // Configuration + ); + +/// Set up the options and input/output files +ppStatsData *ppStatsSetupFromArgs(int *argc, char *argv[], // Command-line arguments + pmConfig *config // Configuration + ); + +bool ppStatsSetupFromRecipe(ppStatsData *data, // Data for running ppStats + pmConfig *config // Configuration + ); + + +/// Return short version information +psString ppStatsVersion(void); + +/// Return source information +psString ppStatsSource(void); + +/// Return long version information +psString ppStatsVersionLong(void); + +/// Populate header with version information +bool ppStatsVersionHeader( + psMetadata *header ///< Header to populate + ); + + +void p_ppStatsGetMetadata(psMetadata *target, // Target for metadata + psMetadata *source, // Source for metadata + psList *list // List containing keywords + ); + +void p_ppStatsGetAnalysis(psMetadata *target, // Output Target for metadata + psList *headers, // List containing desired keywords + psMetadata *source, // Input Source for metadata + psList *list // List containing analysis blocks + ); + +bool p_ppStatsDoThis(psList *toDoList, // List of things to do + const char *this // The name of "this" + ); + +void p_ppStatsAddToHierarchy(psMetadata *source, // Source to add + psMetadata *target, // Target to which to add + const char *name, // Name of source + const char *comment // Comment for source + ); + +psExit ppStatsReadout(psMetadata *cellResults, // Metadata holding the chip results + pmReadout *readout, // Cell for which to get statistics + int nReadout, // readout number + ppStatsData *data, // The data + const pmConfig *config // Configuration + ); + +// used by ppStatsFromMetadata: +psDataType psDataTypeFromString(char *typename); +ppStatsEntry *ppStatsEntryAlloc(void); + +psArray *ppStatsFromMetadataEntries (psMetadata *recipe); +bool ppStatsFromMetadataParse (psMetadata *input, psArray *entries); +bool ppStatsFromMetadataStats (psArray *entries); +bool ppStatsFromMetadataPrint (psArray *entries, char *filename); + +#endif Index: /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsCell.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsCell.c (revision 41082) +++ /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsCell.c (revision 41082) @@ -0,0 +1,123 @@ +# include "ppStatsInternal.h" + +psExit ppStatsCell(psMetadata *chipResults, // Metadata holding the chip results + pmCell *cell, // Cell for which to get statistics + psFits *fits, // FITS file handle + pmFPAview *view, // View for analysis + ppStatsData *data, // The data + pmConfig *config // Configuration + ) +{ + assert(chipResults); + assert(cell); + assert(data); + assert(config); + + const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME"); // Name of cell + + // Check to see if this is a cell of interest + if (!p_ppStatsDoThis(data->cells, cellName)) { + return PS_EXIT_SUCCESS; + } + +#if 0 + // select the header unit for this cell + // XXX why is this needed for Cell but not Chip? + pmHDU *hdu = pmHDUFromCell(cell); // HDU for cell + if (!hdu || hdu->blankPHU) { + if (fits) { + // No HDU means there's no data in this cell + return PS_EXIT_SUCCESS; + } else { + psError(PS_ERR_UNKNOWN, false, "Can't find HDU for cell\n"); + return PS_EXIT_CONFIG_ERROR; + } + } +#endif + + // Extract Header and Concept values from the Cell and Readout->analysis level + + // find or generate a metadata to hold the results + bool mdok; // Status of MD lookup + psMetadata *cellResults = psMemIncrRefCounter(psMetadataLookupMetadata(&mdok, chipResults, cellName)); + if (!mdok || !cellResults) { + cellResults = psMetadataAlloc(); + } + + // Extract Header values + if (psListLength(data->headers)) { + // extract from existing headers + if (cell->hdu) { + if (fits && !pmCellReadHeader(cell, fits, config)) { + psError (PS_ERR_IO, false, "trouble reading cell header\n"); + psFree(cellResults); + return PS_EXIT_DATA_ERROR; + } + pmHDU *hdu = cell->hdu; // HDU for headers + p_ppStatsGetMetadata(cellResults, hdu->header, data->headers); + } + // extract from data->analysis output MD entries + if (psListLength(data->analysis)) { + p_ppStatsGetAnalysis (cellResults, data->headers, cell->analysis, data->analysis); + } + } + + // Extract Concept values + if (psListLength(data->concepts) > 0) { + if (fits && cell->hdu && !pmCellReadHeader(cell, fits, config)) { + psError (PS_ERR_IO, false, "trouble reading cell header\n"); + psFree(cellResults); + return PS_EXIT_DATA_ERROR; + } + pmConceptsReadCell(cell, PM_CONCEPT_SOURCE_ALL, false, config); + p_ppStatsGetMetadata(cellResults, cell->concepts, data->concepts); + } + + // If we want to measure pixel statistics, we must read (or have) the pixel data + if (data->doStats || psListLength(data->summary)) { + // Read the image pixel data + if (fits && !pmCellRead(cell, fits, config)) { + psError (PS_ERR_IO, false, "trouble reading cell data\n"); + psFree(cellResults); + return PS_EXIT_DATA_ERROR; + } + } + + // check if we can legitimately iterate to the readout level + psArray *readouts = cell->readouts; // Array of component readouts + if (view->readout >= readouts->n) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Desired readout view (%d) doesn't match " + "number of readouts (%ld)\n", view->readout, readouts->n); + return PS_EXIT_CONFIG_ERROR; + } + if ((view->readout == -1) && (readouts->n > 1)) { + psWarning("Multiple readouts (%ld) present in cell %s\n", readouts->n, cellName); + if (!data->doFirstReadout3D) { + psWarning ("Skipping Video Cell for ppStatsCell"); + return PS_EXIT_SUCCESS; + } + } + + // Iterate over readouts + for (int i = 0; i < readouts->n; i++) { + if ((view->readout >= 0) && (i != view->readout)) continue; + pmReadout *readout = readouts->data[i]; // Cell of interest + psExit result = ppStatsReadout(cellResults, readout, i, data, config); + if (result != PS_EXIT_SUCCESS) { + psError(PS_ERR_UNKNOWN, false, "trouble with readout stats for %d\n", i); + if (fits) { + pmCellFreeData(cell); + } + psFree(cellResults); + return result; + } + } + + // Add the cell results to the chip + p_ppStatsAddToHierarchy(cellResults, chipResults, cellName, "Results for cell"); + psFree (cellResults); + if (fits) { + pmCellFreeData(cell); + } + return PS_EXIT_SUCCESS; +} Index: /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsChip.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsChip.c (revision 41082) +++ /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsChip.c (revision 41082) @@ -0,0 +1,103 @@ +# include "ppStatsInternal.h" + +psExit ppStatsChip(psMetadata *fpaResults, // Metadata holding the fpa results + pmChip *chip, // Chip for which to get statistics + psFits *fits, // FITS file handle + pmFPAview *view, // View for analysis + ppStatsData *data,// The data + pmConfig *config // Configuration + ) +{ + assert(fpaResults); + assert(chip); + assert(view); + assert(data); + assert(config); + + const char *chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME"); // Name of chip + + // Check to see if this is a chip of interest + if (!p_ppStatsDoThis(data->chips, chipName)) { + return PS_EXIT_SUCCESS; + } + if (!chip->file_exists && !chip->data_exists) { + return PS_EXIT_SUCCESS; + } + + // Extract Header and Concept values from the Chip level + + // find or generate a metadata to hold the results + bool mdok; // Status of MD lookup + psMetadata *chipResults = psMemIncrRefCounter(psMetadataLookupMetadata(&mdok, fpaResults, chipName)); + if (!mdok || !chipResults) { + chipResults = psMetadataAlloc(); + } + + // Extract Header values + if (psListLength(data->headers)) { + // extract from existing headers + if (chip->hdu) { + if (fits && !pmChipReadHeader(chip, fits, config)) { + psError (PS_ERR_IO, false, "trouble reading chip header\n"); + psFree(chipResults); + return PS_EXIT_DATA_ERROR; + } + pmHDU *hdu = chip->hdu; // HDU for headers + p_ppStatsGetMetadata(chipResults, hdu->header, data->headers); + } + // extract from data->analysis output MD entries + if (psListLength(data->analysis)) { + p_ppStatsGetAnalysis (chipResults, data->headers, chip->analysis, data->analysis); + } + } + + // Extract Concept values + if (psListLength(data->concepts) > 0) { + if (fits && chip->hdu && !pmChipReadHeader(chip, fits, config)) { + psError (PS_ERR_IO, false, "trouble reading chip header\n"); + psFree(chipResults); + return PS_EXIT_DATA_ERROR; + } + pmConceptsReadChip(chip, PM_CONCEPT_SOURCE_ALL, false, false, config); + p_ppStatsGetMetadata(chipResults, chip->concepts, data->concepts); + } + + // check if we can legitimately iterate to the cell level + psArray *cells = chip->cells; // Array of cells + if (view->cell >= cells->n) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Desired cell view (%d) doesn't match " + "number of cells (%ld)\n", view->cell, cells->n); + return PS_EXIT_CONFIG_ERROR; + } + + // Iterate over cells (if view->cell is set, skip all others) + for (int i = 0; i < cells->n; i++) { + if ((view->cell >= 0) && (i != view->cell)) continue; + pmCell *cell = cells->data[i]; // Cell of interest + + // XXX for now, skip the video cells (cell->readouts->n > 1) + if ((cell->readouts->n > 1) && !data->doFirstReadout3D){ + psWarning ("Skipping Video Cell for ppStatsCell"); + continue; + } + + psExit result = ppStatsCell(chipResults, cell, fits, view, data, config); + if (result != PS_EXIT_SUCCESS) { + psError(PS_ERR_UNKNOWN, false, "trouble with cell stats for %d\n", i); + if (fits) { + pmChipFreeData(chip); + } + psFree(chipResults); + return result; + } + } + + if (fits) { + pmChipFreeData(chip); + } + + p_ppStatsAddToHierarchy(chipResults, fpaResults, chipName, "Results for chip"); + + psFree (chipResults); + return PS_EXIT_SUCCESS; +} Index: /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsData.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsData.c (revision 41082) +++ /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsData.c (revision 41082) @@ -0,0 +1,53 @@ +#include "ppStatsInternal.h" + +static void statsDataFree(ppStatsData *data // Data to free + ) +{ + if (data->fits) { + psFitsClose(data->fits); + data->fits = NULL; + } + psFree(data->fpa); + psFree(data->view); + psFree(data->headers); + psFree(data->concepts); + psFree(data->analysis); + psFree(data->summary); + + psFree(data->stats); + psFree(data->fileView); + + psFree(data->chips); + psFree(data->cells); + + return; +} + + +ppStatsData *ppStatsDataAlloc(void) +{ + ppStatsData *data = psAlloc(sizeof(ppStatsData)); // Newly allocated data + psMemSetDeallocator(data, (psFreeFunc)statsDataFree); + + data->fits = NULL; + data->fpa = NULL; + data->view = NULL; + + data->headers = psListAlloc(NULL); + data->concepts = psListAlloc(NULL); + data->analysis = psListAlloc(NULL); + data->summary = psListAlloc(NULL); + data->stats = psStatsAlloc(0); + data->doStats = false; + data->fileLevel = false; + data->fileView = NULL; + + data->sample = 0; + data->maskVal = 0; + data->doFirstReadout3D = false; + data->chips = psListAlloc(NULL); + data->cells = psListAlloc(NULL); + + return data; +} + Index: /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsErrorCodes.c.in =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsErrorCodes.c.in (revision 41082) +++ /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsErrorCodes.c.in (revision 41082) @@ -0,0 +1,26 @@ +/* + * The line + { PPSTATS_ERR_$X{ErrorCode}, "$X{ErrorDescription}"}, + * (without the Xs) + * will be replaced by values from errorCodes.dat + */ +#include "pslib.h" +#include "ppStatsErrorCodes.h" + +void ppStatsErrorRegister(void) +{ + static psErrorDescription errors[] = { + { PPSTATS_ERR_BASE, "First value we use; lower values belong to psLib" }, + { PPSTATS_ERR_${ErrorCode}, "${ErrorDescription}"}, + }; + static int nerror = PPSTATS_ERR_NERROR - PPSTATS_ERR_BASE; // number of values in enum + + for (int i = 0; i < nerror; i++) { + psErrorDescription *tmp = psAlloc(sizeof(psErrorDescription)); + p_psMemSetPersistent(tmp, true); + *tmp = errors[i]; + psErrorRegister(tmp, 1); + psFree(tmp); /* it's on the internal list */ + } + nerror = 0; // don't register more than once +} Index: /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsErrorCodes.dat =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsErrorCodes.dat (revision 41082) +++ /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsErrorCodes.dat (revision 41082) @@ -0,0 +1,10 @@ +# +# This file is used to generate ppStatsErrorClasses.h +# +BASE = 2000 First value we use; lower values belong to psLib +# these errors correspond to standard exit conditions +ARGUMENTS Incorrect arguments +SYS System error +CONFIG Problem in configure files +PROG Programming error +DATA invalid data Index: /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsErrorCodes.h.in =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsErrorCodes.h.in (revision 41082) +++ /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsErrorCodes.h.in (revision 41082) @@ -0,0 +1,18 @@ +#if !defined(PPSTATS_ERROR_CODES_H) +#define PPSTATS_ERROR_CODES_H +/* + * The line + * PPSTATS_ERR_$X{ErrorCode}, + * (without the X) + * + * will be replaced by values from errorCodes.dat + */ +typedef enum { + PPSTATS_ERR_BASE = 2000, + PPSTATS_ERR_${ErrorCode}, + PPSTATS_ERR_NERROR +} ppStatsErrorCode; + +void ppStatsErrorRegister(void); + +#endif Index: /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsFPA.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsFPA.c (revision 41082) +++ /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsFPA.c (revision 41082) @@ -0,0 +1,55 @@ +#include "ppStatsInternal.h" + +psMetadata *ppStatsFPA(psMetadata *out, + pmFPA *fpa, // FPA for which to get statistics + pmFPAview *view, // View for analysis + psImageMaskType maskVal, // Value to mask + pmConfig *config // Configuration + ) +{ + PS_ASSERT_PTR_NON_NULL(fpa, NULL) + PS_ASSERT_PTR_NON_NULL(view, NULL); + PS_ASSERT_PTR_NON_NULL(config, NULL); + + ppStatsData *data = ppStatsDataAlloc(); // All the input data + + // Get the options, open the files + if (!ppStatsSetupFromRecipe(data, config)) { + psError(PS_ERR_UNKNOWN, false, "Unable to get ppStats options from recipe."); + psFree(data); + return NULL; + } + + // Override recipe mask value + data->maskVal = maskVal; + + if (data->fpa) { + psFree(data->fpa); + } + data->fpa = psMemIncrRefCounter(fpa); + + if (data->view) { + psFree(data->view); + } + data->view = psMemIncrRefCounter(view); + + // Go through the FPA and do the hard work + psExit status; // Status of statistics loop + psMetadata *result = ppStatsLoop(&status, data, config); + if (status != PS_EXIT_SUCCESS) { + psError (PS_ERR_UNKNOWN, false, "Not able to measure FPA statistics.\n"); + psFree(result); + psFree(data); + return (NULL); + } + + if (out != NULL) { + psMetadataOverlay (out, result); + psFree(result); + psFree(data); + return out; + } + + psFree(data); + return result; +} Index: /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsFringe.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsFringe.c (revision 41082) +++ /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsFringe.c (revision 41082) @@ -0,0 +1,37 @@ +#include "ppStatsInternal.h" + +bool ppStatsFringe(psMetadata *stats, const pmChip *chip, const char *root, const char *fringeName) +{ + PS_ASSERT_PTR_NON_NULL(chip, false); + PS_ASSERT_STRING_NON_EMPTY(fringeName, false); + + bool mdok; // Status of MD lookup + const char *chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME"); // Name of chip + + psMetadata *chipStats = psMetadataLookupMetadata(&mdok, stats, chipName); // Chip statistics + if (!mdok || !chipStats) { + chipStats = psMetadataAlloc(); + psMetadataAddMetadata(stats, PS_LIST_TAIL, chipName, 0, NULL, chipStats); + } else { + psMemIncrRefCounter(chipStats); + } + + pmFringeScale *fringes = psMetadataLookupPtr(NULL, chip->analysis, fringeName); // Solution for fringes + for (int i = 0; i < fringes->nFringeFrames; i++) { + psString name = NULL; // Name of statistic + psStringAppend(&name, "%s_%d", root, i); + psMetadataAddF32(chipStats, PS_LIST_TAIL, name, 0, "Fringe amplitude", + fringes->coeff->data.F32[i + 1]); + psFree(name); + name = NULL; + psStringAppend(&name, "%s_ERR_%d", root, i); + psMetadataAddF32(chipStats, PS_LIST_TAIL, name, 0, "Fringe amplitude error", + fringes->coeffErr->data.F32[i + 1]); + psFree(name); + } + + psFree(chipStats); + + return true; +} + Index: /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsFromMetadata.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsFromMetadata.c (revision 41082) +++ /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsFromMetadata.c (revision 41082) @@ -0,0 +1,68 @@ +# include "ppStatsInternal.h" + +// USAGE: ppStatsFromMetadata (input) (output) [-recipe PPSTATS_METADATA recipe] +int main(int argc, char **argv) +{ + bool status; + + psLibInit(NULL); + + // Parse the configuration and arguments + pmConfig *config = pmConfigRead(&argc, argv, PPSTATS_MD_RECIPE); + if (!config) { + psErrorStackPrint(stderr, "Unable to read configuration.\n"); + exit(PS_EXIT_CONFIG_ERROR); + } + + if (argc != 4) { + psErrorStackPrint(stderr, "USAGE: ppStatsFromMetadata (input) (output) (recipe)\n"); + exit(PS_EXIT_CONFIG_ERROR); + } + + psMetadata *input = NULL; + if (!strcmp (argv[1], "-")) { + psString string = psSlurpFile (stdin); + input = psMetadataConfigParse (NULL, NULL, string, false); + psFree (string); + } else { + input = psMetadataConfigRead (NULL, NULL, argv[1], false); + } + + // parse the recipe to determine the fields of interest + psMetadata *recipes = psMetadataLookupPtr (&status, config->recipes, PPSTATS_MD_RECIPE); + if (!recipes) { + psErrorStackPrint(stderr, "missing recipe set %s\n", PPSTATS_MD_RECIPE); + exit(PS_EXIT_CONFIG_ERROR); + } + + psMetadata *recipe = psMetadataLookupPtr (&status, recipes, argv[3]); + if (!recipes) { + psErrorStackPrint(stderr, "missing ppStatsFromMetadata recipe %s\n", argv[3]); + exit(PS_EXIT_CONFIG_ERROR); + } + + psArray *entries = ppStatsFromMetadataEntries (recipe); + if (!entries) { + psErrorStackPrint(stderr, "problem with recipe %s\n", PPSTATS_MD_RECIPE); + exit(PS_EXIT_CONFIG_ERROR); + } + + if (!ppStatsFromMetadataParse (input, entries)) { + psErrorStackPrint(stderr, "problem with input data\n"); + exit(PS_EXIT_CONFIG_ERROR); + } + + // calculate the stats for the non-constant entries (already calculated) + if (!ppStatsFromMetadataStats (entries)) { + psErrorStackPrint(stderr, "problem calculating statistics\n"); + exit(PS_EXIT_CONFIG_ERROR); + } + + // calculate the stats for the non-constant entries (already calculated) + if (!ppStatsFromMetadataPrint (entries, argv[2])) { + psErrorStackPrint(stderr, "problem calculating statistics\n"); + exit(PS_EXIT_CONFIG_ERROR); + } + + exit (0); +} Index: /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsFromMetadataEntries.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsFromMetadataEntries.c (revision 41082) +++ /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsFromMetadataEntries.c (revision 41082) @@ -0,0 +1,114 @@ +# include "ppStatsInternal.h" + +psArray *ppStatsFromMetadataEntries (psMetadata *recipe) { + + bool status = false; + + psMetadataItem *item = NULL; + + psArray *entries = psArrayAllocEmpty(16); + + // the recipe file should include a collection of psMetadata items with the name ENTRY + + // psMetadataPrint (stderr, recipe, 1); + + // loop over the items, selecting those with the name "ENTRY" + psMetadataIterator *iter = psMetadataIteratorAlloc(recipe, PS_LIST_HEAD, NULL); // Iterator + while ((item = psMetadataGetAndIncrement(iter))) { + if (strcmp (item->name, "ENTRY")) continue; + + if (item->type != PS_DATA_METADATA) { + psError(PS_ERR_UNKNOWN, false, "ppStatsFromMetadata ENTRY not of type METADATA\n"); + return NULL; + } + + psMetadata *entryMD = item->data.md; + + ppStatsEntry *entry = ppStatsEntryAlloc (); + + entry->keyword = psMetadataLookupStr (&status, entryMD, "KEYWORD"); + if (!entry->keyword) { + psError(PS_ERR_UNKNOWN, true, "ppStatsFromMetadata recipe ENTRY missing KEYWORD value\n"); + return NULL; + } + + entry->statistic = psMetadataLookupStr (&status, entryMD, "STATISTIC"); + if (!entry->statistic) { + psError(PS_ERR_UNKNOWN, true, "ppStatsFromMetadata recipe ENTRY %s missing STATISTIC\n", entry->keyword); + return NULL; + } + + entry->flag = psMetadataLookupStr (&status, entryMD, "FLAG"); + if (!entry->flag) { + psError(PS_ERR_UNKNOWN, true, "ppStatsFromMetadata recipe ENTRY %s missing FLAG\n", entry->keyword); + return NULL; + } + + char *typename = psMetadataLookupStr (&status, entryMD, "TYPE"); + if (!typename) { + psError(PS_ERR_UNKNOWN, true, "ppStatsFromMetadata recipe ENTRY %s missing TYPE\n", entry->keyword); + return NULL; + } + + entry->type = psDataTypeFromString (typename); + if (!entry->type) { + psError(PS_ERR_UNKNOWN, true, "ppStatsFromMetadata recipe ENTRY %s has invalid TYPE %s\n", entry->keyword, typename); + return NULL; + } + + // fprintf (stderr, "adding %s to array of %ld\n", entry->keyword, entries->n); + psArrayAdd (entries, 16, entry); + psFree (entry); + } + + return entries; +} + +void ppStatsEntryFree (ppStatsEntry *entry) { + + if (!entry) return; + + psFree (entry->keyword); + psFree (entry->statistic); + psFree (entry->flag); + psFree (entry->value); + psFree (entry->vector); + return; +} + +ppStatsEntry *ppStatsEntryAlloc () { + + ppStatsEntry *entry = (ppStatsEntry *) psAlloc (sizeof(ppStatsEntry)); + psMemSetDeallocator(entry, (psFreeFunc) ppStatsEntryFree); + + entry->keyword = NULL; + entry->statistic = NULL; + entry->flag = NULL; + entry->value = NULL; + entry->vector = NULL; + return entry; +} + +psDataType psDataTypeFromString (char *typename) { + + if (!strcmp (typename, "U8")) return PS_DATA_U8; + if (!strcmp (typename, "U16")) return PS_DATA_U16; + if (!strcmp (typename, "U32")) return PS_DATA_U32; + if (!strcmp (typename, "U64")) return PS_DATA_U64; + + if (!strcmp (typename, "S8")) return PS_DATA_S8; + if (!strcmp (typename, "S16")) return PS_DATA_S16; + if (!strcmp (typename, "S32")) return PS_DATA_S32; + if (!strcmp (typename, "S64")) return PS_DATA_S64; + + if (!strcmp (typename, "F32")) return PS_DATA_F32; + if (!strcmp (typename, "F64")) return PS_DATA_F64; + + if (!strcmp (typename, "BOOL")) return PS_DATA_BOOL; + if (!strcmp (typename, "STR")) return PS_DATA_STRING; + if (!strcmp (typename, "TIME")) return PS_DATA_TIME; + if (!strcmp (typename, "UTC")) return PS_DATA_TIME; + if (!strcmp (typename, "TAI")) return PS_DATA_TIME; + + return 0; +} Index: /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsFromMetadataParse.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsFromMetadataParse.c (revision 41082) +++ /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsFromMetadataParse.c (revision 41082) @@ -0,0 +1,146 @@ +# include "ppStatsInternal.h" + +// loop over the metadata, adding entries of interest to their data arrays +bool ppStatsFromMetadataParse (psMetadata *input, psArray *entries) { + + psMetadataItem *item = NULL; + + // loop over the items, selecting those with the name "ENTRY" + psMetadataIterator *iter = psMetadataIteratorAlloc(input, PS_LIST_HEAD, NULL); // Iterator + + while ((item = psMetadataGetAndIncrement(iter))) { + if (item->type == PS_DATA_METADATA) { + psMetadata *folder = item->data.md; + if (!ppStatsFromMetadataParse (folder, entries)) { + psError(PS_ERR_UNKNOWN, false, "error parse metadata folder %s\n", item->name); + return false; + } + continue; + } + + // find the matching entry, if one exists (if not, this is not an error: we are not + // obliged to determine stats for all entries). we may have more than one entry for a + // given keyword (multiple stats for a single input item are allowed). + for (int i = 0; i < entries->n; i++) { + ppStatsEntry *entry = entries->data[i]; + if (strcmp (entry->keyword, item->name)) continue; + + // check if the types match? + if (entry->type != item->type) { + fprintf (stderr, "WARNING? mismatched type, skipping\n"); + continue; + } + + // save the value + + // if constant, save or compare with existing value + if (!strcasecmp (entry->statistic, "constant")) { + if (entry->value) { + // check that they match + } else { + if ((item->type == PS_DATA_STRING) && (item->data.str == NULL)) { + continue; + } + entry->value = psMemIncrRefCounter (item); + } + continue; + } + + bool useRMS = false; + if (!strcasecmp (entry->statistic, "rms")) { + useRMS = true; + } + + // only numerical values can have stats; all others must be 'constant' + if (entry->type >= PS_DATA_BOOL) { + psError(PS_ERR_UNKNOWN, false, "only numerical types can have non-constant stats: %s\n", entry->keyword); + return false; + } + + if (!entry->vector) { + entry->vector = psVectorAllocEmpty (16, entry->type); + } + + switch (item->type) { + case PS_DATA_U8: + if (useRMS) { + psVectorAppend (entry->vector, PS_SQR(item->data.U8)); + } else { + psVectorAppend (entry->vector, item->data.U8); + } + break; + case PS_DATA_U16: + if (useRMS) { + psVectorAppend (entry->vector, PS_SQR(item->data.U16)); + } else { + psVectorAppend (entry->vector, item->data.U16); + } + break; + case PS_DATA_U32: + if (useRMS) { + psVectorAppend (entry->vector, PS_SQR(item->data.U32)); + } else { + psVectorAppend (entry->vector, item->data.U32); + } + break; + case PS_DATA_U64: + if (useRMS) { + psVectorAppend (entry->vector, PS_SQR(item->data.U64)); + } else { + psVectorAppend (entry->vector, item->data.U64); + } + break; + + case PS_DATA_S8: + if (useRMS) { + psVectorAppend (entry->vector, PS_SQR(item->data.S8)); + } else { + psVectorAppend (entry->vector, item->data.S8); + } + break; + case PS_DATA_S16: + if (useRMS) { + psVectorAppend (entry->vector, PS_SQR(item->data.S16)); + } else { + psVectorAppend (entry->vector, item->data.S16); + } + break; + case PS_DATA_S32: + if (useRMS) { + psVectorAppend (entry->vector, PS_SQR(item->data.S32)); + } else { + psVectorAppend (entry->vector, item->data.S32); + } + break; + case PS_DATA_S64: + if (useRMS) { + psVectorAppend (entry->vector, PS_SQR(item->data.S64)); + } else { + psVectorAppend (entry->vector, item->data.S64); + } + break; + + case PS_DATA_F32: + if (!isfinite(item->data.F32)) continue; + if (useRMS) { + psVectorAppend (entry->vector, PS_SQR(item->data.F32)); + } else { + psVectorAppend (entry->vector, item->data.F32); + } + break; + case PS_DATA_F64: + if (!isfinite(item->data.F64)) continue; + if (useRMS) { + psVectorAppend (entry->vector, PS_SQR(item->data.F64)); + } else { + psVectorAppend (entry->vector, item->data.F64); + } + break; + default: + psError(PS_ERR_UNKNOWN, true, "ppStatsFromMetadata recipe ENTRY %s has invalid type\n", entry->keyword); + return false; + } + } + } + return true; +} Index: /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsFromMetadataPrint.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsFromMetadataPrint.c (revision 41082) +++ /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsFromMetadataPrint.c (revision 41082) @@ -0,0 +1,70 @@ +# include "ppStatsInternal.h" + +// calculate the stats for the non-constant entries (already calculated) +bool ppStatsFromMetadataPrint(psArray *entries, char *filename) +{ + bool status = true; // Status of printing + FILE *f = NULL; + if (!strcmp (filename, "-")) { + f = stdout; + } else { + f = fopen (filename, "w"); + if (f == NULL) { + psError(PS_ERR_UNKNOWN, false, "ppStatsFromMetadata cannot open output file %s\n", filename); + return false; + } + } + + // at this point, we have entries with values (of type STR or F32) or NULL + for (int i = 0; i < entries->n; i++) { + ppStatsEntry *entry = entries->data[i]; + + if (!entry->value) continue; + +// Print a value +#define VALUE_NUMERICAL_CASE(TYPE, FORMAT, NAME) \ + case PS_TYPE_##TYPE: \ + fprintf(f, "%s %" FORMAT " ", entry->flag, entry->value->data.NAME); \ + break; \ + + + switch (entry->value->type) { + VALUE_NUMERICAL_CASE(U8, PRIu8, U8); + VALUE_NUMERICAL_CASE(U16, PRIu16, U16); + VALUE_NUMERICAL_CASE(U32, PRIu32, U32); + VALUE_NUMERICAL_CASE(U64, PRIu64, U64); + VALUE_NUMERICAL_CASE(S8, PRId8, S8); + VALUE_NUMERICAL_CASE(S16, PRId16, S16); + VALUE_NUMERICAL_CASE(S32, PRId32, S32); + VALUE_NUMERICAL_CASE(S64, PRId64, S64); + VALUE_NUMERICAL_CASE(F32, "f", F32); + VALUE_NUMERICAL_CASE(F64, "lf", F64); + case PS_DATA_STRING: + if (entry->value->data.str) { + fprintf(f, "%s '%s' ", entry->flag, entry->value->data.str); + } + break; + case PS_DATA_BOOL: + if (entry->value->data.B) { + fprintf(f, "%s ", entry->flag); + } + break; + case PS_DATA_TIME: { + psTime *t = (psTime*)entry->value->data.V; + if (t) { + psString str = psTimeToISO(t); + fprintf(f, "%s %.19sZ ", entry->flag, str); + psFree(str); + } + break; + } + default: + psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unsupported type: %x", entry->value->type); + status = false; + } + } + fprintf(f, "\n"); + + if (f != stdout) fclose (f); + return status; +} Index: /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsFromMetadataStats.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsFromMetadataStats.c (revision 41082) +++ /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsFromMetadataStats.c (revision 41082) @@ -0,0 +1,101 @@ +# include "ppStatsInternal.h" + +static void computeBitwiseOr(ppStatsEntry *entry); + +// calculate the stats for the non-constant entries (already calculated) +bool ppStatsFromMetadataStats(psArray *entries) +{ + for (int i = 0; i < entries->n; i++) { + ppStatsEntry *entry = entries->data[i]; + + if (!strcasecmp (entry->statistic, "constant")) continue; + + // XXX skip or warn on missing stats? + if (!entry->vector) continue; + + if (!strcasecmp (entry->statistic, "bitwiseor")) { + computeBitwiseOr(entry); + continue; + } + + psStatsOptions option; + if (!strcasecmp (entry->statistic, "RMS")) { + option = psStatsOptionFromString ("SAMPLE_MEAN"); + goto got_stats; + } + if (!strcasecmp (entry->statistic, "SUM")) { + option = psStatsOptionFromString ("SAMPLE_MEAN"); + goto got_stats; + } + if (!strcasecmp (entry->statistic, "UQ")) { + option = psStatsOptionFromString ("ROBUST_QUARTILE"); + goto got_stats; + } + if (!strcasecmp (entry->statistic, "LQ")) { + option = psStatsOptionFromString ("ROBUST_QUARTILE"); + goto got_stats; + } + + option = psStatsOptionFromString (entry->statistic); + + got_stats: + if (entry->vector->n == 0) { + continue; + } + + psStats *stats = psStatsAlloc(option); + if (!psVectorStats(stats, entry->vector, NULL, NULL, 0)) { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats for %s", entry->statistic); + continue; + } + + double value; + if (!strcasecmp (entry->statistic, "RMS")) { + value = sqrt(stats->sampleMean); + goto got_value; + } + if (!strcasecmp (entry->statistic, "SUM")) { + value = stats->sampleMean * entry->vector->n; + goto got_value; + } + if (!strcasecmp (entry->statistic, "UQ")) { + value = stats->robustUQ; + goto got_value; + } + if (!strcasecmp (entry->statistic, "LQ")) { + value = stats->robustLQ; + goto got_value; + } + value = psStatsGetValue(stats, option); + + got_value: + entry->value = psMetadataItemAllocF32(entry->keyword, entry->statistic, value); + psFree(stats); + } + return true; +} + +static void computeBitwiseOr(ppStatsEntry *entry) +{ + psU64 result = 0; + psVector *vector = entry->vector; + for (int j = 0; j < vector->n; j++) { + // XXX: should we handle other types + if (entry->type == PS_DATA_U64) { + result |= vector->data.U64[j]; + } else if (entry->type == PS_DATA_U32) { + result |= vector->data.U32[j]; + } else { + return; + } + } + + if (entry->type == PS_DATA_U64) { + entry->value = psMetadataItemAllocU64(entry->keyword, entry->statistic, result); + } else { + entry->value = psMetadataItemAllocU32(entry->keyword, entry->statistic, result); + } + + return; +} + Index: /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsInternal.h =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsInternal.h (revision 41082) +++ /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsInternal.h (revision 41082) @@ -0,0 +1,17 @@ + +# ifdef HAVE_CONFIG_H +# include +# endif + +#ifndef PP_STATS_STAND_ALONE_H +#define PP_STATS_STAND_ALONE_H + +#include +#include +#include +#include +#include + +#include "ppStats.h" + +#endif Index: /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsLoop.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsLoop.c (revision 41082) +++ /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsLoop.c (revision 41082) @@ -0,0 +1,73 @@ +# include "ppStatsInternal.h" + +psMetadata *ppStatsLoop(psExit *result, + ppStatsData *data, // The data + pmConfig *config // Configuration + ) +{ + PS_ASSERT_PTR_NON_NULL(data, NULL); + pmFPA *fpa = data->fpa; // FPA to analyse + psFits *fits = data->fits; // FITS file handle + PS_ASSERT_PTR_NON_NULL(fpa, NULL); + + *result = PS_EXIT_SUCCESS; + + // allocate or bump the ref counter (so we can just free below) + pmFPAview *view = (data->view) ? psMemIncrRefCounter(data->view) : pmFPAviewAlloc(0); + + // allocate a new one if needed + psMetadata *newResults = psMetadataAlloc(); + + // check if we can legitimately iterate to the chip level + psArray *chips = fpa->chips; // Array of chips + if (view->chip >= chips->n) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Desired chip view (%d) doesn't match " + "number of chips (%ld)\n", view->chip, chips->n); + return NULL; + } + + // Iterate over chips (if view->chip is set, skip all others) + for (int i = 0; i < chips->n; i++) { + if ((view->chip >= 0) && (i != view->chip)) continue; + pmChip *chip = chips->data[i]; // Chip of interest + *result = ppStatsChip(newResults, chip, fits, view, data, config); + if (*result != PS_EXIT_SUCCESS) { + psError(PS_ERR_UNKNOWN, false, "trouble with stats for chip %d\n", i); + psFree (view); + psFree (newResults); + return NULL; + } + } + + // Iterate through the FPA -- Do this after the chip since we in the SPLIT cases, we + // only populate the FPA concepts after chips have been read. + if (psListLength(data->headers) > 0 && fpa->hdu) { + if (fits && !pmFPAReadHeader(fpa, fits, config)) { + psError(PS_ERR_IO, false, "Unable to read header for FPA."); + psFree(view); + psFree(newResults); + *result = PS_EXIT_DATA_ERROR; + return NULL; + } + pmHDU *hdu = fpa->hdu; // HDU for headers + p_ppStatsGetMetadata(newResults, hdu->header, data->headers); + } + if (psListLength(data->concepts) > 0) { + if (fits && fpa->hdu && !pmFPAReadHeader(fpa, fits, config)) { + psError(PS_ERR_IO, false, "Unable to read header for FPA."); + psFree(view); + psFree(newResults); + *result = PS_EXIT_DATA_ERROR; + return NULL; + } + pmConceptsReadFPA(fpa, PM_CONCEPT_SOURCE_ALL, false, config); + p_ppStatsGetMetadata(newResults, fpa->concepts, data->concepts); + } + + if (fits) { + pmFPAFreeData(fpa); + } + + psFree(view); + return newResults; +} Index: /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsMetadata.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsMetadata.c (revision 41082) +++ /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsMetadata.c (revision 41082) @@ -0,0 +1,63 @@ +#include "ppStatsInternal.h" + +// measure only the non-pixel-related statistics (headers, concepts) +psMetadata *ppStatsMetadata(psMetadata *out, + pmFPA *fpa, // FPA for which to get statistics + pmFPAview *view, // View for analysis + psImageMaskType maskVal, // Value to mask + pmConfig *config // Configuration + ) +{ + PS_ASSERT_PTR_NON_NULL(fpa, NULL) + PS_ASSERT_PTR_NON_NULL(view, NULL); + PS_ASSERT_PTR_NON_NULL(config, NULL); + + ppStatsData *data = ppStatsDataAlloc(); // All the input data + + // Get the options, open the files + if (!ppStatsSetupFromRecipe(data, config)) { + psError(PS_ERR_UNKNOWN, false, "Unable to get ppStats options from recipe."); + psFree(data); + return NULL; + } + + // drop the 'concept' and 'header' elements of data: + psFree (data->summary); + psFree (data->stats); + data->summary = psListAlloc(NULL); + data->stats = psStatsAlloc(0); + data->doStats = false; + + // Override recipe mask value + data->maskVal = maskVal; + + if (data->fpa) { + psFree(data->fpa); + } + data->fpa = psMemIncrRefCounter(fpa); + + if (data->view) { + psFree(data->view); + } + data->view = psMemIncrRefCounter(view); + + // Go through the FPA and do the hard work + psExit status; // Status of statistics loop + psMetadata *result = ppStatsLoop(&status, data, config); + if (status != PS_EXIT_SUCCESS) { + psError (PS_ERR_UNKNOWN, false, "Not able to measure FPA statistics.\n"); + psFree(result); + psFree(data); + return (NULL); + } + + if (out != NULL) { + psMetadataOverlay (out, result); + psFree(result); + psFree(data); + return out; + } + + psFree(data); + return result; +} Index: /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsPixels.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsPixels.c (revision 41082) +++ /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsPixels.c (revision 41082) @@ -0,0 +1,62 @@ +#include "ppStatsInternal.h" + +// measure only the pixel-related statistics (stats, summary) +psMetadata *ppStatsPixels(psMetadata *out, + pmFPA *fpa, // FPA for which to get statistics + pmFPAview *view, // View for analysis + psImageMaskType maskVal, // Value to mask + pmConfig *config // Configuration + ) +{ + PS_ASSERT_PTR_NON_NULL(fpa, NULL) + PS_ASSERT_PTR_NON_NULL(view, NULL); + PS_ASSERT_PTR_NON_NULL(config, NULL); + + ppStatsData *data = ppStatsDataAlloc(); // All the input data + + // Get the options, open the files + if (!ppStatsSetupFromRecipe(data, config)) { + psError(PS_ERR_UNKNOWN, false, "Unable to get ppStats options from recipe."); + psFree(data); + return NULL; + } + + // drop the 'concept' and 'header' elements of data: + psFree (data->headers); + psFree (data->concepts); + data->headers = psListAlloc(NULL); + data->concepts = psListAlloc(NULL); + + // Override recipe mask value + data->maskVal = maskVal; + + if (data->fpa) { + psFree(data->fpa); + } + data->fpa = psMemIncrRefCounter(fpa); + + if (data->view) { + psFree(data->view); + } + data->view = psMemIncrRefCounter(view); + + // Go through the FPA and do the hard work + psExit status; // Status of statistics loop + psMetadata *result = ppStatsLoop(&status, data, config); + if (status != PS_EXIT_SUCCESS) { + psError (PS_ERR_UNKNOWN, false, "Not able to measure FPA statistics.\n"); + psFree(result); + psFree(data); + return (NULL); + } + + if (out != NULL) { + psMetadataOverlay (out, result); + psFree(result); + psFree(data); + return out; + } + + psFree(data); + return result; +} Index: /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsReadout.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsReadout.c (revision 41082) +++ /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsReadout.c (revision 41082) @@ -0,0 +1,214 @@ +# include "ppStatsInternal.h" + +static bool warnNonFinite = false; // Have we warned about non-finite values? + +psExit ppStatsReadout(psMetadata *cellResults, // Metadata holding the chip results + pmReadout *readout, // Cell for which to get statistics + int nReadout, // readout number + ppStatsData *data, // The data + const pmConfig *config // Configuration + ) +{ + assert(cellResults); + assert(readout); + assert(data); + assert(config); + + /*** psphot and psastro put their results on the readout->analysis metadata (PSPHOT.HEADER, + PSASTRO.HEADER). we need to pull quantities of interest from those locations. to do + this, we need to select the appropriate readout. ***/ + + // Extract Header and Concept values from the Cell and Readout->analysis level + + psString readoutName = NULL; + psStringAppend (&readoutName, "READOUT.%02d", nReadout); + + // find or generate a metadata to hold the results + bool mdok; // Status of MD lookup + psMetadata *readoutResults = psMemIncrRefCounter(psMetadataLookupMetadata(&mdok, cellResults, readoutName)); + if (!mdok || !readoutResults) { + readoutResults = psMetadataAlloc(); + } + + // Extract Header values + if (psListLength(data->headers)) { + // extract from data->analysis output MD entries + if (psListLength(data->analysis)) { + p_ppStatsGetAnalysis (readoutResults, data->headers, readout->analysis, data->analysis); + } + } + + // Do we want to measure pixel statistics? + if (!data->doStats && !psListLength(data->summary)) { + // Nothing further to do --- don't want to waste our time reading the data + goto readoutDone; + } + + if (!readout->image) { + psLogMsg(__func__, PS_LOG_WARN, "No image associated with readout --- ignored.\n"); + goto readoutDone; + } + + // Measure basic image statistics (means, stdevs, etc). + if (data->sample <= 0.0) { + if (!psImageStats(data->stats, readout->image, readout->mask, data->maskVal)) { + psWarning("Unable to perform statistics on readout %s --- ignored.\n", readoutName); + goto statsDone; + } + } else { + // Apply sampling + psImage *image = readout->image; // The image of interest + psImage *mask = readout->mask; // The mask image + int numSamples = data->sample * image->numCols * image->numRows; // Number of samples + int sampleSpace = 1.0 / data->sample; // Space between samples + psVector *sampleValues = psVectorAlloc(numSamples, PS_TYPE_F32); // Vector of samples + psVector *sampleMask = psVectorAlloc(numSamples, PS_TYPE_VECTOR_MASK); // Corresponding mask + psVectorInit(sampleMask, 0); + + for (int i = 0; i < numSamples; i++) { + int j = i * sampleSpace; + int y = j / image->numCols; + int x = j % image->numCols; + sampleValues->data.F32[i] = image->data.F32[y][x]; + + // ignore the sampleMask if there is no input mask + if (!mask) continue; + + // if this pixel is masked, set the sample mask + if (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & data->maskVal) { + sampleMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 1; + continue; + } + + // mask any unmasked NAN/INF values + if (!isfinite(sampleValues->data.F32[i])) { + // warn for the first unmasked NAN/INF value + if (!warnNonFinite) { + psWarning("Unmasked non-finite value detected at %d,%d; suppressing further warnings", x, y); + warnNonFinite = true; + } + sampleMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 1; + } + } + if (!psVectorStats(data->stats, sampleValues, NULL, sampleMask, 1)) { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats for readout %s", readoutName); + } + // XXX raise if we get a NAN? psWarning("Unable to perform statistics on readout %s.\n", readoutName); + psFree(sampleValues); + psFree(sampleMask); + } + +#define WRITE_STAT(SYMBOL, NAME, SOURCE) \ + if (data->stats->options & SYMBOL) { \ + psMetadataAddF32(readoutResults, PS_LIST_TAIL, NAME, 0, NULL, data->stats->SOURCE); \ + } + + WRITE_STAT(PS_STAT_SAMPLE_MEAN, "SAMPLE_MEAN", sampleMean); + WRITE_STAT(PS_STAT_SAMPLE_MEDIAN, "SAMPLE_MEDIAN", sampleMedian); + WRITE_STAT(PS_STAT_SAMPLE_STDEV, "SAMPLE_STDEV", sampleStdev); + WRITE_STAT(PS_STAT_SAMPLE_SKEWNESS, "SAMPLE_SKEWNESS", sampleSkewness); + WRITE_STAT(PS_STAT_SAMPLE_KURTOSIS, "SAMPLE_KURTOSIS", sampleKurtosis); + WRITE_STAT(PS_STAT_SAMPLE_QUARTILE, "SAMPLE_LQ", sampleLQ); + WRITE_STAT(PS_STAT_SAMPLE_QUARTILE, "SAMPLE_UQ", sampleUQ); + WRITE_STAT(PS_STAT_ROBUST_MEDIAN, "ROBUST_MEDIAN", robustMedian); + WRITE_STAT(PS_STAT_ROBUST_STDEV, "ROBUST_STDEV", robustStdev); + WRITE_STAT(PS_STAT_ROBUST_QUARTILE, "ROBUST_LQ", robustLQ); + WRITE_STAT(PS_STAT_ROBUST_QUARTILE, "ROBUST_UQ", robustUQ); + WRITE_STAT(PS_STAT_ROBUST_QUARTILE, "ROBUST_N50", robustN50); + WRITE_STAT(PS_STAT_FITTED_MEAN, "FITTED_MEAN", fittedMean); + WRITE_STAT(PS_STAT_FITTED_STDEV, "FITTED_STDEV", fittedStdev); + WRITE_STAT(PS_STAT_CLIPPED_MEAN, "CLIPPED_MEAN", clippedMean); + WRITE_STAT(PS_STAT_CLIPPED_STDEV, "CLIPPED_STDEV", clippedStdev); + + // measure other types of statistics tests + +statsDone: + // count saturated pixels + if (psListLength(data->summary) > 0) { + bool get_nSatPixels = false; + bool get_fSatPixels = false; + bool findNumGood = false; // Return the number of good pixels? + bool findFracGood = false; // Return the fraction of good pixels? + + psListIterator *iterator = psListIteratorAlloc(data->summary, PS_LIST_HEAD, false); + psString choice; + while ((choice = psListGetAndIncrement(iterator))) { + if (!strcasecmp(choice, "SAT_PIXEL_NUM")) get_nSatPixels = true; + if (!strcasecmp(choice, "SAT_PIXEL_FRAC")) get_fSatPixels = true; + if (!strcasecmp(choice, "GOOD_PIXEL_NUM")) findNumGood = true; + if (!strcasecmp(choice, "GOOD_PIXEL_FRAC")) findFracGood = true; + } + + if (get_nSatPixels || get_fSatPixels) { + // Find the saturation point + float saturation = psMetadataLookupF32(&mdok, readout->parent->concepts, + "CELL.SATURATION"); // Saturation level + if (!mdok || isnan(saturation)) { + psWarning("CELL.SATURATION is not set --- unable to measure N_SAT_PIXELS.\n"); + if (get_nSatPixels) { + psMetadataAddS32(readoutResults, PS_LIST_TAIL, "SAT_PIXEL_NUM", 0, NULL, 0); + } + if (get_fSatPixels) { + psMetadataAddF32(readoutResults, PS_LIST_TAIL, "SAT_PIXEL_FRAC", 0, NULL, NAN); + } + } else { + int nSatPixels = 0; + for (int j = 0; j < readout->image->numRows; j++) { + for (int i = 0; i < readout->image->numCols; i++) { + if (readout->image->data.F32[j][i] >= saturation) { + nSatPixels ++; + } + } + } + if (get_nSatPixels) { + psMetadataAddS32(readoutResults, PS_LIST_TAIL, "SAT_PIXEL_NUM", 0, + "Number of saturated pixels", nSatPixels); + } + if (get_fSatPixels) { + psMetadataAddF32(readoutResults, PS_LIST_TAIL, "SAT_PIXEL_FRAC", 0, + "Fraction of saturated pixels", + nSatPixels / (float)(readout->image->numRows * readout->image->numCols)); + } + } + } + + if (findNumGood || findFracGood) { + if (!readout->mask || data->maskVal == 0) { + psWarning("Number or fraction of good pixels requested, but no mask provided"); + if (findNumGood) { + psMetadataAddS32(readoutResults, PS_LIST_TAIL, "GOOD_PIXEL_NUM", 0, NULL, 0); + } + if (findFracGood) { + psMetadataAddF32(readoutResults, PS_LIST_TAIL, "GOOD_PIXEL_FRAC", 0, NULL, NAN); + } + } else { + int numBad = 0; // Number of bad pixels + for (int j = 0; j < readout->mask->numRows; j++) { + for (int i = 0; i < readout->mask->numCols; i++) { + if (readout->mask->data.PS_TYPE_IMAGE_MASK_DATA[j][i] & data->maskVal) { + numBad++; + } + } + } + + int numTotal = readout->mask->numRows * readout->mask->numCols; // Total number of pixels + int numGood = numTotal - numBad; // Number of good pixels + if (findNumGood) { + psMetadataAddS32(readoutResults, PS_LIST_TAIL, "GOOD_PIXEL_NUM", 0, + "Number of good pixels", numGood); + } + if (findFracGood) { + psMetadataAddF32(readoutResults, PS_LIST_TAIL, "GOOD_PIXEL_FRAC", 0, + "Fraction of good pixels", numGood / (float)numTotal); + } + } + } + } + +readoutDone: + // Add the readout results to the cell + p_ppStatsAddToHierarchy(readoutResults, cellResults, readoutName, "Results for readout"); + psFree (readoutResults); + psFree (readoutName); + return PS_EXIT_SUCCESS; +} Index: /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsSetupFromArgs.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsSetupFromArgs.c (revision 41082) +++ /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsSetupFromArgs.c (revision 41082) @@ -0,0 +1,230 @@ +#include "ppStatsInternal.h" + +// This file is for setting up the required inputs from the command-line + +// Print usage information and die +static void usageAndDie(char *argv[], // Command-line arguments: only need the first which is always present + pmConfig *config // Configuration + ) +{ + printf("Return headers, concepts and/or image statistics.\n\n" + "Usage:\n" + "\t%s INPUT.fits [OUTPUT_NAME]\n" + "\n", argv[0]); + psArgumentHelp(config->arguments); + psFree(config); + psLibFinalize(); + pmConceptsDone(); + pmConfigDone(); + exit(EXIT_FAILURE); +} + +// Generate a list from the arguments +static void listFromArguments(psMetadata *arguments, // Arguments to parse + const char *name, // Name of the item, for error message + const char *flag, // The flag for the argument of interest + psList *target // The target list + ) +{ + psString regex = NULL; // Regular expression for the flag + psStringAppend(®ex, "^%s$", flag); + psMetadataIterator *iterator = psMetadataIteratorAlloc(arguments, PS_LIST_HEAD, regex); // Iterator + psFree(regex); + psMetadataItem *item; // Item from iteration + while ((item = psMetadataGetAndIncrement(iterator))) { + if (item->type != PS_DATA_STRING) { + psLogMsg(__func__, PS_LOG_WARN, "%s name is not of type STRING (%x) --- ignored.\n", + name, item->type); + continue; + } + if (item->data.V && strlen(item->data.V) > 0) { + psListAdd(target, PS_LIST_TAIL, item->data.V); + } + } + psFree(iterator); + + return; +} + +// Set the statistics option; for arguments +static inline void statsOptionArguments(psMetadata *arguments, // Arguments to parse + const char *name, // Name for option + ppStatsData *data, // Configuration data + psStatsOptions option // Option to check for + ) +{ + if (psMetadataLookupBool(NULL, arguments, name)) { + data->stats->options |= option; + data->doStats = true; + } + return; +} + +// Print out what we're going to do, from the list +static void checkList(psList *list, // List + const char *name // Name of list + ) +{ + psString value; + psListIterator *iterator = psListIteratorAlloc(list, PS_LIST_HEAD, false); + while ((value = psListGetAndIncrement(iterator))) { + printf("%s: %s\n", name, value); + } + psFree(iterator); + return; +} + + +ppStatsData *ppStatsSetupFromArgs(int *argc, char *argv[], // Command-line arguments + pmConfig *config // Configuration + ) +{ + // Setup and parse command-line arguments + psMetadata *arguments = config->arguments; // Arguments + psMetadataAddStr(arguments, PS_LIST_TAIL, "-chip", PS_META_DUPLICATE_OK, "Chip to inspect", NULL); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-cell", PS_META_DUPLICATE_OK, "Cell to inspect", NULL); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-header", PS_META_DUPLICATE_OK, "Header to look up", NULL); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-concept", PS_META_DUPLICATE_OK, "Concept to look up", NULL); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-summary", PS_META_DUPLICATE_OK, "Summary statistic to calculate", NULL); + psMetadataAddBool(arguments, PS_LIST_TAIL, "-mean", 0, "Calculate sample mean", false); + psMetadataAddBool(arguments, PS_LIST_TAIL, "-stdev", 0, "Calculate sample standard deviation", false); + psMetadataAddBool(arguments, PS_LIST_TAIL, "-median", 0, "Calculate sample median", false); + psMetadataAddBool(arguments, PS_LIST_TAIL, "-quartile", 0, "Calculate sample quartiles", false); + psMetadataAddBool(arguments, PS_LIST_TAIL, "-skewness", 0, "Calculate sample skewness", false); + psMetadataAddBool(arguments, PS_LIST_TAIL, "-kurtosis", 0, "Calculate sample kurtosis", false); + psMetadataAddBool(arguments, PS_LIST_TAIL, "-robust-median", 0, "Calculate robust median", false); + psMetadataAddBool(arguments, PS_LIST_TAIL, "-robust-stdev", 0, "Calculate robust standard deviation", false); + psMetadataAddBool(arguments, PS_LIST_TAIL, "-robust-quartile", 0, "Calculate robust quartile range", false); + psMetadataAddBool(arguments, PS_LIST_TAIL, "-fitted-mean", 0, "Calculate fitted mean", false); + psMetadataAddBool(arguments, PS_LIST_TAIL, "-fitted-stdev", 0, "Calculate fitted standard deviation", false); + psMetadataAddBool(arguments, PS_LIST_TAIL, "-clipped-mean", 0, "Calculate clipped median", false); + psMetadataAddBool(arguments, PS_LIST_TAIL, "-clipped-stdev", 0, "Calculate clipped standard deviation", false); + psMetadataAddS32(arguments, PS_LIST_TAIL, "-iter", 0, "Clipping iterations", 0); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-rej", 0, "Clipping level", 0.0); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-sample", 0, "Sampling fraction", 0.0); + psMetadataAddBool(arguments, PS_LIST_TAIL, "-level", 0, "File level", 0); + psMetadataAddBool(arguments, PS_LIST_TAIL, "-format", 0, "Show File format", 0); + psMetadataAddBool(arguments, PS_LIST_TAIL, "-show-camera", 0, "Show Camera name", 0); + + if (*argc == 1) { + // No command-line arguments: print the help + usageAndDie(argv, config); + } + if (!psArgumentParse(arguments, argc, argv) || + (*argc != 2 && *argc != 3)) { + printf("Unable to parse command-line arguments.\n\n"); + usageAndDie(argv, config); + } + + // Parse the command-line options + ppStatsData *data = ppStatsDataAlloc(); // The data + const char *inName = argv[1]; // Input file name + psArgumentRemove(1, argc, argv); + + listFromArguments(arguments, "Chip", "-chip", data->chips); + listFromArguments(arguments, "Cell", "-cell", data->cells); + listFromArguments(arguments, "Header", "-header", data->headers); + listFromArguments(arguments, "Concept", "-concept", data->concepts); + listFromArguments(arguments, "Summary", "-summary", data->summary); + + // Set the statistics options + statsOptionArguments(arguments, "-mean", data, PS_STAT_SAMPLE_MEAN); + statsOptionArguments(arguments, "-stdev", data, PS_STAT_SAMPLE_STDEV); + statsOptionArguments(arguments, "-median", data, PS_STAT_SAMPLE_MEDIAN); + statsOptionArguments(arguments, "-quartile", data, PS_STAT_SAMPLE_QUARTILE); + statsOptionArguments(arguments, "-skewness", data, PS_STAT_SAMPLE_SKEWNESS); + statsOptionArguments(arguments, "-kurtosis", data, PS_STAT_SAMPLE_KURTOSIS); + statsOptionArguments(arguments, "-robust-median", data, PS_STAT_ROBUST_MEDIAN); + statsOptionArguments(arguments, "-robust-stdev", data, PS_STAT_ROBUST_STDEV); + statsOptionArguments(arguments, "-robust-quartile", data, PS_STAT_ROBUST_QUARTILE); + statsOptionArguments(arguments, "-fitted-mean", data, PS_STAT_FITTED_MEAN); + statsOptionArguments(arguments, "-fitted-stdev", data, PS_STAT_FITTED_STDEV); + statsOptionArguments(arguments, "-clipped-mean", data, PS_STAT_CLIPPED_MEAN); + statsOptionArguments(arguments, "-clipped-stdev", data, PS_STAT_CLIPPED_STDEV); + data->stats->clipSigma = psMetadataLookupF32(NULL, arguments, "-rej"); + data->stats->clipIter = psMetadataLookupS32(NULL, arguments, "-iter"); + data->sample = psMetadataLookupF32(NULL, arguments, "-sample"); + data->fileLevel = psMetadataLookupBool(NULL, arguments, "-level"); + data->showFormat = psMetadataLookupBool(NULL, arguments, "-format"); + data->showCamera = psMetadataLookupBool(NULL, arguments, "-show-camera"); + data->fileView = NULL; + + // Open the input file, determine the camera + int result = PS_EXIT_UNKNOWN_ERROR; + { + psString resolved = pmConfigConvertFilename(inName, config, false, false); // Resolved filename + data->fits = psFitsOpen(resolved, "r"); + if (!data->fits) { + psError(PS_ERR_IO, false, "Unable to open input file %s\n", resolved); + psFree(resolved); + result = PS_EXIT_DATA_ERROR; + goto die; + } + psFree(resolved); + + psMetadata *header = psFitsReadHeader(NULL, data->fits); // The FITS (primary) header + psMetadata *format = pmConfigCameraFormatFromHeader(NULL, NULL, NULL, config, header, true); + if (!format) { + psError(PS_ERR_UNKNOWN, false, "Unable to determine camera format for %s\n", inName); + psFree(header); + result = PS_EXIT_DATA_ERROR; + goto die; + } + + // need to handle optional alternate EXTWORD value + psMetadata *fileMenu = psMetadataLookupMetadata (NULL, format, "FILE"); + if (!fileMenu) { + psError (PS_ERR_IO, true, "FILE METADATA missing from camera format\n"); + return false; + } + bool status = false; + char *extword = psMetadataLookupStr (&status, fileMenu, "EXTWORD"); + if (status) { + psFitsSetExtnameWord (data->fits, extword); + } + + data->fpa = pmFPAConstruct(config->camera, config->cameraName); + if (!data->fpa) { + psError(PS_ERR_UNKNOWN, false, "Unable to construct FPA for %s\n", resolved); + psFree(header); + psFree(format); + result = PS_EXIT_CONFIG_ERROR; + goto die; + } + data->fileView = pmFPAAddSourceFromHeader(data->fpa, header, format); + psFree(header); + psFree(format); + if (!data->fileView) { + psError(PS_ERR_UNKNOWN, false, "Unable to add input file %s to FPA.\n", inName); + result = PS_EXIT_CONFIG_ERROR; + goto die; + } + } + + // Get the rest from the recipe + if (!ppStatsSetupFromRecipe(data, config)) { + psError(PS_ERR_UNKNOWN, false, "Unable to read ppStats options from recipe."); + psFree(data); + return NULL; + } + + // Print out what we're going to do + if (psTraceGetLevel("ppStats") > 9) { + checkList(data->chips, "CHIP"); + checkList(data->cells, "CELL"); + checkList(data->headers, "HEADER"); + checkList(data->concepts, "CONCEPT"); + } + + return data; + + // Common path for error conditions: clean up and exit. +die: + psErrorStackPrint(stderr, "Unable to set ppStats parameters from command-line arguments"); + psFree(config); + psFree(data); + pmConceptsDone(); + pmConfigDone(); + psLibFinalize(); + exit(result); +} Index: /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsSetupFromRecipe.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsSetupFromRecipe.c (revision 41082) +++ /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsSetupFromRecipe.c (revision 41082) @@ -0,0 +1,145 @@ +#include "ppStatsInternal.h" + +// Strings in a recipe may be defined multiply (with MULTI) or listed on a single line +static void listFromRecipe(psList *target, // The target list + const psMetadata *recipe, // Recipe to search + const char *name // Name for item within recipe + ) +{ + // If the list already has entries, don't read anything else + if (psListLength(target) > 0) { + return; + } + + // First check that at least one item of interest exists + psMetadataItem *checkItem = psMetadataLookup(recipe, name); + if (!checkItem) { + // Nothing to see here + return; + } + + psString regex = NULL; // Regular expression for the flag + psStringAppend(®ex, "^%s$", name); + psMetadataIterator *iterator = psMetadataIteratorAlloc(recipe, PS_LIST_HEAD, regex); + psFree(regex); + psMetadataItem *item; + int numItem = 0; // Occurrence of the item in the recipe; to help the user in case of trouble + while ((item = psMetadataGetAndIncrement(iterator))) { + numItem++; + if (item->type != PS_DATA_STRING) { + psLogMsg(__func__, PS_LOG_WARN, "Occurrence %d of %s in the recipe is " + "not of type STRING (%x) --- ignored.\n", numItem, name, item->type); + continue; + } + // Parse into a list of independent values + psList *values = psStringSplit(item->data.V, " ,;", false); + // Copy into the target + psListIterator *valuesIter = psListIteratorAlloc(values, PS_LIST_HEAD, false); + psString valueString; + while ((valueString = psListGetAndIncrement(valuesIter))) { + // Remove the default _UNDEF entries in the recipe, which come from the limitation of the metadata + // config language to have an empty MULTI, combined with the constraint that higher-level recipes + // aren't allowed to add new values but only change what already exists. + if (strcmp(valueString, "_UNDEF") != 0) { + psListAdd(target, PS_LIST_TAIL, valueString); + } + } + psFree(valuesIter); + psFree(values); + } + psFree(iterator); + + return; +} + + +bool ppStatsSetupFromRecipe(ppStatsData *data, // Data for running ppStats + pmConfig *config // Configuration + ) +{ + PS_ASSERT_PTR_NON_NULL(data, false); + PS_ASSERT_PTR_NON_NULL(config, false); + + // Determine recipe parameters + bool mdok; // Status of MD lookup + psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSTATS_RECIPE); + if (!mdok || !recipe) { + psLogMsg(__func__, PS_LOG_WARN, "Unable to find recipe %s.\n", PPSTATS_RECIPE); + return data; + } + + listFromRecipe(data->chips, recipe, "CHIP"); + listFromRecipe(data->cells, recipe, "CELL"); + listFromRecipe(data->headers, recipe, "HEADER"); + listFromRecipe(data->concepts, recipe, "CONCEPT"); + listFromRecipe(data->analysis, recipe, "ANALYSIS"); + + // Parse SUMMARY statistics information + listFromRecipe(data->summary, recipe, "SUMMARY"); + + // Parse the statistics options + psList *recipeStats = psListAlloc(NULL); // List of statistics options + listFromRecipe(recipeStats, recipe, "STAT"); + + // validate STATs choices + if (psListLength(recipeStats) > 0) { + psListIterator *iterator = psListIteratorAlloc(recipeStats, PS_LIST_HEAD, false); + psString statString; // Statistic string, from iteration + while ((statString = psListGetAndIncrement(iterator))) { + psStatsOptions stat = psStatsOptionFromString(statString); + if (stat == 0) { + psLogMsg(__func__, PS_LOG_WARN, "Can't interpret STATS entry in recipe: " + "%s --- ignored.\n", statString); + continue; + } + data->stats->options |= stat; + data->doStats = true; + } + } + psFree(recipeStats); + + // Clipping options + if (data->stats->clipIter == 0 && isnan(data->stats->clipSigma)) { + int iter = psMetadataLookupS32(&mdok, recipe, "ITER"); // Number of clipping iterations + if (mdok && iter > 0) { + data->stats->clipIter = iter; + } else { + psLogMsg(__func__, PS_LOG_WARN, "ITER in recipe is not of type S32 and positive --- " + "retaining default.\n"); + } + float rej = psMetadataLookupF32(&mdok, recipe, "REJ"); // Clipping level + if (mdok && rej > 0) { + data->stats->clipSigma = rej; + } else { + psLogMsg(__func__, PS_LOG_WARN, "REJ in recipe is not of type F32 and positive --- " + "retaining default.\n"); + } + } + + if (data->sample == 0) { + float sample = psMetadataLookupF32(&mdok, recipe, "SAMPLE"); // Sample fraction + if (mdok && sample > 0) { + data->sample = sample; + } else { + psLogMsg(__func__, PS_LOG_WARN, "SAMPLE in recipe is not of type F32 and positive --- " + "retaining default.\n"); + } + } + + bool doFirst = psMetadataLookupBool(&mdok, recipe, "DO.FIRST.READOUT.3D"); // Sample fraction + if (mdok) { + data->doFirstReadout3D = doFirst; + } + + // set the mask value used for stand-alone analyses. + if (data->maskVal == 0) { + const char *names = psMetadataLookupStr(&mdok, recipe, "MASKVAL"); // Names for mask value + if (mdok) { + data->maskVal = pmConfigMaskGet(names, config); + } else { + psWarning("MASKVAL in recipe is not of type STR --- retaining default.\n"); + } + } + + return data; +} Index: /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsUtils.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsUtils.c (revision 41082) +++ /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsUtils.c (revision 41082) @@ -0,0 +1,83 @@ +# include "ppStatsInternal.h" + +void p_ppStatsGetMetadata(psMetadata *target, // Target for metadata + psMetadata *source, // Source for metadata + psList *list // List containing keywords + ) +{ + assert(target); + assert(list); + if (!source) { + // Nothing to get from! + return; + } + + psListIterator *iterator = psListIteratorAlloc(list, PS_LIST_HEAD, false); // Iterator + psString name; // Name from iteration + while ((name = psListGetAndIncrement(iterator))) { + psMetadataItem *item = psMetadataLookup(source, name); // Item of interest, or NULL + if (item) { + psMetadataAddItem(target, item, PS_LIST_TAIL, PS_META_REPLACE); + } + } + psFree(iterator); + return; +} + +void p_ppStatsGetAnalysis(psMetadata *target, // Output Target for metadata + psList *headers, // List containing desired keywords + psMetadata *source, // Input Source for metadata + psList *list // List containing analysis blocks + ) +{ + bool status; + + psListIterator *iterator = psListIteratorAlloc(list, PS_LIST_HEAD, false); // Iterator + psString name; // Name from iteration + while ((name = psListGetAndIncrement(iterator))) { + psMetadata *folder = psMetadataLookupMetadata(&status, source, name); // Item of interest, or NULL + if (folder) { + p_ppStatsGetMetadata (target, folder, headers); + } + } + psFree(iterator); + return; +} + +bool p_ppStatsDoThis(psList *toDoList, // List of things to do + const char *this // The name of "this" + ) +{ + if (psListLength(toDoList) == 0) { + // No list --- do everything + return true; + } + + psListIterator *iterator = psListIteratorAlloc(toDoList, PS_LIST_HEAD, false); // Iterator + psString test; // Test string, from iteration + while ((test = psListGetAndIncrement(iterator))) { + if (strcmp(this, test) == 0) { + // It's in the list --- do it + psFree(iterator); + return true; + } + } + psFree(iterator); + // Couldn't find it --- don't do it + return false; +} + +void p_ppStatsAddToHierarchy(psMetadata *source, // Source to add + psMetadata *target, // Target to which to add + const char *name, // Name of source + const char *comment // Comment for source + ) +{ + if (psListLength(source->list) > 0 && !psMetadataLookup(target, name)) { + psMetadataAdd(target, PS_LIST_TAIL, name, PS_DATA_METADATA | PS_META_REPLACE, + comment, source); + } + return; +} + + Index: /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsVersion.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsVersion.c (revision 41082) +++ /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsVersion.c (revision 41082) @@ -0,0 +1,63 @@ +#include "ppStatsInternal.h" +#include "ppStatsVersionDefinitions.h" + +#ifndef PPSTATS_VERSION +#error "PPSTATS_VERSION is not set" +#endif +#ifndef PPSTATS_BRANCH +#error "PPSTATS_BRANCH is not set" +#endif +#ifndef PPSTATS_SOURCE +#error "PPSTATS_SOURCE is not set" +#endif + +psString ppStatsVersion(void) +{ + char *value = NULL; + psStringAppend(&value, "%s@%s", PPSTATS_BRANCH, PPSTATS_VERSION); + return value; +} + +psString ppStatsSource(void) +{ + return psStringCopy(PPSTATS_SOURCE); +} + +psString ppStatsVersionLong(void) +{ + psString version = ppStatsVersion(); // Version, to return + psString source = ppStatsSource(); // Source + + psStringPrepend(&version, "ppStats "); + psStringAppend(&version, " from %s, built %s, %s", source, __DATE__, __TIME__); + psFree(source); + +#ifdef __OPTIMIZE__ + psStringAppend(&version, " optimised"); +#else + psStringAppend(&version, " unoptimised"); +#endif + + return version; +}; + +bool ppStatsVersionHeader(psMetadata *header) +{ + PS_ASSERT_METADATA_NON_NULL(header, false); + + psString version = ppStatsVersion(); // Software version + psString source = ppStatsSource(); // Software source + + psMetadataAddStr(header, PS_LIST_TAIL, "STATS_V", PS_META_REPLACE, NULL, PPSTATS_VERSION); + + psStringPrepend(&version, "ppStats version: "); + psStringPrepend(&source, "ppStats source: "); + + psMetadataAddStr(header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, NULL, version); + psMetadataAddStr(header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, NULL, source); + + psFree(version); + psFree(source); + + return true; +} Index: /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsVersionDefinitions.h.in =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsVersionDefinitions.h.in (revision 41082) +++ /branches/ccl_branches/ipponly-20191108/ppStats/src/ppStatsVersionDefinitions.h.in (revision 41082) @@ -0,0 +1,8 @@ +#ifndef PPSTATS_VERSION_DEFINITIONS_H +#define PPSTATS_VERSION_DEFINITIONS_H + +#define PPSTATS_VERSION @PPSTATS_VERSION@ // SVN version +#define PPSTATS_BRANCH @PPSTATS_BRANCH@ // SVN branch +#define PPSTATS_SOURCE @PPSTATS_SOURCE@ // SVN source + +#endif