Index: branches/eam_branches/20090522/psLib/configure.ac
===================================================================
--- branches/eam_branches/20090522/psLib/configure.ac	(revision 24238)
+++ branches/eam_branches/20090522/psLib/configure.ac	(revision 24557)
@@ -392,4 +392,5 @@
   test/sys/Makefile
   test/types/Makefile
+  test/optime/Makefile
   utils/Makefile
 ])
Index: branches/eam_branches/20090522/psLib/src/fits/psFits.c
===================================================================
--- branches/eam_branches/20090522/psLib/src/fits/psFits.c	(revision 24238)
+++ branches/eam_branches/20090522/psLib/src/fits/psFits.c	(revision 24557)
@@ -931,9 +931,11 @@
     }
 
+    // use strncmp so that we can resolve the string valuves that cfitsio puts into headers (RICE_1 GZIP_1)
+    // into psFitsCompressionType
     if (strcmp(string, "NONE") == 0) return PS_FITS_COMPRESS_NONE;
-    if (strcmp(string, "GZIP") == 0) return PS_FITS_COMPRESS_GZIP;
-    if (strcmp(string, "RICE") == 0) return PS_FITS_COMPRESS_RICE;
-    if (strcmp(string, "HCOMPRESS") == 0) return PS_FITS_COMPRESS_HCOMPRESS;
-    if (strcmp(string, "PLIO") == 0) return PS_FITS_COMPRESS_PLIO;
+    if (strncmp(string, "GZIP", 4) == 0) return PS_FITS_COMPRESS_GZIP;
+    if (strncmp(string, "RICE", 4) == 0) return PS_FITS_COMPRESS_RICE;
+    if (strncmp(string, "HCOMPRESS", 9) == 0) return PS_FITS_COMPRESS_HCOMPRESS;
+    if (strncmp(string, "PLIO", 4) == 0) return PS_FITS_COMPRESS_PLIO;
 
     psWarning("Unable to identify compression type (%s) --- none set.", string);
Index: branches/eam_branches/20090522/psLib/src/fits/psFitsTable.c
===================================================================
--- branches/eam_branches/20090522/psLib/src/fits/psFitsTable.c	(revision 24238)
+++ branches/eam_branches/20090522/psLib/src/fits/psFitsTable.c	(revision 24557)
@@ -390,9 +390,11 @@
 
 
-bool psFitsInsertTable(psFits* fits,
-                       const psMetadata* header,
-                       const psArray* table,
-                       const char *extname,
-                       bool after)
+static bool fitsInsertTable(psFits* fits,             // FITS file
+                            const psMetadata* header, // FITS header to write
+                            const psArray* table,     // Table to write
+                            const char *extname,      // Extension name to give table
+                            bool after,               // Write table after current extension?
+                            bool writeData            // Write data?
+                            )
 {
     PS_ASSERT_FITS_NON_NULL(fits, false);
@@ -403,5 +405,5 @@
 
     long numRows = table->n;
-    if (numRows < 1) {
+    if (writeData && numRows < 1) {
         // no table data, what can I do?
         psError(PS_ERR_BAD_PARAMETER_SIZE, true,
@@ -504,5 +506,5 @@
         fits_create_tbl(fits->fd,
                         BINARY_TBL,
-                        table->n, // number of rows in table
+                        writeData ? numRows : 0, // number of rows in table
                         numColumns, // number of columns in table
                         (char**)columnNames->data, // names of the columns
@@ -524,5 +526,5 @@
         // Insert the table
         fits_insert_btbl(fits->fd,
-                         table->n, // number of rows in table
+                         writeData ? numRows : 0, // number of rows in table
                          numColumns, // number of columns in table
                          (char**)columnNames->data, // names of the columns
@@ -562,76 +564,78 @@
 
     // cfitsio requires that we write the data by columns --- urgh!
-    psMetadataIteratorSet(colSpecsIter, PS_LIST_HEAD);
-    for (long colNum = 1; (colSpecItem = psMetadataGetAndIncrement(colSpecsIter)); colNum++) {
-        // Note: colNum is unit-indexed, because it's for cfitsio
-        colSpec *spec = colSpecItem->data.V; // The specification
-        if (PS_DATA_IS_PRIMITIVE(spec->type)) {
-            size_t dataSize = PSELEMTYPE_SIZEOF(spec->type); // Size (in bytes) of this type
-            psVector *columnData = psVectorAlloc(table->n, spec->type); // The raw row data, to be written
-            psVectorInit(columnData, 0);
-            for (long i = 0; i < table->n; i++) {
-                psMetadata *row = table->data[i]; // The row of interest
-                psMetadataItem *dataItem = psMetadataLookup(row, colSpecItem->name); // The value of interest
-                memcpy(&columnData->data.U8[i * dataSize], &dataItem->data, dataSize);
+    if (writeData) {
+        psMetadataIteratorSet(colSpecsIter, PS_LIST_HEAD);
+        for (long colNum = 1; (colSpecItem = psMetadataGetAndIncrement(colSpecsIter)); colNum++) {
+            // Note: colNum is unit-indexed, because it's for cfitsio
+            colSpec *spec = colSpecItem->data.V; // The specification
+            if (PS_DATA_IS_PRIMITIVE(spec->type)) {
+                size_t dataSize = PSELEMTYPE_SIZEOF(spec->type); // Size (in bytes) of this type
+                psVector *columnData = psVectorAlloc(table->n, spec->type); // The raw row data, to be written
+                psVectorInit(columnData, 0);
+                for (long i = 0; i < table->n; i++) {
+                    psMetadata *row = table->data[i]; // The row of interest
+                    psMetadataItem *dataItem = psMetadataLookup(row, colSpecItem->name); // Value of interest
+                    memcpy(&columnData->data.U8[i * dataSize], &dataItem->data, dataSize);
+                }
+
+                int fitsDataType;           // Data type for cfitsio
+                p_psFitsTypeToCfitsio(spec->type, NULL, NULL, &fitsDataType);
+                fits_write_col(fits->fd,
+                               fitsDataType,
+                               colNum, // column number
+                               1, // first row
+                               1, // first element
+                               table->n, // number of rows
+                               columnData->data.U8, // the data
+                               &status);
+                psFree(columnData);
+            } else {
+                switch (spec->type) {
+                  case PS_DATA_STRING: {
+                      psArray *strings = psArrayAlloc(table->n); // Array of strings
+                      for (long i = 0; i < table->n; i++) {
+                          psMetadata *row = table->data[i]; // The row of interest
+                          strings->data[i] = psMemIncrRefCounter(psMetadataLookupStr(NULL, row,
+                                                                                     colSpecItem->name));
+                      }
+                      fits_write_col_str(fits->fd, colNum, 1, 1, table->n, (char**)strings->data, &status);
+                      psFree(strings);
+                      break;
+                  }
+                  case PS_DATA_VECTOR: {
+                      size_t dataSize = PSELEMTYPE_SIZEOF(spec->vectorType); // Size of data, in bytes
+                      psVector *columnData = psVectorAlloc(spec->size * table->n * dataSize, PS_TYPE_U8);
+                      psVectorInit(columnData, 0);
+                      for (long i = 0; i < table->n; i++) {
+                          psMetadata *row = table->data[i]; // The row of interest
+                          psMetadataItem* dataItem = psMetadataLookup(row, colSpecItem->name);
+                          if (dataItem->type != PS_DATA_VECTOR) {
+                              // Just in case --- get a zero instead of some weird result
+                              continue;
+                          }
+                          psVector *vector = dataItem->data.V;
+                          memcpy(&columnData->data.U8[i * dataSize * spec->size], vector->data.U8,
+                                 vector->n * dataSize);
+                      }
+
+                      int fitsDataType;           // Data type for cfitsio
+                      p_psFitsTypeToCfitsio(spec->vectorType, NULL, NULL, &fitsDataType);
+                      fits_write_col(fits->fd, fitsDataType, colNum, 1, 1, table->n * spec->size,
+                                     columnData->data.U8, &status);
+                      psFree(columnData);
+                      break;
+                  }
+                  default:
+                    psAbort("Should never get here.\n");
+                }
             }
 
-            int fitsDataType;           // Data type for cfitsio
-            p_psFitsTypeToCfitsio(spec->type, NULL, NULL, &fitsDataType);
-            fits_write_col(fits->fd,
-                           fitsDataType,
-                           colNum, // column number
-                           1, // first row
-                           1, // first element
-                           table->n, // number of rows
-                           columnData->data.U8, // the data
-                           &status);
-            psFree(columnData);
-        } else {
-            switch (spec->type) {
-            case PS_DATA_STRING: {
-                    psArray *strings = psArrayAlloc(table->n); // Array of strings
-                    for (long i = 0; i < table->n; i++) {
-                        psMetadata *row = table->data[i]; // The row of interest
-                        strings->data[i] = psMemIncrRefCounter(psMetadataLookupStr(NULL, row,
-                                                               colSpecItem->name));
-                    }
-                    fits_write_col_str(fits->fd, colNum, 1, 1, table->n, (char**)strings->data, &status);
-                    psFree(strings);
-                    break;
-                }
-            case PS_DATA_VECTOR: {
-                    size_t dataSize = PSELEMTYPE_SIZEOF(spec->vectorType); // Size of data, in bytes
-                    psVector *columnData = psVectorAlloc(spec->size * table->n * dataSize, PS_TYPE_U8);
-                    psVectorInit(columnData, 0);
-                    for (long i = 0; i < table->n; i++) {
-                        psMetadata *row = table->data[i]; // The row of interest
-                        psMetadataItem* dataItem = psMetadataLookup(row, colSpecItem->name);
-                        if (dataItem->type != PS_DATA_VECTOR) {
-                            // Just in case --- get a zero instead of some weird result
-                            continue;
-                        }
-                        psVector *vector = dataItem->data.V;
-                        memcpy(&columnData->data.U8[i * dataSize * spec->size], vector->data.U8,
-                               vector->n * dataSize);
-                    }
-
-                    int fitsDataType;           // Data type for cfitsio
-                    p_psFitsTypeToCfitsio(spec->vectorType, NULL, NULL, &fitsDataType);
-                    fits_write_col(fits->fd, fitsDataType, colNum, 1, 1, table->n * spec->size,
-                                   columnData->data.U8, &status);
-                    psFree(columnData);
-                    break;
-                }
-            default:
-                psAbort("Should never get here.\n");
+            // Check error status from writing column
+            if (status != 0) {
+                psFitsError(status, true, "Unable to write column %ld of FITS table", colNum);
+                psFree(colSpecsIter);
+                psFree(colSpecs);
+                return false;
             }
-        }
-
-        // Check error status from writing column
-        if (status != 0) {
-            psFitsError(status, true, "Unable to write column %ld of FITS table", colNum);
-            psFree(colSpecsIter);
-            psFree(colSpecs);
-            return false;
         }
     }
@@ -650,4 +654,42 @@
     return true;
 }
+
+
+bool psFitsInsertTable(psFits* fits, const psMetadata* header, const psArray* table, const char *extname,
+                       bool after)
+{
+    PS_ASSERT_FITS_NON_NULL(fits, false);
+    PS_ASSERT_FITS_WRITABLE(fits, false);
+    return fitsInsertTable(fits, header, table, extname, after, true);
+}
+
+bool psFitsWriteTableEmpty(psFits *fits, const psMetadata *header, const psMetadata *columns,
+                           const char *extname)
+{
+    PS_ASSERT_FITS_NON_NULL(fits, false);
+    PS_ASSERT_FITS_WRITABLE(fits, false);
+    if (!psFitsMoveLast(fits)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to move to last extension to write table");
+        return false;
+    }
+    psArray *table = psArrayAlloc(1);   // Dummy table carrying column definitions
+    table->data[0] = psMemIncrRefCounter((psPtr)columns); // Casting away const
+    bool status = fitsInsertTable(fits, header, table, extname, true, false); // Status of insertion
+    psFree(table);
+    return status;
+}
+
+bool psFitsInsertTableEmpty(psFits *fits, const psMetadata *header, const psMetadata *columns,
+                            const char *extname, bool after)
+{
+    PS_ASSERT_FITS_NON_NULL(fits, false);
+    PS_ASSERT_FITS_WRITABLE(fits, false);
+    psArray *table = psArrayAlloc(1);   // Dummy table carrying column definitions
+    table->data[0] = psMemIncrRefCounter((psPtr)columns); // Casting away const
+    bool status = fitsInsertTable(fits, header, table, extname, after, false); // Status of insertion
+    psFree(table);
+    return status;
+}
+
 
 bool psFitsUpdateTable(psFits* fits,
Index: branches/eam_branches/20090522/psLib/src/fits/psFitsTable.h
===================================================================
--- branches/eam_branches/20090522/psLib/src/fits/psFitsTable.h	(revision 24238)
+++ branches/eam_branches/20090522/psLib/src/fits/psFitsTable.h	(revision 24557)
@@ -90,4 +90,12 @@
 );
 
+/// Write an empty table
+bool psFitsWriteTableEmpty(
+    psFits *fits,                       ///< FITS file pointer
+    const psMetadata *header,           ///< Header to write
+    const psMetadata *columns,          ///< Column definitions; no data used except name,type
+    const char *extname                 ///< Extension name for table
+    );
+
 /** Inserts a whole FITS table. A new HDU of the type BINTABLE is inserted either
  *  before or after, depending on the AFTER parameter, the current HDU.
@@ -104,4 +112,14 @@
     bool after    ///< TRUE if insert is done after CHDU, otherwise table is inserted before CHDU
 );
+
+/// Insert an empty table
+bool psFitsInsertTableEmpty(
+    psFits *fits,              ///< FITS file pointer
+    const psMetadata *header,  ///< Header to write
+    const psMetadata *columns, ///< Column definitions; no data used except name,type
+    const char *extname,       ///< Extension name for table
+    bool after                 ///< Insert after current HDU?
+    );
+
 
 /** Updates a FITS table.  The current HDU type must be either
Index: branches/eam_branches/20090522/psLib/src/imageops/psImageBackground.c
===================================================================
--- branches/eam_branches/20090522/psLib/src/imageops/psImageBackground.c	(revision 24238)
+++ branches/eam_branches/20090522/psLib/src/imageops/psImageBackground.c	(revision 24557)
@@ -52,24 +52,41 @@
 
     // Minimum and maximum values
-    float min = values->data.F32[0];
-    float max = values->data.F32[0];
+    float min = +PS_MAX_F32;
+    float max = -PS_MAX_F32;
 
     // select a subset of the image pixels to measure the stats
     long n = 0;                         // Number of actual pixels in subset
-    for (long i = 0; i < Nsubset; i++) {
-        double frnd = psRandomUniform(rng);
-        int pixel = Npixels * frnd;
-        int ix = pixel % nx;
-        int iy = pixel / nx;
+    if (Nsubset >= Npixels) {
+	// if we have an image smaller than Nsubset, just loop over the image pixels
+	for (int iy = 0; iy < ny; iy++) {
+	    for (int ix = 0; ix < nx; ix++) {
+		if (!isfinite(image->data.F32[iy][ix]) || (mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskValue)) {
+		    continue;
+		}
 
-        if (!isfinite(image->data.F32[iy][ix]) || (mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskValue)) {
-            continue;
-        }
+		float value = image->data.F32[iy][ix];
+		min = PS_MIN(value, min);
+		max = PS_MAX(value, max);
+		values->data.F32[n] = value;
+		n++;
+	    }
+	}
+    } else {
+	for (long i = 0; i < Nsubset; i++) {
+	    double frnd = psRandomUniform(rng);
+	    int pixel = Npixels * frnd;
+	    int ix = pixel % nx;
+	    int iy = pixel / nx;
 
-        float value = image->data.F32[iy][ix];
-        min = PS_MIN(value, min);
-        max = PS_MIN(value, max);
-        values->data.F32[n] = value;
-        n++;
+	    if (!isfinite(image->data.F32[iy][ix]) || (mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskValue)) {
+		continue;
+	    }
+
+	    float value = image->data.F32[iy][ix];
+	    min = PS_MIN(value, min);
+	    max = PS_MAX(value, max);
+	    values->data.F32[n] = value;
+	    n++;
+	}
     }
     if (n < 0.01*Nsubset) {
Index: branches/eam_branches/20090522/psLib/src/imageops/psImageConvolve.c
===================================================================
--- branches/eam_branches/20090522/psLib/src/imageops/psImageConvolve.c	(revision 24238)
+++ branches/eam_branches/20090522/psLib/src/imageops/psImageConvolve.c	(revision 24557)
@@ -122,4 +122,22 @@
 
     return kernel;
+}
+
+psKernel *psKernelCopy(const psKernel *in)
+{
+    PS_ASSERT_KERNEL_NON_NULL(in, NULL);
+
+    psKernel *out = psAlloc(sizeof(psKernel)); // The copied kernel, to be returned
+    psMemSetDeallocator(out,(psFreeFunc)kernelFree);
+
+    out->image = psImageCopy(NULL, in->image, PS_TYPE_KERNEL);
+    out->xMin = in->xMin;
+    out->xMax = in->xMax;
+    out->yMin = in->yMin;
+    out->yMax = in->yMax;
+
+    kernelRedirects(out, out->image->numRows);
+
+    return out;
 }
 
Index: branches/eam_branches/20090522/psLib/src/imageops/psImageConvolve.h
===================================================================
--- branches/eam_branches/20090522/psLib/src/imageops/psImageConvolve.h	(revision 24238)
+++ branches/eam_branches/20090522/psLib/src/imageops/psImageConvolve.h	(revision 24557)
@@ -96,4 +96,11 @@
     );
 
+/// Copy a kernel
+///
+/// Performs a deep copy of the input kernel
+psKernel *psKernelCopy(
+    const psKernel *in                  ///< Kernel to be copied
+    );
+
 /// Checks the type of a particular pointer.
 ///
Index: branches/eam_branches/20090522/psLib/src/imageops/psImageUnbin.c
===================================================================
--- branches/eam_branches/20090522/psLib/src/imageops/psImageUnbin.c	(revision 24238)
+++ branches/eam_branches/20090522/psLib/src/imageops/psImageUnbin.c	(revision 24557)
@@ -30,5 +30,5 @@
 // DX, DY are the binning factor
 // dx, dy is the distance is high-res pixels to the 0,0 corner of the first
-// binned pixel. 
+// binned pixel.
 // XXX check that this is still consistent with psphotImageMedian...
 psImage *psImageUnbin(psImage *out, const psImage *in, const psImageBinning *binning)
@@ -69,5 +69,5 @@
             // (Xs,Ys) : (Xe,Ye) : binned pixel centers in unbinned coords
             // corresponding to (Ix,Iy), (Ix+1,Iy+1)
-	    // XXX should this be "+ dx" and + dy?
+            // XXX should this be "+ dx" and + dy?
             int Xs = PS_MAX (0, PS_MIN (Nx, psImageBinningGetFineX(binning, Ix + 0.5)));
             int Ys = PS_MAX (0, PS_MIN (Ny, psImageBinningGetFineY(binning, Iy + 0.5)));
@@ -76,28 +76,28 @@
 
             for (int iy = Ys; (iy < Ye) && (iy < Ny); iy++) {
-		float dY = (iy - Ys) / (float) DY;
-		float rY = 1.0 - dY;
+                float dY = (iy - Ys) / (float) DY;
+                float rY = 1.0 - dY;
                 float Vxs = V10*dY + V00*rY;
                 float Vxe = V11*dY + V01*rY;
 
-                // Vxs = (V10 - V00)*(iy - Ys) / DY + V00;                                                                                                       
-                // Vxe = (V11 - V01)*(iy - Ys) / DY + V01;                                                                                                       
-
-		// dVx = Vxs_1 - Vxs_1
-		// dVx = (V10*dY_1 + V00*rY_1) - (V10*dY_0 + V00*rY_0);
-		// dY_0 = (iy - Ys)/DY     = iy/DY - Ys/DY;
-		// dY_1 = (iy + 1 - Ys)/DY = iy/DY - Ys/DY + 1/DY;
-		// rY_0 = 1 - dY_0;
-		// rY_1 = 1 - dY_1;
-		// dVx = V10*(ddY) + V00*(drY);
-		// ddY = 1/DY;
-		// drY = -1/DY;
-
-		// dVxs = (V10 - V00)/DY;
-		// dVxe = (V11 - V01)/DY;
-
-		// ddV = (Vxe_1 - Vxs_1)/DX - (Vxe_0 - Vxs_0)/DX;
-		// ddV = (Vxe_1 - Vxe_0)/DX - (Vxs_1 - Vxs_0)/DX;
-		// ddV = (V11 - V01 - V10 + V00)/(DX*DY);
+                // Vxs = (V10 - V00)*(iy - Ys) / DY + V00;
+                // Vxe = (V11 - V01)*(iy - Ys) / DY + V01;
+
+                // dVx = Vxs_1 - Vxs_1
+                // dVx = (V10*dY_1 + V00*rY_1) - (V10*dY_0 + V00*rY_0);
+                // dY_0 = (iy - Ys)/DY     = iy/DY - Ys/DY;
+                // dY_1 = (iy + 1 - Ys)/DY = iy/DY - Ys/DY + 1/DY;
+                // rY_0 = 1 - dY_0;
+                // rY_1 = 1 - dY_1;
+                // dVx = V10*(ddY) + V00*(drY);
+                // ddY = 1/DY;
+                // drY = -1/DY;
+
+                // dVxs = (V10 - V00)/DY;
+                // dVxe = (V11 - V01)/DY;
+
+                // ddV = (Vxe_1 - Vxs_1)/DX - (Vxe_0 - Vxs_0)/DX;
+                // ddV = (Vxe_1 - Vxe_0)/DX - (Vxs_1 - Vxs_0)/DX;
+                // ddV = (V11 - V01 - V10 + V00)/(DX*DY);
 
                 float dV = (Vxe - Vxs) / DX;
@@ -128,5 +128,5 @@
             for (int ix = 0; ix < Xs; ix++) {
                 vOut[iy][ix] = V;
-		// assert (fabs(V - psImageUnbinPixel(ix, iy, in, binning)) < UNBIN_TOL*fabs(V));
+                // assert (fabs(V - psImageUnbinPixel(ix, iy, in, binning)) < UNBIN_TOL*fabs(V));
             }
             V += dV;
@@ -141,5 +141,5 @@
             for (int ix = Xe; ix < Nx; ix++) {
                 vOut[iy][ix] = V;
-		// assert (fabs(V - psImageUnbinPixel(ix, iy, in, binning)) < UNBIN_TOL*fabs(V));
+                // assert (fabs(V - psImageUnbinPixel(ix, iy, in, binning)) < UNBIN_TOL*fabs(V));
             }
             V += dV;
@@ -163,5 +163,5 @@
             for (int ix = Xs; (ix < Xe) && (ix < Nx); ix++) {
                 vOut[iy][ix] = V;
-		// assert (fabs(V - psImageUnbinPixel(ix, iy, in, binning)) < UNBIN_TOL*fabs(V));
+                // assert (fabs(V - psImageUnbinPixel(ix, iy, in, binning)) < UNBIN_TOL*fabs(V));
                 V += dV;
             }
@@ -176,5 +176,5 @@
             for (int ix = Xs; (ix < Xe) && (ix < Nx); ix++) {
                 vOut[iy][ix] = V;
-		// assert (fabs(V - psImageUnbinPixel(ix, iy, in, binning)) < UNBIN_TOL*fabs(V));
+                // assert (fabs(V - psImageUnbinPixel(ix, iy, in, binning)) < UNBIN_TOL*fabs(V));
                 V += dV;
             }
@@ -186,13 +186,13 @@
     {
         float V;
-	// center of last pixel 
-	int Xs = PS_MAX (0, PS_MIN (Nx, psImageBinningGetFineX(binning, 0  + 0.5)));
-	int Xe = PS_MAX (0, PS_MIN (Nx, psImageBinningGetFineX(binning, nx - 0.5)));
-	int Ys = PS_MAX (0, PS_MIN (Ny, psImageBinningGetFineY(binning, 0  + 0.5)));
-	int Ye = PS_MAX (0, PS_MIN (Ny, psImageBinningGetFineY(binning, ny - 0.5)));
+        // center of last pixel
+        int Xs = PS_MAX (0, PS_MIN (Nx, psImageBinningGetFineX(binning, 0  + 0.5)));
+        int Xe = PS_MAX (0, PS_MIN (Nx, psImageBinningGetFineX(binning, nx - 0.5)));
+        int Ys = PS_MAX (0, PS_MIN (Ny, psImageBinningGetFineY(binning, 0  + 0.5)));
+        int Ye = PS_MAX (0, PS_MIN (Ny, psImageBinningGetFineY(binning, ny - 0.5)));
 
         // 0,0
         V = vIn[0][0];
-	// assert (fabs(V - psImageUnbinPixel(0, 0, in, binning)) < UNBIN_TOL*fabs(V));
+        // assert (fabs(V - psImageUnbinPixel(0, 0, in, binning)) < UNBIN_TOL*fabs(V));
 
         for (int iy = 0; iy < Ys; iy++)
@@ -204,5 +204,5 @@
         // Nx,0
         V = vIn[0][nx-1];
-	// assert (fabs(V - psImageUnbinPixel(Nx-1, 0, in, binning)) < UNBIN_TOL*fabs(V));
+        // assert (fabs(V - psImageUnbinPixel(Nx-1, 0, in, binning)) < UNBIN_TOL*fabs(V));
 
         for (int iy = 0; iy < Ys; iy++)
@@ -214,5 +214,5 @@
         // 0,Ny
         V = vIn[ny-1][0];
-	// assert (fabs(V - psImageUnbinPixel(0, Ny-1, in, binning)) < UNBIN_TOL*fabs(V));
+        // assert (fabs(V - psImageUnbinPixel(0, Ny-1, in, binning)) < UNBIN_TOL*fabs(V));
 
         for (int iy = Ye; iy < Ny; iy++)
@@ -224,5 +224,5 @@
         // Nx,Ny
         V = vIn[ny-1][nx-1];
-	// assert (fabs(V - psImageUnbinPixel(Nx-1, Ny-1, in, binning)) < UNBIN_TOL*fabs(V));
+        // assert (fabs(V - psImageUnbinPixel(Nx-1, Ny-1, in, binning)) < UNBIN_TOL*fabs(V));
 
         for (int iy = Ye; iy < Ny; iy++)
@@ -243,6 +243,6 @@
 
 double psImageUnbinPixel(const double xFine, const double yFine, // desired Unbinned point (parent coords)
-			 const psImage *in, // binned image
-			 const psImageBinning *binning)   //!< Overhang
+                         const psImage *in, // binned image
+                         const psImageBinning *binning)   //!< Overhang
 {
     PS_ASSERT_IMAGE_NON_NULL(in, NAN);
@@ -273,5 +273,5 @@
 
     if ((x < -nXedge) || (x > in->numCols + nXedge) || (y < -nYedge) || (y > in->numRows + nYedge)) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Point (%f,%f) lies outside binned image", x, y);
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Point (%lf,%lf) lies outside binned image", x, y);
         return NAN;
     }
@@ -281,6 +281,6 @@
     // if we have a single pixel, there is no spatial information
     if ((in->numCols == 1) && (in->numRows == 1)) {
-	const double value = in->data.F32[0][0];
-	return value;
+        const double value = in->data.F32[0][0];
+        return value;
     }
 
@@ -306,21 +306,21 @@
     // if Nx == 1, we have no x-dir spatial information
     if (in->numCols == 1) {
-	double V0 = in->data.F32[Ys][Xs];
-	double V1 = in->data.F32[Ye][Xs];
-
-	const double value = V0*ry + V1*dy;
-	return value;
-    }	
+        double V0 = in->data.F32[Ys][Xs];
+        double V1 = in->data.F32[Ye][Xs];
+
+        const double value = V0*ry + V1*dy;
+        return value;
+    }
 
     // if Ny == 1, we have no y-dir spatial information
     if (in->numRows == 1) {
-	double V0 = in->data.F32[Ys][Xs];
-	double V1 = in->data.F32[Ys][Xe];
-
-	const double value = V0*rx + V1*dx;
-	return value;
-    }	
-
-    // Vxy 
+        double V0 = in->data.F32[Ys][Xs];
+        double V1 = in->data.F32[Ys][Xe];
+
+        const double value = V0*rx + V1*dx;
+        return value;
+    }
+
+    // Vxy
     double V00 = in->data.F32[Ys][Xs];
     double V01 = in->data.F32[Ye][Xs];
@@ -332,33 +332,33 @@
     // corners
     if ((dx < 0.0) && (dy < 0.0)) {
-	return V00;
-    }	
+        return V00;
+    }
     if ((dx > 1.0) && (dy < 0.0)) {
-	return V10;
-    }	
+        return V10;
+    }
     if ((dx < 0.0) && (dy > 1.0)) {
-	return V01;
-    }	
+        return V01;
+    }
     if ((dx > 1.0) && (dy > 1.0)) {
-	return V11;
-    }	
+        return V11;
+    }
 
     // sides
     if (dx < 0.0) {
-	value = V00*ry + V01*dy;
-	return value;
-    }	
+        value = V00*ry + V01*dy;
+        return value;
+    }
     if (dy < 0.0) {
-	value = V00*rx + V10*dx;
-	return value;
-    }	
+        value = V00*rx + V10*dx;
+        return value;
+    }
     if (dx > 1.0) {
-	value = V10*ry + V11*dy;
-	return value;
-    }	
+        value = V10*ry + V11*dy;
+        return value;
+    }
     if (dy > 1.0) {
-	value = V01*rx + V11*dx;
-	return value;
-    }	
+        value = V01*rx + V11*dx;
+        return value;
+    }
 
     // bilinear interpolation
Index: branches/eam_branches/20090522/psLib/src/sys/psAbort.c
===================================================================
--- branches/eam_branches/20090522/psLib/src/sys/psAbort.c	(revision 24238)
+++ branches/eam_branches/20090522/psLib/src/sys/psAbort.c	(revision 24557)
@@ -10,5 +10,5 @@
  *  @author Eric Van Alst, MHPCC
  *  @author Joshua Hoblitt, University of Hawaii
- *   
+ *
  *  @version $Revision: 1.16 $ $Name: not supported by cvs2svn $
  *  @date $Date: 2008-04-13 08:18:27 $
@@ -34,5 +34,5 @@
                ...)
 {
-    psErrorStackPrint(stderr, "Aborting. Error stack:");
+    psErrorStackPrint(stderr, "Aborting in function %s at %s:%d. Error stack:", func, file, lineno);
 
     va_list argPtr;             // variable list arguement pointer
@@ -51,12 +51,12 @@
 
 void p_psAssert(const char *file,
-		unsigned int lineno,
-		const char *func,
-		const bool value,
-		const char *format,
-		...)
+                unsigned int lineno,
+                const char *func,
+                const bool value,
+                const char *format,
+                ...)
 {
     if (value) return;
-    psErrorStackPrint(stderr, "Aborting. Error stack:");
+    psErrorStackPrint(stderr, "Assertion failed in function %s at %s:%d. Error stack:", func, file, lineno);
 
     va_list argPtr;             // variable list arguement pointer
Index: branches/eam_branches/20090522/psLib/test/Makefile.am
===================================================================
--- branches/eam_branches/20090522/psLib/test/Makefile.am	(revision 24238)
+++ branches/eam_branches/20090522/psLib/test/Makefile.am	(revision 24557)
@@ -1,3 +1,3 @@
-SUBDIRS = tap pstap $(SRCDIRS)
+SUBDIRS = tap pstap optime $(SRCDIRS)
 
 TESTS = test.pl
Index: branches/eam_branches/20090522/psLib/test/astro/tap_psTime_01.c
===================================================================
--- branches/eam_branches/20090522/psLib/test/astro/tap_psTime_01.c	(revision 24238)
+++ branches/eam_branches/20090522/psLib/test/astro/tap_psTime_01.c	(revision 24557)
@@ -842,10 +842,8 @@
     {
         psMemId id = psMemGetId();
-        bool status = psTimeConvert(NULL, PS_TIME_TAI);
-
-        ok(status == false, "psTimeConvert(NULL, PS_TIME_TAI) returned NULL");
-        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-
-        psFree(time2);
+        bool status= psTimeConvert(NULL, PS_TIME_TAI);
+
+        ok(status== false, "psTimeConvert(NULL, PS_TIME_TAI) returned NULL");
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
     }
 
@@ -858,8 +856,7 @@
         psMemId id = psMemGetId();
         psTime *time1 = psTimeAlloc(PS_TIME_TAI);
-        psTime *time2 = psTimeConvert(time1,-100);
-        ok(time2 == NULL, "psTimeConvert(time1, -100) returned NULL");
-        psFree(time1);
-        psFree(time2);
+        bool status = psTimeConvert(time1,-100);
+        ok(status == false, "psTimeConvert(time1, -100) returned false");
+        psFree(time1);
         ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
     }
@@ -873,8 +870,7 @@
         psMemId id = psMemGetId();
         psTime *time1 = psTimeAlloc(PS_TIME_UTC);
-        psTime *time2 = psTimeConvert(time1,-100);
-        ok(time2 == NULL, "psTimeConvert(time1, -100) returned NULL");
-        psFree(time1);
-        psFree(time2);
+        bool status = psTimeConvert(time1,-100);
+        ok(status == false, "psTimeConvert(time1, -100) returned false");
+        psFree(time1);
         ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
     }
@@ -888,8 +884,7 @@
         psMemId id = psMemGetId();
         psTime *time1 = psTimeAlloc(PS_TIME_TT);
-        psTime *time2 = psTimeConvert(time1,-100);
-        ok(time2 == NULL, "psTimeConvert(time1, -100) returned NULL");
-        psFree(time1);
-        psFree(time2);
+        bool status = psTimeConvert(time1,-100);
+        ok(status == false, "psTimeConvert(time1, -100) returned false");
+        psFree(time1);
         ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
     }
@@ -903,8 +898,7 @@
         psTime *time1 = psTimeAlloc(PS_TIME_TAI);
         time1->type = -100;
-        psTime *time2 = psTimeConvert(time1,PS_TIME_TAI);
-        ok(time2 == NULL, "psTimeConvert(time1, PS_TIME_TAI) returned NULL");
-        psFree(time1);
-        psFree(time2);
+        bool status = psTimeConvert(time1,PS_TIME_TAI);
+        ok(status == false, "psTimeConvert(time1, PS_TIME_TAI) returned false");
+        psFree(time1);
         ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
     }
@@ -919,8 +913,7 @@
         psTime *time1 = psTimeAlloc(PS_TIME_TAI);
         time1->nsec = 2e9;
-        psTime *time2 = psTimeConvert(time1, PS_TIME_TAI);
-        ok(time2 == NULL, "psTimeConvert(time1, PS_TIME_TAI) returns NULL for incorrect psTime object");
-        psFree(time1);
-        psFree(time2);
+        bool status = psTimeConvert(time1, PS_TIME_TAI);
+        ok(status == false, "psTimeConvert(time1, PS_TIME_TAI) returns NULL for incorrect psTime object");
+        psFree(time1);
         ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
     }
@@ -929,4 +922,5 @@
     // psTimeConvert()
     //Attempt to convert a time to the same type
+    //Should return true because time->type == type
     {
         psMemId id = psMemGetId();
@@ -936,10 +930,6 @@
         time1->type = PS_TIME_TAI;
         time1->leapsecond = false;
-        psTime *time2 = psTimeConvert(time1, PS_TIME_TAI);
-        ok(time2 == time1, "psTimeConvert(time, ...) returns time for conversion to same type");
-        is_long(time2->sec, 1, "time->sec");
-        is_long(time2->nsec, 2, "time->nsec");
-        ok(time2->type == PS_TIME_TAI, "time->type");
-        is_bool(time2->leapsecond, false, "time->leapsecond");
+        bool status = psTimeConvert(time1, PS_TIME_TAI);
+        ok(status == true, "psTimeConvert(time, ...) returns true for conversion to same type");
 
         psFree(time1);
@@ -957,10 +947,10 @@
         time1->leapsecond = false;
 
-        psTime *time2 = psTimeConvert(time1,PS_TIME_TAI);
-        ok(time2 == time1, "psTimeConvert(time, ...) returns time after conversion to a different type");
-        is_long(time2->sec, testTime1SecondsTAI, "psTimeConvert() returned the correct ->sec");
-        is_long(time2->nsec, testTime1NanosecondsTAI, "psTimeConvert() returned the correct ->nsec");
-        ok(time2->type == PS_TIME_TAI, "psTimeConvert() returned the correct type");
-        is_bool(time2->leapsecond, false, "time->leapsecond");
+        bool status = psTimeConvert(time1,PS_TIME_TAI);
+        ok(status == true, "psTimeConvert(time, ...) returns true after conversion to a different type");
+        is_long(time1->sec, testTime1SecondsTAI, "psTimeConvert() returned the correct ->sec");
+        is_long(time1->nsec, testTime1NanosecondsTAI, "psTimeConvert() returned the correct ->nsec");
+        ok(time1->type == PS_TIME_TAI, "psTimeConvert() returned the correct type");
+        is_bool(time1->leapsecond, false, "time->leapsecond");
 
         psFree(time1);
@@ -978,10 +968,10 @@
         time1->leapsecond = false;
 
-        psTime *time2 = psTimeConvert(time1,PS_TIME_TT);
-        ok(time2 == time1, "psTimeConvert(time, ...) returns time after conversion to a different type");
-        is_long(time2->sec, testTime1SecondsTT, "psTimeConvert() returned the correct ->sec");
-        is_long(time2->nsec, testTime1NanosecondsTT, "psTimeConvert() returned the correct ->nsec");
-        ok(time2->type == PS_TIME_TT, "psTimeConvert() returned the correct type");
-        is_bool(time2->leapsecond, false, "time->leapsecond");
+        bool status = psTimeConvert(time1,PS_TIME_TT);
+        ok(status == true, "psTimeConvert(time, ...) returns true after conversion to a different type");
+        is_long(time1->sec, testTime1SecondsTT, "psTimeConvert() returned the correct ->sec");
+        is_long(time1->nsec, testTime1NanosecondsTT, "psTimeConvert() returned the correct ->nsec");
+        ok(time1->type == PS_TIME_TT, "psTimeConvert() returned the correct type");
+        is_bool(time1->leapsecond, false, "time->leapsecond");
 
         psFree(time1);
@@ -999,10 +989,10 @@
         time1->leapsecond = false;
 
-        psTime *time2 = psTimeConvert(time1,PS_TIME_UT1);
-        ok(time2 == time1, "psTimeConvert(time, ...) returns time after conversion to a different type");
-        is_long(time2->sec, testTime1SecondsUT1, "psTimeConvert() returned the correct ->sec");
-        is_long(time2->nsec, testTime1NanosecondsUT1, "psTimeConvert() returned the correct ->nsec");
-        ok(time2->type == PS_TIME_UT1, "psTimeConvert() returned the correct type");
-        is_bool(time2->leapsecond, false, "time->leapsecond");
+        bool status = psTimeConvert(time1,PS_TIME_UT1);
+        ok(status == true, "psTimeConvert(time, ...) returns true after conversion to a different type");
+        is_long(time1->sec, testTime1SecondsUT1, "psTimeConvert() returned the correct ->sec");
+        is_long(time1->nsec, testTime1NanosecondsUT1, "psTimeConvert() returned the correct ->nsec");
+        ok(time1->type == PS_TIME_UT1, "psTimeConvert() returned the correct type");
+        is_bool(time1->leapsecond, false, "time->leapsecond");
 
         psFree(time1);
@@ -1019,10 +1009,10 @@
         time1->leapsecond = false;
 
-        psTime *time2 = psTimeConvert(time1,PS_TIME_UTC);
-        ok(time2 == time1, "psTimeConvert(time, ...) returns time after conversion to a different type");
-        is_long(time2->sec, testTime1SecondsUTC, "psTimeConvert() returned the correct ->sec");
-        is_long(time2->nsec, testTime1NanosecondsUTC, "psTimeConvert() returned the correct ->nsec");
-        ok(time2->type == PS_TIME_UTC, "psTimeConvert() returned the correct type");
-        is_bool(time2->leapsecond, false, "time->leapsecond");
+        bool status = psTimeConvert(time1,PS_TIME_UTC);
+        ok(status == true, "psTimeConvert(time, ...) returns true after conversion to a different type");
+        is_long(time1->sec, testTime1SecondsUTC, "psTimeConvert() returned the correct ->sec");
+        is_long(time1->nsec, testTime1NanosecondsUTC, "psTimeConvert() returned the correct ->nsec");
+        ok(time1->type == PS_TIME_UTC, "psTimeConvert() returned the correct type");
+        is_bool(time1->leapsecond, false, "time->leapsecond");
 
         psFree(time1);
@@ -1041,10 +1031,10 @@
         time1->leapsecond = false;
 
-        psTime *time2 = psTimeConvert(time1,PS_TIME_TT);
-        ok(time2 == time1, "psTimeConvert() returned time for conversion to same type");
-        is_long(time2->sec, testTime1SecondsTT, "psTimeConvert() returned the correct ->sec");
-        is_long(time2->nsec, testTime1NanosecondsTT, "psTimeConvert() returned the correct ->nsec");
-        ok(time2->type == PS_TIME_TT, "psTimeConvert() returned the correct type");
-        is_bool(time2->leapsecond, false, "time->leapsecond");
+        bool status = psTimeConvert(time1,PS_TIME_TT);
+        ok(status == true, "psTimeConvert() returned true for conversion to same type");
+        is_long(time1->sec, testTime1SecondsTT, "psTimeConvert() returned the correct ->sec");
+        is_long(time1->nsec, testTime1NanosecondsTT, "psTimeConvert() returned the correct ->nsec");
+        ok(time1->type == PS_TIME_TT, "psTimeConvert() returned the correct type");
+        is_bool(time1->leapsecond, false, "time->leapsecond");
 
         psFree(time1);
@@ -1062,10 +1052,10 @@
         time1->leapsecond = false;
 
-        psTime *time2 = psTimeConvert(time1,PS_TIME_UT1);
-        ok(time1 == time2, "psTimeConvert() returned time for conversion to same type");
-        is_long(time2->sec, testTime1SecondsUT1, "psTimeConvert() returned the correct ->sec");
-        is_long(time2->nsec, testTime1NanosecondsUT1, "psTimeConvert() returned the correct ->nsec");
-        ok(time2->type == PS_TIME_UT1, "psTimeConvert() returned the correct type");
-        is_bool(time2->leapsecond, false, "time->leapsecond");
+        bool status = psTimeConvert(time1,PS_TIME_UT1);
+        ok(status == true, "psTimeConvert() returned true for conversion to same type");
+        is_long(time1->sec, testTime1SecondsUT1, "psTimeConvert() returned the correct ->sec");
+        is_long(time1->nsec, testTime1NanosecondsUT1, "psTimeConvert() returned the correct ->nsec");
+        ok(time1->type == PS_TIME_UT1, "psTimeConvert() returned the correct type");
+        is_bool(time1->leapsecond, false, "time->leapsecond");
 
         psFree(time1);
@@ -1083,10 +1073,10 @@
         time1->leapsecond = false;
 
-        psTime *time2 = psTimeConvert(time1,PS_TIME_UTC);
-        ok(time2 == time1, "psTimeConvert() returned time for conversion to same type");
-        is_long(time2->sec, testTime1SecondsUTC, "psTimeConvert() returned the correct ->sec");
-        is_long(time2->nsec, testTime1NanosecondsUTC, "psTimeConvert() returned the correct ->nsec");
-        ok(time2->type == PS_TIME_UTC, "psTimeConvert() returned the correct type");
-        is_bool(time2->leapsecond, false, "time->leapsecond");
+        bool status = psTimeConvert(time1,PS_TIME_UTC);
+        ok(status == true, "psTimeConvert() returned true for conversion to same type");
+        is_long(time1->sec, testTime1SecondsUTC, "psTimeConvert() returned the correct ->sec");
+        is_long(time1->nsec, testTime1NanosecondsUTC, "psTimeConvert() returned the correct ->nsec");
+        ok(time1->type == PS_TIME_UTC, "psTimeConvert() returned the correct type");
+        is_bool(time1->leapsecond, false, "time->leapsecond");
 
         psFree(time1);
@@ -1104,10 +1094,10 @@
         time1->leapsecond = false;
 
-        psTime *time2 = psTimeConvert(time1,PS_TIME_TAI);
-        ok(time2 == time1, "psTimeConvert() returned time for conversion to same type");
-        is_long(time2->sec, testTime1SecondsTAI, "psTimeConvert() returned the correct ->sec");
-        is_long(time2->nsec, testTime1NanosecondsTAI, "psTimeConvert() returned the correct ->nsec");
-        ok(time2->type == PS_TIME_TAI, "psTimeConvert() returned the correct type");
-        is_bool(time2->leapsecond, false, "time->leapsecond");
+        bool status = psTimeConvert(time1,PS_TIME_TAI);
+        ok(status == true, "psTimeConvert() returned true for conversion to same type");
+        is_long(time1->sec, testTime1SecondsTAI, "psTimeConvert() returned the correct ->sec");
+        is_long(time1->nsec, testTime1NanosecondsTAI, "psTimeConvert() returned the correct ->nsec");
+        ok(time1->type == PS_TIME_TAI, "psTimeConvert() returned the correct type");
+        is_bool(time1->leapsecond, false, "time->leapsecond");
 
         psFree(time1);
@@ -1125,10 +1115,10 @@
         time1->leapsecond = false;
 
-        psTime *time2 = psTimeConvert(time1,PS_TIME_UT1);
-        ok(time1 == time2, "psTimeConvert() returned time for conversion to same type");
-        is_long(time2->sec, testTime1SecondsUT1, "psTimeConvert() returned the correct ->sec");
-        is_long(time2->nsec, testTime1NanosecondsUT1, "psTimeConvert() returned the correct ->nsec");
-        ok(time2->type == PS_TIME_UT1, "psTimeConvert() returned the correct type");
-        is_bool(time2->leapsecond, false, "time->leapsecond");
+        bool status = psTimeConvert(time1,PS_TIME_UT1);
+        ok(status == true, "psTimeConvert() returned true for conversion to same type");
+        is_long(time1->sec, testTime1SecondsUT1, "psTimeConvert() returned the correct ->sec");
+        is_long(time1->nsec, testTime1NanosecondsUT1, "psTimeConvert() returned the correct ->nsec");
+        ok(time1->type == PS_TIME_UT1, "psTimeConvert() returned the correct type");
+        is_bool(time1->leapsecond, false, "time->leapsecond");
 
         psFree(time1);
@@ -1147,8 +1137,7 @@
         time1->leapsecond = false;
 
-        psTime *time2 = psTimeConvert(time1, PS_TIME_UTC);
-        ok(time2 == NULL, "psTimeConvert() returned NULL for conversion from UT1");
-        psFree(time1);
-        psFree(time2);
+        bool status = psTimeConvert(time1, PS_TIME_UTC);
+        ok(status == false, "psTimeConvert() returned false for conversion from UT1");
+        psFree(time1);
         ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
     }
Index: branches/eam_branches/20090522/psLib/test/astro/tap_psTime_03.c
===================================================================
--- branches/eam_branches/20090522/psLib/test/astro/tap_psTime_03.c	(revision 24238)
+++ branches/eam_branches/20090522/psLib/test/astro/tap_psTime_03.c	(revision 24557)
@@ -494,5 +494,5 @@
             psFree(timeStr);
 
-            time = psTimeConvert(time, PS_TIME_TT);
+            bool status = psTimeConvert(time, PS_TIME_TT);
             timeStr = psTimeToISO(time);
             is_str(timeStr, testTimeBStrTT[i], "TT ISO string");
@@ -500,5 +500,5 @@
 
             // Verify UTC ISO string
-            time = psTimeConvert(time, PS_TIME_UTC);
+            status = psTimeConvert(time, PS_TIME_UTC);
             time->leapsecond = testTimeBLeapsecond[i];
             timeStr = psTimeToISO(time);
@@ -506,5 +506,5 @@
             psFree(timeStr);
 
-            time = psTimeConvert(time, PS_TIME_UT1);
+            status = psTimeConvert(time, PS_TIME_UT1);
             timeStr = psTimeToISO(time);
             is_str(timeStr, testTimeBStrUT1[i], "UT1 ISO string");
Index: branches/eam_branches/20090522/psLib/test/math/Makefile.am
===================================================================
--- branches/eam_branches/20090522/psLib/test/math/Makefile.am	(revision 24238)
+++ branches/eam_branches/20090522/psLib/test/math/Makefile.am	(revision 24557)
@@ -43,5 +43,4 @@
 	tap_psStats08 \
 	tap_psStats09 \
-	tap_psStatsTiming \
 	tap_psFunc01 \
 	tap_psStats_Sample_01 \
Index: branches/eam_branches/20090522/psLib/test/math/tap_psStatsTiming.c
===================================================================
--- branches/eam_branches/20090522/psLib/test/math/tap_psStatsTiming.c	(revision 24238)
+++ 	(revision )
@@ -1,828 +1,0 @@
-#include <stdio.h>
-#include <string.h>
-#include <pslib.h>
-
-#include "tap.h"
-#include "pstap.h"
-
-// example tap lines:
-// ok(condition, "condition succeeded");
-// skip_start(condition, Nskip, "Skipping tests because of failure");
-
-# define DTIME(A,B) ((A.tv_sec - B.tv_sec) + 1e-6*(A.tv_usec - B.tv_usec))
-struct timeval start, mark;
-
-int main (void)
-{
-    plan_tests(68);
-
-//    diag("psStats timing tests");
-
-    // build a gauss-deviate vector (mean = 0.0, sigma = 1.0) for tests
-    psRandom *seed = psRandomAllocSpecific (PS_RANDOM_TAUS, 0);
-    psVector *rnd = psVectorAlloc (1000, PS_TYPE_F32);
-    for (int i = 0; i < rnd->n; i++) {
-        rnd->data.F32[i] = psRandomGaussian (seed);
-    }
-
-//    diag ("timing for sample mean");
-    /********** SAMPLE MEAN ***********/
-    // test stat sample mean (no mask, no range)
-    {
-        psMemId id = psMemGetId();
-
-        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN);
-
-        gettimeofday (&start, NULL);
-        for (int i = 0; i < 10000; i++)
-        {
-            psVectorStats (stats, rnd, NULL, NULL, 0);
-        }
-        gettimeofday (&mark, NULL);
-        psF64 delta = DTIME(mark, start);
-        ok (delta < 0.1, "sample mean %f (mask: 0, range: 0): %.3f sec", stats->sampleMean, delta);
-        psFree (stats);
-
-        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-
-    // test stat sample mean (mask, no range)
-    {
-        psMemId id = psMemGetId();
-
-        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN);
-        psVector *mask = psVectorAlloc (1000, PS_TYPE_U8);
-        psVectorInit (mask, 0);
-        mask->data.U8[100] = 1;
-        mask->data.U8[200] = 1;
-        mask->data.U8[300] = 1;
-
-        gettimeofday (&start, NULL);
-        for (int i = 0; i < 10000; i++)
-        {
-            psVectorStats (stats, rnd, NULL, mask, 1);
-        }
-        gettimeofday (&mark, NULL);
-        psF64 delta = DTIME(mark, start);
-        ok (delta < 0.12, "sample mean %f (mask: 1, range: 0): %.3f sec", stats->sampleMean, delta);
-        psFree (stats);
-        psFree (mask);
-
-        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-
-    // test stat sample mean (no mask, range)
-    {
-        psMemId id = psMemGetId();
-
-        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_USE_RANGE);
-        stats->min = -10;
-        stats->max = +10;
-
-        gettimeofday (&start, NULL);
-        for (int i = 0; i < 10000; i++)
-        {
-            psVectorStats (stats, rnd, NULL, NULL, 0);
-        }
-        gettimeofday (&mark, NULL);
-        psF64 delta = DTIME(mark, start);
-        ok (delta < 0.18, "sample mean %f (mask: 0, range: 1): %.3f sec", stats->sampleMean, delta);
-        psFree (stats);
-
-        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-
-    // test stat sample mean (mask, range)
-    {
-        psMemId id = psMemGetId();
-
-        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_USE_RANGE);
-        stats->min = -10;
-        stats->max = +10;
-        psVector *mask = psVectorAlloc (1000, PS_TYPE_U8);
-        psVectorInit (mask, 0);
-        mask->data.U8[100] = 1;
-        mask->data.U8[200] = 1;
-        mask->data.U8[300] = 1;
-
-        gettimeofday (&start, NULL);
-        for (int i = 0; i < 10000; i++)
-        {
-            psVectorStats (stats, rnd, NULL, mask, 1);
-        }
-        gettimeofday (&mark, NULL);
-        psF64 delta = DTIME(mark, start);
-        ok (delta < 0.2, "sample mean %f (mask: 1, range: 1): %.3f sec", stats->sampleMean, delta);
-        psFree (stats);
-        psFree (mask);
-
-        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-
-    // test stat sample mean (mask, range : small sample)
-    {
-        psMemId id = psMemGetId();
-
-        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_USE_RANGE);
-        stats->min = -10;
-        stats->max = +10;
-        psVector *mask = psVectorAlloc (10, PS_TYPE_U8);
-        psVectorInit (mask, 0);
-        mask->data.U8[3] = 1;
-        int nOld = rnd->n;
-
-        rnd->n = 10;
-        gettimeofday (&start, NULL);
-        for (int i = 0; i < 1000000; i++)
-        {
-            psVectorStats (stats, rnd, NULL, mask, 1);
-        }
-        gettimeofday (&mark, NULL);
-        rnd->n = nOld;
-
-        psF64 delta = DTIME(mark, start);
-        ok (delta < 0.2, "sample mean %f (mask: 1, range: 1): %.3f sec", stats->sampleMean, delta);
-        psFree (stats);
-        psFree (mask);
-
-        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-
-//    diag ("timing for sample median");
-    /********** SAMPLE MEDIAN ***********/
-    // test stat sample median (no mask, no range)
-    {
-        psMemId id = psMemGetId();
-
-        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
-
-        gettimeofday (&start, NULL);
-        for (int i = 0; i < 10000; i++)
-        {
-            psVectorStats (stats, rnd, NULL, NULL, 0);
-        }
-        gettimeofday (&mark, NULL);
-        psF64 delta = DTIME(mark, start);
-        ok (delta < 2.8, "sample median %f (mask: 0, range: 0): %.3f sec", stats->sampleMedian, delta);
-        psFree (stats);
-
-        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-
-    // test stat sample median (mask, no range)
-    {
-        psMemId id = psMemGetId();
-
-        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
-        psVector *mask = psVectorAlloc (1000, PS_TYPE_U8);
-        psVectorInit (mask, 0);
-        mask->data.U8[100] = 1;
-        mask->data.U8[200] = 1;
-        mask->data.U8[300] = 1;
-
-        gettimeofday (&start, NULL);
-        for (int i = 0; i < 10000; i++)
-        {
-            psVectorStats (stats, rnd, NULL, mask, 1);
-        }
-        gettimeofday (&mark, NULL);
-        psF64 delta = DTIME(mark, start);
-        ok (delta < 2.8, "sample median %f (mask: 1, range: 0): %.3f sec", stats->sampleMedian, delta);
-        psFree (stats);
-        psFree (mask);
-
-        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-
-    // test stat sample median (no mask, range)
-    {
-        psMemId id = psMemGetId();
-
-        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_USE_RANGE);
-        stats->min = -10;
-        stats->max = +10;
-
-        gettimeofday (&start, NULL);
-        for (int i = 0; i < 10000; i++)
-        {
-            psVectorStats (stats, rnd, NULL, NULL, 0);
-        }
-        gettimeofday (&mark, NULL);
-        psF64 delta = DTIME(mark, start);
-        ok (delta < 2.8, "sample median %f (mask: 0, range: 1): %.3f sec", stats->sampleMedian, delta);
-        psFree (stats);
-
-        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-
-    // test stat sample median (mask, range)
-    {
-        psMemId id = psMemGetId();
-
-        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_USE_RANGE);
-        stats->min = -10;
-        stats->max = +10;
-        psVector *mask = psVectorAlloc (1000, PS_TYPE_U8);
-        psVectorInit (mask, 0);
-        mask->data.U8[100] = 1;
-        mask->data.U8[200] = 1;
-        mask->data.U8[300] = 1;
-
-        gettimeofday (&start, NULL);
-        for (int i = 0; i < 10000; i++)
-        {
-            psVectorStats (stats, rnd, NULL, mask, 1);
-        }
-        gettimeofday (&mark, NULL);
-        psF64 delta = DTIME(mark, start);
-        ok (delta < 2.8, "sample median %f (mask: 1, range: 1): %.3f sec", stats->sampleMedian, delta);
-        psFree (stats);
-        psFree (mask);
-
-        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-
-//    diag ("timing for sample stdev");
-    /********** SAMPLE STDEV ***********/
-    // test stat sample stdev (no mask, no range)
-    {
-        psMemId id = psMemGetId();
-
-        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_STDEV);
-
-        gettimeofday (&start, NULL);
-        for (int i = 0; i < 10000; i++)
-        {
-            psVectorStats (stats, rnd, NULL, NULL, 0);
-        }
-        gettimeofday (&mark, NULL);
-        psF64 delta = DTIME(mark, start);
-        ok (delta < 0.2, "sample stdev %f (mask: 0, range: 0): %.3f sec", stats->sampleStdev, delta);
-        psFree (stats);
-
-        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-
-    // test stat sample stdev (mask, no range)
-    {
-        psMemId id = psMemGetId();
-
-        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_STDEV);
-        psVector *mask = psVectorAlloc (1000, PS_TYPE_U8);
-        psVectorInit (mask, 0);
-        mask->data.U8[100] = 1;
-        mask->data.U8[200] = 1;
-        mask->data.U8[300] = 1;
-
-        gettimeofday (&start, NULL);
-        for (int i = 0; i < 10000; i++)
-        {
-            psVectorStats (stats, rnd, NULL, mask, 1);
-        }
-        gettimeofday (&mark, NULL);
-        psF64 delta = DTIME(mark, start);
-        ok (delta < 0.27, "sample stdev %f (mask: 1, range: 0): %.3f sec", stats->sampleStdev, delta);
-        psFree (stats);
-        psFree (mask);
-
-        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-
-    // test stat sample stdev (no mask, range)
-    {
-        psMemId id = psMemGetId();
-
-        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_STDEV | PS_STAT_USE_RANGE);
-        stats->min = -10;
-        stats->max = +10;
-
-        gettimeofday (&start, NULL);
-        for (int i = 0; i < 10000; i++)
-        {
-            psVectorStats (stats, rnd, NULL, NULL, 0);
-        }
-        gettimeofday (&mark, NULL);
-        psF64 delta = DTIME(mark, start);
-        ok (delta < 0.36, "sample stdev %f (mask: 0, range: 1): %.3f sec", stats->sampleStdev, delta);
-        psFree (stats);
-
-        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-
-    // test stat sample stdev (mask, range)
-    {
-        psMemId id = psMemGetId();
-
-        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_STDEV | PS_STAT_USE_RANGE);
-        stats->min = -10;
-        stats->max = +10;
-        psVector *mask = psVectorAlloc (1000, PS_TYPE_U8);
-        psVectorInit (mask, 0);
-        mask->data.U8[100] = 1;
-        mask->data.U8[200] = 1;
-        mask->data.U8[300] = 1;
-
-        gettimeofday (&start, NULL);
-        for (int i = 0; i < 10000; i++)
-        {
-            psVectorStats (stats, rnd, NULL, mask, 1);
-        }
-        gettimeofday (&mark, NULL);
-        psF64 delta = DTIME(mark, start);
-        ok (delta < 0.42, "sample stdev %f (mask: 1, range: 1): %.3f sec", stats->sampleStdev, delta);
-        psFree (stats);
-        psFree (mask);
-
-        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-
-    // test stat sample stdev (mask, range)
-    {
-        psMemId id = psMemGetId();
-
-        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_STDEV | PS_STAT_USE_RANGE);
-        stats->min = -10;
-        stats->max = +10;
-        psVector *mask = psVectorAlloc (10, PS_TYPE_U8);
-        psVectorInit (mask, 0);
-        mask->data.U8[1] = 1;
-        int nOld = rnd->n;
-
-        rnd->n = 10;
-        gettimeofday (&start, NULL);
-        for (int i = 0; i < 1000000; i++)
-        {
-            psVectorStats (stats, rnd, NULL, mask, 1);
-        }
-        gettimeofday (&mark, NULL);
-        rnd->n = nOld;
-
-        psF64 delta = DTIME(mark, start);
-        ok (delta < 0.42, "sample stdev %f (mask: 1, range: 1): %.3f sec", stats->sampleStdev, delta);
-        psFree (stats);
-        psFree (mask);
-
-        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-
-//    diag ("timing for sample min,max");
-    /*************** MIN,MAX ******************/
-    // test stat min,max (no mask, no range)
-    {
-        psMemId id = psMemGetId();
-
-        psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
-        gettimeofday (&start, NULL);
-        for (int i = 0; i < 10000; i++)
-        {
-            psVectorStats (stats, rnd, NULL, NULL, 1);
-        }
-        gettimeofday (&mark, NULL);
-        psF64 delta = DTIME(mark, start);
-        ok (delta < 0.17, "sample min,max %f,%f (mask: 0, range: 0): %.3f sec", stats->min, stats->max, delta);
-        psFree (stats);
-
-        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-    // test stat min,max (no mask, no range)
-    {
-        psMemId id = psMemGetId();
-
-        psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
-        psVector *mask = psVectorAlloc (1000, PS_TYPE_U8);
-        psVectorInit (mask, 0);
-        mask->data.U8[100] = 1;
-        mask->data.U8[200] = 1;
-        mask->data.U8[300] = 1;
-
-        gettimeofday (&start, NULL);
-        for (int i = 0; i < 10000; i++)
-        {
-            psVectorStats (stats, rnd, NULL, mask, 1);
-        }
-        gettimeofday (&mark, NULL);
-        psF64 delta = DTIME(mark, start);
-        ok (delta < 0.18, "sample min,max %f,%f (mask: 1, range: 0): %.3f sec", stats->min, stats->max, delta);
-        psFree (stats);
-        psFree (mask);
-
-        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-    // test stat min,max (no mask, no range)
-    {
-        psMemId id = psMemGetId();
-
-        psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX | PS_STAT_USE_RANGE);
-        stats->min = -10;
-        stats->max = +10;
-
-        gettimeofday (&start, NULL);
-        for (int i = 0; i < 10000; i++)
-        {
-            psVectorStats (stats, rnd, NULL, NULL, 1);
-        }
-        gettimeofday (&mark, NULL);
-        psF64 delta = DTIME(mark, start);
-        ok (delta < 0.22, "sample min,max %f,%f (mask: 0, range: 1): %.3f sec", stats->min, stats->max, delta);
-        psFree (stats);
-
-        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-    // test stat min,max (no mask, no range)
-    {
-        psMemId id = psMemGetId();
-
-        psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX | PS_STAT_USE_RANGE);
-        stats->min = -10;
-        stats->max = +10;
-        psVector *mask = psVectorAlloc (1000, PS_TYPE_U8);
-        psVectorInit (mask, 0);
-        mask->data.U8[100] = 1;
-        mask->data.U8[200] = 1;
-        mask->data.U8[300] = 1;
-
-        gettimeofday (&start, NULL);
-        for (int i = 0; i < 10000; i++)
-        {
-            psVectorStats (stats, rnd, NULL, mask, 1);
-        }
-        gettimeofday (&mark, NULL);
-        psF64 delta = DTIME(mark, start);
-        ok (delta < 0.26, "sample min,max %f,%f (mask: 1, range: 1): %.3f sec", stats->min, stats->max, delta);
-        psFree (stats);
-        psFree (mask);
-
-        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-
-//    diag ("timing for clipped stats");
-    /********** CLIPPED STATS ***********/
-    {
-        psMemId id = psMemGetId();
-
-        psStats *stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
-        psVector *rnd2 = psVectorAlloc (1000, PS_TYPE_F32);
-        for (int i = 0; i < rnd2->n; i++)
-        {
-            rnd2->data.F32[i] = psRandomGaussian (seed);
-        }
-
-        gettimeofday (&start, NULL);
-        for (int i = 0; i < 1000; i++)
-        {
-            psVectorStats (stats, rnd2, NULL, NULL, 1);
-        }
-        gettimeofday (&mark, NULL);
-        psF64 delta = DTIME(mark, start);
-        ok (delta < 0.3, "clipped mean %f, stdev %f (mask: 0, range: 0): %.3f sec (1000 pts / 1000 loops)", stats->clippedMean, stats->clippedStdev, delta);
-        psFree (stats);
-        psFree (rnd2);
-
-        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-    {
-        psMemId id = psMemGetId();
-
-        psStats *stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
-        psVector *rnd2 = psVectorAlloc (3000, PS_TYPE_F32);
-        for (int i = 0; i < rnd2->n; i++)
-        {
-            rnd2->data.F32[i] = psRandomGaussian (seed);
-        }
-
-        gettimeofday (&start, NULL);
-        for (int i = 0; i < 1000; i++)
-        {
-            psVectorStats (stats, rnd2, NULL, NULL, 1);
-        }
-        gettimeofday (&mark, NULL);
-        psF64 delta = DTIME(mark, start);
-        ok (delta < 0.5, "clipped mean %f, stdev %f (mask: 0, range: 0): %.3f sec (3000 pts / 1000 loops)", stats->clippedMean, stats->clippedStdev, delta);
-        psFree (stats);
-        psFree (rnd2);
-
-        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-    {
-        psMemId id = psMemGetId();
-
-        psStats *stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
-        psVector *rnd2 = psVectorAlloc (10000, PS_TYPE_F32);
-        for (int i = 0; i < rnd2->n; i++)
-        {
-            rnd2->data.F32[i] = psRandomGaussian (seed);
-        }
-
-        gettimeofday (&start, NULL);
-        for (int i = 0; i < 1000; i++)
-        {
-            psVectorStats (stats, rnd2, NULL, NULL, 1);
-        }
-        gettimeofday (&mark, NULL);
-        psF64 delta = DTIME(mark, start);
-        ok (delta < 1.2, "clipped mean %f, stdev %f (mask: 0, range: 0): %.3f sec (10000 pts / 1000 loops)", stats->clippedMean, stats->clippedStdev, delta);
-        psFree (stats);
-        psFree (rnd2);
-
-        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-
-//    diag ("timing for robust stats");
-    /********** ROBUST STATS ***********/
-    {
-        psMemId id = psMemGetId();
-
-        psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_ROBUST_QUARTILE);
-        psVector *rnd2 = psVectorAlloc (1000, PS_TYPE_F32);
-        for (int i = 0; i < rnd2->n; i++)
-        {
-            rnd2->data.F32[i] = psRandomGaussian (seed);
-        }
-
-        gettimeofday (&start, NULL);
-        for (int i = 0; i < 1000; i++)
-        {
-            psVectorStats (stats, rnd2, NULL, NULL, 1);
-        }
-        gettimeofday (&mark, NULL);
-        psF64 delta = DTIME(mark, start);
-        ok (delta < 0.3, "robust mean %f, stdev %f (mask: 0, range: 0): %.3f sec (1000 pts / 1000 loops)", stats->robustMedian, stats->robustStdev, delta);
-        psFree (stats);
-        psFree (rnd2);
-
-        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-    {
-        psMemId id = psMemGetId();
-
-        psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_ROBUST_QUARTILE);
-        psVector *rnd2 = psVectorAlloc (3000, PS_TYPE_F32);
-        for (int i = 0; i < rnd2->n; i++)
-        {
-            rnd2->data.F32[i] = psRandomGaussian (seed);
-        }
-
-        gettimeofday (&start, NULL);
-        for (int i = 0; i < 1000; i++)
-        {
-            psVectorStats (stats, rnd2, NULL, NULL, 1);
-        }
-        gettimeofday (&mark, NULL);
-        psF64 delta = DTIME(mark, start);
-        ok (delta < 0.5, "robust mean %f, stdev %f (mask: 0, range: 0): %.3f sec (3000 pts / 1000 loops)", stats->robustMedian, stats->robustStdev, delta);
-        psFree (stats);
-        psFree (rnd2);
-
-        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-    {
-        psMemId id = psMemGetId();
-
-        psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_ROBUST_QUARTILE);
-        psVector *rnd2 = psVectorAlloc (10000, PS_TYPE_F32);
-        for (int i = 0; i < rnd2->n; i++)
-        {
-            rnd2->data.F32[i] = psRandomGaussian (seed);
-        }
-
-        gettimeofday (&start, NULL);
-        for (int i = 0; i < 1000; i++)
-        {
-            psVectorStats (stats, rnd2, NULL, NULL, 1);
-        }
-        gettimeofday (&mark, NULL);
-        psF64 delta = DTIME(mark, start);
-        ok (delta < 1.2, "robust mean %f, stdev %f (mask: 0, range: 0): %.3f sec (10000 pts / 1000 loops)", stats->robustMedian, stats->robustStdev, delta);
-        psFree (stats);
-        psFree (rnd2);
-
-        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-
-//    diag ("timing for fitted stats");
-    /********** FITTED TIMING ***********/
-    {
-        psMemId id = psMemGetId();
-
-        psStats *stats = psStatsAlloc (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV);
-        psVector *rnd2 = psVectorAlloc (1000, PS_TYPE_F32);
-        for (int i = 0; i < rnd2->n; i++)
-        {
-            rnd2->data.F32[i] = psRandomGaussian (seed);
-        }
-
-        gettimeofday (&start, NULL);
-        for (int i = 0; i < 1000; i++)
-        {
-            psVectorStats (stats, rnd2, NULL, NULL, 1);
-
-        }
-        gettimeofday (&mark, NULL);
-        psF64 delta = DTIME(mark, start);
-        ok (delta < 0.7, "fitted mean %f, stdev %f (mask: 0, range: 0): %.3f sec (1000 pts / 1000 loops)", stats->fittedMean, stats->fittedStdev, delta);
-        psFree (stats);
-        psFree (rnd2);
-
-        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-    {
-        psMemId id = psMemGetId();
-
-        psStats *stats = psStatsAlloc (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV);
-        psVector *rnd2 = psVectorAlloc (3000, PS_TYPE_F32);
-        for (int i = 0; i < rnd2->n; i++)
-        {
-            rnd2->data.F32[i] = psRandomGaussian (seed);
-        }
-
-        gettimeofday (&start, NULL);
-        for (int i = 0; i < 1000; i++)
-        {
-            psVectorStats (stats, rnd2, NULL, NULL, 1);
-        }
-        gettimeofday (&mark, NULL);
-        psF64 delta = DTIME(mark, start);
-        ok (delta < 0.8, "fitted mean %f, stdev %f (mask: 0, range: 0): %.3f sec (3000 pts / 1000 loops)", stats->fittedMean, stats->fittedStdev, delta);
-        psFree (stats);
-        psFree (rnd2);
-
-        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-    {
-        psMemId id = psMemGetId();
-
-        psStats *stats = psStatsAlloc (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV);
-        psVector *rnd2 = psVectorAlloc (10000, PS_TYPE_F32);
-        for (int i = 0; i < rnd2->n; i++)
-        {
-            rnd2->data.F32[i] = psRandomGaussian (seed);
-        }
-
-        gettimeofday (&start, NULL);
-        for (int i = 0; i < 1000; i++)
-        {
-            psVectorStats (stats, rnd2, NULL, NULL, 1);
-        }
-        gettimeofday (&mark, NULL);
-        psF64 delta = DTIME(mark, start);
-        ok (delta < 2.2, "fitted mean %f, stdev %f (mask: 0, range: 0): %.3f sec (10000 pts / 1000 loops)", stats->fittedMean, stats->fittedStdev, delta);
-        psFree (stats);
-        psFree (rnd2);
-
-        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-
-//    diag ("timing for fitted (v2) stats");
-    /********** FITTED (v2) TIMING ***********/
-    {
-        psMemId id = psMemGetId();
-
-        psStats *stats = psStatsAlloc (PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2);
-        psVector *rnd2 = psVectorAlloc (1000, PS_TYPE_F32);
-        for (int i = 0; i < rnd2->n; i++)
-        {
-            rnd2->data.F32[i] = psRandomGaussian (seed);
-        }
-
-        gettimeofday (&start, NULL);
-        for (int i = 0; i < 1000; i++)
-        {
-            psVectorStats (stats, rnd2, NULL, NULL, 1);
-
-        }
-        gettimeofday (&mark, NULL);
-        psF64 delta = DTIME(mark, start);
-        ok (delta < 0.7, "fitted mean %f, stdev %f (mask: 0, range: 0): %.3f sec (1000 pts / 1000 loops)", stats->fittedMean, stats->fittedStdev, delta);
-        psFree (stats);
-        psFree (rnd2);
-
-        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-    {
-        psMemId id = psMemGetId();
-
-        psStats *stats = psStatsAlloc (PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2);
-        psVector *rnd2 = psVectorAlloc (3000, PS_TYPE_F32);
-        for (int i = 0; i < rnd2->n; i++)
-        {
-            rnd2->data.F32[i] = psRandomGaussian (seed);
-        }
-
-        gettimeofday (&start, NULL);
-        for (int i = 0; i < 1000; i++)
-        {
-            psVectorStats (stats, rnd2, NULL, NULL, 1);
-        }
-        gettimeofday (&mark, NULL);
-        psF64 delta = DTIME(mark, start);
-        ok (delta < 0.8, "fitted mean %f, stdev %f (mask: 0, range: 0): %.3f sec (3000 pts / 1000 loops)", stats->fittedMean, stats->fittedStdev, delta);
-        psFree (stats);
-        psFree (rnd2);
-
-        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-    {
-        psMemId id = psMemGetId();
-
-        psStats *stats = psStatsAlloc (PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2);
-        psVector *rnd2 = psVectorAlloc (10000, PS_TYPE_F32);
-        for (int i = 0; i < rnd2->n; i++)
-        {
-            rnd2->data.F32[i] = psRandomGaussian (seed);
-        }
-
-        gettimeofday (&start, NULL);
-        for (int i = 0; i < 1000; i++)
-        {
-            psVectorStats (stats, rnd2, NULL, NULL, 1);
-        }
-        gettimeofday (&mark, NULL);
-        psF64 delta = DTIME(mark, start);
-        ok (delta < 2.2, "fitted mean %f, stdev %f (mask: 0, range: 0): %.3f sec (10000 pts / 1000 loops)", stats->fittedMean, stats->fittedStdev, delta);
-        psFree (stats);
-        psFree (rnd2);
-
-        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-
-//    diag ("compare sample, robust, and fitted mean and stdev to theoretical");
-    // compare SAMPLE, FITTED, ROBUST mean to theoretical
-    {
-        psMemId id = psMemGetId();
-
-        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV);
-        psVector *sample = psVectorAlloc (1000, PS_TYPE_F32);
-        psVector *robust = psVectorAlloc (1000, PS_TYPE_F32);
-        psVector *fitted = psVectorAlloc (1000, PS_TYPE_F32);
-
-        for (int i = 0; i < 1000; i++)
-        {
-            // generate a new sample
-            for (int j = 0; j < rnd->n; j++) {
-                rnd->data.F32[j] = psRandomGaussian (seed);
-            }
-            // measure the stats
-            psVectorStats (stats, rnd, NULL, NULL, 1);
-            sample->data.F32[i] = stats->sampleMean;
-            robust->data.F32[i] = stats->robustMedian;
-            fitted->data.F32[i] = stats->fittedMean;
-        }
-        psFree (stats);
-
-        stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
-        psVectorStats (stats, sample, NULL, NULL, 1);
-        ok (stats->sampleStdev < 2/sqrt(1000), "sample mean %f, stdev %f (1000 tries)", stats->sampleMean, stats->sampleStdev);
-        psVectorStats (stats, robust, NULL, NULL, 1);
-        ok (stats->sampleStdev < 2/sqrt(1000), "robust mean %f, stdev %f (1000 tries)", stats->sampleMean, stats->sampleStdev);
-        psVectorStats (stats, fitted, NULL, NULL, 1);
-        ok (stats->sampleStdev < 2/sqrt(1000), "fitted mean %f, stdev %f (1000 tries)", stats->sampleMean, stats->sampleStdev);
-        psFree (stats);
-        psFree (sample);
-        psFree (robust);
-        psFree (fitted);
-
-        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-
-//    diag ("compare sample, robust, and fitted mean and stdev to theoretical");
-    // compare SAMPLE, FITTED_V2, ROBUST mean to theoretical
-    {
-        psMemId id = psMemGetId();
-
-        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2);
-        psVector *sample = psVectorAlloc (1000, PS_TYPE_F32);
-        psVector *robust = psVectorAlloc (1000, PS_TYPE_F32);
-        psVector *fitted = psVectorAlloc (1000, PS_TYPE_F32);
-
-        for (int i = 0; i < 1000; i++)
-        {
-            // generate a new sample
-            for (int j = 0; j < rnd->n; j++) {
-                rnd->data.F32[j] = psRandomGaussian (seed);
-            }
-            // measure the stats
-            psVectorStats (stats, rnd, NULL, NULL, 1);
-            sample->data.F32[i] = stats->sampleMean;
-            robust->data.F32[i] = stats->robustMedian;
-            fitted->data.F32[i] = stats->fittedMean;
-        }
-        psFree (stats);
-
-        stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
-        psVectorStats (stats, sample, NULL, NULL, 1);
-        ok (stats->sampleStdev < 2/sqrt(1000), "sample mean %f, stdev %f (1000 tries)", stats->sampleMean, stats->sampleStdev);
-        psVectorStats (stats, robust, NULL, NULL, 1);
-        ok (stats->sampleStdev < 2/sqrt(1000), "robust mean %f, stdev %f (1000 tries)", stats->sampleMean, stats->sampleStdev);
-        psVectorStats (stats, fitted, NULL, NULL, 1);
-        ok (stats->sampleStdev < 2/sqrt(1000), "fitted mean %f, stdev %f (1000 tries)", stats->sampleMean, stats->sampleStdev);
-        psFree (stats);
-        psFree (sample);
-        psFree (robust);
-        psFree (fitted);
-
-        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-
-    return exit_status();
-}
-
Index: branches/eam_branches/20090522/psLib/test/math/tap_psStatsTiming.txt
===================================================================
--- branches/eam_branches/20090522/psLib/test/math/tap_psStatsTiming.txt	(revision 24238)
+++ 	(revision )
@@ -1,47 +1,0 @@
-
-Running tap_psStatsTiming on alala (dual AMD Opteron, 64bit, 2.2GHz)
-yields the following timing results:
-
-# timing for sample mean (1000 loops of 10000 pts)
-ok 1 - sample mean 0.009719 (mask: 0, range: 0): 0.072 sec
-ok 3 - sample mean 0.011060 (mask: 1, range: 0): 0.119 sec
-ok 5 - sample mean 0.009719 (mask: 0, range: 1): 0.170 sec
-ok 7 - sample mean 0.011060 (mask: 1, range: 1): 0.198 sec
-
-# timing for sample median (1000 loops of 10000 pts)
-ok 9 - sample median 0.021781 (mask: 0, range: 0): 2.625 sec
-ok 11 - sample median 0.023795 (mask: 1, range: 0): 2.646 sec
-ok 13 - sample median 0.021781 (mask: 0, range: 1): 2.703 sec
-ok 15 - sample median 0.023795 (mask: 1, range: 1): 2.716 sec
-
-# timing for sample stdev (1000 loops of 10000 pts)
-ok 17 - sample stdev 0.964753 (mask: 0, range: 0): 0.193 sec
-ok 19 - sample stdev 0.965887 (mask: 1, range: 0): 0.257 sec
-ok 21 - sample stdev 0.964753 (mask: 0, range: 1): 0.353 sec
-ok 23 - sample stdev 0.965887 (mask: 1, range: 1): 0.401 sec
-
-# timing for sample min,max (1000 loops of 10000 pts)
-ok 25 - sample min,max -3.205688,2.706797 (mask: 0, range: 0): 0.125 sec
-ok 27 - sample min,max -3.205688,2.706797 (mask: 1, range: 0): 0.152 sec
-ok 29 - sample min,max -3.205688,2.706797 (mask: 0, range: 1): 0.201 sec
-ok 31 - sample min,max -3.205688,2.706797 (mask: 1, range: 1): 0.238 sec
-
-# timing for clipped stats
-not ok 33 - clipped mean -0.047714, stdev 0.991979 (mask: 0, range: 0): 0.369 sec (1000 pts / 1000 loops)
-not ok 35 - clipped mean 0.023963, stdev 0.972186 (mask: 0, range: 0): 1.219 sec (3000 pts / 1000 loops)
-not ok 37 - clipped mean -0.007020, stdev 0.985410 (mask: 0, range: 0): 4.883 sec (10000 pts / 1000 loops)
-
-NOTE: these fail because they are being compared to the 'robust' stats
-limits below.  The clipped mean algorithm should not be so slow (and
-apparently non-linear in npts).
-
-# timing for robust stats
-ok 39 - robust mean 0.123348, stdev 1.014896 (mask: 0, range: 0): 0.187 sec (1000 pts / 1000 loops)
-ok 41 - robust mean -0.006812, stdev 0.974468 (mask: 0, range: 0): 0.382 sec (3000 pts / 1000 loops)
-ok 43 - robust mean -0.013591, stdev 1.001539 (mask: 0, range: 0): 1.076 sec (10000 pts / 1000 loops)
-
-# timing for fitted stats
-ok 45 - fitted mean -0.029859, stdev 0.982947 (mask: 0, range: 0): 0.381 sec (1000 pts / 1000 loops)
-ok 47 - fitted mean 0.014660, stdev 0.956168 (mask: 0, range: 0): 0.727 sec (3000 pts / 1000 loops)
-ok 49 - fitted mean -0.008402, stdev 1.001366 (mask: 0, range: 0): 1.914 sec (10000 pts / 1000 loops)
-
Index: branches/eam_branches/20090522/psLib/test/optime/Makefile.am
===================================================================
--- branches/eam_branches/20090522/psLib/test/optime/Makefile.am	(revision 24557)
+++ branches/eam_branches/20090522/psLib/test/optime/Makefile.am	(revision 24557)
@@ -0,0 +1,27 @@
+AM_CPPFLAGS = \
+	$(SRCINC) \
+	-I$(top_srcdir)/test/tap/src \
+	-I$(top_srcdir)/test/pstap/src \
+	$(PSLIB_CFLAGS)
+AM_LDFLAGS = \
+	$(top_builddir)/src/libpslib.la  \
+	$(top_builddir)/test/tap/src/libtap.la \
+	$(top_builddir)/test/pstap/src/libpstap.la \
+	$(PSLIB_LIBS)
+
+TEST_PROGS = \
+	tap_psStatsTiming
+
+if BUILD_TESTS
+bin_PROGRAMS = $(TEST_PROGS)
+TESTS = $(TEST_PROGS)
+else
+check_PROGRAMS = $(TEST_PROGS)
+endif
+
+CLEANFILES = $(tmp_files) core core.* *~ *.bb *.bbg *.da gmon.out
+
+tests: $(check_PROGRAMS)
+	$(top_srcdir)/test/test.pl
+
+test: check
Index: branches/eam_branches/20090522/psLib/test/optime/tap_psStatsTiming.c
===================================================================
--- branches/eam_branches/20090522/psLib/test/optime/tap_psStatsTiming.c	(revision 24557)
+++ branches/eam_branches/20090522/psLib/test/optime/tap_psStatsTiming.c	(revision 24557)
@@ -0,0 +1,828 @@
+#include <stdio.h>
+#include <string.h>
+#include <pslib.h>
+
+#include "tap.h"
+#include "pstap.h"
+
+// example tap lines:
+// ok(condition, "condition succeeded");
+// skip_start(condition, Nskip, "Skipping tests because of failure");
+
+# define DTIME(A,B) ((A.tv_sec - B.tv_sec) + 1e-6*(A.tv_usec - B.tv_usec))
+struct timeval start, mark;
+
+int main (void)
+{
+    plan_tests(68);
+
+//    diag("psStats timing tests");
+
+    // build a gauss-deviate vector (mean = 0.0, sigma = 1.0) for tests
+    psRandom *seed = psRandomAllocSpecific (PS_RANDOM_TAUS, 0);
+    psVector *rnd = psVectorAlloc (1000, PS_TYPE_F32);
+    for (int i = 0; i < rnd->n; i++) {
+        rnd->data.F32[i] = psRandomGaussian (seed);
+    }
+
+//    diag ("timing for sample mean");
+    /********** SAMPLE MEAN ***********/
+    // test stat sample mean (no mask, no range)
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN);
+
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 10000; i++)
+        {
+            psVectorStats (stats, rnd, NULL, NULL, 0);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 0.1, "sample mean %f (mask: 0, range: 0): %.3f sec", stats->sampleMean, delta);
+        psFree (stats);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+
+    // test stat sample mean (mask, no range)
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN);
+        psVector *mask = psVectorAlloc (1000, PS_TYPE_U8);
+        psVectorInit (mask, 0);
+        mask->data.U8[100] = 1;
+        mask->data.U8[200] = 1;
+        mask->data.U8[300] = 1;
+
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 10000; i++)
+        {
+            psVectorStats (stats, rnd, NULL, mask, 1);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 0.12, "sample mean %f (mask: 1, range: 0): %.3f sec", stats->sampleMean, delta);
+        psFree (stats);
+        psFree (mask);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+
+    // test stat sample mean (no mask, range)
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_USE_RANGE);
+        stats->min = -10;
+        stats->max = +10;
+
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 10000; i++)
+        {
+            psVectorStats (stats, rnd, NULL, NULL, 0);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 0.18, "sample mean %f (mask: 0, range: 1): %.3f sec", stats->sampleMean, delta);
+        psFree (stats);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+
+    // test stat sample mean (mask, range)
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_USE_RANGE);
+        stats->min = -10;
+        stats->max = +10;
+        psVector *mask = psVectorAlloc (1000, PS_TYPE_U8);
+        psVectorInit (mask, 0);
+        mask->data.U8[100] = 1;
+        mask->data.U8[200] = 1;
+        mask->data.U8[300] = 1;
+
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 10000; i++)
+        {
+            psVectorStats (stats, rnd, NULL, mask, 1);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 0.2, "sample mean %f (mask: 1, range: 1): %.3f sec", stats->sampleMean, delta);
+        psFree (stats);
+        psFree (mask);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+
+    // test stat sample mean (mask, range : small sample)
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_USE_RANGE);
+        stats->min = -10;
+        stats->max = +10;
+        psVector *mask = psVectorAlloc (10, PS_TYPE_U8);
+        psVectorInit (mask, 0);
+        mask->data.U8[3] = 1;
+        int nOld = rnd->n;
+
+        rnd->n = 10;
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 1000000; i++)
+        {
+            psVectorStats (stats, rnd, NULL, mask, 1);
+        }
+        gettimeofday (&mark, NULL);
+        rnd->n = nOld;
+
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 0.2, "sample mean %f (mask: 1, range: 1): %.3f sec", stats->sampleMean, delta);
+        psFree (stats);
+        psFree (mask);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+
+//    diag ("timing for sample median");
+    /********** SAMPLE MEDIAN ***********/
+    // test stat sample median (no mask, no range)
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
+
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 10000; i++)
+        {
+            psVectorStats (stats, rnd, NULL, NULL, 0);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 2.8, "sample median %f (mask: 0, range: 0): %.3f sec", stats->sampleMedian, delta);
+        psFree (stats);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+
+    // test stat sample median (mask, no range)
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
+        psVector *mask = psVectorAlloc (1000, PS_TYPE_U8);
+        psVectorInit (mask, 0);
+        mask->data.U8[100] = 1;
+        mask->data.U8[200] = 1;
+        mask->data.U8[300] = 1;
+
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 10000; i++)
+        {
+            psVectorStats (stats, rnd, NULL, mask, 1);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 2.8, "sample median %f (mask: 1, range: 0): %.3f sec", stats->sampleMedian, delta);
+        psFree (stats);
+        psFree (mask);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+
+    // test stat sample median (no mask, range)
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_USE_RANGE);
+        stats->min = -10;
+        stats->max = +10;
+
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 10000; i++)
+        {
+            psVectorStats (stats, rnd, NULL, NULL, 0);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 2.8, "sample median %f (mask: 0, range: 1): %.3f sec", stats->sampleMedian, delta);
+        psFree (stats);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+
+    // test stat sample median (mask, range)
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_USE_RANGE);
+        stats->min = -10;
+        stats->max = +10;
+        psVector *mask = psVectorAlloc (1000, PS_TYPE_U8);
+        psVectorInit (mask, 0);
+        mask->data.U8[100] = 1;
+        mask->data.U8[200] = 1;
+        mask->data.U8[300] = 1;
+
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 10000; i++)
+        {
+            psVectorStats (stats, rnd, NULL, mask, 1);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 2.8, "sample median %f (mask: 1, range: 1): %.3f sec", stats->sampleMedian, delta);
+        psFree (stats);
+        psFree (mask);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+
+//    diag ("timing for sample stdev");
+    /********** SAMPLE STDEV ***********/
+    // test stat sample stdev (no mask, no range)
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_STDEV);
+
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 10000; i++)
+        {
+            psVectorStats (stats, rnd, NULL, NULL, 0);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 0.2, "sample stdev %f (mask: 0, range: 0): %.3f sec", stats->sampleStdev, delta);
+        psFree (stats);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+
+    // test stat sample stdev (mask, no range)
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_STDEV);
+        psVector *mask = psVectorAlloc (1000, PS_TYPE_U8);
+        psVectorInit (mask, 0);
+        mask->data.U8[100] = 1;
+        mask->data.U8[200] = 1;
+        mask->data.U8[300] = 1;
+
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 10000; i++)
+        {
+            psVectorStats (stats, rnd, NULL, mask, 1);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 0.27, "sample stdev %f (mask: 1, range: 0): %.3f sec", stats->sampleStdev, delta);
+        psFree (stats);
+        psFree (mask);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+
+    // test stat sample stdev (no mask, range)
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_STDEV | PS_STAT_USE_RANGE);
+        stats->min = -10;
+        stats->max = +10;
+
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 10000; i++)
+        {
+            psVectorStats (stats, rnd, NULL, NULL, 0);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 0.36, "sample stdev %f (mask: 0, range: 1): %.3f sec", stats->sampleStdev, delta);
+        psFree (stats);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+
+    // test stat sample stdev (mask, range)
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_STDEV | PS_STAT_USE_RANGE);
+        stats->min = -10;
+        stats->max = +10;
+        psVector *mask = psVectorAlloc (1000, PS_TYPE_U8);
+        psVectorInit (mask, 0);
+        mask->data.U8[100] = 1;
+        mask->data.U8[200] = 1;
+        mask->data.U8[300] = 1;
+
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 10000; i++)
+        {
+            psVectorStats (stats, rnd, NULL, mask, 1);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 0.42, "sample stdev %f (mask: 1, range: 1): %.3f sec", stats->sampleStdev, delta);
+        psFree (stats);
+        psFree (mask);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+
+    // test stat sample stdev (mask, range)
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_STDEV | PS_STAT_USE_RANGE);
+        stats->min = -10;
+        stats->max = +10;
+        psVector *mask = psVectorAlloc (10, PS_TYPE_U8);
+        psVectorInit (mask, 0);
+        mask->data.U8[1] = 1;
+        int nOld = rnd->n;
+
+        rnd->n = 10;
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 1000000; i++)
+        {
+            psVectorStats (stats, rnd, NULL, mask, 1);
+        }
+        gettimeofday (&mark, NULL);
+        rnd->n = nOld;
+
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 0.42, "sample stdev %f (mask: 1, range: 1): %.3f sec", stats->sampleStdev, delta);
+        psFree (stats);
+        psFree (mask);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+
+//    diag ("timing for sample min,max");
+    /*************** MIN,MAX ******************/
+    // test stat min,max (no mask, no range)
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 10000; i++)
+        {
+            psVectorStats (stats, rnd, NULL, NULL, 1);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 0.17, "sample min,max %f,%f (mask: 0, range: 0): %.3f sec", stats->min, stats->max, delta);
+        psFree (stats);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+    // test stat min,max (no mask, no range)
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
+        psVector *mask = psVectorAlloc (1000, PS_TYPE_U8);
+        psVectorInit (mask, 0);
+        mask->data.U8[100] = 1;
+        mask->data.U8[200] = 1;
+        mask->data.U8[300] = 1;
+
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 10000; i++)
+        {
+            psVectorStats (stats, rnd, NULL, mask, 1);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 0.18, "sample min,max %f,%f (mask: 1, range: 0): %.3f sec", stats->min, stats->max, delta);
+        psFree (stats);
+        psFree (mask);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+    // test stat min,max (no mask, no range)
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX | PS_STAT_USE_RANGE);
+        stats->min = -10;
+        stats->max = +10;
+
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 10000; i++)
+        {
+            psVectorStats (stats, rnd, NULL, NULL, 1);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 0.22, "sample min,max %f,%f (mask: 0, range: 1): %.3f sec", stats->min, stats->max, delta);
+        psFree (stats);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+    // test stat min,max (no mask, no range)
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX | PS_STAT_USE_RANGE);
+        stats->min = -10;
+        stats->max = +10;
+        psVector *mask = psVectorAlloc (1000, PS_TYPE_U8);
+        psVectorInit (mask, 0);
+        mask->data.U8[100] = 1;
+        mask->data.U8[200] = 1;
+        mask->data.U8[300] = 1;
+
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 10000; i++)
+        {
+            psVectorStats (stats, rnd, NULL, mask, 1);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 0.26, "sample min,max %f,%f (mask: 1, range: 1): %.3f sec", stats->min, stats->max, delta);
+        psFree (stats);
+        psFree (mask);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+
+//    diag ("timing for clipped stats");
+    /********** CLIPPED STATS ***********/
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
+        psVector *rnd2 = psVectorAlloc (1000, PS_TYPE_F32);
+        for (int i = 0; i < rnd2->n; i++)
+        {
+            rnd2->data.F32[i] = psRandomGaussian (seed);
+        }
+
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 1000; i++)
+        {
+            psVectorStats (stats, rnd2, NULL, NULL, 1);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 0.3, "clipped mean %f, stdev %f (mask: 0, range: 0): %.3f sec (1000 pts / 1000 loops)", stats->clippedMean, stats->clippedStdev, delta);
+        psFree (stats);
+        psFree (rnd2);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
+        psVector *rnd2 = psVectorAlloc (3000, PS_TYPE_F32);
+        for (int i = 0; i < rnd2->n; i++)
+        {
+            rnd2->data.F32[i] = psRandomGaussian (seed);
+        }
+
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 1000; i++)
+        {
+            psVectorStats (stats, rnd2, NULL, NULL, 1);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 0.5, "clipped mean %f, stdev %f (mask: 0, range: 0): %.3f sec (3000 pts / 1000 loops)", stats->clippedMean, stats->clippedStdev, delta);
+        psFree (stats);
+        psFree (rnd2);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
+        psVector *rnd2 = psVectorAlloc (10000, PS_TYPE_F32);
+        for (int i = 0; i < rnd2->n; i++)
+        {
+            rnd2->data.F32[i] = psRandomGaussian (seed);
+        }
+
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 1000; i++)
+        {
+            psVectorStats (stats, rnd2, NULL, NULL, 1);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 1.2, "clipped mean %f, stdev %f (mask: 0, range: 0): %.3f sec (10000 pts / 1000 loops)", stats->clippedMean, stats->clippedStdev, delta);
+        psFree (stats);
+        psFree (rnd2);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+
+//    diag ("timing for robust stats");
+    /********** ROBUST STATS ***********/
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_ROBUST_QUARTILE);
+        psVector *rnd2 = psVectorAlloc (1000, PS_TYPE_F32);
+        for (int i = 0; i < rnd2->n; i++)
+        {
+            rnd2->data.F32[i] = psRandomGaussian (seed);
+        }
+
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 1000; i++)
+        {
+            psVectorStats (stats, rnd2, NULL, NULL, 1);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 0.3, "robust mean %f, stdev %f (mask: 0, range: 0): %.3f sec (1000 pts / 1000 loops)", stats->robustMedian, stats->robustStdev, delta);
+        psFree (stats);
+        psFree (rnd2);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_ROBUST_QUARTILE);
+        psVector *rnd2 = psVectorAlloc (3000, PS_TYPE_F32);
+        for (int i = 0; i < rnd2->n; i++)
+        {
+            rnd2->data.F32[i] = psRandomGaussian (seed);
+        }
+
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 1000; i++)
+        {
+            psVectorStats (stats, rnd2, NULL, NULL, 1);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 0.5, "robust mean %f, stdev %f (mask: 0, range: 0): %.3f sec (3000 pts / 1000 loops)", stats->robustMedian, stats->robustStdev, delta);
+        psFree (stats);
+        psFree (rnd2);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_ROBUST_QUARTILE);
+        psVector *rnd2 = psVectorAlloc (10000, PS_TYPE_F32);
+        for (int i = 0; i < rnd2->n; i++)
+        {
+            rnd2->data.F32[i] = psRandomGaussian (seed);
+        }
+
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 1000; i++)
+        {
+            psVectorStats (stats, rnd2, NULL, NULL, 1);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 1.2, "robust mean %f, stdev %f (mask: 0, range: 0): %.3f sec (10000 pts / 1000 loops)", stats->robustMedian, stats->robustStdev, delta);
+        psFree (stats);
+        psFree (rnd2);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+
+//    diag ("timing for fitted stats");
+    /********** FITTED TIMING ***********/
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV);
+        psVector *rnd2 = psVectorAlloc (1000, PS_TYPE_F32);
+        for (int i = 0; i < rnd2->n; i++)
+        {
+            rnd2->data.F32[i] = psRandomGaussian (seed);
+        }
+
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 1000; i++)
+        {
+            psVectorStats (stats, rnd2, NULL, NULL, 1);
+
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 0.7, "fitted mean %f, stdev %f (mask: 0, range: 0): %.3f sec (1000 pts / 1000 loops)", stats->fittedMean, stats->fittedStdev, delta);
+        psFree (stats);
+        psFree (rnd2);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV);
+        psVector *rnd2 = psVectorAlloc (3000, PS_TYPE_F32);
+        for (int i = 0; i < rnd2->n; i++)
+        {
+            rnd2->data.F32[i] = psRandomGaussian (seed);
+        }
+
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 1000; i++)
+        {
+            psVectorStats (stats, rnd2, NULL, NULL, 1);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 0.8, "fitted mean %f, stdev %f (mask: 0, range: 0): %.3f sec (3000 pts / 1000 loops)", stats->fittedMean, stats->fittedStdev, delta);
+        psFree (stats);
+        psFree (rnd2);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV);
+        psVector *rnd2 = psVectorAlloc (10000, PS_TYPE_F32);
+        for (int i = 0; i < rnd2->n; i++)
+        {
+            rnd2->data.F32[i] = psRandomGaussian (seed);
+        }
+
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 1000; i++)
+        {
+            psVectorStats (stats, rnd2, NULL, NULL, 1);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 2.2, "fitted mean %f, stdev %f (mask: 0, range: 0): %.3f sec (10000 pts / 1000 loops)", stats->fittedMean, stats->fittedStdev, delta);
+        psFree (stats);
+        psFree (rnd2);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+
+//    diag ("timing for fitted (v2) stats");
+    /********** FITTED (v2) TIMING ***********/
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2);
+        psVector *rnd2 = psVectorAlloc (1000, PS_TYPE_F32);
+        for (int i = 0; i < rnd2->n; i++)
+        {
+            rnd2->data.F32[i] = psRandomGaussian (seed);
+        }
+
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 1000; i++)
+        {
+            psVectorStats (stats, rnd2, NULL, NULL, 1);
+
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 0.7, "fitted mean %f, stdev %f (mask: 0, range: 0): %.3f sec (1000 pts / 1000 loops)", stats->fittedMean, stats->fittedStdev, delta);
+        psFree (stats);
+        psFree (rnd2);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2);
+        psVector *rnd2 = psVectorAlloc (3000, PS_TYPE_F32);
+        for (int i = 0; i < rnd2->n; i++)
+        {
+            rnd2->data.F32[i] = psRandomGaussian (seed);
+        }
+
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 1000; i++)
+        {
+            psVectorStats (stats, rnd2, NULL, NULL, 1);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 0.8, "fitted mean %f, stdev %f (mask: 0, range: 0): %.3f sec (3000 pts / 1000 loops)", stats->fittedMean, stats->fittedStdev, delta);
+        psFree (stats);
+        psFree (rnd2);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2);
+        psVector *rnd2 = psVectorAlloc (10000, PS_TYPE_F32);
+        for (int i = 0; i < rnd2->n; i++)
+        {
+            rnd2->data.F32[i] = psRandomGaussian (seed);
+        }
+
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 1000; i++)
+        {
+            psVectorStats (stats, rnd2, NULL, NULL, 1);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 2.2, "fitted mean %f, stdev %f (mask: 0, range: 0): %.3f sec (10000 pts / 1000 loops)", stats->fittedMean, stats->fittedStdev, delta);
+        psFree (stats);
+        psFree (rnd2);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+
+//    diag ("compare sample, robust, and fitted mean and stdev to theoretical");
+    // compare SAMPLE, FITTED, ROBUST mean to theoretical
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV);
+        psVector *sample = psVectorAlloc (1000, PS_TYPE_F32);
+        psVector *robust = psVectorAlloc (1000, PS_TYPE_F32);
+        psVector *fitted = psVectorAlloc (1000, PS_TYPE_F32);
+
+        for (int i = 0; i < 1000; i++)
+        {
+            // generate a new sample
+            for (int j = 0; j < rnd->n; j++) {
+                rnd->data.F32[j] = psRandomGaussian (seed);
+            }
+            // measure the stats
+            psVectorStats (stats, rnd, NULL, NULL, 1);
+            sample->data.F32[i] = stats->sampleMean;
+            robust->data.F32[i] = stats->robustMedian;
+            fitted->data.F32[i] = stats->fittedMean;
+        }
+        psFree (stats);
+
+        stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
+        psVectorStats (stats, sample, NULL, NULL, 1);
+        ok (stats->sampleStdev < 2/sqrt(1000), "sample mean %f, stdev %f (1000 tries)", stats->sampleMean, stats->sampleStdev);
+        psVectorStats (stats, robust, NULL, NULL, 1);
+        ok (stats->sampleStdev < 2/sqrt(1000), "robust mean %f, stdev %f (1000 tries)", stats->sampleMean, stats->sampleStdev);
+        psVectorStats (stats, fitted, NULL, NULL, 1);
+        ok (stats->sampleStdev < 2/sqrt(1000), "fitted mean %f, stdev %f (1000 tries)", stats->sampleMean, stats->sampleStdev);
+        psFree (stats);
+        psFree (sample);
+        psFree (robust);
+        psFree (fitted);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+
+//    diag ("compare sample, robust, and fitted mean and stdev to theoretical");
+    // compare SAMPLE, FITTED_V2, ROBUST mean to theoretical
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2);
+        psVector *sample = psVectorAlloc (1000, PS_TYPE_F32);
+        psVector *robust = psVectorAlloc (1000, PS_TYPE_F32);
+        psVector *fitted = psVectorAlloc (1000, PS_TYPE_F32);
+
+        for (int i = 0; i < 1000; i++)
+        {
+            // generate a new sample
+            for (int j = 0; j < rnd->n; j++) {
+                rnd->data.F32[j] = psRandomGaussian (seed);
+            }
+            // measure the stats
+            psVectorStats (stats, rnd, NULL, NULL, 1);
+            sample->data.F32[i] = stats->sampleMean;
+            robust->data.F32[i] = stats->robustMedian;
+            fitted->data.F32[i] = stats->fittedMean;
+        }
+        psFree (stats);
+
+        stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
+        psVectorStats (stats, sample, NULL, NULL, 1);
+        ok (stats->sampleStdev < 2/sqrt(1000), "sample mean %f, stdev %f (1000 tries)", stats->sampleMean, stats->sampleStdev);
+        psVectorStats (stats, robust, NULL, NULL, 1);
+        ok (stats->sampleStdev < 2/sqrt(1000), "robust mean %f, stdev %f (1000 tries)", stats->sampleMean, stats->sampleStdev);
+        psVectorStats (stats, fitted, NULL, NULL, 1);
+        ok (stats->sampleStdev < 2/sqrt(1000), "fitted mean %f, stdev %f (1000 tries)", stats->sampleMean, stats->sampleStdev);
+        psFree (stats);
+        psFree (sample);
+        psFree (robust);
+        psFree (fitted);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+
+    return exit_status();
+}
+
Index: branches/eam_branches/20090522/psLib/test/optime/tap_psStatsTiming.txt
===================================================================
--- branches/eam_branches/20090522/psLib/test/optime/tap_psStatsTiming.txt	(revision 24557)
+++ branches/eam_branches/20090522/psLib/test/optime/tap_psStatsTiming.txt	(revision 24557)
@@ -0,0 +1,47 @@
+
+Running tap_psStatsTiming on alala (dual AMD Opteron, 64bit, 2.2GHz)
+yields the following timing results:
+
+# timing for sample mean (1000 loops of 10000 pts)
+ok 1 - sample mean 0.009719 (mask: 0, range: 0): 0.072 sec
+ok 3 - sample mean 0.011060 (mask: 1, range: 0): 0.119 sec
+ok 5 - sample mean 0.009719 (mask: 0, range: 1): 0.170 sec
+ok 7 - sample mean 0.011060 (mask: 1, range: 1): 0.198 sec
+
+# timing for sample median (1000 loops of 10000 pts)
+ok 9 - sample median 0.021781 (mask: 0, range: 0): 2.625 sec
+ok 11 - sample median 0.023795 (mask: 1, range: 0): 2.646 sec
+ok 13 - sample median 0.021781 (mask: 0, range: 1): 2.703 sec
+ok 15 - sample median 0.023795 (mask: 1, range: 1): 2.716 sec
+
+# timing for sample stdev (1000 loops of 10000 pts)
+ok 17 - sample stdev 0.964753 (mask: 0, range: 0): 0.193 sec
+ok 19 - sample stdev 0.965887 (mask: 1, range: 0): 0.257 sec
+ok 21 - sample stdev 0.964753 (mask: 0, range: 1): 0.353 sec
+ok 23 - sample stdev 0.965887 (mask: 1, range: 1): 0.401 sec
+
+# timing for sample min,max (1000 loops of 10000 pts)
+ok 25 - sample min,max -3.205688,2.706797 (mask: 0, range: 0): 0.125 sec
+ok 27 - sample min,max -3.205688,2.706797 (mask: 1, range: 0): 0.152 sec
+ok 29 - sample min,max -3.205688,2.706797 (mask: 0, range: 1): 0.201 sec
+ok 31 - sample min,max -3.205688,2.706797 (mask: 1, range: 1): 0.238 sec
+
+# timing for clipped stats
+not ok 33 - clipped mean -0.047714, stdev 0.991979 (mask: 0, range: 0): 0.369 sec (1000 pts / 1000 loops)
+not ok 35 - clipped mean 0.023963, stdev 0.972186 (mask: 0, range: 0): 1.219 sec (3000 pts / 1000 loops)
+not ok 37 - clipped mean -0.007020, stdev 0.985410 (mask: 0, range: 0): 4.883 sec (10000 pts / 1000 loops)
+
+NOTE: these fail because they are being compared to the 'robust' stats
+limits below.  The clipped mean algorithm should not be so slow (and
+apparently non-linear in npts).
+
+# timing for robust stats
+ok 39 - robust mean 0.123348, stdev 1.014896 (mask: 0, range: 0): 0.187 sec (1000 pts / 1000 loops)
+ok 41 - robust mean -0.006812, stdev 0.974468 (mask: 0, range: 0): 0.382 sec (3000 pts / 1000 loops)
+ok 43 - robust mean -0.013591, stdev 1.001539 (mask: 0, range: 0): 1.076 sec (10000 pts / 1000 loops)
+
+# timing for fitted stats
+ok 45 - fitted mean -0.029859, stdev 0.982947 (mask: 0, range: 0): 0.381 sec (1000 pts / 1000 loops)
+ok 47 - fitted mean 0.014660, stdev 0.956168 (mask: 0, range: 0): 0.727 sec (3000 pts / 1000 loops)
+ok 49 - fitted mean -0.008402, stdev 1.001366 (mask: 0, range: 0): 1.914 sec (10000 pts / 1000 loops)
+