Index: branches/tap_branches/psLib/share/Makefile.am =================================================================== --- branches/tap_branches/psLib/share/Makefile.am (revision 25900) +++ branches/tap_branches/psLib/share/Makefile.am (revision 27838) @@ -7,5 +7,4 @@ tab5.2b.dat \ tab5.2c.dat \ - finals2000A.dat \ - finals_all.dat + finals2000A.dat Index: branches/tap_branches/psLib/src/fits/psFits.c =================================================================== --- branches/tap_branches/psLib/src/fits/psFits.c (revision 25900) +++ branches/tap_branches/psLib/src/fits/psFits.c (revision 27838) @@ -149,5 +149,5 @@ int thisErrno = errno; // Error number char errorBuf[MAX_STRING_LENGTH], *errorMsg; -#if ((_POSIX_C_SOURCE >= 200112L || _XOPEN_SOURCE >= 600) && ! _GNU_SOURCE) +#if (((_POSIX_C_SOURCE >= 200112L || _XOPEN_SOURCE >= 600) && ! _GNU_SOURCE) || __APPLE__) strerror_r(thisErrno, errorBuf, MAX_STRING_LENGTH); errorMsg = errorBuf; @@ -208,5 +208,5 @@ int thisErrno = errno; char errorBuf[64], *errorMsg; -# if ((_POSIX_C_SOURCE >= 200112L || _XOPEN_SOURCE >= 600) && ! _GNU_SOURCE) +#if (((_POSIX_C_SOURCE >= 200112L || _XOPEN_SOURCE >= 600) && ! _GNU_SOURCE) || __APPLE__) strerror_r (errno, errorBuf, 64); errorMsg = errorBuf; Index: branches/tap_branches/psLib/src/fits/psFitsHeader.c =================================================================== --- branches/tap_branches/psLib/src/fits/psFitsHeader.c (revision 25900) +++ branches/tap_branches/psLib/src/fits/psFitsHeader.c (revision 27838) @@ -570,4 +570,10 @@ // to preserve NAXISn etc for reference, so we don't do this. + if (keywordInList(name, noWriteFitsKeys) || + (keyStarts && keywordStartsWith(name, noWriteFitsKeyStarts))) { + psTrace("psLib.fits", 3, "Not writing FITS keyword %s", name); + continue; + } + // Options for compression if (compressing) { @@ -582,10 +588,4 @@ } else if (keywordInList(name, noWriteCompressedKeys) || (keyStarts && keywordStartsWith(name, noWriteCompressedKeyStarts))) { - psTrace("psLib.fits", 3, "Not writing FITS keyword %s", name); - continue; - } - - if (keywordInList(name, noWriteFitsKeys) || - (keyStarts && keywordStartsWith(name, noWriteFitsKeyStarts))) { psTrace("psLib.fits", 3, "Not writing FITS keyword %s", name); continue; Index: branches/tap_branches/psLib/src/fits/psFitsImage.c =================================================================== --- branches/tap_branches/psLib/src/fits/psFitsImage.c (revision 25900) +++ branches/tap_branches/psLib/src/fits/psFitsImage.c (revision 27838) @@ -498,6 +498,11 @@ psImage *inImage; // Image to read in if (floatType == PS_FITS_FLOAT_NONE) { - inImage = psImageRecycle(outImage, numCols, numRows, info->psDatatype); - outImage = psMemIncrRefCounter(inImage); + if (!outImage || outImage->type.type == info->psDatatype) { + outImage = psImageRecycle(outImage, numCols, numRows, info->psDatatype); + inImage = psMemIncrRefCounter(outImage); + } else { + outImage = psImageRecycle(outImage, numCols, numRows, outImage->type.type); + inImage = psImageAlloc(numCols, numRows, info->psDatatype); + } } else { inImage = psImageAlloc(numCols, numRows, info->psDatatype); @@ -511,9 +516,11 @@ return NULL; } - psFree(info); if (floatType != PS_FITS_FLOAT_NONE) { outImage = psFitsFloatImageFromDisk(outImage, inImage, floatType); - } + } else if (outImage->type.type != info->psDatatype) { + outImage = psImageCopy(outImage, inImage, outImage->type.type); + } + psFree(info); psFree(inImage); Index: branches/tap_branches/psLib/src/fits/psFitsScale.c =================================================================== --- branches/tap_branches/psLib/src/fits/psFitsScale.c (revision 25900) +++ branches/tap_branches/psLib/src/fits/psFitsScale.c (revision 27838) @@ -15,6 +15,7 @@ #include "psImage.h" #include "psFits.h" +#include "psStats.h" +#include "psImageStats.h" #include "psImageBackground.h" -#include "psStats.h" #include "psImageStructManip.h" @@ -29,4 +30,7 @@ #define MEAN_STAT PS_STAT_ROBUST_MEDIAN // Statistic to use for mean #define STDEV_STAT PS_STAT_ROBUST_STDEV // Statistic to use for stdev + +#define DESPERATE_MEAN_STAT PS_STAT_SAMPLE_MEDIAN // Statistic to use for mean when deperate +#define DESPERATE_STDEV_STAT PS_STAT_SAMPLE_QUARTILE // Statistic to use for stdev when desperate @@ -118,14 +122,76 @@ psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); psStats *stats = psStatsAlloc(MEAN_STAT | STDEV_STAT); // Statistics object + double mean, stdev; // Mean and standard deviation if (!psImageBackground(stats, NULL, image, mask, maskVal, rng)) { - psError(PS_ERR_UNKNOWN, false, "Unable to perform statistics on image"); - psFree(rng); - psFree(stats); - return false; + // It could be because the image is entirely masked, in which case we don't want to error + bool good = false; // Any good pixels? + + +// Find good pixels in an image, by image type +#define GOOD_PIXELS_CASE(TYPE) \ + case PS_TYPE_##TYPE: \ + for (int y = 0; y < image->numRows && !good; y++) { \ + for (int x = 0; x < image->numCols && !good; x++) { \ + if (mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal)) { \ + continue; \ + } \ + if (!isfinite(image->data.TYPE[y][x])) { \ + continue; \ + } \ + good = true; \ + } \ + } \ + break; + + switch (image->type.type) { + GOOD_PIXELS_CASE(F32); + GOOD_PIXELS_CASE(F64); + default: + psAbort("Unsupported case: %x", image->type.type); + } + + if (!good) { + psLogMsg("psLib.fits", PS_LOG_DETAIL, "Image has no good pixels, setting BSCALE = 1, BZERO = 0"); + psErrorClear(); + *bscale = 1.0; + *bzero = 0.0; + psFree(rng); + psFree(stats); + return true; + } + + // There are some good pixels in there somewhere; psImageBackground just didn't find them + psLogMsg("psLib.fits", PS_LOG_DETAIL, + "Couldn't measure background statistics for image quantisation; retrying."); + psErrorClear(); + // Retry using all the available pixels + stats->nSubsample = image->numCols * image->numRows + 1; + if (!psImageStats(stats, image, mask, maskVal)) { + psLogMsg("psLib.fits", PS_LOG_DETAIL, + "Couldn't measure background statistics for image quantisation (attempt 2); retrying."); + psErrorClear(); + // Retry with desperate statistic + stats->options = DESPERATE_MEAN_STAT | DESPERATE_STDEV_STAT; + if (!psImageStats(stats, image, mask, maskVal)) { + psError(PS_ERR_UNKNOWN, false, "Unable to measure background statistics for image"); + psFree(rng); + psFree(stats); + return false; + } else { + // Desperate retry + mean = psStatsGetValue(stats, DESPERATE_MEAN_STAT); + stdev = psStatsGetValue(stats, DESPERATE_STDEV_STAT); + } + } else { + // Retry with all available pixels + mean = psStatsGetValue(stats, MEAN_STAT); + stdev = psStatsGetValue(stats, STDEV_STAT); + } + } else { + // First attempt + mean = psStatsGetValue(stats, MEAN_STAT); + stdev = psStatsGetValue(stats, STDEV_STAT); } psFree(rng); - - double mean = psStatsGetValue(stats, MEAN_STAT); // Mean - double stdev = psStatsGetValue(stats, STDEV_STAT); // Standard deviation psFree(stats); if (!isfinite(mean) || !isfinite(stdev)) { @@ -318,5 +384,5 @@ } else { \ value = (value - zero) * scale; \ - if (options->fuzz) { \ + if (options->fuzz && (value - (int)value != 0.0)) { \ /* Add random factor [-0.5,0.5): adds a variance of 1/12, */ \ /* but preserves the expectation value */ \ Index: branches/tap_branches/psLib/src/fits/psFitsTable.c =================================================================== --- branches/tap_branches/psLib/src/fits/psFitsTable.c (revision 25900) +++ branches/tap_branches/psLib/src/fits/psFitsTable.c (revision 27838) @@ -162,8 +162,9 @@ switch (typecode) { + // TBYTE and TSHORT fall though to read into psS32 case TBYTE: case TSHORT: - case TLONGLONG: READ_TABLE_ROW_CASE(TLONG, long, S32, S32); + READ_TABLE_ROW_CASE(TLONGLONG, long, S64, S64); READ_TABLE_ROW_CASE(TFLOAT, float, F32, F32); READ_TABLE_ROW_CASE(TDOUBLE, double, F64, F64); Index: branches/tap_branches/psLib/src/imageops/psImageBackground.c =================================================================== --- branches/tap_branches/psLib/src/imageops/psImageBackground.c (revision 25900) +++ branches/tap_branches/psLib/src/imageops/psImageBackground.c (revision 27838) @@ -15,4 +15,11 @@ #include "psRandom.h" #include "psError.h" + +static int nFailures = 0; + +void psImageBackgroundInit() { + nFailures = 0; + return; +} // XXX allow the user to choose the stats method? @@ -39,5 +46,18 @@ long ny = image->numRows; - const int Npixels = nx*ny; // Total number of pixels + psImage *bad = psImageAlloc(nx, ny, PS_TYPE_U8); // Image with bad pixels + psImageInit(bad, 0); + + int Npixels = 0; // Total number of pixels + for (int y = 0; y < ny; y++) { + for (int x = 0; x < nx; x++) { + if (!isfinite(image->data.F32[y][x]) || + (mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskValue)) { + bad->data.U8[y][x] = 0xFF; + } + Npixels++; + } + } + const int Nsubset = (stats->nSubsample == 0) ? Npixels : PS_MIN(stats->nSubsample, Npixels); // Number of pixels in subset @@ -58,39 +78,48 @@ long n = 0; // Number of actual pixels in subset 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 we have an image smaller than Nsubset, just loop over the (good) image pixels + for (int iy = 0; iy < ny; iy++) { + for (int ix = 0; ix < nx; ix++) { + if (bad->data.U8[iy][ix]) { + continue; + } - float value = image->data.F32[iy][ix]; - min = PS_MIN(value, min); - max = PS_MAX(value, max); - values->data.F32[n] = value; - n++; - } - } + 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; + // Subsample all pixels + // This is not optimal, since there may be a large masked fraction that leaves us with few good pixels. + // In that case, you should have set Nsubset....... + Npixels = nx * ny; + for (long i = 0; i < Nsubset; i++) { + double frnd = psRandomUniform(rng); + int pixel = Npixels * frnd; + int ix = pixel % nx; + int iy = pixel / nx; - if (!isfinite(image->data.F32[iy][ix]) || (mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskValue)) { - continue; - } + if (bad->data.U8[iy][ix]) { + continue; + } - float value = image->data.F32[iy][ix]; - min = PS_MIN(value, min); - max = PS_MAX(value, max); - values->data.F32[n] = value; - n++; - } + float value = image->data.F32[iy][ix]; + min = PS_MIN(value, min); + max = PS_MAX(value, max); + values->data.F32[n] = value; + n++; + } } + + psFree(bad); + if (n < 0.01*Nsubset) { - psLogMsg("psLib.psImageBackground", PS_LOG_DETAIL, - "Unable to measure image background: too few data points (%ld)", n); + if ((nFailures < 3) || (nFailures % 100 == 0)) { + psLogMsg("psLib.psImageBackground", PS_LOG_DETAIL, "Unable to measure image background: too few data points (%ld) (%d failures)", n, nFailures); + } + nFailures ++; psFree(values); return false; @@ -108,5 +137,8 @@ if (!psVectorSort(values, values)) { - psError(PS_ERR_UNKNOWN, false, "Unable to sort values.\n"); + if ((nFailures < 3) || (nFailures % 100 == 0)) { + psError(PS_ERR_UNKNOWN, false, "Unable to sort values.(%d failures)\n", nFailures); + } + nFailures ++; psFree(values); return false; @@ -132,7 +164,10 @@ fclose (f); } - psError(PS_ERR_UNKNOWN, false, "Unable to measure statistics for image background " - "(%dx%d, (row0,col0) = (%d,%d)", - image->numRows, image->numCols, image->row0, image->col0); + if ((nFailures < 3) || (nFailures % 100 == 0)) { + psError(PS_ERR_UNKNOWN, false, "Unable to measure statistics for image background " + "(%dx%d, (row0,col0) = (%d,%d) (%d failures)\n", + image->numRows, image->numCols, image->row0, image->col0, nFailures); + } + nFailures ++; psFree(values); return false; Index: branches/tap_branches/psLib/src/imageops/psImageBackground.h =================================================================== --- branches/tap_branches/psLib/src/imageops/psImageBackground.h (revision 25900) +++ branches/tap_branches/psLib/src/imageops/psImageBackground.h (revision 27838) @@ -22,4 +22,6 @@ #include +void psImageBackgroundInit(); + // Get the background for an image bool psImageBackground(psStats *stats, // desired measurement and options Index: branches/tap_branches/psLib/src/imageops/psImageConvolve.c =================================================================== --- branches/tap_branches/psLib/src/imageops/psImageConvolve.c (revision 25900) +++ branches/tap_branches/psLib/src/imageops/psImageConvolve.c (revision 27838) @@ -246,4 +246,84 @@ return kernel; } + +bool psKernelTruncate(psKernel *kernel, float frac) +{ + PS_ASSERT_KERNEL_NON_NULL(kernel, false); + PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(frac, 0.0, false); + PS_ASSERT_FLOAT_LESS_THAN(frac, 1.0, false); + + if (frac == 0.0) { + // Nothing to do + return true; + } + + int xMin = kernel->xMin, xMax = kernel->xMax, yMin = kernel->yMin, yMax = kernel->yMax; // Bounds + int maxSize = PS_MAX(PS_MAX(PS_MAX(-xMin, xMax), -yMin), yMax); // Maximum size + + // Determine the threshold + // Summing absolute values because large negative deviations have power as well + double sumKernel = 0.0; // Sum of the kernel + for (int y = yMin; y <= yMax; y++) { + for (int x = xMin; x <= xMax; x++) { + sumKernel += fabsf(kernel->kernel[y][x]); + } + } + + float threshold = sumKernel * (1.0 - frac); // Threshold for truncation + + // Find truncation size + int truncate = 0; // Truncation radius + for (int radius = 1; truncate == 0 && radius < maxSize; radius++) { + int uMin = PS_MAX(-radius, xMin); + int uMax = PS_MIN(radius, xMax); + int vMin = PS_MAX(-radius, yMin); + int vMax = PS_MIN(radius, yMax); + int r2 = PS_SQR(radius); + double sum = 0.0; + for (int v = vMin; v <= vMax; v++) { + int v2 = PS_SQR(v); + for (int u = uMin; u <= uMax; u++) { + int u2 = PS_SQR(u); + if (u2 + v2 <= r2) { + sum += fabsf(kernel->kernel[v][u]); + } + } + } + if (sum > threshold) { + // This is the truncation radius + truncate = radius; + } + } + if (truncate == maxSize) { + // No truncation possible + return true; + } + + // Truncate the kernel + { + int uMin = PS_MAX(-truncate, xMin); + int uMax = PS_MIN(truncate, xMax); + int vMin = PS_MAX(-truncate, yMin); + int vMax = PS_MIN(truncate, yMax); + int r2 = PS_SQR(truncate); + for (int v = vMin; v <= vMax; v++) { + int v2 = PS_SQR(v); + for (int u = uMin; u <= uMax; u++) { + int u2 = PS_SQR(u); + if (u2 + v2 > r2) { + kernel->kernel[v][u] = 0.0; + } + } + } + } + kernel->xMin = PS_MAX(-truncate, kernel->xMin); + kernel->xMax = PS_MIN(truncate, kernel->xMax); + kernel->yMin = PS_MAX(-truncate, kernel->yMin); + kernel->yMax = PS_MIN(truncate, kernel->yMax); + + return true; + } + + psImage *psImageConvolveDirect(psImage *out, const psImage *in, const psKernel *kernel) Index: branches/tap_branches/psLib/src/imageops/psImageConvolve.h =================================================================== --- branches/tap_branches/psLib/src/imageops/psImageConvolve.h (revision 25900) +++ branches/tap_branches/psLib/src/imageops/psImageConvolve.h (revision 27838) @@ -138,4 +138,14 @@ ); +/// Truncate a kernel +/// +/// Truncates the outer parts of the kernel where the contribution is below the nominated fraction of the +/// total kernel. +bool psKernelTruncate( + psKernel *in, ///< Kernel to be truncated + float frac ///< Fraction for truncation threshold + ); + + /// Convolve an image with a kernel, using a direct convolution /// Index: branches/tap_branches/psLib/src/imageops/psImageCovariance.c =================================================================== --- branches/tap_branches/psLib/src/imageops/psImageCovariance.c (revision 25900) +++ branches/tap_branches/psLib/src/imageops/psImageCovariance.c (revision 27838) @@ -14,6 +14,11 @@ #include "psTrace.h" #include "psBinaryOp.h" +#include "psScalar.h" +#include "psThread.h" #include "psImageCovariance.h" + +static bool threaded = false; // Run threaded? + psKernel *psImageCovarianceNone(void) @@ -24,54 +29,12 @@ } - -psKernel *psImageCovarianceCalculate(const psKernel *kernel, const psKernel *covariance) -{ - PS_ASSERT_KERNEL_NON_NULL(kernel, NULL); - - // See http://en.wikipedia.org/wiki/Error_propagation - // - // If - // f_k = sum_i A_ik x_i - // is a set of functions, then the covariance matrix for f is given by: - // M^f_ij = sum_k sum_l A_ik M^x_kl A_lj - // where M^x is the covariance matrix for x. - // Note that the errors in f are correlated (covariance) even if the errors in x are not. - - psKernel *covar; // Covariance matrix to use - if (covariance) { - covar = psMemIncrRefCounter((psKernel*)covariance); // Casting away const - } else { - covar = psImageCovarianceNone(); - } - - // Check for non-finite elements - for (int y = kernel->yMin; y <= kernel->yMax; y++) { - for (int x = kernel->xMin; x <= kernel->xMax; x++) { - if (!isfinite(kernel->kernel[y][x])) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, - "Non-finite covariance matrix element in kernel at %d,%d", x, y); - psFree(covar); - return NULL; - } - } - } - for (int y = covar->yMin; y <= covar->yMax; y++) { - for (int x = covar->xMin; x <= covar->xMax; x++) { - if (!isfinite(covar->kernel[y][x])) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, - "Non-finite covariance matrix element in covariance matrix at %d,%d", x, y); - psFree(covar); - return NULL; - } - } - } - - // The above (double) sum for the covariance matrix means that, for each point in the output covariance - // matrix, we need to work out all combinations of getting to the central point via a kernel, input - // covariance matrix and another kernel. This means that the resultant covariance matrix has twice the - // size of the kernel plus the size of the input covariance matrix. - int xMin = kernel->xMin - kernel->xMax + covar->xMin, xMax = kernel->xMax - kernel->xMin + covar->xMax; - int yMin = kernel->yMin - kernel->yMax + covar->yMin, yMax = kernel->yMax - kernel->yMin + covar->yMax; - psKernel *out = psKernelAlloc(xMin, xMax, yMin, yMax); // Covariance matrix for output +/// Calculation of covariance matrix element when convolving +static float imageCovarianceCalculate(const psKernel *covar, // Original covariance matrix + const psKernel *kernel, // Convolution kernel + int x, int y // Coordinates in output covariance matrix + ) +{ + psAssert(covar, "Require covariance matrix"); + psAssert(kernel, "Require kernel"); // Need to go: @@ -89,44 +52,234 @@ // from the source coordinate), we take the smallest possible (because everything else is zero outside). - double total = 0.0; // Total covariance + // Range for v + int vMin = PS_MAX(kernel->yMin + covar->yMin, y + kernel->yMin); + int vMax = PS_MIN(kernel->yMax + covar->yMax, y + kernel->yMax); + // Range for u + int uMin = PS_MAX(kernel->xMin + covar->xMin, x + kernel->xMin); + int uMax = PS_MIN(kernel->xMax + covar->xMax, x + kernel->xMax); + + double sum = 0.0; // Sum for value of covariance matrix at (x,y) + for (int v = vMin; v <= vMax; v++) { + // Range for q + int qMin = PS_MAX(v + covar->yMin, kernel->yMin); + int qMax = PS_MIN(v + covar->yMax, kernel->yMax); + for (int u = uMin; u <= uMax; u++) { + // Range for p + int pMin = PS_MAX(u + covar->xMin, kernel->xMin); + int pMax = PS_MIN(u + covar->xMax, kernel->xMax); + + double xyuvValue = kernel->kernel[v-y][u-x]; // Value for (x,y) --> (u,v) + + double uvpqValue = 0.0; // Value for (u,v) --> (p,q) --> (0,0) + for (int q = qMin; q <= qMax; q++) { + for (int p = pMin; p <= pMax; p++) { + uvpqValue += (double)covar->kernel[q-v][p-u] * (double)kernel->kernel[q][p]; + } + } + sum += xyuvValue * uvpqValue; + } + } + + return sum; +} + +/// Thread entry point for calculation of covariance matrix element when convolving +static bool imageCovarianceCalculateThread(psThreadJob *job) +{ + PS_ASSERT_THREAD_JOB_NON_NULL(job, false); + psAssert(job->args, "No job arguments"); + psAssert(job->args->n == 5, "Wrong number of job arguments: %ld", job->args->n); + + psKernel *out = job->args->data[0]; // Output covariance matrix + const psKernel *covar = job->args->data[1]; // Input covariance matrix + const psKernel *kernel = job->args->data[2]; // Convolution kernel + int x = PS_SCALAR_VALUE(job->args->data[3], S32); // x coordinate in output covariance matrix + int y = PS_SCALAR_VALUE(job->args->data[4], S32); // y coordinate in output covariance matrix + + out->kernel[y][x] = imageCovarianceCalculate(covar, kernel, x, y); + + return true; +} + + + +psKernel *psImageCovarianceCalculate(const psKernel *kernel, const psKernel *covariance) +{ + PS_ASSERT_KERNEL_NON_NULL(kernel, NULL); + + // See http://en.wikipedia.org/wiki/Error_propagation + // + // If + // f_k = sum_i A_ik x_i + // is a set of functions, then the covariance matrix for f is given by: + // M^f_ij = sum_k sum_l A_ik M^x_kl A_lj + // where M^x is the covariance matrix for x. + // Note that the errors in f are correlated (covariance) even if the errors in x are not. + + psKernel *covar; // Covariance matrix to use + if (covariance) { + covar = psMemIncrRefCounter((psKernel*)covariance); // Casting away const + } else { + covar = psImageCovarianceNone(); + } + + // Check for non-finite elements + for (int y = kernel->yMin; y <= kernel->yMax; y++) { + for (int x = kernel->xMin; x <= kernel->xMax; x++) { + if (!isfinite(kernel->kernel[y][x])) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + "Non-finite covariance matrix element in kernel at %d,%d", x, y); + psFree(covar); + return NULL; + } + } + } + for (int y = covar->yMin; y <= covar->yMax; y++) { + for (int x = covar->xMin; x <= covar->xMax; x++) { + if (!isfinite(covar->kernel[y][x])) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + "Non-finite covariance matrix element in covariance matrix at %d,%d", x, y); + psFree(covar); + return NULL; + } + } + } + + // The above (double) sum for the covariance matrix means that, for each point in the output covariance + // matrix, we need to work out all combinations of getting to the central point via a kernel, input + // covariance matrix and another kernel. This means that the resultant covariance matrix has twice the + // size of the kernel plus the size of the input covariance matrix. + int xMin = kernel->xMin - kernel->xMax + covar->xMin, xMax = kernel->xMax - kernel->xMin + covar->xMax; + int yMin = kernel->yMin - kernel->yMax + covar->yMin, yMax = kernel->yMax - kernel->yMin + covar->yMax; + psKernel *out = psKernelAlloc(xMin, xMax, yMin, yMax); // Covariance matrix for output + for (int y = yMin; y <= yMax; y++) { - // Range for v - int vMin = PS_MAX(kernel->yMin + covar->yMin, y + kernel->yMin); - int vMax = PS_MIN(kernel->yMax + covar->yMax, y + kernel->yMax); for (int x = xMin; x <= xMax; x++) { - // Range for u - int uMin = PS_MAX(kernel->xMin + covar->xMin, x + kernel->xMin); - int uMax = PS_MIN(kernel->xMax + covar->xMax, x + kernel->xMax); - - double sum = 0.0; // Sum for value of covariance matrix at (x,y) - for (int v = vMin; v <= vMax; v++) { - // Range for q - int qMin = PS_MAX(v + covar->yMin, kernel->yMin); - int qMax = PS_MIN(v + covar->yMax, kernel->yMax); - for (int u = uMin; u <= uMax; u++) { - // Range for p - int pMin = PS_MAX(u + covar->xMin, kernel->xMin); - int pMax = PS_MIN(u + covar->xMax, kernel->xMax); - - double xyuvValue = kernel->kernel[v-y][u-x]; // Value for (x,y) --> (u,v) - - double uvpqValue = 0.0; // Value for (u,v) --> (p,q) --> (0,0) - for (int q = qMin; q <= qMax; q++) { - for (int p = pMin; p <= pMax; p++) { - uvpqValue += (double)covar->kernel[q-v][p-u] * (double)kernel->kernel[q][p]; - } - } - sum += xyuvValue * uvpqValue; + if (threaded) { + psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_COVARIANCE_CALCULATE"); + psArrayAdd(job->args, 1, out); + psArrayAdd(job->args, 1, covar); + psArrayAdd(job->args, 1, (psKernel*)kernel); // Casting away const + PS_ARRAY_ADD_SCALAR(job->args, x, PS_TYPE_S32); + PS_ARRAY_ADD_SCALAR(job->args, y, PS_TYPE_S32); + if (!psThreadJobAddPending(job)) { + psFree(covar); + return NULL; } - } - out->kernel[y][x] = sum; - total += sum; - } - } - psTrace("psLib.imageops", 3, "Total covariance: %lf ; Central variance: %f\n", total, out->kernel[0][0]); - + psFree(job); + } else { + out->kernel[y][x] = imageCovarianceCalculate(covar, kernel, x, y); + } + } + } psFree(covar); + + if (threaded && !psThreadPoolWait(true)) { + psError(PS_ERR_UNKNOWN, false, "Error waiting for threads."); + return false; + } + return out; } + +float psImageCovarianceCalculateFactor(const psKernel *kernel, const psKernel *covariance) +{ + psKernel *covar; // Covariance matrix to use + if (covariance) { + covar = psMemIncrRefCounter((psKernel*)covariance); // Casting away const + } else { + covar = psImageCovarianceNone(); + } + + // Check for non-finite elements + for (int y = kernel->yMin; y <= kernel->yMax; y++) { + for (int x = kernel->xMin; x <= kernel->xMax; x++) { + if (!isfinite(kernel->kernel[y][x])) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + "Non-finite covariance matrix element in kernel at %d,%d", x, y); + psFree(covar); + return NAN; + } + } + } + for (int y = covar->yMin; y <= covar->yMax; y++) { + for (int x = covar->xMin; x <= covar->xMax; x++) { + if (!isfinite(covar->kernel[y][x])) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + "Non-finite covariance matrix element in covariance matrix at %d,%d", x, y); + psFree(covar); + return NAN; + } + } + } + + float factor = imageCovarianceCalculate(covar, kernel, 0, 0); // Covariance factor + psFree(covar); + return factor; +} + +// Calculation of covariance matrix element when binning +static float imageCovarianceBin(const psKernel *covar, // Original covariance matrix + int bin, // Binning factor + float binVal, // Convolution kernel value for binning + int x, int y // Coordinates in output covariance matrix + ) +{ + psAssert(covar, "Require covariance matrix"); + psAssert(bin > 0 && binVal > 0, "Require binning: %d %f", bin, binVal); + + int binMin = -(bin - 1) / 2, binMax = bin / 2; // Range of "kernel" + + // Range for v + int vMin = PS_MAX(binMin + covar->yMin, bin * y + binMin); + int vMax = PS_MIN(binMax + covar->yMax, bin * y + binMax); + // Range for u + int uMin = PS_MAX(binMin + covar->xMin, bin * x + binMin); + int uMax = PS_MIN(binMax + covar->xMax, bin * x + binMax); + + double sum = 0.0; // Sum for value of covariance matrix at (x,y) + for (int v = vMin; v <= vMax; v++) { + // Range for q + int qMin = PS_MAX(v + covar->yMin, binMin); + int qMax = PS_MIN(v + covar->yMax, binMax); + for (int u = uMin; u <= uMax; u++) { + // Range for p + int pMin = PS_MAX(u + covar->xMin, binMin); + int pMax = PS_MIN(u + covar->xMax, binMax); + + double xyuvValue = binVal; // Value for (x,y) --> (u,v) + + double uvpqValue = 0.0; // Value for (u,v) --> (p,q) --> (0,0) + for (int q = qMin; q <= qMax; q++) { + for (int p = pMin; p <= pMax; p++) { + uvpqValue += (double)covar->kernel[q-v][p-u] * (double)binVal; + } + } + sum += xyuvValue * uvpqValue; + } + } + + return sum; +} + +/// Thread entry point for calculation of covariance matrix element when binning +static bool imageCovarianceBinThread(psThreadJob *job) +{ + PS_ASSERT_THREAD_JOB_NON_NULL(job, false); + psAssert(job->args, "No job arguments"); + psAssert(job->args->n == 6, "Wrong number of job arguments: %ld", job->args->n); + + psKernel *out = job->args->data[0]; // Output covariance matrix + const psKernel *covar = job->args->data[1]; // Input covariance matrix + int bin = PS_SCALAR_VALUE(job->args->data[2], S32); // Binning factor + float binVal = PS_SCALAR_VALUE(job->args->data[3], F32); // Convolution kernel value for binning + int x = PS_SCALAR_VALUE(job->args->data[4], S32); // x coordinate in output covariance matrix + int y = PS_SCALAR_VALUE(job->args->data[5], S32); // y coordinate in output covariance matrix + + out->kernel[y][x] = imageCovarianceBin(covar, bin, binVal, x, y); + + return true; +} + psKernel *psImageCovarianceBin(int bin, const psKernel *covariance, bool average) @@ -165,43 +318,31 @@ psKernel *out = psKernelAlloc(xMin, xMax, yMin, yMax); // Covariance matrix for output - double total = 0.0; // Total covariance for (int y = yMin; y <= yMax; y++) { - // Range for v - int vMin = PS_MAX(binMin + covar->yMin, bin * y + binMin); - int vMax = PS_MIN(binMax + covar->yMax, bin * y + binMax); for (int x = xMin; x <= xMax; x++) { - // Range for u - int uMin = PS_MAX(binMin + covar->xMin, bin * x + binMin); - int uMax = PS_MIN(binMax + covar->xMax, bin * x + binMax); - - double sum = 0.0; // Sum for value of covariance matrix at (x,y) - for (int v = vMin; v <= vMax; v++) { - // Range for q - int qMin = PS_MAX(v + covar->yMin, binMin); - int qMax = PS_MIN(v + covar->yMax, binMax); - for (int u = uMin; u <= uMax; u++) { - // Range for p - int pMin = PS_MAX(u + covar->xMin, binMin); - int pMax = PS_MIN(u + covar->xMax, binMax); - - double xyuvValue = binVal; // Value for (x,y) --> (u,v) - - double uvpqValue = 0.0; // Value for (u,v) --> (p,q) --> (0,0) - for (int q = qMin; q <= qMax; q++) { - for (int p = pMin; p <= pMax; p++) { - uvpqValue += (double)covar->kernel[q-v][p-u] * (double)binVal; - } - } - sum += xyuvValue * uvpqValue; + if (threaded) { + psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_COVARIANCE_BIN"); + psArrayAdd(job->args, 1, out); + psArrayAdd(job->args, 1, covar); + PS_ARRAY_ADD_SCALAR(job->args, bin, PS_TYPE_S32); + PS_ARRAY_ADD_SCALAR(job->args, binVal, PS_TYPE_F32); + PS_ARRAY_ADD_SCALAR(job->args, x, PS_TYPE_S32); + PS_ARRAY_ADD_SCALAR(job->args, y, PS_TYPE_S32); + if (!psThreadJobAddPending(job)) { + psFree(covar); + return NULL; } - } - out->kernel[y][x] = sum; - total += sum; - } - } - psTrace("psLib.imageops", 3, "Total covariance: %lf ; Central variance: %f\n", total, out->kernel[0][0]); - + psFree(job); + } else { + out->kernel[y][x] = imageCovarianceBin(covar, bin, binVal, x, y); + } + } + } psFree(covar); + if (threaded && !psThreadPoolWait(true)) { + psError(PS_ERR_UNKNOWN, false, "Error waiting for threads."); + return false; + } + return out; } @@ -212,72 +353,53 @@ } -psKernel *psImageCovarianceSum(const psArray *array) +float psImageCovarianceFactorForAperture(const psKernel *covar, float radius) +{ + if (!covar) return 1.0; + + float Sum = 0.0; + + for (int y = covar->yMin; y <= covar->yMax; y++) { + if (y < -radius) continue; + if (y > +radius) continue; + for (int x = covar->xMin; x <= covar->xMax; x++) { + if (x < -radius) continue; + if (x > +radius) continue; + + if (hypot(x, y) > radius) continue; + + psAssert (isfinite(covar->kernel[y][x]), "invalid NAN in covariance matrix"); + Sum += covar->kernel[y][x]; + } + } + + return Sum; +} + +float psImageCovarianceSum(const psKernel *covar) +{ + PS_ASSERT_KERNEL_NON_NULL(covar, NAN); + + int xMin = covar->xMin, xMax = covar->xMax, yMin = covar->yMin, yMax = covar->yMax; // Range for covariance + double sum = 0.0; // Sum of covariance + for (int y = yMin; y <= yMax; y++) { + for (int x = xMin; x <= xMax; x++) { + sum += covar->kernel[y][x]; + } + } + + return sum; +} + + +psKernel *psImageCovarianceAverage(const psArray *array) { PS_ASSERT_ARRAY_NON_NULL(array, NULL); PS_ASSERT_ARRAY_NON_EMPTY(array, NULL); - int xMin = INT_MAX, xMax = INT_MIN, yMin = INT_MAX, yMax = INT_MIN; // Range for covariance - int num = 0; // Number of good matrices to sum - for (int i = 0; i < array->n; i++) { - psKernel *covar = array->data[i]; // Covariance matrix - if (!covar) { - continue; - } - xMin = PS_MIN(xMin, covar->xMin); - xMax = PS_MAX(xMax, covar->xMax); - yMin = PS_MIN(yMin, covar->yMin); - yMax = PS_MAX(yMax, covar->yMax); - num++; - } - if (num == 0) { - psError(PS_ERR_BAD_PARAMETER_SIZE, true, "No covariance matrices supplied for summation"); - return NULL; - } - - psKernel *sum = psKernelAlloc(xMin, xMax, yMin, yMax); // Summed covariance - for (int i = 0; i < array->n; i++) { - psKernel *covar = array->data[i]; // Covariance matrix - if (!covar) { - continue; - } - for (int y = covar->yMin; y <= covar->yMax; y++) { - for (int x = covar->xMin; x <= covar->xMax; x++) { - if (!isfinite(covar->kernel[y][x])) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, - "Non-finite covariance matrix element at %d,%d for input %d", - x, y, i); - psFree(sum); - return NULL; - } - sum->kernel[y][x] += covar->kernel[y][x]; - } - } - } - - return sum; -} - - -psKernel *psImageCovarianceAverage(const psArray *array) -{ - PS_ASSERT_ARRAY_NON_NULL(array, NULL); - PS_ASSERT_ARRAY_NON_EMPTY(array, NULL); - - int num = 0; // Number of good matrices to average - for (int i = 0; i < array->n; i++) { - psKernel *covar = array->data[i]; // Covariance matrix - if (covar) { - num++; - } - } - if (num == 0) { - psError(PS_ERR_BAD_PARAMETER_SIZE, true, "No covariance matrices supplied for averaging."); - return NULL; - } - - psKernel *sum = psImageCovarianceSum(array); // Sum of covariances - psBinaryOp(sum->image, sum->image, "/", psScalarAlloc(num, PS_TYPE_F32)); - - return sum; + psVector *weights = psVectorAlloc(array->n, PS_TYPE_F32); // Weights to apply + psVectorInit(weights, 1.0); + psKernel *out = psImageCovarianceAverageWeighted(array, weights); + psFree(weights); + return out; } @@ -310,4 +432,5 @@ psKernel *sum = psKernelAlloc(xMin, xMax, yMin, yMax); // Summed covariance + psImageInit(sum->image, 0.0); for (int i = 0; i < array->n; i++) { psKernel *covar = array->data[i]; // Covariance matrix @@ -405,2 +528,32 @@ return true; } + + +bool psImageCovarianceSetThreads(bool set) +{ + bool old = threaded; // Old value + if (set && !threaded) { + { + psThreadTask *task = psThreadTaskAlloc("PSLIB_IMAGE_COVARIANCE_CALCULATE", 5); + task->function = &imageCovarianceCalculateThread; + psThreadTaskAdd(task); + psFree(task); + } + { + psThreadTask *task = psThreadTaskAlloc("PSLIB_IMAGE_COVARIANCE_BIN", 6); + task->function = &imageCovarianceBinThread; + psThreadTaskAdd(task); + psFree(task); + } + } else if (!set && threaded) { + psThreadTaskRemove("PSLIB_IMAGE_COVARIANCE_CALCULATE"); + psThreadTaskRemove("PSLIB_IMAGE_COVARIANCE_BIN"); + } + threaded = set; + return old; +} + +bool psImageCovarianceGetThreads(void) +{ + return threaded; +} Index: branches/tap_branches/psLib/src/imageops/psImageCovariance.h =================================================================== --- branches/tap_branches/psLib/src/imageops/psImageCovariance.h (revision 25900) +++ branches/tap_branches/psLib/src/imageops/psImageCovariance.h (revision 27838) @@ -47,7 +47,19 @@ ); -/// Sum multiple covariance pseudo-matrices -psKernel *psImageCovarianceSum( - const psArray *array ///< Array of covariance pseudo-matrices +/// Return the pixel-to-pixel covariance factor following calculation +/// +/// This doesn't require calculation of the entire covariance matrix, so is much faster. +float psImageCovarianceCalculateFactor( + const psKernel *kernel, ///< Convolution kernel + const psKernel *covariance ///< Current covariance pseudo-matrix + ); + + +/// Return the covariance factor for an aperture of a given radius +float psImageCovarianceFactorForAperture(const psKernel *covar, float radius); + +/// Return the sum of the covariance pseudo-matrix +float psImageCovarianceSum( + const psKernel *covariance ///< Covariance pseudo-matrix ); @@ -78,4 +90,12 @@ ); +/// Control threading for image covariance functions +/// +/// Returns old threading status +bool psImageCovarianceSetThreads(bool threaded ///< Run image covariance functions threaded? + ); + +/// Return whether image covariance functions are threaded +bool psImageCovarianceGetThreads(void); /// @} Index: branches/tap_branches/psLib/src/imageops/psImageInterpolate.c =================================================================== --- branches/tap_branches/psLib/src/imageops/psImageInterpolate.c (revision 25900) +++ branches/tap_branches/psLib/src/imageops/psImageInterpolate.c (revision 27838) @@ -197,10 +197,10 @@ { // Casting away const - psFree((psImage*)interp->image); - psFree((psImage*)interp->mask); - psFree((psImage*)interp->variance); - psFree((psImage*)interp->kernel); - psFree((psImage*)interp->kernel2); - psFree((psVector*)interp->sumKernel2); + psFree(interp->image); + psFree(interp->mask); + psFree(interp->variance); + psFree(interp->kernel); + psFree(interp->kernel2); + psFree(interp->sumKernel2); } Index: branches/tap_branches/psLib/src/jpeg/psImageJpeg.c =================================================================== --- branches/tap_branches/psLib/src/jpeg/psImageJpeg.c (revision 25900) +++ branches/tap_branches/psLib/src/jpeg/psImageJpeg.c (revision 27838) @@ -202,19 +202,27 @@ if (isfinite(row[i])) { pixel = PS_JPEG_SCALEVALUE(row[i],zero,scale); - outPix[0] = Rpix[pixel]; - outPix[1] = Gpix[pixel]; - outPix[2] = Bpix[pixel]; + outPix[0] = Rpix[pixel]; + outPix[1] = Gpix[pixel]; + outPix[2] = Bpix[pixel]; } else { - // XXX NAN value should be set per-color map - outPix[0] = 0xff; - outPix[1] = 0x00; - outPix[2] = 0xff; - } - } - jpeg_write_scanlines(&cinfo, jpegLineList, 1); + // XXX NAN value should be set per-color map + outPix[0] = 0xff; + outPix[1] = 0x00; + outPix[2] = 0xff; + } + } + if (jpeg_write_scanlines(&cinfo, jpegLineList, 1) == 0) { + psError(PS_ERR_IO, true, "Unable to write line %d to JPEG", j); + psFree(jpegLine); + return false; + } } jpeg_finish_compress(&cinfo); - fclose(f); + if (fclose(f) == EOF) { + psError(PS_ERR_IO, true, "Failed to close %s", filename); + psFree(jpegLine); + return false; + } jpeg_destroy_compress(&cinfo); Index: branches/tap_branches/psLib/src/math/psBinaryOp.c =================================================================== --- branches/tap_branches/psLib/src/math/psBinaryOp.c (revision 25900) +++ branches/tap_branches/psLib/src/math/psBinaryOp.c (revision 27838) @@ -379,5 +379,5 @@ } -psMathType* psBinaryOp(psPtr out, const psPtr in1, const char *op, const psPtr in2) +psMathType* psBinaryOp(psPtr out, psPtr in1, const char *op, psPtr in2) { Index: branches/tap_branches/psLib/src/math/psBinaryOp.h =================================================================== --- branches/tap_branches/psLib/src/math/psBinaryOp.h (revision 25900) +++ branches/tap_branches/psLib/src/math/psBinaryOp.h (revision 27838) @@ -55,7 +55,7 @@ psMathType* psBinaryOp( psPtr out, ///< Output type, either psImage or psVector. - const psPtr in1, ///< First input, either psImage or psVector. + psPtr in1, ///< First input, either psImage or psVector. const char *op, ///< Operator. - const psPtr in2 ///< Second input, either psImage or psVector. + psPtr in2 ///< Second input, either psImage or psVector. ); Index: branches/tap_branches/psLib/src/math/psHistogram.c =================================================================== --- branches/tap_branches/psLib/src/math/psHistogram.c (revision 25900) +++ branches/tap_branches/psLib/src/math/psHistogram.c (revision 27838) @@ -58,6 +58,6 @@ psTrace("psLib.math", 3, "---- %s() begin ----\n", __func__); psTrace("psLib.math", 5, "(lower, upper, n) is (%f, %f, %d)\n", lower, upper, n); - PS_ASSERT_INT_POSITIVE(n, NULL); - PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(upper, lower, NULL); + psAssert(n > 0, "Number of bins must be positive"); + psAssert(upper >= lower, "Bounds must be sensical"); // Allocate memory for the new histogram structure. If there are N bins, then there are N+1 bounds to @@ -127,5 +127,5 @@ static void histogramFree(psHistogram* myHist) { - psFree((void *)myHist->bounds); + psFree(myHist->bounds); psFree(myHist->nums); } @@ -287,23 +287,23 @@ double binSize = (out->bounds->data.F32[out->nums->n] - out->bounds->data.F32[0]) / (float) out->nums->n; // Histogram bin size binNum = (inF32->data.F32[i] - out->bounds->data.F32[0]) / binSize; - binNum = PS_MAX (binNum, 0); - binNum = PS_MIN (binNum, numBins - 1); - - // value is in bin 'i' if bound[i] <= value < bound[i] - - // we may slightly overshoot or undershoot. creep up or down on the true bin - if (inF32->data.F32[i] < out->bounds->data.F32[binNum]) { - psTrace("psLib.math", 6, "missed target bin, adjusting: %f vs %f to %f\n", inF32->data.F32[i], out->bounds->data.F32[binNum], out->bounds->data.F32[binNum+1]); - while ((inF32->data.F32[i] < out->bounds->data.F32[binNum]) && (binNum > 0)) { - binNum --; - } - - } - if (inF32->data.F32[i] >= out->bounds->data.F32[binNum+1]) { - psTrace("psLib.math", 6, "missed target bin, adjusting: %f vs %f to %f\n", inF32->data.F32[i], out->bounds->data.F32[binNum], out->bounds->data.F32[binNum+1]); - while ((inF32->data.F32[i] >= out->bounds->data.F32[binNum+1]) && (binNum < numBins - 1)) { - binNum ++; - } - } + binNum = PS_MAX (binNum, 0); + binNum = PS_MIN (binNum, numBins - 1); + + // value is in bin 'i' if bound[i] <= value < bound[i] + + // we may slightly overshoot or undershoot. creep up or down on the true bin + if (inF32->data.F32[i] < out->bounds->data.F32[binNum]) { + psTrace("psLib.math", 6, "missed target bin, adjusting: %f vs %f to %f\n", inF32->data.F32[i], out->bounds->data.F32[binNum], out->bounds->data.F32[binNum+1]); + while ((inF32->data.F32[i] < out->bounds->data.F32[binNum]) && (binNum > 0)) { + binNum --; + } + + } + if (inF32->data.F32[i] >= out->bounds->data.F32[binNum+1]) { + psTrace("psLib.math", 6, "missed target bin, adjusting: %f vs %f to %f\n", inF32->data.F32[i], out->bounds->data.F32[binNum], out->bounds->data.F32[binNum+1]); + while ((inF32->data.F32[i] >= out->bounds->data.F32[binNum+1]) && (binNum < numBins - 1)) { + binNum ++; + } + } if (errorsF32) { @@ -326,17 +326,17 @@ // correct bin number requires a bit more work. tmpScalar.data.F32 = inF32->data.F32[i]; - psVectorBinaryDisectResult result; + psVectorBinaryDisectResult result; binNum = psVectorBinaryDisect(&result, out->bounds, &tmpScalar); - if (result != PS_BINARY_DISECT_PASS) { - continue; - } - if (errorsF32 != NULL) { - if (!UpdateHistogramBins(binNum, out, inF32->data.F32[i], errors->data.F32[i])) { - psLogMsg(__func__, PS_LOG_WARN, "WARNING: Failed to update the histogram " - "bins with the errors vector.\n"); - } - } else { - out->nums->data.F32[binNum] += 1.0; - } + if (result != PS_BINARY_DISECT_PASS) { + continue; + } + if (errorsF32 != NULL) { + if (!UpdateHistogramBins(binNum, out, inF32->data.F32[i], errors->data.F32[i])) { + psLogMsg(__func__, PS_LOG_WARN, "WARNING: Failed to update the histogram " + "bins with the errors vector.\n"); + } + } else { + out->nums->data.F32[binNum] += 1.0; + } } } Index: branches/tap_branches/psLib/src/math/psMatrix.c =================================================================== --- branches/tap_branches/psLib/src/math/psMatrix.c (revision 25900) +++ branches/tap_branches/psLib/src/math/psMatrix.c (revision 27838) @@ -94,102 +94,109 @@ LHS_NAME.data = RHS_NAME; - -/*****************************************************************************/ -/* FILE STATIC FUNCTIONS */ -/*****************************************************************************/ - -static void psVectorToGslVector(gsl_vector *outGslVector, const psVector *inVector); -static void gslVectorToPsVector(psVector *outVector, gsl_vector *inGslVector); -static void psImageToGslMatrix(gsl_matrix *outGslMatrix, const psImage *inImage); -static void gslMatrixToPsImage(psImage *outImage, gsl_matrix *inGslMatrix); +//////////////////////////////////////////////////////////////////////////////// +// Conversion functions +//////////////////////////////////////////////////////////////////////////////// + +// gsl_vector holds *doubles*, so we can directly copy F64, but need to convert F32 /** Static function to copy psF32 or psF64 vector data to a GSL vector */ -static void psVectorToGslVector(gsl_vector *outGslVector, - const psVector *inVector) -{ - psU32 i = 0; - psU32 n = 0; - - - n = inVector->n; - for(i=0; itype.type == PS_TYPE_F32) { - outGslVector->data[i] = (psF64)inVector->data.F32[i]; +static void vectorPStoGSL(gsl_vector *out, const psVector *in) +{ + psAssert(out->size == in->n, "Sizes don't match!"); + + long n = in->n; // Size of input + switch (in->type.type) { + case PS_TYPE_F32: + for (long i = 0; i < n; i++) { + out->data[i] = in->data.F32[i]; + } + break; + case PS_TYPE_F64: + memcpy(out->data, in->data.F64, n * PSELEMTYPE_SIZEOF(PS_TYPE_F64)); + break; + default: + psAbort("Unsupported vector type: %x\n", in->type.type); + } + return; +} + +/** Static function to copy GSL vector data to a psF32 or psF64 vector */ +static void vectorGSLtoPS(psVector *out, const gsl_vector *in) +{ + psAssert(in->size == out->n, "Sizes don't match!"); + + long n = out->n; // Size of output + switch (out->type.type) { + case PS_TYPE_F32: + for (long i = 0; i < n; i++) { + out->data.F32[i] = in->data[i]; + } + break; + case PS_TYPE_F64: + memcpy(out->data.F64, in->data, n * PSELEMTYPE_SIZEOF(PS_TYPE_F64)); + break; + default: + psAbort("Unsupported vector type: %x\n", out->type.type); + } + return; +} + + +/** Static function to copy psF32 or psF64 image data to a GSL matrix */ +static void matrixPStoGSL(gsl_matrix *out, const psImage *in) +{ + psAssert(out->size1 == in->numRows && out->size2 == in->numCols, "Sizes don't match!"); + + int numCols = in->numCols, numRows = in->numRows; // Size of matrix + switch (in->type.type) { + case PS_TYPE_F32: + for (int y = 0, i = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++, i++) { + out->data[i] = in->data.F32[y][x]; + } + } + break; + case PS_TYPE_F64: + if (in->parent|| out->tda != out->size1) { + for (int y = 0, i = 0; y < numRows; y++, i += numCols) { + memcpy(&out->data[i], in->data.F64[y], numCols * PSELEMTYPE_SIZEOF(PS_TYPE_F64)); + } } else { - outGslVector->data[i] = inVector->data.F64[i]; - } - } -} - -/** Static function to copy GSL vector data to a psF32 or psF64 vector */ -static void gslVectorToPsVector(psVector *outVector, - gsl_vector *inGslVector) -{ - psU32 i = 0; - psU32 n = 0; - - - n = outVector->n; - for(i=0; itype.type == PS_TYPE_F32) { - outVector->data.F32[i] = (psF32)inGslVector->data[i]; + memcpy(out->data, in->p_rawDataBuffer, numCols * numRows * PSELEMTYPE_SIZEOF(PS_TYPE_F64)); + } + break; + default: + psAbort("Unsupported vector type: %x\n", in->type.type); + } + return; +} + +/** Static function to copy GSL matrix data to a psF32 or psF64 image */ +static void matrixGSLtoPS(psImage *out, const gsl_matrix *in) +{ + psAssert(in->size1 == out->numRows && in->size2 == out->numCols, "Sizes don't match!"); + + int numCols = out->numCols, numRows = out->numRows; // Size of matrix + switch (out->type.type) { + case PS_TYPE_F32: + for (int y = 0, i = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++, i++) { + out->data.F32[y][x] = in->data[i]; + } + } + break; + case PS_TYPE_F64: + if (out->parent || in->tda != in->size1) { + for (int y = 0, i = 0; y < numRows; y++, i += numCols) { + memcpy(out->data.F64[y], &in->data[i], numCols * PSELEMTYPE_SIZEOF(PS_TYPE_F64)); + } } else { - outVector->data.F64[i] = inGslVector->data[i]; - } - } -} - -/** Static function to copy psF32 or psF64 image data to a GSL matrix */ -static void psImageToGslMatrix(gsl_matrix *outGslMatrix, - const psImage *inImage) -{ - psU32 i = 0; - psU32 j = 0; - psU32 numRows = 0; - psU32 numCols = 0; - - - numRows = inImage->numRows; - numCols = inImage->numCols; - if(inImage->type.type == PS_TYPE_F32) { - for(i=0; idata[i*numCols+j] = inImage->data.F32[i][j]; - } - } - } else { - for(i=0; idata[i*numCols+j] = inImage->data.F64[i][j]; - } - } - } -} - -/** Static function to copy GSL matrix data to a psF32 or psF64 image */ -static void gslMatrixToPsImage(psImage *outImage, - gsl_matrix *inGslMatrix) -{ - psU32 i = 0; - psU32 j = 0; - psU32 numRows = 0; - psU32 numCols = 0; - - - numRows = outImage->numRows; - numCols = outImage->numCols; - if(outImage->type.type == PS_TYPE_F32) { - for(i=0; idata.F32[i][j] = inGslMatrix->data[i*numCols+j]; - } - } - } else { - for(i=0; idata.F64[i][j] = inGslMatrix->data[i*numCols+j]; - } - } - } + memcpy(out->p_rawDataBuffer, in->data, numCols * numRows * PSELEMTYPE_SIZEOF(PS_TYPE_F64)); + } + break; + default: + psAbort("Unsupported vector type: %x\n", out->type.type); + } + return; } @@ -245,5 +252,5 @@ // Copy psImage data into GSL matrix data - psImageToGslMatrix(lu, in); + matrixPStoGSL(lu, in); // Calculate LU decomposition @@ -251,5 +258,5 @@ // Copy GSL matrix data to psImage data - gslMatrixToPsImage(out, lu); + matrixGSLtoPS(out, lu); // Free GSL data @@ -293,7 +300,7 @@ // Initialize GSL data lu = gsl_matrix_alloc(numRows, numCols); - psImageToGslMatrix(lu, LU); + matrixPStoGSL(lu, LU); b = gsl_vector_alloc(RHS->n); - psVectorToGslVector(b, RHS); + vectorPStoGSL(b, RHS); x = gsl_vector_alloc(RHS->n); @@ -306,5 +313,5 @@ // Copy GSL vector data to psVector data - gslVectorToPsVector(out, x); + vectorGSLtoPS(out, x); // Free GSL data @@ -341,5 +348,5 @@ // Initialize GSL data lu = gsl_matrix_alloc(numRows, numCols); - psImageToGslMatrix(lu, LU); + matrixPStoGSL(lu, LU); permGSL.size = perm->n; @@ -352,5 +359,5 @@ // Copy GSL vector data to psVector data - gslMatrixToPsImage(out, inverse); + matrixGSLtoPS(out, inverse); // Free GSL data @@ -379,37 +386,37 @@ switch (a->type.type) { case PS_TYPE_F32: { - psF32 **values = a->data.F32; /* Dereference */ - int numCols = a->numCols, numRows = a->numRows; /* Size of matrix */ - for (int i = 0; i < numRows; i++) { - for (int j = 0; j < numCols; j++) { - if (!isfinite(values[i][j])) { - // psError(PS_ERR_BAD_PARAMETER_VALUE, 3, - // "Input matrix contains non-finite elements: matrix[%d][%d] is %.2f\n", - // i, j, values[i][j]); - return false; - } - } - } - break; + psF32 **values = a->data.F32; /* Dereference */ + int numCols = a->numCols, numRows = a->numRows; /* Size of matrix */ + for (int i = 0; i < numRows; i++) { + for (int j = 0; j < numCols; j++) { + if (!isfinite(values[i][j])) { + // psError(PS_ERR_BAD_PARAMETER_VALUE, 3, + // "Input matrix contains non-finite elements: matrix[%d][%d] is %.2f\n", + // i, j, values[i][j]); + return false; + } + } + } + break; } case PS_TYPE_F64: { - psF64 **values = a->data.F64; /* Dereference */ - int numCols = a->numCols, numRows = a->numRows; /* Size of matrix */ - for (int i = 0; i < numRows; i++) { - for (int j = 0; j < numCols; j++) { - if (!isfinite(values[i][j])) { - // psError(PS_ERR_BAD_PARAMETER_VALUE, 3, - // "Input matrix contains non-finite elements: matrix[%d][%d] is %.2f\n", - // i, j, values[i][j]); - return false; - } - } - } - break; + psF64 **values = a->data.F64; /* Dereference */ + int numCols = a->numCols, numRows = a->numRows; /* Size of matrix */ + for (int i = 0; i < numRows; i++) { + for (int j = 0; j < numCols; j++) { + if (!isfinite(values[i][j])) { + // psError(PS_ERR_BAD_PARAMETER_VALUE, 3, + // "Input matrix contains non-finite elements: matrix[%d][%d] is %.2f\n", + // i, j, values[i][j]); + return false; + } + } + } + break; } - // MATRIX_CHECK_NONFINITE_CASE(F32, a); - // MATRIX_CHECK_NONFINITE_CASE(F64, a); + // MATRIX_CHECK_NONFINITE_CASE(F32, a); + // MATRIX_CHECK_NONFINITE_CASE(F64, a); default: - psAbort("Should never get here."); + psAbort("Should never get here."); } @@ -471,31 +478,31 @@ switch (a->type.type) { case PS_TYPE_F32: { - psF32 **values = a->data.F32; /* Dereference */ - int numCols = a->numCols, numRows = a->numRows; /* Size of matrix */ - for (int i = 0; i < numRows; i++) { - for (int j = 0; j < numCols; j++) { - if (!isfinite(values[i][j])) { - return false; - } - } - } - break; + psF32 **values = a->data.F32; /* Dereference */ + int numCols = a->numCols, numRows = a->numRows; /* Size of matrix */ + for (int i = 0; i < numRows; i++) { + for (int j = 0; j < numCols; j++) { + if (!isfinite(values[i][j])) { + return false; + } + } + } + break; } case PS_TYPE_F64: { - psF64 **values = a->data.F64; /* Dereference */ - int numCols = a->numCols, numRows = a->numRows; /* Size of matrix */ - for (int i = 0; i < numRows; i++) { - for (int j = 0; j < numCols; j++) { - if (!isfinite(values[i][j])) { - return false; - } - } - } - break; + psF64 **values = a->data.F64; /* Dereference */ + int numCols = a->numCols, numRows = a->numRows; /* Size of matrix */ + for (int i = 0; i < numRows; i++) { + for (int j = 0; j < numCols; j++) { + if (!isfinite(values[i][j])) { + return false; + } + } + } + break; } default: - psAbort("Should never get here."); - } - + psAbort("Should never get here."); + } + // Following the algorithm laid out by Press et al., we loop along the matrix diagonal, but // we do not operate on the diagonal elements in order. Instead, we are looking for the @@ -511,149 +518,149 @@ if (a->type.type == PS_TYPE_F32) { - psF32 **A = a->data.F32; - psF32 *B = b->data.F32; - int *colIndex = colIndexV->data.S32; - int *rowIndex = rowIndexV->data.S32; - int *pivot = pivotV->data.S32; - psF32 growth = 1.0; - - for (int diag = 0; diag < nSquare; diag++) { - - psF32 maxval = 0.0; - int maxrow = 0; - int maxcol = 0; - - // search for the next pivot - for (int row = 0; row < nSquare; row++) { - if (!isfinite(A[row][diag])) goto escape; - - // if we have already operated on this row (pivot[row] is true), skip it - if (pivot[row]) continue; - - // if we have not yet operated on this row (pivot[row] is false), look for pivot for this row - for (int col = 0; col < nSquare; col++) { - if (pivot[col]) continue; - if (fabs (A[row][col]) < maxval) continue; - maxval = fabs (A[row][col]); - maxrow = row; - maxcol = col; - } - } - - // if pivot[maxcol] is set, we have already done this row: this implies a singular matrix - if (pivot[maxcol]) goto escape; - pivot[maxcol] = 1; - - // if the selected pivot is off the diagonal, do a row swap - if (maxrow != maxcol) { - for (int col = 0; col < nSquare; col++) PS_SWAP (A[maxrow][col], A[maxcol][col]); - PS_SWAP (B[maxrow], B[maxcol]); - } - rowIndex[diag] = maxrow; - colIndex[diag] = maxcol; - if (A[maxcol][maxcol] == 0.0) goto escape; - // Kahan replaces the 0.0 pivot with epsilon*(largest element in column) + underFlow. - // Here we are going to raise an error if the dynamic range is too large. - - /* rescale by pivot reciprocal */ - psF32 tmpval = 1.0 / A[maxcol][maxcol]; - A[maxcol][maxcol] = 1.0; - for (int col = 0; col < nSquare; col++) A[maxcol][col] *= tmpval; - B[maxcol] *= tmpval; - - // check for ill-conditioned matrix: measure the pivot growth and trigger on over/under flow - growth *= tmpval; - psTrace ("psLib.math", 4, "growth : %e\n", growth); - if (fabs(growth) > MAX_RANGE) goto escape; - - /* adjust the elements above the pivot */ - for (int row = 0; row < nSquare; row++) { - if (row == maxcol) continue; - tmpval = A[row][maxcol]; - A[row][maxcol] = 0.0; - for (int col = 0; col < nSquare; col++) A[row][col] -= A[maxcol][col]*tmpval; - B[row] -= B[maxcol]*tmpval; - } - } - - // swap back the inverse matrix based on the row swaps above - for (int col = nSquare - 1; col >= 0; col--) { - if (rowIndex[col] != colIndex[col]) { - for (int row = 0; row < nSquare; row++) PS_SWAP (A[row][rowIndex[col]], A[row][colIndex[col]]); - } - } + psF32 **A = a->data.F32; + psF32 *B = b->data.F32; + int *colIndex = colIndexV->data.S32; + int *rowIndex = rowIndexV->data.S32; + int *pivot = pivotV->data.S32; + psF32 growth = 1.0; + + for (int diag = 0; diag < nSquare; diag++) { + + psF32 maxval = 0.0; + int maxrow = 0; + int maxcol = 0; + + // search for the next pivot + for (int row = 0; row < nSquare; row++) { + if (!isfinite(A[row][diag])) goto escape; + + // if we have already operated on this row (pivot[row] is true), skip it + if (pivot[row]) continue; + + // if we have not yet operated on this row (pivot[row] is false), look for pivot for this row + for (int col = 0; col < nSquare; col++) { + if (pivot[col]) continue; + if (fabs (A[row][col]) < maxval) continue; + maxval = fabs (A[row][col]); + maxrow = row; + maxcol = col; + } + } + + // if pivot[maxcol] is set, we have already done this row: this implies a singular matrix + if (pivot[maxcol]) goto escape; + pivot[maxcol] = 1; + + // if the selected pivot is off the diagonal, do a row swap + if (maxrow != maxcol) { + for (int col = 0; col < nSquare; col++) PS_SWAP (A[maxrow][col], A[maxcol][col]); + PS_SWAP (B[maxrow], B[maxcol]); + } + rowIndex[diag] = maxrow; + colIndex[diag] = maxcol; + if (A[maxcol][maxcol] == 0.0) goto escape; + // Kahan replaces the 0.0 pivot with epsilon*(largest element in column) + underFlow. + // Here we are going to raise an error if the dynamic range is too large. + + /* rescale by pivot reciprocal */ + psF32 tmpval = 1.0 / A[maxcol][maxcol]; + A[maxcol][maxcol] = 1.0; + for (int col = 0; col < nSquare; col++) A[maxcol][col] *= tmpval; + B[maxcol] *= tmpval; + + // check for ill-conditioned matrix: measure the pivot growth and trigger on over/under flow + growth *= tmpval; + psTrace ("psLib.math", 4, "growth : %e\n", growth); + if (fabs(growth) > MAX_RANGE) goto escape; + + /* adjust the elements above the pivot */ + for (int row = 0; row < nSquare; row++) { + if (row == maxcol) continue; + tmpval = A[row][maxcol]; + A[row][maxcol] = 0.0; + for (int col = 0; col < nSquare; col++) A[row][col] -= A[maxcol][col]*tmpval; + B[row] -= B[maxcol]*tmpval; + } + } + + // swap back the inverse matrix based on the row swaps above + for (int col = nSquare - 1; col >= 0; col--) { + if (rowIndex[col] != colIndex[col]) { + for (int row = 0; row < nSquare; row++) PS_SWAP (A[row][rowIndex[col]], A[row][colIndex[col]]); + } + } } else { - psF64 **A = a->data.F64; - psF64 *B = b->data.F64; - int *colIndex = colIndexV->data.S32; - int *rowIndex = rowIndexV->data.S32; - int *pivot = pivotV->data.S32; - psF64 growth = 1.0; - - for (int diag = 0; diag < nSquare; diag++) { - - psF64 maxval = 0.0; - int maxrow = 0; - int maxcol = 0; - - // search for the next pivot - for (int row = 0; row < nSquare; row++) { - if (!isfinite(A[row][diag])) goto escape; - - // if we have already operated on this row (pivot[row] is true), skip it - if (pivot[row]) continue; - - // if we have not yet operated on this row (pivot[row] is false), look for pivot for this row - for (int col = 0; col < nSquare; col++) { - if (pivot[col]) continue; - if (fabs (A[row][col]) < maxval) continue; - maxval = fabs (A[row][col]); - maxrow = row; - maxcol = col; - } - } - - // if pivot[maxcol] is set, we have already done this row: this implies a singular matrix - if (pivot[maxcol]) goto escape; - pivot[maxcol] = 1; - - // if the selected pivot is off the diagonal, do a row swap - if (maxrow != maxcol) { - for (int col = 0; col < nSquare; col++) PS_SWAP (A[maxrow][col], A[maxcol][col]); - PS_SWAP (B[maxrow], B[maxcol]); - } - rowIndex[diag] = maxrow; - colIndex[diag] = maxcol; - if (A[maxcol][maxcol] == 0.0) goto escape; - // Kahan replaces the 0.0 pivot with epsilon*(largest element in column) + underFlow. - // Here we are going to raise an error if the dynamic range is too large. - - /* rescale by pivot reciprocal */ - psF64 tmpval = 1.0 / A[maxcol][maxcol]; - A[maxcol][maxcol] = 1.0; - for (int col = 0; col < nSquare; col++) A[maxcol][col] *= tmpval; - B[maxcol] *= tmpval; - - // check for ill-conditioned matrix: measure the pivot growth and trigger on over/under flow - growth *= tmpval; - psTrace ("psLib.math", 4, "growth : %e\n", growth); - if (fabs(growth) > MAX_RANGE) goto escape; - - /* adjust the elements above the pivot */ - for (int row = 0; row < nSquare; row++) { - if (row == maxcol) continue; - tmpval = A[row][maxcol]; - A[row][maxcol] = 0.0; - for (int col = 0; col < nSquare; col++) A[row][col] -= A[maxcol][col]*tmpval; - B[row] -= B[maxcol]*tmpval; - } - } - - // swap back the inverse matrix based on the row swaps above - for (int col = nSquare - 1; col >= 0; col--) { - if (rowIndex[col] != colIndex[col]) { - for (int row = 0; row < nSquare; row++) PS_SWAP (A[row][rowIndex[col]], A[row][colIndex[col]]); - } - } + psF64 **A = a->data.F64; + psF64 *B = b->data.F64; + int *colIndex = colIndexV->data.S32; + int *rowIndex = rowIndexV->data.S32; + int *pivot = pivotV->data.S32; + psF64 growth = 1.0; + + for (int diag = 0; diag < nSquare; diag++) { + + psF64 maxval = 0.0; + int maxrow = 0; + int maxcol = 0; + + // search for the next pivot + for (int row = 0; row < nSquare; row++) { + if (!isfinite(A[row][diag])) goto escape; + + // if we have already operated on this row (pivot[row] is true), skip it + if (pivot[row]) continue; + + // if we have not yet operated on this row (pivot[row] is false), look for pivot for this row + for (int col = 0; col < nSquare; col++) { + if (pivot[col]) continue; + if (fabs (A[row][col]) < maxval) continue; + maxval = fabs (A[row][col]); + maxrow = row; + maxcol = col; + } + } + + // if pivot[maxcol] is set, we have already done this row: this implies a singular matrix + if (pivot[maxcol]) goto escape; + pivot[maxcol] = 1; + + // if the selected pivot is off the diagonal, do a row swap + if (maxrow != maxcol) { + for (int col = 0; col < nSquare; col++) PS_SWAP (A[maxrow][col], A[maxcol][col]); + PS_SWAP (B[maxrow], B[maxcol]); + } + rowIndex[diag] = maxrow; + colIndex[diag] = maxcol; + if (A[maxcol][maxcol] == 0.0) goto escape; + // Kahan replaces the 0.0 pivot with epsilon*(largest element in column) + underFlow. + // Here we are going to raise an error if the dynamic range is too large. + + /* rescale by pivot reciprocal */ + psF64 tmpval = 1.0 / A[maxcol][maxcol]; + A[maxcol][maxcol] = 1.0; + for (int col = 0; col < nSquare; col++) A[maxcol][col] *= tmpval; + B[maxcol] *= tmpval; + + // check for ill-conditioned matrix: measure the pivot growth and trigger on over/under flow + growth *= tmpval; + psTrace ("psLib.math", 4, "growth : %e\n", growth); + if (fabs(growth) > MAX_RANGE) goto escape; + + /* adjust the elements above the pivot */ + for (int row = 0; row < nSquare; row++) { + if (row == maxcol) continue; + tmpval = A[row][maxcol]; + A[row][maxcol] = 0.0; + for (int col = 0; col < nSquare; col++) A[row][col] -= A[maxcol][col]*tmpval; + B[row] -= B[maxcol]*tmpval; + } + } + + // swap back the inverse matrix based on the row swaps above + for (int col = nSquare - 1; col >= 0; col--) { + if (rowIndex[col] != colIndex[col]) { + for (int row = 0; row < nSquare; row++) PS_SWAP (A[row][rowIndex[col]], A[row][colIndex[col]]); + } + } } @@ -701,5 +708,5 @@ lu = gsl_matrix_alloc(numRows, numCols); inv = gsl_matrix_alloc(numRows, numCols); - psImageToGslMatrix(lu, in); + matrixPStoGSL(lu, in); // Invert data and calculate determinant @@ -708,5 +715,5 @@ if (determinant) { // XXX this is getting the wrong value: is it the wrong calculation? - // it disagrees with the results of + // it disagrees with the results of // det = (psF32)gsl_linalg_LU_det(lu, signum); // used in psMatrixDeterminatn @@ -716,5 +723,5 @@ // Copy GSL matrix data to psImage data - gslMatrixToPsImage(out, inv); + matrixGSLtoPS(out, inv); // Free GSL structs @@ -749,9 +756,9 @@ perm = gsl_permutation_alloc(numRows); lu = gsl_matrix_alloc(numRows, numCols); - psImageToGslMatrix(lu, in); + matrixPStoGSL(lu, in); // Calculate determinant gsl_linalg_LU_decomp(lu, perm, &signum); - det = (psF32)gsl_linalg_LU_det(lu, signum); + det = (psF32)gsl_linalg_LU_lndet(lu); // Free GSL structs @@ -877,5 +884,5 @@ inGSL = gsl_matrix_alloc(numRows, numCols); - psImageToGslMatrix(inGSL, in); + matrixPStoGSL(inGSL, in); outGSL = gsl_matrix_alloc(numRows, numCols); @@ -892,5 +899,5 @@ // Copy GSL matrix data to psImage data - gslMatrixToPsImage(out, outGSL); + matrixGSLtoPS(out, outGSL); // Free GSL structs @@ -1026,6 +1033,91 @@ } +psVector *psMatrixSolveSVD(psVector *out, const psImage *matrix, const psVector *vector, float thresh) +{ + #define psMatrixSolveSVD_EXIT {psFree(out); return NULL; } + PS_ASSERT_GENERAL_IMAGE_NON_NULL(matrix, psMatrixSolveSVD_EXIT); + PS_CHECK_DIMEN_AND_TYPE(matrix, PS_DIMEN_IMAGE, psMatrixSolveSVD_EXIT); + PS_ASSERT_GENERAL_VECTOR_NON_NULL(vector, psMatrixSolveSVD_EXIT); + PS_CHECK_DIMEN_AND_TYPE(vector, PS_DIMEN_VECTOR, psMatrixSolveSVD_EXIT); + + int numCols = matrix->numCols, numRows = matrix->numRows; // Size of matrix + + // Decompose matrix: A = U S V^T + gsl_matrix *A = gsl_matrix_alloc(numRows, numCols); // Input matrix in GSL-speak; becomes matrix U + gsl_matrix *V = gsl_matrix_alloc(numCols, numCols); // Untransposed matrix V + gsl_vector *S = gsl_vector_alloc(numCols); // Singular values + gsl_vector *work = gsl_vector_alloc(numCols); // Work space for GSL + + matrixPStoGSL(A, matrix); + + int gslStatus = 0; // Status of GSL + if ((gslStatus = gsl_linalg_SV_decomp(A, V, S, work))) { + const char *err = gsl_strerror(gslStatus); + psError(PS_ERR_UNKNOWN, true, "Unable to decompose matrix: %s", err); + gsl_matrix_free(A); + gsl_matrix_free(V); + gsl_vector_free(S); + gsl_vector_free(work); + return NULL; + } + gsl_vector_free(work); + + if (isfinite(thresh) && thresh > 0.0) { + // Trim the singular values + double total = 0.0; // Total of singular values + for (int i = 0; i < numCols; i++) { + total += gsl_vector_get(S, i); + } + thresh *= total; + for (int i = 0; i < numCols; i++) { + double value = gsl_vector_get(S, i); // Singular value + if (value < thresh) { + psTrace("psLib.math", 5, "Trimming singular value %d: %lg", i, value); + gsl_vector_set(S, i, 0.0); +#if 0 + for (int j = 0; j < numCols; j++) { + // Being thorough; probably unnecessary + gsl_matrix_set(V, j, i, 0.0); + gsl_matrix_set(A, j, i, 0.0); + } +#endif + } else { + psTrace("psLib.math", 5, "Singular value %d: %lg", i, value); + } + } + } + + // Solve system (or minimise least-squares if overconstrained): Ax = b + gsl_vector *b = gsl_vector_alloc(numCols); // Vector b + gsl_vector *x = gsl_vector_alloc(numCols); // Solution + + vectorPStoGSL(b, vector); + + if ((gslStatus = gsl_linalg_SV_solve(A, V, S, b, x))) { + const char *err = gsl_strerror(gslStatus); + psError(PS_ERR_UNKNOWN, true, "Unable to solve matrix equation: %s", err); + gsl_matrix_free(A); + gsl_matrix_free(V); + gsl_vector_free(S); + gsl_vector_free(b); + gsl_vector_free(x); + return NULL; + } + + gsl_matrix_free(A); + gsl_matrix_free(V); + gsl_vector_free(S); + gsl_vector_free(b); + + out = psVectorRecycle(out, numCols, PS_TYPE_F64); + + vectorGSLtoPS(out, x); + gsl_vector_free(x); + + return out; +} + // This code supplied by Andy Becker (becker@astro.washington.edu) -psImage *psMatrixSVD(psImage* evec, psVector* eval, const psImage* in) +psImage *psMatrixSVD_old(psImage* evec, psVector* eval, const psImage* in) { #define psMatrixSVD_EXIT {psFree(evec); psFree(eval); return NULL;} @@ -1050,5 +1142,5 @@ // Copy psImage data into GSL matrix data - psImageToGslMatrix(A, in); + matrixPStoGSL(A, in); // Calculate SVD decomposition @@ -1056,6 +1148,6 @@ // Copy GSL matrix data to psImage data - gslMatrixToPsImage(evec, V); - gslVectorToPsVector(eval, S); + matrixGSLtoPS(evec, V); + vectorGSLtoPS(eval, S); // Take the square root of eval @@ -1076,2 +1168,56 @@ return evec; } + +// this is basically a wrapper for the gsl function: gsl_linalg_SV_decomp() SVD decomposes +// matrix A based on the following equation: A = U w V^T . This function (as usual for SVD +// implementations) returns V not V^T. U and V are returned to images; w is returned to a +// vector representing the diagonal of w. The input image A is not modified. U, V, and w may +// be supplied as NULL or may be allocated; their lengths are set here to match the +// dimensionality of A. XXX there is no error handling for the gsl functions (anywhere in +// psMatrix.c) +bool psMatrixSVD(psImage **U, psVector **w, psImage **V, const psImage *A) +{ + // Error checks Missing one for eval + PS_ASSERT_PTR_NON_NULL(U, false); + PS_ASSERT_PTR_NON_NULL(w, false); + PS_ASSERT_PTR_NON_NULL(V, false); + PS_ASSERT_PTR_NON_NULL(A, false); + + // A is provided with size Nx,Ny = numCols,numRows + // U has size Nx,Ny + // V has size Nx,Nx + // w has size Nx + + // Initialize data + int numRows = A->numRows; + int numCols = A->numCols; + + *U = psImageRecycle(*U, numCols, numRows, A->type.type); + *V = psImageRecycle(*V, numCols, numCols, A->type.type); + *w = psVectorRecycle(*w, numCols, A->type.type); + + gsl_matrix *Agsl = gsl_matrix_alloc(numRows, numCols); + gsl_matrix *Vgsl = gsl_matrix_alloc(numCols, numCols); + gsl_vector *Sgsl = gsl_vector_alloc(numCols); + gsl_vector *work = gsl_vector_alloc(numCols); + + // Copy psImage data into GSL matrix data + matrixPStoGSL(Agsl, A); + + // Calculate SVD decomposition + gsl_linalg_SV_decomp(Agsl, Vgsl, Sgsl, work); + + // Copy GSL matrix data to psImage data + matrixGSLtoPS(*V, Vgsl); + matrixGSLtoPS(*U, Agsl); // gsl_linalg_SV_decomp replaces A with U + vectorGSLtoPS(*w, Sgsl); + + // Free GSL data + gsl_matrix_free(Agsl); + gsl_matrix_free(Vgsl); + gsl_vector_free(Sgsl); + gsl_vector_free(work); + + return true; +} + Index: branches/tap_branches/psLib/src/math/psMatrix.h =================================================================== --- branches/tap_branches/psLib/src/math/psMatrix.h (revision 25900) +++ branches/tap_branches/psLib/src/math/psMatrix.h (revision 27838) @@ -66,7 +66,7 @@ */ psImage *psMatrixLUInvert( - psImage *out, ///< place result here if not NULL - const psImage* LU, ///< LU-decomposed matrix. - const psVector* perm ///< Permutation vector resulting from psMatrixLUD function. + psImage *out, ///< place result here if not NULL + const psImage* LU, ///< LU-decomposed matrix. + const psVector* perm ///< Permutation vector resulting from psMatrixLUD function. ); @@ -186,6 +186,16 @@ ); -/// Single value decomposition, provided by Andy Becker -psImage *psMatrixSVD(psImage* evec, psVector* eval, const psImage* in); +/// Solve a matrix equation using Singular Value Decomposition +/// +/// Solves Ax = b for x +psVector *psMatrixSolveSVD( + psVector *solution, ///< Solution to output, or NULL + const psImage *matrix, ///< Matrix to be solved + const psVector *vector, ///< Vector of values + float thresh ///< Threshold relative to maximum for trimming singular values + ); + +/// Single value decomposition (original by Andy Becker, updated by EAM) +bool psMatrixSVD(psImage **U, psVector **w, psImage **V, const psImage *A); /// @} Index: branches/tap_branches/psLib/src/math/psMinimizeLMM.c =================================================================== --- branches/tap_branches/psLib/src/math/psMinimizeLMM.c (revision 25900) +++ branches/tap_branches/psLib/src/math/psMinimizeLMM.c (revision 27838) @@ -99,25 +99,25 @@ // XXX check that the GJ solver works: -# if (TESTGJ) +# if (TESTGJ) psImage *out = psImageAlloc (alpha->numRows, alpha->numCols, PS_TYPE_F32); for (int oy = 0; oy < out->numRows; oy++) { - for (int ox = 0; ox < out->numCols; ox++) { - float value = 0; - for (int i = 0; i < alpha->numCols; i++) { - value += alpha->data.F32[i][ox]*Alpha->data.F32[oy][i]; - } - out->data.F32[oy][ox] = value; - } + for (int ox = 0; ox < out->numCols; ox++) { + float value = 0; + for (int i = 0; i < alpha->numCols; i++) { + value += alpha->data.F32[i][ox]*Alpha->data.F32[oy][i]; + } + out->data.F32[oy][ox] = value; + } } psVector *vect = psVectorAlloc (beta->n, PS_TYPE_F32); for (int oy = 0; oy < vect->n; oy++) { - float value = 0; - for (int i = 0; i < alpha->numCols; i++) { - value += alpha->data.F32[oy][i]*Beta->data.F32[i]; - } - vect->data.F32[oy] = value; - } - + float value = 0; + for (int i = 0; i < alpha->numCols; i++) { + value += alpha->data.F32[oy][i]*Beta->data.F32[i]; + } + vect->data.F32[oy] = value; + } + psFree (out); psFree (vect); @@ -223,5 +223,5 @@ if (isnan(chisq)) { psTrace ("psLib.math", 5, "psMinLM_SetABX() returned a NAN chisq.\n"); - psVectorInit (delta, NAN); + psVectorInit (delta, NAN); retValue = false; } @@ -238,5 +238,5 @@ if (!status) { psTrace ("psLib.math", 5, "psMinLM_GuessABP() returned FALSE.\n"); - psVectorInit (delta, NAN); + psVectorInit (delta, NAN); retValue = false; } @@ -301,7 +301,7 @@ PS_ASSERT_VECTOR_NON_NULL(dy, NAN); - PS_ASSERT_VECTOR_TYPE(params, PS_TYPE_F32, false); + PS_ASSERT_VECTOR_TYPE(params, PS_TYPE_F32, NAN); if (paramMask) { - PS_ASSERT_VECTOR_TYPE(paramMask, PS_TYPE_VECTOR_MASK, false); + PS_ASSERT_VECTOR_TYPE(paramMask, PS_TYPE_VECTOR_MASK, NAN); } @@ -325,7 +325,7 @@ chisq += PS_SQR(delta) * dy->data.F32[i]; - assert (!isnan(dy->data.F32[i])); - assert (!isnan(delta)); - assert (!isnan(chisq)); + if (isnan(dy->data.F32[i])) return NAN; + if (isnan(delta)) return NAN; + if (isnan(chisq)) return NAN; // we track alpha,beta and params,deriv separately Index: branches/tap_branches/psLib/src/math/psStats.c =================================================================== --- branches/tap_branches/psLib/src/math/psStats.c (revision 25900) +++ branches/tap_branches/psLib/src/math/psStats.c (revision 27838) @@ -749,7 +749,13 @@ // Iterate to get the best bin size; an iteration limit is enforced at the bottom of the loop. for (int iterate = 1; iterate > 0; iterate++) { - psTrace(TRACE, 6, - "-------------------- Iterating on Bin size. Iteration number %d --------------------\n", - iterate); + psTrace(TRACE, 6, "-------------------- Iterating on Bin size. Iteration number %d --------------------\n", iterate); + + if (iterate >= PS_ROBUST_MAX_ITERATIONS) { + // This occurs when a large number of the values are identical --- a bin size cannot be found + // that will spread out the distribution. Therefore, set what we can, and fall over + // gracefully. + COUNT_WARNING(10, 100, "Maximum number of iterations (%d) exceeded.", PS_ROBUST_MAX_ITERATIONS); + goto escape; + } // Get the minimum and maximum values @@ -791,12 +797,15 @@ psTrace(TRACE, 6, "Initial robust bin size is %.2f\n", binSize); - // ADD step 0: Construct the histogram with the specified bin size. NOTE: we can not specify the bin - // size precisely since the argument to psHistogramAlloc() is the number of bins, not the binSize. If - // we get here, we know that binSize != 0.0. - long numBins = (max - min) / binSize; // Number of bins + // ADD step 0: Construct the histogram with the specified bin size. NOTE: we can + // not specify the bin size precisely since the argument to psHistogramAlloc() is + // the number of bins, not the binSize. If we get here, we know that binSize != + // 0.0. We can also have a floating-point round-off error such that the last bin + // of the histogram does not correspond exactly with the value of 'max'. Let's be + // a bit generous and extend the histogram by two bins in either direction + long numBins = 4 + (max - min) / binSize; // Number of bins psTrace(TRACE, 6, "Numbins is %ld\n", numBins); psTrace(TRACE, 6, "Creating a robust histogram from data range (%.2f - %.2f)\n", min, max); // Generate the histogram - histogram = psHistogramAlloc(min, max, numBins); + histogram = psHistogramAlloc(min - 2.0*binSize, max + 2.0*binSize, numBins); // XXXXX we need to consider this step if errors -> variance if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal)) { @@ -807,5 +816,4 @@ psFree(statsMinMax); psFree(mask); - return false; } @@ -813,4 +821,43 @@ PS_VECTOR_PRINT_F32(histogram->bounds); PS_VECTOR_PRINT_F32(histogram->nums); + } + + // perversity check: if most of the values land in a single bin, then we probably + // have a perverse case (eg, small number of points at extremely large / small + // values; nearly bi-modal distribution). if so, keep only points within 5? 10? + // bins of that excess bin: + int nMaxBin = 0; + int iMaxBin = 0; + for (long i = 1; i < histogram->nums->n; i++) { + if (histogram->nums->data.F32[i] > nMaxBin) { + nMaxBin = histogram->nums->data.F32[i]; + iMaxBin = i; + } + } + if (nMaxBin > numValid / 2) { + float minKeep = histogram->bounds->data.F32[iMaxBin] - 10*binSize; + float maxKeep = histogram->bounds->data.F32[iMaxBin + 1] + 10*binSize; + int nInvalid = 0; + for (long i = 0; i < myVector->n; i++) { + // skip the already-masked values + if (mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & maskVal) continue; + bool invalid = false; + invalid |= (myVector->data.F32[i] <= minKeep); + invalid |= (myVector->data.F32[i] >= maxKeep); + invalid |= (!isfinite(myVector->data.F32[i])); + if (!invalid) continue; + mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = maskVal; + nInvalid ++; + } + + if (nInvalid) { + psTrace(TRACE, 6, "data is concentrated in a single bin, masking %d extreme outliers and retrying\n", nInvalid); + psFree(histogram); + psFree(cumulative); + histogram = NULL; + cumulative = NULL; + continue; + } + // if we did not mask anything, give up. } @@ -1007,4 +1054,5 @@ stats->robustN50 = N50; psTrace(TRACE, 6, "The robustN50 is %ld.\n", N50); + psTrace(TRACE, 6, "The robust median and stdev are %f, %f\n", stats->robustMedian, stats->robustStdev); // Clean up @@ -1825,4 +1873,5 @@ COUNT_WARNING(10, 100, "Failed to calculate the min/max of the input vector.\n"); psFree(statsMinMax); + psFree(histogram); goto escape; } @@ -1894,5 +1943,5 @@ if (!status) { - psErrorClear(); + psErrorClear(); COUNT_WARNING(10, 100, "Failed to fit a gaussian to the robust histogram.\n"); psFree(poly); @@ -1984,5 +2033,5 @@ if (!status) { - psErrorClear(); + psErrorClear(); COUNT_WARNING(10, 100, "Failed to fit a gaussian to the robust histogram.\n"); psFree(poly); Index: branches/tap_branches/psLib/src/math/psUnaryOp.c =================================================================== --- branches/tap_branches/psLib/src/math/psUnaryOp.c (revision 25900) +++ branches/tap_branches/psLib/src/math/psUnaryOp.c (revision 27838) @@ -211,5 +211,5 @@ } -psMathType* psUnaryOp(psPtr out, const psPtr in, const char *op) +psMathType* psUnaryOp(psPtr out, psPtr in, const char *op) { #define psUnaryOp_EXIT { \ Index: branches/tap_branches/psLib/src/math/psUnaryOp.h =================================================================== --- branches/tap_branches/psLib/src/math/psUnaryOp.h (revision 25900) +++ branches/tap_branches/psLib/src/math/psUnaryOp.h (revision 27838) @@ -59,5 +59,5 @@ psMathType* psUnaryOp( psPtr out, ///< Output type, either psImage or psVector. - const psPtr in, ///< Input, either psImage or psVector. + psPtr in, ///< Input, either psImage or psVector. const char *op ///< Operator. ); Index: branches/tap_branches/psLib/src/mathtypes/psImage.h =================================================================== --- branches/tap_branches/psLib/src/mathtypes/psImage.h (revision 25900) +++ branches/tap_branches/psLib/src/mathtypes/psImage.h (revision 27838) @@ -66,4 +66,5 @@ #define P_PSIMAGE_SET_ROW0(img,r0) {*(int*)&img->row0 = r0;} #define P_PSIMAGE_SET_TYPE(img,t) {*(psMathType*)&img->type = t;} +#define P_PSIMAGE_GET_TYPE(img) ((img)->type->type) /** Create an image of the specified size and type. Index: branches/tap_branches/psLib/src/mathtypes/psVector.c =================================================================== --- branches/tap_branches/psLib/src/mathtypes/psVector.c (revision 25900) +++ branches/tap_branches/psLib/src/mathtypes/psVector.c (revision 27838) @@ -729,5 +729,5 @@ char line[1024]; - sprintf (line, "vector: %s\n", name); + sprintf (line, "# vector: %s\n", name); if (write(fd, line, strlen(line))) {;} //ignore return value Index: branches/tap_branches/psLib/src/pslib_strict.h =================================================================== --- branches/tap_branches/psLib/src/pslib_strict.h (revision 25900) +++ branches/tap_branches/psLib/src/pslib_strict.h (revision 27838) @@ -102,5 +102,5 @@ #include "psType.h" #include "psArray.h" -#include "psBitSet.h" +#include "psBits.h" #include "psHash.h" #include "psList.h" Index: branches/tap_branches/psLib/src/sys/psErrorCodes.c.in =================================================================== --- branches/tap_branches/psLib/src/sys/psErrorCodes.c.in (revision 25900) +++ branches/tap_branches/psLib/src/sys/psErrorCodes.c.in (revision 27838) @@ -67,5 +67,5 @@ static void freeErrorDescription(psErrorDescription* err) { - psFree((void *)err->description); + psFree(err->description); } Index: branches/tap_branches/psLib/src/sys/psMemory.h =================================================================== --- branches/tap_branches/psLib/src/sys/psMemory.h (revision 25900) +++ branches/tap_branches/psLib/src/sys/psMemory.h (revision 27838) @@ -326,6 +326,7 @@ /** Free memory. This operates much like free(). - * + * * @see psAlloc, psRealloc + * note: we cast ptr to (void *) in case we are supplied a const pointer. */ #ifdef DOXYGEN @@ -336,5 +337,5 @@ #ifndef SWIG #define psFree(ptr) \ - psMemDecrRefCounter(ptr) + ptr = psMemDecrRefCounter((void *)ptr); #endif // ifndef SWIG #endif // ifdef DOXYGEN Index: branches/tap_branches/psLib/src/sys/psTrace.c =================================================================== --- branches/tap_branches/psLib/src/sys/psTrace.c (revision 25900) +++ branches/tap_branches/psLib/src/sys/psTrace.c (revision 27838) @@ -119,5 +119,5 @@ psMemSetPersistent((psPtr)comp->name,false); - psFree((void *)comp->name); + psFree(comp->name); } Index: branches/tap_branches/psLib/src/sys/psType.c =================================================================== --- branches/tap_branches/psLib/src/sys/psType.c (revision 25900) +++ branches/tap_branches/psLib/src/sys/psType.c (revision 27838) @@ -20,5 +20,5 @@ #include "psType.h" -#include "psBitSet.h" +#include "psBits.h" #include "psFits.h" #include "psPixels.h" @@ -45,6 +45,6 @@ } break; - case PS_DATA_BITSET: - if (psMemCheckBitSet(ptr)) { + case PS_DATA_BITS: + if (psMemCheckBits(ptr)) { return true; } Index: branches/tap_branches/psLib/src/sys/psType.h =================================================================== --- branches/tap_branches/psLib/src/sys/psType.h (revision 25900) +++ branches/tap_branches/psLib/src/sys/psType.h (revision 27838) @@ -107,5 +107,5 @@ PS_DATA_STRING = 0x10000, ///< psString (char *) PS_DATA_ARRAY, ///< psArray - PS_DATA_BITSET, ///< psBitSet + PS_DATA_BITS, ///< psBits PS_DATA_CUBE, ///< psCube PS_DATA_FITS, ///< psFits Index: branches/tap_branches/psLib/src/types/Makefile.am =================================================================== --- branches/tap_branches/psLib/src/types/Makefile.am (revision 25900) +++ branches/tap_branches/psLib/src/types/Makefile.am (revision 27838) @@ -6,5 +6,5 @@ libpslibtypes_la_SOURCES = \ psArray.c \ - psBitSet.c \ + psBits.c \ psHash.c \ psList.c \ @@ -23,5 +23,5 @@ pkginclude_HEADERS = \ psArray.h \ - psBitSet.h \ + psBits.h \ psHash.h \ psList.h \ Index: branches/tap_branches/psLib/src/types/psArray.c =================================================================== --- branches/tap_branches/psLib/src/types/psArray.c (revision 25900) +++ branches/tap_branches/psLib/src/types/psArray.c (revision 27838) @@ -166,5 +166,5 @@ // drop an item from the array and free it bool psArrayRemoveData(psArray* array, - const psPtr data) + psPtr data) { PS_ASSERT_ARRAY_NON_NULL(array, false); Index: branches/tap_branches/psLib/src/types/psArray.h =================================================================== --- branches/tap_branches/psLib/src/types/psArray.h (revision 25900) +++ branches/tap_branches/psLib/src/types/psArray.h (revision 27838) @@ -186,5 +186,5 @@ bool psArrayRemoveData( psArray* array, ///< array to operate on - const psPtr data ///< the data pointer to remove from psArray + psPtr data ///< the data pointer to remove from psArray ); Index: branches/tap_branches/psLib/src/types/psBitSet.c =================================================================== --- branches/tap_branches/psLib/src/types/psBitSet.c (revision 25900) +++ (revision ) @@ -1,303 +1,0 @@ -/** @file psBitSet.c - * - * @brief Creates an array of bytes of arbitrary length for storing individual bits. - * - * Bit masks are useful tools for toggling various flags and options. This set of functions module provides - * a mechanism to create an array of bits of arbitrary length and manipulate them with basic binary - * operations. A print function is also provided to display the entire set of bits in binary format as a - * string. - * - * @author Ross Harman, MHPCC - * @author Robert DeSonia, MHPCC - * - * @version $Revision: 1.42 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-08-09 03:30:16 $ - * - * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii - */ - -#ifdef HAVE_CONFIG_H -# include "config.h" -#endif - -#include -#include -#include -#include - -#include "psBitSet.h" -#include "psMemory.h" -#include "psError.h" -#include "psAbort.h" -#include "psString.h" -#include "psAssert.h" - - - -enum { - UNKNOWN_OP, - AND_OP, - OR_OP, - XOR_OP, - NOT_OP -}; - -static void bitSetFree(psBitSet* inBitSet); - -/** Private function to create a mask. - * - * Creates an eight bit mask with the given bit set. All other bits in the byte are zero. The input bit uses - * zero-based indexing, and is the cumulitive index within the array, not the localized byte's bit position. - * - * @return char*: Pointer to byte in which bit is contained. - */ -PS_ATTR_PURE static char mask(psS32 bit) -{ - char mask = (char)0x01; - - // Ignore splint warning about negative bit shifts - /* @i@ */ - mask = mask << (bit % 8); - - return mask; -} - -static void bitSetFree(psBitSet* inBitSet) -{ - psFree(inBitSet->bits); -} - -bool psMemCheckBitSet(psPtr ptr) -{ - PS_ASSERT_PTR(ptr, false); - return ( psMemGetDeallocator(ptr) == (psFreeFunc)bitSetFree ); -} - - -psBitSet* p_psBitSetAlloc(const char *file, - unsigned int lineno, - const char *func, - long nalloc) -{ - if (nalloc < 0) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, - _("The number of bit in a psBitSet (%ld) must be greater than zero."), - nalloc); - return NULL; - } - - psS32 numBytes = 0; - psBitSet* newObj = NULL; - - numBytes = ceil(nalloc / 8.0); - newObj = p_psAlloc(file, lineno, func, sizeof(psBitSet)); - psMemSetDeallocator(newObj, (psFreeFunc) bitSetFree); - newObj->n = numBytes; - - // Ignore splint warning about releasing pointer members, since they've not been allocated yet - /* @i@ */ - newObj->bits = psAlloc(sizeof(char) * numBytes); - - memset(newObj->bits, 0, numBytes); - - return newObj; -} - - -bool psBitSetSet(psBitSet* bitSet, - long bit) -{ - unsigned char *byte = NULL; - - if (bitSet == NULL) { - psError(PS_ERR_BAD_PARAMETER_NULL, true, - _("Can not operate on a NULL psBitSet.")); - return false; - } else if ( (bit < 0) || - (bit > bitSet->n * 8 - 1) ) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, - _("The specified bit position (%ld) is invalid. Position must be between 0 and %ld."), - bit, bitSet->n * 8 - 1); - return false; - } - // Variable byte is the byte in the array that contains the bit to be set - byte = bitSet->bits + bit / 8; - *byte |= mask(bit); - - return true; -} - -bool psBitSetClear(psBitSet* bitSet, - long bit) -{ - unsigned char *byte = NULL; - - if (bitSet == NULL) { - psError(PS_ERR_BAD_PARAMETER_NULL, true, - _("Can not operate on a NULL psBitSet.")); - return false; - } else if ( (bit < 0) || - (bit > bitSet->n * 8 - 1) ) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, - _("The specified bit position (%ld) is invalid. Position must be between 0 and %ld."), - bit, bitSet->n * 8 - 1); - return false; - } - // Variable byte is the byte in the array that contains the bit to be set - byte = bitSet->bits + bit / 8; - *byte &= ! mask(bit); - - return true; -} - -bool psBitSetTest(const psBitSet* bitSet, - long bit) -{ - unsigned char *byte = NULL; - - if (bitSet == NULL) { - psError(PS_ERR_BAD_PARAMETER_NULL, true, - _("Can not operate on a NULL psBitSet.")); - return false; - } else if ( (bit < 0) || - (bit > bitSet->n * 8 - 1) ) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, - _("The specified bit position (%ld) is invalid. Position must be between 0 and %ld."), - bit,bitSet->n * 8 - 1); - return false; - } - - // Variable byte is the byte in the array that contains the bit to be tested - byte = bitSet->bits + bit / 8; - return ((*byte & mask(bit)) != 0); -} - -psBitSet* psBitSetOp(psBitSet* outBitSet, - const psBitSet* inBitSet1, - const char *operator, - const psBitSet* inBitSet2) -{ - if (inBitSet1 == NULL) { - psError(PS_ERR_BAD_PARAMETER_NULL, true, - _("First psBitSet operand can not be NULL.")); - psFree(outBitSet); - return NULL; - } - if (operator == NULL) { - psError(PS_ERR_BAD_PARAMETER_NULL, true, - _("Specified operator is NULL. Must specify desired operator.")); - psFree(outBitSet); - return NULL; - } - psS32 i = 0; - psS32 n = 0; - unsigned char* outBits = NULL; - unsigned char* inBits1 = NULL; - unsigned char* inBits2 = NULL; - psS32 op = UNKNOWN_OP; - - inBits1 = inBitSet1->bits; - - - // parse the operator - if (strcmp(operator,"AND")==0) { - op = AND_OP; - } else if (strcmp(operator,"OR")==0) { - op = OR_OP; - } else if (strcmp(operator,"XOR")==0) { - op = XOR_OP; - } else if (strcmp(operator,"NOT")==0) { - op = NOT_OP; - } else { - psFree(outBitSet); - psError(PS_ERR_BAD_PARAMETER_VALUE, true, - _("Specified operator, %s, is invalid. Valid operators are AND, OR, and XOR."), - operator); - return NULL; - } - - if (op != NOT_OP) { - if (inBitSet2 == NULL) { - psError(PS_ERR_BAD_PARAMETER_NULL, true, - _("Second psBitSet operand can not be NULL.")); - psFree(outBitSet); - return NULL; - } - - if (inBitSet1->n != inBitSet2->n) { - psError(PS_ERR_BAD_PARAMETER_SIZE, true, - _("The psBitSet operand must be the same size.")); - psFree(outBitSet); - return NULL; - } - inBits2 = inBitSet2->bits; - } - - if (outBitSet == NULL) { - outBitSet = psBitSetAlloc(inBitSet1->n*8); - } else if (outBitSet->n != inBitSet1->n) { - outBitSet->n = inBitSet1->n; - outBitSet->bits = psRealloc(outBitSet->bits, inBitSet1->n); - } - - n = outBitSet->n; - outBits = outBitSet->bits; - - switch (op) { - case AND_OP: - for (i = 0; i < n; i++) { - outBits[i] = inBits1[i] & inBits2[i]; - } - break; - case OR_OP: - for (i = 0; i < n; i++) { - outBits[i] = inBits1[i] | inBits2[i]; - } - break; - case XOR_OP: - for (i = 0; i < n; i++) { - outBits[i] = inBits1[i] ^ inBits2[i]; - } - break; - case NOT_OP: - default: - for (i = 0; i < n; i++) { - outBits[i] = ~inBits1[i]; - } - break; - } - - return outBitSet; -} - -psBitSet* psBitSetNot(psBitSet* outBitSet, - const psBitSet* inBitSet) -{ - if (inBitSet == NULL) { - psError(PS_ERR_BAD_PARAMETER_NULL, true, - _("Operand can not be NULL.")); - psFree(outBitSet); - return NULL; - } - - outBitSet = psBitSetOp(outBitSet,inBitSet,"NOT",NULL); - - return outBitSet; -} - -psString psBitSetToString(const psBitSet* bitSet) -{ - PS_ASSERT_PTR_NON_NULL(bitSet, NULL); - psS32 i = 0; - psS32 numBits = bitSet->n * 8; - // char *outString = psAlloc((size_t) numBits + 1); - psString outString = psStringAlloc(numBits + 1); - - for (i = 0; i < numBits; i++) { - outString[numBits - i - 1] = psBitSetTest(bitSet, i) ? '1' : '0'; - } - - outString[numBits] = 0; - - return outString; -} Index: branches/tap_branches/psLib/src/types/psBitSet.h =================================================================== --- branches/tap_branches/psLib/src/types/psBitSet.h (revision 25900) +++ (revision ) @@ -1,175 +1,0 @@ -/** @file psBitSet.h - * - * @brief Creates an array of bytes of arbitrary length for storing individual bits. - * - * Bit masks are useful tools for toggling various flags and options. This set of functions module provides - * a mechanism to create an array of bits of arbitrary length and manipulate them with basic binary - * operations. A print function is also provided to display the entire set of bits in binary format as a - * string. - * - * @author PAP, EAM, IfA - * @author Ross Harman, MHPCC - * - * @version $Revision: 1.32 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-08-14 03:18:41 $ - * - * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii - */ - -#ifndef PSBITSET_H -#define PSBITSET_H - -#include "psType.h" -#include "psMutex.h" - -/// @addtogroup DataContainer Data Containers -/// @{ - -/******************************************************************************/ -/* TYPE DEFINITIONS */ -/******************************************************************************/ - -/** Struct containing array of bytes to hold bit data and corresponding array length. - * - * The bits in the struct are assembled in as an array of bytes with eight bits per - * byte. The bits are arranged with the LSB in first (right most) position of the - * first array element. - */ -typedef struct -{ - long n; ///< Number of bytes in the array - psU8 *bits; ///< Aray of bytes holding bits - psMutex lock; ///< Optional lock for thread safety -} -psBitSet; - -/*****************************************************************************/ -/* FUNCTION PROTOTYPES */ -/*****************************************************************************/ - -/** Checks the type of a particular pointer. - * - * Uses the appropriate deallocation function in psMemBlock to check the ptr - * datatype. - * - * @return bool: True if the pointer matches a psBitSet structure, false otherwise. - */ -bool psMemCheckBitSet( - psPtr ptr ///< the pointer whose type to check -); - - -/** Allocate a psBitSet. - * - * Create a psBitSet with the number of bits specified by the user. All bits - * are set to zero upon allocation. - * - * @return psBitSet* : Pointer to struct containing array of bits and size of array. - */ -#ifdef DOXYGEN -psBitSet* psBitSetAlloc( - long nalloc ///< Number of bits in psBitSet array -); -#else // ifdef DOXYGEN -psBitSet* p_psBitSetAlloc( - const char *file, ///< File of caller - unsigned int lineno, ///< Line number of caller - const char *func, ///< Function name of caller - long nalloc ///< Number of bits in psBitSet array -) PS_ATTR_MALLOC; -#define psBitSetAlloc(nalloc) \ - p_psBitSetAlloc(__FILE__, __LINE__, __func__, nalloc) -#endif // ifdef DOXYGEN - - - - -/** Set a bit. - * - * Sets a bit at a given bit location. The bit is set based on a zero index - * with the first bit set in the zero bit slot of the zero element of the byte - * array. As an example, setting bit 3 in an array with two elements would - * result in an psBitSet that looks like 00000000 00001000. - * - * @return bool : Successful operation? - */ -bool psBitSetSet( - psBitSet* bitSet, ///< Pointer to psBitSet to be set. - long bit ///< Bit to be set. -); - - -/** Clear a bit. - * - * Clear a bit at a given bit location. The bit is cleared based on a zero - * index with the first bit set in the zero bit slot of the zero element of - * the byte array. - * - * @return bool : Successful operation? - */ -bool psBitSetClear( - psBitSet* bitSet, ///< Pointer to psBitSet to be cleared. - long bit ///< Bit to be cleared. -); - - -/** Test the value of a bit. - * - * Prints the value of a bit at a given bit location, either one or zero. The - * resulting bit is based on a zero index format with the first bit set in the - * zero bit slot of the zero element of the byte array. As an example, - * testing bit 3 in a psBitSet with two bytes that looks like 00000000 - * 00001000 would return a value of one, since that is the value that was set. - * - * @return bool: True if successful, otherwise false - */ -bool psBitSetTest( - const psBitSet* bitSet, ///< Pointer psBitSet to be tested. - long bit ///< Bit to be tested. -); - - -/** Perform a binary operation on two psBitSets - * - * Perform an AND, OR, or XOR on two psBitSets. If the BitMasks are not the - * same size, the operation will not be performed and an error message will be - * logged. - * - * @return psBitSet* : Pointer to struct containing result of binary operation. - */ -psBitSet* psBitSetOp( - psBitSet* outBitSet, ///< Resulting psBitSet from binary operation - const psBitSet* inBitSet1, ///< First psBitSet on which to operate - const char *operator, ///< Bit operation - const psBitSet* inBitSet2 ///< Second psBitSet on which to operate -); - - -/** Perform a not operation on a psBitSet - * - * Toggles bits in a psBitset. All zero bits are set to one and all one bits - * are set to zero. - * - * @return psBitSet* : Pointer to struct containing result of operation. - */ -psBitSet* psBitSetNot( - psBitSet* outBitSet, ///< Resulting psBitSet from operation - const psBitSet* inBitSet ///< Input psBitSet -); - - -/** Convert the psBitSet to a string of ones and zeros. - * - * Converts the contents of a psBitSet to a string representation of its - * binary form of ones and zeros. The LSB is the right-most chracter. Each set - * of eight characters represents one byte. - * - * @return psString: Pointer to character array containing string data. - */ -psString psBitSetToString( - const psBitSet* bitSet ///< psBitSet to convert -); - - -/// @} -#endif // #ifndef PSBITSET_H Index: branches/tap_branches/psLib/src/types/psBits.c =================================================================== --- branches/tap_branches/psLib/src/types/psBits.c (revision 27838) +++ branches/tap_branches/psLib/src/types/psBits.c (revision 27838) @@ -0,0 +1,219 @@ +/** @file psBits.c + * + * @brief Creates an array of bytes of arbitrary length for storing individual bits. + * + * Bit masks are useful tools for toggling various flags and options. This set of functions module provides + * a mechanism to create an array of bits of arbitrary length and manipulate them with basic binary + * operations. A print function is also provided to display the entire set of bits in binary format as a + * string. + * + * @author Ross Harman, MHPCC + * @author Robert DeSonia, MHPCC + * + * @version $Revision: 1.42 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-08-09 03:30:16 $ + * + * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii + */ + +#ifdef HAVE_CONFIG_H +# include "config.h" +#endif + +#include +#include +#include +#include + +#include "psBits.h" +#include "psMemory.h" +#include "psError.h" +#include "psAbort.h" +#include "psString.h" +#include "psAssert.h" + + +static void bitsFree(psBits *inBits); + +/** Private function to create a mask. + * + * Creates an eight bit mask with the given bit set. All other bits in the byte are zero. The input bit uses + * zero-based indexing, and is the cumulitive index within the array, not the localized byte's bit position. + * + * @return char*: Pointer to byte in which bit is contained. + */ +PS_ATTR_PURE static inline char mask(long bit) +{ + psAssert(bit < 8, "Bad bit: %ld", bit); + return (char)0x01 << (bit % 8); +} + +// Return the byte with the bit of interest +static inline psU8 *bitsGetByte(const psBits *bits, long bit) +{ + return bits->bits + bit / 8; +} + + +static void bitsFree(psBits *inBits) +{ + psFree(inBits->bits); +} + +bool psMemCheckBits(psPtr ptr) +{ + PS_ASSERT_PTR(ptr, false); + return (psMemGetDeallocator(ptr) == (psFreeFunc)bitsFree); +} + + +psBits *p_psBitsAlloc(const char *file, unsigned int lineno, const char *func, long nalloc) +{ + psAssert(nalloc >= 0, "The number of bits in a psBits (%ld) must be greater than zero.", nalloc); + + int numBytes = ceil((float)nalloc / 8.0); // Number of bytes to use + psBits *bits = p_psAlloc(file, lineno, func, sizeof(psBits)); + psMemSetDeallocator(bits, (psFreeFunc)bitsFree); + bits->n = nalloc; + + bits->bits = p_psAlloc(file, lineno, func, numBytes); + memset(bits->bits, 0, numBytes); + + return bits; +} + + +bool psBitsSet(psBits *bits, long bit) +{ + PS_ASSERT_BITS_NON_NULL(bits, false); + PS_ASSERT_BITS_VALID_BIT(bits, bit, false); + + psU8 *byte = bitsGetByte(bits, bit); // Byte with the bit of interest + *byte |= mask(bit); + + return true; +} + +bool psBitsClear(psBits *bits, long bit) +{ + PS_ASSERT_BITS_NON_NULL(bits, false); + PS_ASSERT_BITS_VALID_BIT(bits, bit, false); + + psU8 *byte = bitsGetByte(bits, bit); // Byte with the bit of interest + *byte &= !mask(bit); + + return true; +} + + +bool psBitsTest(const psBits *bits, long bit) +{ + // XXX These errors probably cannot be caught + PS_ASSERT_BITS_NON_NULL(bits, false); + PS_ASSERT_BITS_VALID_BIT(bits, bit, false); + + psU8 *byte = bitsGetByte(bits, bit); // Byte with the bit of interest + return ((*byte & mask(bit)) != 0); +} + +psBits *psBitsOp(psBits *outBits, const psBits *bits1, const char *operator, const psBits *bits2) +{ + PS_ASSERT_BITS_NON_NULL(bits1, NULL); + PS_ASSERT_STRING_NON_EMPTY(operator, NULL); + + enum { + UNKNOWN_OP, + AND_OP, + OR_OP, + XOR_OP, + NOT_OP + } op = UNKNOWN_OP; + + // Parse the operator + if (strcmp(operator, "AND") == 0) { + op = AND_OP; + } else if (strcmp(operator, "OR") == 0) { + op = OR_OP; + } else if (strcmp(operator, "XOR") == 0) { + op = XOR_OP; + } else if (strcmp(operator, "NOT") == 0) { + op = NOT_OP; + } else { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + _("Specified operator, %s, is invalid. Valid operators are AND, OR, XOR, NOT."), + operator); + return NULL; + } + + if (op != NOT_OP) { + PS_ASSERT_BITS_NON_NULL(bits2, false); + if (bits1->n != bits2->n) { + psError(PS_ERR_BAD_PARAMETER_SIZE, true, + "The psBits operands must be the same size: %ld vs %ld", + bits1->n, bits2->n); + return NULL; + } + } + + int numBits = bits1->n; // Number of bits + int numBytes = ceil((float)numBits / 8.0); // Number of bytes + + if (!outBits) { + outBits = psBitsAlloc(numBits); + } else if (outBits->n != numBits) { + outBits->n = numBits; + outBits->bits = psRealloc(outBits->bits, numBytes); + } + + psU8 *in1 = bits1->bits; // Bits for input 1 + psU8 *in2 = (op == NOT_OP) ? NULL : bits2->bits; // Bits for input 2 + psU8 *out = outBits->bits; // Bits for output + + switch (op) { + case AND_OP: + for (int i = 0; i < numBytes; i++) { + out[i] = in1[i] & in2[i]; + } + break; + case OR_OP: + for (int i = 0; i < numBytes; i++) { + out[i] = in1[i] | in2[i]; + } + break; + case XOR_OP: + for (int i = 0; i < numBytes; i++) { + out[i] = in1[i] ^ in2[i]; + } + break; + case NOT_OP: + default: + for (int i = 0; i < numBytes; i++) { + out[i] = ~in1[i]; + } + break; + } + + return outBits; +} + +psBits *psBitsNot(psBits *outBits, const psBits *inBits) +{ + return psBitsOp(outBits, inBits, "NOT", NULL); +} + +psString psBitsToString(const psBits *bits) +{ + PS_ASSERT_BITS_NON_NULL(bits, NULL); + + psS32 numBits = bits->n; + psString string = psStringAlloc(numBits + 1); + + char *out = &string[numBits]; + *out = '\0'; + out--; + for (long i = 0; i < numBits; i++, out--) { + *out = psBitsTest(bits, i) ? '1' : '0'; + } + + return string; +} Index: branches/tap_branches/psLib/src/types/psBits.h =================================================================== --- branches/tap_branches/psLib/src/types/psBits.h (revision 27838) +++ branches/tap_branches/psLib/src/types/psBits.h (revision 27838) @@ -0,0 +1,191 @@ +/** @file psBits.h + * + * @brief Creates an array of bytes of arbitrary length for storing individual bits. + * + * Bit masks are useful tools for toggling various flags and options. This set of functions module provides + * a mechanism to create an array of bits of arbitrary length and manipulate them with basic binary + * operations. A print function is also provided to display the entire set of bits in binary format as a + * string. + * + * @author PAP, EAM, IfA + * @author Ross Harman, MHPCC + * + * @version $Revision: 1.32 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-08-14 03:18:41 $ + * + * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii + */ + +#ifndef PSBITS_H +#define PSBITS_H + +#include "psType.h" +#include "psMutex.h" + +/// @addtogroup DataContainer Data Containers +/// @{ + +/******************************************************************************/ +/* TYPE DEFINITIONS */ +/******************************************************************************/ + +/** Struct containing array of bytes to hold bit data and corresponding array length. + * + * The bits in the struct are assembled in as an array of bytes with eight bits per + * byte. The bits are arranged with the LSB in first (right most) position of the + * first array element. + */ +typedef struct +{ + long n; ///< Number of bits in the array + psU8 *bits; ///< Aray of bytes holding bits + psMutex lock; ///< Optional lock for thread safety +} psBits; + +/*****************************************************************************/ +/* FUNCTION PROTOTYPES */ +/*****************************************************************************/ + +/** Checks the type of a particular pointer. + * + * Uses the appropriate deallocation function in psMemBlock to check the ptr + * datatype. + * + * @return bool: True if the pointer matches a psBits structure, false otherwise. + */ +bool psMemCheckBits( + psPtr ptr ///< the pointer whose type to check +); + + +/** Allocate a psBits. + * + * Create a psBits with the number of bits specified by the user. All bits + * are set to zero upon allocation. + * + * @return psBits* : Pointer to struct containing array of bits and size of array. + */ +#ifdef DOXYGEN +psBits* psBitsAlloc( + long nalloc ///< Number of bits in psBits array +); +#else // ifdef DOXYGEN +psBits* p_psBitsAlloc( + const char *file, ///< File of caller + unsigned int lineno, ///< Line number of caller + const char *func, ///< Function name of caller + long nalloc ///< Number of bits in psBits array +) PS_ATTR_MALLOC; +#define psBitsAlloc(nalloc) \ + p_psBitsAlloc(__FILE__, __LINE__, __func__, nalloc) +#endif // ifdef DOXYGEN + + + + +/** Set a bit. + * + * Sets a bit at a given bit location. The bit is set based on a zero index + * with the first bit set in the zero bit slot of the zero element of the byte + * array. As an example, setting bit 3 in an array with two elements would + * result in an psBits that looks like 00000000 00001000. + * + * @return bool : Successful operation? + */ +bool psBitsSet( + psBits* bits, ///< Pointer to psBits to be set. + long bit ///< Bit to be set. +); + + +/** Clear a bit. + * + * Clear a bit at a given bit location. The bit is cleared based on a zero + * index with the first bit set in the zero bit slot of the zero element of + * the byte array. + * + * @return bool : Successful operation? + */ +bool psBitsClear( + psBits* bits, ///< Pointer to psBits to be cleared. + long bit ///< Bit to be cleared. +); + + +/** Test the value of a bit. + * + * Prints the value of a bit at a given bit location, either one or zero. The + * resulting bit is based on a zero index format with the first bit set in the + * zero bit slot of the zero element of the byte array. As an example, + * testing bit 3 in a psBits with two bytes that looks like 00000000 + * 00001000 would return a value of one, since that is the value that was set. + * + * @return bool: True if successful, otherwise false + */ +bool psBitsTest( + const psBits* bits, ///< Pointer psBits to be tested. + long bit ///< Bit to be tested. +); + + +/** Perform a binary operation on two psBitss + * + * Perform an AND, OR, or XOR on two psBitss. If the BitMasks are not the + * same size, the operation will not be performed and an error message will be + * logged. + * + * @return psBits* : Pointer to struct containing result of binary operation. + */ +psBits* psBitsOp( + psBits* outBits, ///< Resulting psBits from binary operation + const psBits* inBits1, ///< First psBits on which to operate + const char *operator, ///< Bit operation + const psBits* inBits2 ///< Second psBits on which to operate +); + + +/** Perform a not operation on a psBits + * + * Toggles bits in a psBits. All zero bits are set to one and all one bits + * are set to zero. + * + * @return psBits* : Pointer to struct containing result of operation. + */ +psBits* psBitsNot( + psBits* outBits, ///< Resulting psBits from operation + const psBits* inBits ///< Input psBits +); + + +/** Convert the psBits to a string of ones and zeros. + * + * Converts the contents of a psBits to a string representation of its + * binary form of ones and zeros. The LSB is the right-most chracter. Each set + * of eight characters represents one byte. + * + * @return psString: Pointer to character array containing string data. + */ +psString psBitsToString( + const psBits* bits ///< psBits to convert +); + +#define PS_ASSERT_BITS_NON_NULL(NAME, RVAL) \ +if ((NAME) == NULL || (NAME)->bits == NULL || (NAME)->n < 0) { \ + psError(PS_ERR_BAD_PARAMETER_NULL, true, \ + "Unallowable operation: psBits %s or its data is NULL.", \ + #NAME); \ + return RVAL; \ +} \ + +#define PS_ASSERT_BITS_VALID_BIT(NAME, BIT, RVAL) \ + if (BIT < 0 || BIT >= (NAME)->n) { \ + psError(PS_ERR_BAD_PARAMETER_VALUE, true, \ + "Unallowable operation: psBits %s or its data is NULL.", \ + #NAME); \ + return RVAL; \ +} \ + + + +/// @} +#endif // #ifndef PSBITS_H Index: branches/tap_branches/psLib/src/types/psList.c =================================================================== --- branches/tap_branches/psLib/src/types/psList.c (revision 25900) +++ branches/tap_branches/psLib/src/types/psList.c (revision 27838) @@ -210,4 +210,5 @@ } + // XXX remove this as an error if (location < 0 || location >= (int)list->n) { psError(PS_ERR_BAD_PARAMETER_VALUE, true, @@ -446,4 +447,5 @@ psListIterator* iterator = list->iterators->data[0]; + // XXX remove this as an eror if (! psListIteratorSet(iterator,location)) { psError(PS_ERR_BAD_PARAMETER_VALUE, true, Index: branches/tap_branches/psLib/src/types/psLookupTable.c =================================================================== --- branches/tap_branches/psLib/src/types/psLookupTable.c (revision 25900) +++ branches/tap_branches/psLib/src/types/psLookupTable.c (revision 27838) @@ -31,4 +31,6 @@ #include "psString.h" #include "psError.h" +#include "psString.h" +#include "psSlurp.h" #include "psLookupTable.h" @@ -153,5 +155,7 @@ char *end = NULL; \ ps##TYPE value = FUNC(strValue, &end, 0); \ - if (*end != '\0' && !isspace(*end)) { \ + if (*end != '\0' && *end != '\n' && !isspace(*end)) { \ + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Characters left over after parsing %s: %s", \ + strValue, end); \ *status = PS_PARSE_ERROR_VALUE; \ } \ @@ -164,5 +168,7 @@ char *end = NULL; \ ps##TYPE value = FUNC(strValue, &end); \ - if (*end != '\0' && !isspace(*end)) { \ + if (*end != '\0' && *end != '\n' && !isspace(*end)) { \ + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Characters left over after parsing %s: %s", \ + strValue, end); \ *status = PS_PARSE_ERROR_VALUE; \ } \ @@ -244,175 +250,132 @@ } -psArray *psVectorsReadFromFile(const char *filename, - const char *format) +psArray *psVectorsReadFromFile(const char *filename, const char *format) { PS_ASSERT_STRING_NON_EMPTY(filename, NULL); PS_ASSERT_STRING_NON_EMPTY(format, NULL); - psArray* outputArray = NULL; - psVector* colVector = NULL; - char* strValue = NULL; - char* strNum = NULL; - char* line = NULL; - char* linePtr = NULL; - int numCols = 0; - int numRows = 0; - FILE* fp = NULL; - const char* tempFormat = NULL; - psParseErrorType parseStatus = PS_PARSE_SUCCESS; - - // Create temp pointer which can then be used several times - tempFormat = format; - - // Create output array and set array elements to zero - outputArray = psArrayAllocEmpty(INITIAL_NUM); - - // Parse the format string to determine how many vectors - // and whether the format string is valid - while ((strValue = getToken((char**)&tempFormat, " \t", &parseStatus))) { - - // Check for %d format sub string + psArray *outputArray = psArrayAllocEmpty(INITIAL_NUM); // Array of vectors to return + psParseErrorType parseStatus = PS_PARSE_SUCCESS; // Status of parsing + + // Parse the format string to determine how many vectors and whether the format string is valid + const char *tempFormat = format; // Pointer into format + psString strValue; // Format of interest + int numCols = 0; // Number of columns found in format + while ((strValue = getToken((char**)&tempFormat, " \t", &parseStatus)) && + parseStatus == PS_PARSE_SUCCESS) { + if (strstr(strValue,"\%*") != 0) { + // Don't increase number of columns + continue; + } + psElemType type; // Type specified if (strcmp(strValue,"\%d") == 0 ) { - numCols++; - colVector = psVectorAlloc(1,PS_TYPE_S32); - outputArray = psArrayAdd(outputArray, ARRAY_STRIDE, colVector); - psFree(colVector); + type = PS_TYPE_S32; } else if (strcmp(strValue,"\%ld") == 0) { - numCols++; - colVector = psVectorAlloc(1,PS_TYPE_S64); - outputArray = psArrayAdd(outputArray, ARRAY_STRIDE, colVector); - psFree(colVector); + type = PS_TYPE_S64; } else if (strcmp(strValue,"\%f") == 0) { - numCols++; - colVector = psVectorAlloc(1,PS_TYPE_F32); - outputArray = psArrayAdd(outputArray, ARRAY_STRIDE, colVector); - psFree(colVector); + type = PS_TYPE_F32; } else if (strcmp(strValue,"\%lf") == 0) { - numCols++; - colVector = psVectorAlloc(1,PS_TYPE_F64); - outputArray = psArrayAdd(outputArray, ARRAY_STRIDE, colVector); - psFree(colVector); - } else if (strstr(strValue,"\%*") != 0) { - // Don't increase number of columns + type = PS_TYPE_F64; } else { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, - "Invalid format specifier"); + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Invalid format specifier: %s", strValue); psFree(strValue); - numCols = 0; - break; - } + psFree(outputArray); + return NULL; + } + psVector *colVector = psVectorAllocEmpty(1, type); // Vector for type + outputArray = psArrayAdd(outputArray, ARRAY_STRIDE, colVector); + psFree(colVector); + numCols++; psFree(strValue); + } + if (parseStatus != PS_PARSE_SUCCESS) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Failed to parse format at column %d: %s", + numCols, strValue); + psFree(strValue); + psFree(outputArray); + return NULL; + } + + if (numCols == 0) { + // Format string parse error detected + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Format string was not parsed sucessfully"); + psFree(outputArray); + return NULL; } // If the format string was parsed successfully and return numCols the // prepare to open file and read values - if (numCols > 0) { - - // Open specified file - if ((fp=fopen(filename, "r")) == NULL) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, _("Failed to open file %s."), - filename); - psFree(outputArray); - return NULL; - } else { - // Initialize array index - int arrayIndex = 0; - - // Create reusable line for continuous read - line = (char*)psAlloc(MAX_STRING_LENGTH*sizeof(char)); - - // Loop through file to get numRows, numCols, and column data types - while ((fgets(line, MAX_STRING_LENGTH, fp) != NULL) && - (parseStatus == PS_PARSE_SUCCESS)) { - - // Copy pointer to line for parsing - linePtr = line; - - // If line is not a comment or blank, then extract data - if (!ignoreLine(linePtr)) { - numRows++; - - // Copy format pointer for parsing - tempFormat = format; - arrayIndex = 0; - parseStatus = PS_PARSE_SUCCESS; - - // Loop through format and line strings to get values in text table file - while ((strValue=getToken((char**)&tempFormat," \t",&parseStatus)) - && (strNum=getToken((char**)&linePtr," \t",&parseStatus)) ) { - // Set column vector - colVector = outputArray->data[arrayIndex]; - - // Set column entries based on format string defining the type - if (strcmp(strValue,"\%d") == 0 ) { - colVector = psVectorRecycle(colVector, numRows, - colVector->type.type); - parseValue(colVector,numRows-1,strNum,&parseStatus); - arrayIndex++; - } else if (strcmp(strValue,"\%ld") == 0) { - colVector = psVectorRecycle(colVector, numRows, - colVector->type.type); - parseValue(colVector,numRows-1,strNum,&parseStatus); - arrayIndex++; - } else if (strcmp(strValue,"\%f") == 0) { - colVector = psVectorRecycle(colVector, numRows, - colVector->type.type); - parseValue(colVector,numRows-1,strNum,&parseStatus); - arrayIndex++; - } else if (strcmp(strValue,"\%lf") == 0) { - colVector = psVectorRecycle(colVector, numRows, - colVector->type.type); - parseValue(colVector,numRows-1,strNum,&parseStatus); - arrayIndex++; - } else if (strstr(strValue,"\%*") != 0) { - // Don't increase number of columns - } - psFree(strValue); - psFree(strNum); - - // If the file line was not parsed successful report - // error and return NULL - if (parseStatus != PS_PARSE_SUCCESS) { - psError(PS_ERR_UNKNOWN, true, - "Parsing text file failed."); - fclose(fp); - psFree(outputArray); - psFree(line); - return NULL; - } - } - if (strValue != NULL && strNum == NULL) { - psError(PS_ERR_UNKNOWN, true, - "Parsing text file failed - missing table value(s)."); - fclose(fp); - psFree(outputArray); - psFree(line); - psFree(strValue); - return NULL; - } - } // ignore line + + psString file = psSlurpFilename(filename); // Contents of file + if (!file) { + psError(psErrorCodeLast(), false, "Unable to read file of vectors"); + psFree(outputArray); + return NULL; + } + + psArray *lines = psStringSplitArray(file, "\n", false); // Lines of file + psFree(file); + long numRows = 0; // Number of rows + for (long i = 0; i < lines->n; i++) { + psString line = lines->data[i]; // Line of interest + if (ignoreLine(line)) { + continue; + } + numRows++; + + char *linePtr = line; // Pointer into line + + // Copy format pointer for parsing + const char *tempFormat = format; // Pointer into format + long arrayIndex = 0; // Index in array + parseStatus = PS_PARSE_SUCCESS; + + // Loop through format and line strings to get values in text table file + char *strNum; // Number within line + while ((strValue=getToken((char**)&tempFormat," \t",&parseStatus)) && + (strNum=getToken((char**)&linePtr," \t",&parseStatus)) && + parseStatus == PS_PARSE_SUCCESS) { + if (strstr(strValue,"\%*") != 0) { + continue; } - //Return NULL for an empty table - if (numRows == 0) { - psError(PS_ERR_UNKNOWN, true, - "Parsing text file failed - input table is empty."); - fclose(fp); + // Set column vector + psVector *colVector = outputArray->data[arrayIndex]; // Column vector of interest + + outputArray->data[arrayIndex] = colVector = psVectorRecycle(colVector, numRows, + colVector->type.type); + parseValue(colVector, numRows - 1, strNum, &parseStatus); + arrayIndex++; + + if (parseStatus != PS_PARSE_SUCCESS) { + psError(PS_ERR_UNKNOWN, false, "Parsing text file failed: %s as %s", strNum, strValue); psFree(outputArray); - psFree(line); + psFree(lines); + psFree(strNum); + psFree(strValue); return NULL; } - - // Read on the lines in the file - close file pointer - fclose(fp); - psFree(line); - } - } else { - // Format string parse error detected - psError(PS_ERR_UNKNOWN, true, - "Format string was not parsed sucessfully"); + psFree(strValue); + psFree(strNum); + + } + if (strValue != NULL && strNum == NULL) { + psError(PS_ERR_UNKNOWN, true, + "Parsing text file failed - missing table value(s)."); + psFree(outputArray); + psFree(lines); + psFree(strValue); + return NULL; + } + } + if (parseStatus != PS_PARSE_SUCCESS) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Failed to parse format at column %d: %s", + numCols, strValue); + psFree(strValue); psFree(outputArray); return NULL; } + + psFree(lines); // Return populated array Index: branches/tap_branches/psLib/src/types/psMetadata.c =================================================================== --- branches/tap_branches/psLib/src/types/psMetadata.c (revision 25900) +++ branches/tap_branches/psLib/src/types/psMetadata.c (revision 27838) @@ -280,5 +280,5 @@ } case PS_DATA_ARRAY: // psArray - case PS_DATA_BITSET: // psBitSet + case PS_DATA_BITS: // psBits case PS_DATA_CUBE: // psCube case PS_DATA_FITS: // psFits @@ -622,5 +622,6 @@ // may need to extend this to change the keyname in the copy -bool psMetadataItemSupplement(psMetadata *out, +bool psMetadataItemSupplement(bool *status, + psMetadata *out, const psMetadata *in, const char *key) @@ -632,5 +633,9 @@ psMetadataItem *item = psMetadataLookup(in, key); if (!item) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Could not find '%s' in metadata.\n", key); + if (status) { + *status = false; + } else { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Could not find '%s' in metadata.\n", key); + } return false; } @@ -968,5 +973,12 @@ } -psMetadataItem* psMetadataLookup(const psMetadata *md, +// Get entry by index +psMetadataItem *psMetadataGetIndex(psMetadata *md, long location) +{ + PS_ASSERT_METADATA_NON_NULL(md, false); + return psListGet(md->list, location); +} + +psMetadataItem *psMetadataLookup(const psMetadata *md, const char *key) { @@ -1147,4 +1159,5 @@ PS_ASSERT_METADATA_NON_NULL(md,NULL); + // XXX remove this as an error entry = (psMetadataItem*) psListGet(md->list, location); if (entry == NULL) { Index: branches/tap_branches/psLib/src/types/psMetadata.h =================================================================== --- branches/tap_branches/psLib/src/types/psMetadata.h (revision 25900) +++ branches/tap_branches/psLib/src/types/psMetadata.h (revision 27838) @@ -545,4 +545,5 @@ */ bool psMetadataItemSupplement( + bool *status, ///< if supplied, returns true/false if key is found (suppresses the error) psMetadata *out, ///< output Metadata container for copying. const psMetadata *in, ///< Metadata collection from which to copy. Index: branches/tap_branches/psLib/src/types/psMetadataConfig.c =================================================================== --- branches/tap_branches/psLib/src/types/psMetadataConfig.c (revision 25900) +++ branches/tap_branches/psLib/src/types/psMetadataConfig.c (revision 27838) @@ -1649,5 +1649,10 @@ return false; } - fprintf(file, "%s", fileString); + if (fprintf(file, "%s", fileString) != strlen(fileString)) { + psError(PS_ERR_IO, true, "Failed to write contents of configuration file %s", filename); + psFree(fileString); + fclose(file); + return false; + } psFree(fileString); if (fclose(file) == EOF) { Index: branches/tap_branches/psLib/test/types/Makefile.am =================================================================== --- branches/tap_branches/psLib/test/types/Makefile.am (revision 25900) +++ branches/tap_branches/psLib/test/types/Makefile.am (revision 27838) @@ -29,5 +29,5 @@ tap_psPixels_all \ tap_psHash_all \ - tap_psBitSet_all \ + tap_psBits_all \ tap_psList_all \ tap_psLookupTable_all \ Index: branches/tap_branches/psLib/test/types/tap_psBitSet_all.c =================================================================== --- branches/tap_branches/psLib/test/types/tap_psBitSet_all.c (revision 25900) +++ (revision ) @@ -1,279 +1,0 @@ -/** - * C Implementation: tap_psBitSet_all - * - * Description: Tests for psBitSetAlloc, psBitSetSet, psMemCheckBitSet, psBitSetClear, - * psBitSetTest, psBitSetOp, psBitSetNot, psBitSetToString - * - * Author: dRob , (C) 2006 - * - * Copyright: See COPYING file that comes with this distribution - * - */ -#include -#include -#include "tap.h" -#include "pstap.h" - - -int main(void) -{ - psLogSetFormat("HLNM"); - psLogSetLevel(PS_LOG_INFO); - plan_tests(31); - - - // testBitSetBasics() - { - psMemId id = psMemGetId(); - psBitSet *noBits = NULL; - psBitSet *bs = NULL; - - //Return NULL for attempt to Allocate BitSet of negative size. - { - noBits = psBitSetAlloc(-1); - ok( noBits == NULL, - "psBitSetAlloc: return NULL for negative BitSet-size input."); - } - //Return properly allocated 0-size psBitSet - { - bs = psBitSetAlloc(0); - ok( bs != NULL && psMemCheckBitSet(bs) && bs->n == 0, - "psBitSetAlloc: return properly allocated psBitSet."); - } - //Return properly allocated psBitSet - { - psFree(bs); - bs = psBitSetAlloc(8); - ok( bs != NULL && psMemCheckBitSet(bs) && bs->n == 1, - "psBitSetAlloc: return properly allocated psBitSet."); - } - //Make sure psMemCheckBitSet works correctly - return false - if (0) { - int j = 2; - ok( !psMemCheckBitSet(&j), - "psMemCheckBitSet: return false for non-BitSet input."); - } - - //BitSetSet Tests - //Return FALSE for NULL input psBitSet - { - bool rc = psBitSetSet(NULL, 0); - ok(rc == false, - "psBitSetSet: return FALSE for NULL BitSet input."); - } - //Return FALSE for negative bit input - { - bool rc = psBitSetSet(bs, -1); - ok( rc == false, - "psBitSetSet: return TRUE for negative bits input."); - noBits = NULL; - } - //Return FALSE for out-of-range bits - { - psFree(bs); - bs = psBitSetAlloc(8); - bool rc = psBitSetSet(bs, 8); - ok( rc == false, - "psBitSetSet: return TRUE for out-of-range bits input."); - noBits = NULL; - } - - //Return set BitSet for valid inputs - { - psBitSetSet(bs, 2); - ok( bs->bits[0] == 4, - "psBitSetSet: return properly set BitSet for valid inputs."); - } - - //BitSetClear Tests - //Return FALSE for NULL input psBitSet - { - bool rc = psBitSetClear(noBits, 0); - ok( rc == false, - "psBitSetClear: return FALSE for NULL BitSet input."); - } - //Return FALSE for negative bit input - { - bool rc = psBitSetClear(bs, -1); - ok( rc == false, - "psBitSetClear: return TRUE for negative bits input."); - noBits = NULL; - } - //Return FALSE for out-of-range bits - { - bool rc = psBitSetClear(bs, 8); - ok( rc == false, - "psBitSetClear: return FALSE for out-of-range bits input."); - noBits = NULL; - } - - //Return cleared BitSet for valid inputs - { - psBitSetClear(bs, 2); - ok( bs->bits[0] == 0, - "psBitSetClear: return properly cleared BitSet for valid inputs."); - } - - //BitSetTest Tests - //Return false for NULL input psBitSet - { - ok( !psBitSetTest(noBits, 0), - "psBitSetTest: return false for NULL BitSet input."); - } - //Return false for negative bit input - { - ok( !psBitSetTest(bs, -1), - "psBitSetTest: return false for negative bits input."); - } - //Return false for out-of-range bits - { - ok( !psBitSetTest(bs, 8), - "psBitSetTest: return false for out-of-range bits input."); - } - //Return false for non-matching bit in BitSet - { - ok( !psBitSetTest(bs, 2), - "psBitSetTest: return false for non-matching bit in BitSet."); - } - //Return false for non-matching bit in BitSet - { - psBitSetSet(bs, 2); - ok( psBitSetTest(bs, 2), - "psBitSetTest: return true for matching bit in BitSet."); - } - - psFree(bs); - ok(!psMemCheckLeaks(id, NULL, NULL, false), "no memory leaks"); - } - - - // testBitSetOps() - { - psMemId id = psMemGetId(); - psBitSet *noBits = NULL; - psBitSet *bs = NULL; - bs = psBitSetAlloc(8); - psBitSetSet(bs, 2); // 0000 0100 == 4 - psBitSetSet(bs, 3); // 0000 1100 == 12 - psBitSetSet(bs, 4); // 0001 1100 == 28 - psBitSetSet(bs, 5); // 0011 1100 == 60 - psBitSetSet(bs, 6); // 0111 1100 == 124 - psBitSetSet(bs, 7); // 1111 1100 == 252 - psBitSet *out = NULL; - - //psBitSetNot Tests - //Return NULL for NULL BitSet input - { - out = psBitSetNot(out, noBits); - ok( out == NULL, - "psBitSetNot: return NULL for NULL BitSet input."); - } - //Return correct BitSet for valid BitSet input - { - out = psBitSetNot(out, bs); //bs = 1111 1100 so out should = 0000 0011 = 3 - ok( out->bits[0] == 3, - "psBitSetNot: return correct BitSet for valid BitSet input."); - } - - //psBitSetOp Tests out = psBitSetOp(out, bs1, "op", bs2); - //Return NULL for NULL BitSet input - { - psFree(out); - out = NULL; - out = psBitSetOp(out, noBits, "op", noBits); - ok( out == NULL, - "psBitSetOp: return NULL for NULL BitSet input."); - } - //Return NULL for NULL operator input - { - out = psBitSetOp(out, bs, NULL, noBits); - ok( out == NULL, - "psBitSetOp: return NULL for NULL operator input."); - } - //Return NULL for invalid operator input - { - out = psBitSetOp(out, bs, "XAND", noBits); - ok( out == NULL, - "psBitSetOp: return NULL for invalid operator input."); - } - //Return NULL for AND operator with NULL second BitSet input - { - out = psBitSetOp(out, bs, "AND", noBits); - ok( out == NULL, - "psBitSetOp: return NULL for AND operator with NULL second BitSet input."); - } - //Return NULL for AND operator with BitSet inputs of differing size. - psBitSet *bs2 = psBitSetAlloc(16); - { - out = psBitSetOp(out, bs, "AND", bs2); - ok( out == NULL, - "psBitSetOp: return NULL for AND operator with BitSet inputs of" - " differing size."); - } - psFree(bs); - bs = psBitSetAlloc(16); - psBitSetSet(bs, 1); // 0000 0010 == 2 - psBitSetSet(bs2, 2); // 0000 0100 == 4 - //Return correct psBitSet output for valid inputs with AND operator - { - out = psBitSetOp(out, bs, "AND", bs2); - ok( out->bits[0] == 0, - "psBitSetOp: return correct psBitSet output for valid inputs" - " with AND operator."); - } - //Return correct psBitSet output for valid inputs with OR operator - { - out = psBitSetOp(out, bs, "OR", bs2); - ok( out->bits[0] == 6, - "psBitSetOp: return correct psBitSet output for valid inputs" - " with OR operator."); - } - //Return correct psBitSet output for valid inputs with XOR operator - psBitSetSet(bs2, 1); // 0000 0110 == 6 - { - out = psBitSetOp(out, bs, "XOR", bs2); - ok( out->bits[0] == 4, - "psBitSetOp: return correct psBitSet output for valid inputs" - " with XOR operator."); - } - //Return correct psBitSet output for valid inputs with NOT operator - { - psFree(out); - out = psBitSetAlloc(0); - psBitSetSet(bs, 2); // 0000 0110 == 4 - psBitSetSet(bs, 3); // 0000 1110 == 12 - psBitSetSet(bs, 4); // 0001 1110 == 28 - psBitSetSet(bs, 5); // 0011 1110 == 60 - psBitSetSet(bs, 6); // 0111 1110 == 124 - psBitSetSet(bs, 7); // 1111 1110 == 252 - out = psBitSetOp(out, bs, "NOT", bs2); - ok( out->bits[0] == 1, - "psBitSetOp: return correct psBitSet output for valid inputs" - " with NOT operator."); - } - - //psBitSetToString Tests - //Return NULL for NULL BitSet input - psString bitStr = NULL; - { - bitStr = psBitSetToString(noBits); - ok( bitStr == NULL, - "psBitSetToString: return NULL for NULL BitSet input."); - } - //Return correct string for valid BitSet input - { - psFree(bs); - bs = psBitSetAlloc(8); - psBitSetSet(bs, 2); // 0000 0100 == 4 - bitStr = psBitSetToString(bs); - ok( !strncmp(bitStr, "00000100", 10), - "psBitSetToString: return correct string for valid BitSet input."); - } - - psFree(bitStr); - psFree(out); - psFree(bs); - psFree(bs2); - ok(!psMemCheckLeaks(id, NULL, NULL, false), "no memory leaks"); - } -} Index: branches/tap_branches/psLib/test/types/tap_psBits_all.c =================================================================== --- branches/tap_branches/psLib/test/types/tap_psBits_all.c (revision 27838) +++ branches/tap_branches/psLib/test/types/tap_psBits_all.c (revision 27838) @@ -0,0 +1,273 @@ +/** + * C Implementation: tap_psBits_all + * + * Description: Tests for psBitsAlloc, psBitsSet, psMemCheckBits, psBitsClear, + * psBitsTest, psBitsOp, psBitsNot, psBitsToString + * + * Author: dRob , (C) 2006 + * + * Copyright: See COPYING file that comes with this distribution + * + */ +#include +#include +#include "tap.h" +#include "pstap.h" + + +int main(void) +{ + psLogSetFormat("HLNM"); + psLogSetLevel(PS_LOG_INFO); + plan_tests(30); + + + // testBitsBasics() + { + psMemId id = psMemGetId(); + psBits *noBits = NULL; + psBits *bs = NULL; + + //Return properly allocated 0-size psBits + { + bs = psBitsAlloc(0); + ok( bs != NULL && psMemCheckBits(bs) && bs->n == 0, + "psBitsAlloc: return properly allocated psBits."); + } + //Return properly allocated psBits + { + psFree(bs); + bs = psBitsAlloc(8); + ok( bs != NULL && psMemCheckBits(bs) && bs->n == 8, + "psBitsAlloc: return properly allocated psBits."); + } + //Make sure psMemCheckBits works correctly - return false + if (0) { + int j = 2; + ok( !psMemCheckBits(&j), + "psMemCheckBits: return false for non-Bits input."); + } + + //BitsSet Tests + //Return FALSE for NULL input psBits + { + bool rc = psBitsSet(NULL, 0); + ok(rc == false, + "psBitsSet: return FALSE for NULL Bits input."); + } + //Return FALSE for negative bit input + { + bool rc = psBitsSet(bs, -1); + ok( rc == false, + "psBitsSet: return TRUE for negative bits input."); + noBits = NULL; + } + //Return FALSE for out-of-range bits + { + psFree(bs); + bs = psBitsAlloc(8); + bool rc = psBitsSet(bs, 8); + ok( rc == false, + "psBitsSet: return TRUE for out-of-range bits input."); + noBits = NULL; + } + + //Return set Bits for valid inputs + { + psBitsSet(bs, 2); + ok( bs->bits[0] == 4, + "psBitsSet: return properly set Bits for valid inputs."); + } + + //BitsClear Tests + //Return FALSE for NULL input psBits + { + bool rc = psBitsClear(noBits, 0); + ok( rc == false, + "psBitsClear: return FALSE for NULL Bits input."); + } + //Return FALSE for negative bit input + { + bool rc = psBitsClear(bs, -1); + ok( rc == false, + "psBitsClear: return TRUE for negative bits input."); + noBits = NULL; + } + //Return FALSE for out-of-range bits + { + bool rc = psBitsClear(bs, 8); + ok( rc == false, + "psBitsClear: return FALSE for out-of-range bits input."); + noBits = NULL; + } + + //Return cleared Bits for valid inputs + { + psBitsClear(bs, 2); + ok( bs->bits[0] == 0, + "psBitsClear: return properly cleared Bits for valid inputs."); + } + + //BitsTest Tests + //Return false for NULL input psBits + { + ok( !psBitsTest(noBits, 0), + "psBitsTest: return false for NULL Bits input."); + } + //Return false for negative bit input + { + ok( !psBitsTest(bs, -1), + "psBitsTest: return false for negative bits input."); + } + //Return false for out-of-range bits + { + ok( !psBitsTest(bs, 8), + "psBitsTest: return false for out-of-range bits input."); + } + //Return false for non-matching bit in Bits + { + ok( !psBitsTest(bs, 2), + "psBitsTest: return false for non-matching bit in Bits."); + } + //Return false for non-matching bit in Bits + { + psBitsSet(bs, 2); + ok( psBitsTest(bs, 2), + "psBitsTest: return true for matching bit in Bits."); + } + + psFree(bs); + ok(!psMemCheckLeaks(id, NULL, NULL, false), "no memory leaks"); + } + + + // testBitsOps() + { + psMemId id = psMemGetId(); + psBits *noBits = NULL; + psBits *bs = NULL; + bs = psBitsAlloc(8); + psBitsSet(bs, 2); // 0000 0100 == 4 + psBitsSet(bs, 3); // 0000 1100 == 12 + psBitsSet(bs, 4); // 0001 1100 == 28 + psBitsSet(bs, 5); // 0011 1100 == 60 + psBitsSet(bs, 6); // 0111 1100 == 124 + psBitsSet(bs, 7); // 1111 1100 == 252 + psBits *out = NULL; + + //psBitsNot Tests + //Return NULL for NULL Bits input + { + out = psBitsNot(out, noBits); + ok( out == NULL, + "psBitsNot: return NULL for NULL Bits input."); + } + //Return correct Bits for valid Bits input + { + out = psBitsNot(out, bs); //bs = 1111 1100 so out should = 0000 0011 = 3 + ok( out->bits[0] == 3, + "psBitsNot: return correct Bits for valid Bits input."); + } + + //psBitsOp Tests out = psBitsOp(out, bs1, "op", bs2); + //Return NULL for NULL Bits input + { + psFree(out); + out = NULL; + out = psBitsOp(out, noBits, "op", noBits); + ok( out == NULL, + "psBitsOp: return NULL for NULL Bits input."); + } + //Return NULL for NULL operator input + { + out = psBitsOp(out, bs, NULL, noBits); + ok( out == NULL, + "psBitsOp: return NULL for NULL operator input."); + } + //Return NULL for invalid operator input + { + out = psBitsOp(out, bs, "XAND", noBits); + ok( out == NULL, + "psBitsOp: return NULL for invalid operator input."); + } + //Return NULL for AND operator with NULL second Bits input + { + out = psBitsOp(out, bs, "AND", noBits); + ok( out == NULL, + "psBitsOp: return NULL for AND operator with NULL second Bits input."); + } + //Return NULL for AND operator with Bits inputs of differing size. + psBits *bs2 = psBitsAlloc(16); + { + out = psBitsOp(out, bs, "AND", bs2); + ok( out == NULL, + "psBitsOp: return NULL for AND operator with Bits inputs of" + " differing size."); + } + psFree(bs); + bs = psBitsAlloc(16); + psBitsSet(bs, 1); // 0000 0010 == 2 + psBitsSet(bs2, 2); // 0000 0100 == 4 + //Return correct psBits output for valid inputs with AND operator + { + out = psBitsOp(out, bs, "AND", bs2); + ok( out->bits[0] == 0, + "psBitsOp: return correct psBits output for valid inputs" + " with AND operator."); + } + //Return correct psBits output for valid inputs with OR operator + { + out = psBitsOp(out, bs, "OR", bs2); + ok( out->bits[0] == 6, + "psBitsOp: return correct psBits output for valid inputs" + " with OR operator."); + } + //Return correct psBits output for valid inputs with XOR operator + psBitsSet(bs2, 1); // 0000 0110 == 6 + { + out = psBitsOp(out, bs, "XOR", bs2); + ok( out->bits[0] == 4, + "psBitsOp: return correct psBits output for valid inputs" + " with XOR operator."); + } + //Return correct psBits output for valid inputs with NOT operator + { + psFree(out); + out = psBitsAlloc(0); + psBitsSet(bs, 2); // 0000 0110 == 4 + psBitsSet(bs, 3); // 0000 1110 == 12 + psBitsSet(bs, 4); // 0001 1110 == 28 + psBitsSet(bs, 5); // 0011 1110 == 60 + psBitsSet(bs, 6); // 0111 1110 == 124 + psBitsSet(bs, 7); // 1111 1110 == 252 + out = psBitsOp(out, bs, "NOT", bs2); + ok( out->bits[0] == 1, + "psBitsOp: return correct psBits output for valid inputs" + " with NOT operator."); + } + + //psBitsToString Tests + //Return NULL for NULL Bits input + psString bitStr = NULL; + { + bitStr = psBitsToString(noBits); + ok( bitStr == NULL, + "psBitsToString: return NULL for NULL Bits input."); + } + //Return correct string for valid Bits input + { + psFree(bs); + bs = psBitsAlloc(8); + psBitsSet(bs, 2); // 0000 0100 == 4 + bitStr = psBitsToString(bs); + ok( !strncmp(bitStr, "00000100", 10), + "psBitsToString: return correct string for valid Bits input (%s).", bitStr); + } + + psFree(bitStr); + psFree(out); + psFree(bs); + psFree(bs2); + ok(!psMemCheckLeaks(id, NULL, NULL, false), "no memory leaks"); + } +}