/* $Id: cpl_fit-test.c,v 1.15 2008/02/07 10:54:01 llundin Exp $ * * This file is part of the ESO Common Pipeline Library * Copyright (C) 2001-2004 European Southern Observatory * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation; either version 2 of the License, or * (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program; if not, write to the Free Software * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ /* * $Author: llundin $ * $Date: 2008/02/07 10:54:01 $ * $Revision: 1.15 $ * $Name: $ */ #ifdef HAVE_CONFIG_H # include #endif /*----------------------------------------------------------------------------- Includes -----------------------------------------------------------------------------*/ #include #include #include "cpl_fit.h" #include "cpl_memory.h" #include "cpl_stats.h" #include "cpl_tools.h" /*----------------------------------------------------------------------------- Defines -----------------------------------------------------------------------------*/ #ifndef IMAGESZ #define IMAGESZ 10 #endif #ifndef IMAGESZFIT #define IMAGESZFIT 256 #endif /*----------------------------------------------------------------------------*/ /** * @defgroup cpl_fit_test Testing of the CPL utilities */ /*----------------------------------------------------------------------------*/ /*----------------------------------------------------------------------------- Defines -----------------------------------------------------------------------------*/ #define cpl_fit_imagelist_is_zero(A, B) \ cpl_fit_imagelist_is_zero_macro(A, B, __LINE__) #define cpl_fit_image_is_zero(A, B) \ cpl_fit_image_is_zero_macro(A, B, __LINE__) /*----------------------------------------------------------------------------- Private Function prototypes -----------------------------------------------------------------------------*/ static void cpl_fit_imagelist_polynomial_tests(void); static void cpl_fit_imagelist_is_zero_macro(const cpl_imagelist *, double, int); static void cpl_fit_image_is_zero_macro(const cpl_image *, double, int); /*----------------------------------------------------------------------------*/ /** @brief Unit tests of fit module **/ /*----------------------------------------------------------------------------*/ int main(void) { cpl_test_init(PACKAGE_BUGREPORT, CPL_MSG_WARNING); /* Insert tests below */ cpl_fit_imagelist_polynomial_tests(); /* End of tests */ return cpl_test_end(0); } static void cpl_fit_imagelist_polynomial_tests(void) { const double ditval[] = {0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0}; cpl_imagelist * fit; cpl_imagelist * input; cpl_image * dfiterror = cpl_image_new(IMAGESZFIT, IMAGESZFIT, CPL_TYPE_DOUBLE); cpl_image * ffiterror = cpl_image_new(IMAGESZFIT, IMAGESZFIT, CPL_TYPE_FLOAT); cpl_image * ifiterror = cpl_image_new(IMAGESZFIT, IMAGESZFIT, CPL_TYPE_INT); const int ndits = (int)(sizeof(ditval)/sizeof(double)); cpl_vector * vdit = cpl_vector_wrap(ndits, (double*)ditval); cpl_polynomial * vfiller = cpl_polynomial_new(1); const double sqsum = 204.0; /* Sum of squares of ditvals */ const double mytol = 5.52 * FLT_EPSILON; int i; const cpl_type pixel_type[] = {CPL_TYPE_DOUBLE, CPL_TYPE_FLOAT, CPL_TYPE_INT}; const cpl_type pixeltype = CPL_TYPE_FLOAT; int ntest; cpl_flops flopsum = 0; cpl_flops flops; double secs; double demax; cpl_msg_info(cpl_func, "Testing with %d %d X %d images", ndits, IMAGESZFIT, IMAGESZFIT); fit = cpl_fit_imagelist_polynomial(NULL, NULL, 0, 0, CPL_FALSE, CPL_TYPE_FLOAT, NULL); cpl_test_error(CPL_ERROR_NULL_INPUT); cpl_test_null( fit ); cpl_imagelist_delete(fit); input = cpl_imagelist_new(); fit = cpl_fit_imagelist_polynomial(vdit, input, 0, 0, CPL_FALSE, CPL_TYPE_FLOAT, NULL); cpl_test_error(CPL_ERROR_INCOMPATIBLE_INPUT); cpl_test_null( fit ); cpl_imagelist_delete(fit); fit = cpl_fit_imagelist_polynomial(vdit, input, 1, 1, CPL_FALSE, CPL_TYPE_FLOAT, NULL); cpl_test_error(CPL_ERROR_INCOMPATIBLE_INPUT); cpl_test_null( fit ); cpl_imagelist_delete(fit); cpl_tools_get_cputime(CPL_CLOCK_STOP); secs = cpl_tools_get_cputime(CPL_CLOCK_TIME); /* Test with all types of pixels */ for (ntest = 0; ntest < 3; ntest++) { const cpl_type test_type = pixel_type[ntest]; cpl_image * image = cpl_image_new(IMAGESZFIT, IMAGESZ, test_type); cpl_msg_info(cpl_func, "Fitting with pixel type %u", (unsigned)test_type); cpl_test_zero(cpl_imagelist_set(input, image, 0)); image = cpl_image_duplicate(image); cpl_test_zero(cpl_image_fill_noise_uniform(image, 1.0, 20.0)); cpl_test_zero(cpl_image_multiply_scalar(image, ditval[1])); cpl_test_zero(cpl_imagelist_set(input, image, 1)); /* A perfectly linear set */ for (i=2; i < ndits; i++) { image = cpl_image_multiply_scalar_create(cpl_imagelist_get(input, 1), ditval[i]); cpl_test_zero(cpl_imagelist_set(input, image, i)); } flops = cpl_tools_get_flops(); cpl_tools_get_cputime(CPL_CLOCK_START); fit = cpl_fit_imagelist_polynomial(vdit, input, 1, ndits-1, CPL_FALSE, test_type, NULL); cpl_tools_get_cputime(CPL_CLOCK_STOP); flopsum += cpl_tools_get_flops() - flops; cpl_test_error(CPL_ERROR_NONE ); cpl_test( cpl_imagelist_get_size(fit) == ndits - 1 ); /* The linarity must be equal to the values in image 1 - normalize */ cpl_test_zero(cpl_image_divide(cpl_imagelist_get(fit, 0), cpl_imagelist_get(input, 1))); /* Subtract the expected value in the 1st image */ cpl_test_zero(cpl_image_subtract_scalar(cpl_imagelist_get(fit, 0), 1.0)); cpl_fit_imagelist_is_zero(fit, IMAGESZFIT * mytol); cpl_imagelist_delete(fit); cpl_imagelist_delete(input); input = cpl_imagelist_new(); } /* Create a list of images with a 2nd order function */ for (i=0; i < ndits; i++) { cpl_image * image = cpl_image_new(IMAGESZFIT, IMAGESZFIT, CPL_TYPE_DOUBLE); cpl_test_zero(cpl_image_add_scalar(image, ditval[i]*ditval[i])); cpl_test_zero(cpl_imagelist_set(input, image, i)); cpl_msg_debug(cpl_func, "Dit and mean of input image no. %d: %g %g", i, ditval[i], cpl_image_get_mean(image)); } fit = cpl_fit_imagelist_polynomial(vdit, input, 1, ndits, CPL_FALSE, pixeltype, NULL); if (cpl_error_get_code() != CPL_ERROR_NONE ) { /* Fails on 32-bit intel, but not on others */ cpl_test_error(CPL_ERROR_SINGULAR_MATRIX ); cpl_test_null( fit ); } cpl_imagelist_delete(fit); /* Illegal max-degree */ fit = cpl_fit_imagelist_polynomial(vdit, input, 1, 0, CPL_FALSE, pixeltype, NULL); cpl_test_error(CPL_ERROR_ILLEGAL_INPUT); cpl_test_null( fit ); cpl_imagelist_delete(fit); /* Illegal min-degree */ fit = cpl_fit_imagelist_polynomial(vdit, input, -1, 0, CPL_FALSE, pixeltype, NULL); cpl_test_error(CPL_ERROR_ILLEGAL_INPUT); cpl_test_null( fit ); cpl_imagelist_delete(fit); /* Illegal pixeltype */ fit = cpl_fit_imagelist_polynomial(vdit, input, 0, 0, CPL_FALSE, CPL_TYPE_INVALID, NULL); cpl_test_error(CPL_ERROR_UNSUPPORTED_MODE); cpl_test_null( fit ); cpl_imagelist_delete(fit); flops = cpl_tools_get_flops(); cpl_tools_get_cputime(CPL_CLOCK_START); /* Fit with zero-order term */ /* Also, try to use an integer-type image for fitting error */ fit = cpl_fit_imagelist_polynomial(vdit, input, 0, 2, CPL_TRUE, pixeltype, ifiterror); cpl_tools_get_cputime(CPL_CLOCK_STOP); flopsum += cpl_tools_get_flops() - flops; cpl_test_error(CPL_ERROR_NONE ); cpl_test( cpl_imagelist_get_size(fit) == 3 ); cpl_fit_image_is_zero(ifiterror, mytol); cpl_test_zero(cpl_image_subtract_scalar(cpl_imagelist_get(fit, 2), 1.0)); cpl_fit_imagelist_is_zero(fit, mytol); cpl_imagelist_delete(fit); flops = cpl_tools_get_flops(); cpl_tools_get_cputime(CPL_CLOCK_START); /* Fit with zero-order term */ /* Also, try to use an integer-type image for fitting error */ fit = cpl_fit_imagelist_polynomial(vdit, input, 0, ndits-1, CPL_TRUE, pixeltype, ifiterror); cpl_tools_get_cputime(CPL_CLOCK_STOP); flopsum += cpl_tools_get_flops() - flops; cpl_test_error(CPL_ERROR_NONE ); cpl_test( cpl_imagelist_get_size(fit) == ndits ); cpl_fit_image_is_zero(ifiterror, mytol); cpl_test_zero(cpl_image_subtract_scalar(cpl_imagelist_get(fit, 2), 1.0)); cpl_fit_imagelist_is_zero(fit, mytol); cpl_imagelist_delete(fit); flops = cpl_tools_get_flops(); cpl_tools_get_cputime(CPL_CLOCK_START); /* Fit without zero-order term */ fit = cpl_fit_imagelist_polynomial(vdit, input, 1, ndits-1, CPL_FALSE, pixeltype, dfiterror); cpl_tools_get_cputime(CPL_CLOCK_STOP); flopsum += cpl_tools_get_flops() - flops; cpl_test_error(CPL_ERROR_NONE ); cpl_test( cpl_imagelist_get_size(fit) == ndits-1 ); cpl_test_zero(cpl_image_subtract_scalar(cpl_imagelist_get(fit, 1), 1.0)); cpl_fit_imagelist_is_zero(fit, mytol); cpl_fit_image_is_zero(dfiterror, mytol); cpl_imagelist_delete(fit); flops = cpl_tools_get_flops(); cpl_tools_get_cputime(CPL_CLOCK_START); /* Fit with no zero- and 1st-order terms */ fit = cpl_fit_imagelist_polynomial(vdit, input, 2, ndits, CPL_TRUE, pixeltype, ffiterror); cpl_tools_get_cputime(CPL_CLOCK_STOP); flopsum += cpl_tools_get_flops() - flops; cpl_test_error(CPL_ERROR_NONE ); cpl_test( cpl_imagelist_get_size(fit) == ndits-1 ); cpl_test_zero(cpl_image_subtract_scalar(cpl_imagelist_get(fit, 0), 1.0)); cpl_fit_imagelist_is_zero(fit, mytol); cpl_fit_image_is_zero(ffiterror, mytol); cpl_imagelist_delete(fit); flops = cpl_tools_get_flops(); cpl_tools_get_cputime(CPL_CLOCK_START); /* Fit with one zero-term */ fit = cpl_fit_imagelist_polynomial(vdit, input, 0, 0, CPL_TRUE, pixeltype, dfiterror); cpl_tools_get_cputime(CPL_CLOCK_STOP); flopsum += cpl_tools_get_flops() - flops; cpl_test_error(CPL_ERROR_NONE ); cpl_test( cpl_imagelist_get_size(fit) == 1 ); cpl_test_zero(cpl_image_subtract_scalar(cpl_imagelist_get(fit, 0), sqsum/(double)ndits)); cpl_fit_imagelist_is_zero(fit, mytol); cpl_imagelist_delete(fit); cpl_imagelist_delete(input); (void)cpl_vector_unwrap(vdit); /* Try to fit as many coefficients are there are data points */ /* Also, use floats this time */ /* Use a polynomial to filler vdit */ i = 0; cpl_polynomial_set_coeff(vfiller, &i, 0.0); i = 1; cpl_polynomial_set_coeff(vfiller, &i, 1.0); input = cpl_imagelist_new(); vdit = cpl_vector_new(1); for (ntest = 1; ntest <= ndits; ntest++) { const double gain = 4.0; /* Some random number */ cpl_msg_info(cpl_func, "Fitting %d coefficients to as many points", ntest); cpl_vector_set_size(vdit, ntest); cpl_vector_fill_polynomial(vdit, vfiller, 0.0, 1.0); /* Create a list of images with a 2nd order function */ for (i = ntest - 1; i < ntest; i++) { cpl_image * image = cpl_image_new(IMAGESZFIT, IMAGESZFIT, pixeltype); cpl_test_zero(cpl_image_add_scalar(image, gain * ditval[i]*ditval[i])); cpl_test_zero(cpl_imagelist_set(input, image, i)); cpl_msg_debug(cpl_func, "Dit and mean of input image no. %d: %g %g", i, ditval[i], cpl_image_get_mean(image)); } /* Ready for fitting */ flops = cpl_tools_get_flops(); cpl_tools_get_cputime(CPL_CLOCK_START); /* Fit with zero-order term */ fit = cpl_fit_imagelist_polynomial(vdit, input, 0, ntest-1, CPL_TRUE, pixeltype, ffiterror); cpl_tools_get_cputime(CPL_CLOCK_STOP); flopsum += cpl_tools_get_flops() - flops; if (cpl_error_get_code() != CPL_ERROR_NONE) { cpl_test_null( fit ); cpl_msg_warning(cpl_func, "Could not fit %d coefficients to as many " "points", ntest); cpl_test_error(CPL_ERROR_SINGULAR_MATRIX ); break; } cpl_test( cpl_imagelist_get_size(fit) == ntest ); if (ntest == 2) { cpl_test_zero(cpl_image_subtract_scalar(cpl_imagelist_get(fit, 1), gain)); } else if (ntest > 2) { cpl_test_zero(cpl_image_subtract_scalar(cpl_imagelist_get(fit, 2), gain)); } cpl_fit_imagelist_is_zero(fit, ntest * mytol); cpl_fit_image_is_zero(ffiterror, ntest * mytol); cpl_imagelist_delete(fit); } cpl_imagelist_delete(input); cpl_vector_delete(vdit); /* Multiple samples at three different (equidistant) points */ input = cpl_imagelist_new(); vdit = cpl_vector_new(9); cpl_vector_set(vdit, 0, 5.0); cpl_vector_set(vdit, 1, 5.0); cpl_vector_set(vdit, 2, 5.0); cpl_vector_set(vdit, 3, 3.0); cpl_vector_set(vdit, 4, 3.0); cpl_vector_set(vdit, 5, 3.0); cpl_vector_set(vdit, 6, 7.0); cpl_vector_set(vdit, 7, 7.0); cpl_vector_set(vdit, 8, 7.0); cpl_imagelist_set(input, cpl_image_fill_test_create(IMAGESZFIT, IMAGESZFIT), 0); cpl_imagelist_set(input, cpl_image_fill_test_create(IMAGESZFIT, IMAGESZFIT), 1); cpl_imagelist_set(input, cpl_image_fill_test_create(IMAGESZFIT, IMAGESZFIT), 2); cpl_imagelist_set(input, cpl_image_fill_test_create(IMAGESZFIT, IMAGESZFIT), 3); cpl_imagelist_set(input, cpl_image_fill_test_create(IMAGESZFIT, IMAGESZFIT), 4); cpl_imagelist_set(input, cpl_image_fill_test_create(IMAGESZFIT, IMAGESZFIT), 5); cpl_imagelist_set(input, cpl_image_fill_test_create(IMAGESZFIT, IMAGESZFIT), 6); cpl_imagelist_set(input, cpl_image_fill_test_create(IMAGESZFIT, IMAGESZFIT), 7); cpl_imagelist_set(input, cpl_image_fill_test_create(IMAGESZFIT, IMAGESZFIT), 8); /* Fit a non-over-determined set without constant and linear terms */ flops = cpl_tools_get_flops(); cpl_tools_get_cputime(CPL_CLOCK_START); fit = cpl_fit_imagelist_polynomial(vdit, input, 2, 4, CPL_FALSE, CPL_TYPE_DOUBLE, dfiterror); cpl_tools_get_cputime(CPL_CLOCK_STOP); flopsum += cpl_tools_get_flops() - flops; cpl_test_error(CPL_ERROR_NONE ); cpl_test( cpl_imagelist_get_size(fit) == 3 ); cpl_fit_image_is_zero(dfiterror, DBL_MAX); cpl_imagelist_delete(fit); /* Fit a non-over-determined set with a constant term */ flops = cpl_tools_get_flops(); cpl_tools_get_cputime(CPL_CLOCK_START); fit = cpl_fit_imagelist_polynomial(vdit, input, 0, 2, CPL_FALSE, CPL_TYPE_DOUBLE, dfiterror); cpl_tools_get_cputime(CPL_CLOCK_STOP); flopsum += cpl_tools_get_flops() - flops; cpl_test_error(CPL_ERROR_NONE ); cpl_test( cpl_imagelist_get_size(fit) == 3 ); cpl_fit_image_is_zero(dfiterror, DBL_MAX); demax = cpl_image_get_sqflux(dfiterror); cpl_imagelist_delete(fit); /* Repeat the previous test - this time exploiting that the sampling points are equidistant */ flops = cpl_tools_get_flops(); cpl_tools_get_cputime(CPL_CLOCK_START); fit = cpl_fit_imagelist_polynomial(vdit, input, 0, 2, CPL_TRUE, CPL_TYPE_DOUBLE, dfiterror); cpl_tools_get_cputime(CPL_CLOCK_STOP); flopsum += cpl_tools_get_flops() - flops; cpl_test_error(CPL_ERROR_NONE ); cpl_test( cpl_imagelist_get_size(fit) == 3 ); cpl_fit_image_is_zero(dfiterror, DBL_MAX); cpl_test( cpl_image_get_sqflux(dfiterror) <= demax ); cpl_imagelist_delete(fit); /* Repeat the previous test - after sorting the sampling points */ cpl_vector_sort(vdit, 1); flops = cpl_tools_get_flops(); cpl_tools_get_cputime(CPL_CLOCK_START); fit = cpl_fit_imagelist_polynomial(vdit, input, 0, 2, CPL_TRUE, CPL_TYPE_DOUBLE, dfiterror); cpl_tools_get_cputime(CPL_CLOCK_STOP); flopsum += cpl_tools_get_flops() - flops; cpl_test_error(CPL_ERROR_NONE ); cpl_test( cpl_imagelist_get_size(fit) == 3 ); cpl_fit_image_is_zero(dfiterror, DBL_MAX); cpl_imagelist_delete(fit); /* Fit an under-determined set without a constant term */ fit = cpl_fit_imagelist_polynomial(vdit, input, 1, 4, CPL_FALSE, CPL_TYPE_DOUBLE, dfiterror); cpl_test_error(CPL_ERROR_SINGULAR_MATRIX ); cpl_test_null( fit ); cpl_imagelist_delete(input); cpl_vector_delete(vdit); secs = cpl_tools_get_cputime(CPL_CLOCK_TIME) - secs; cpl_msg_info("","Speed while fitting with image size %d in " "%g secs [Mflop/s]: %g", IMAGESZFIT, secs, flopsum/secs/1e6); cpl_tools_get_cputime(CPL_CLOCK_START); /* Done testing */ cpl_image_delete(dfiterror); cpl_image_delete(ffiterror); cpl_image_delete(ifiterror); cpl_polynomial_delete(vfiller); } /*----------------------------------------------------------------------------*/ /** @brief Verify that all elements in an imagelist are zero (within a tolerance) @param self The list of images to check @param tol The non-negative tolerance param line The line number of the caller @return void */ /*----------------------------------------------------------------------------*/ static void cpl_fit_imagelist_is_zero_macro(const cpl_imagelist * self, double tol, int line) { const int n = cpl_imagelist_get_size(self); int i; for (i = 0; i < n; i++) { /* FIXME: Need traceback of failure */ cpl_fit_image_is_zero_macro(cpl_imagelist_get_const(self, i), tol, line); } } /*----------------------------------------------------------------------------*/ /** @brief Verify that all elements in an image are zero (within a tolerance) @param self The image to check @param tol The non-negative tolerance param line The line number of the caller @return void */ /*----------------------------------------------------------------------------*/ static void cpl_fit_image_is_zero_macro(const cpl_image * self, double tol, int line) { cpl_stats * stats = cpl_stats_new_from_image(self, CPL_STATS_MIN | CPL_STATS_MAX | CPL_STATS_MEAN); const double mymin = cpl_stats_get_min(stats); const double mymax = cpl_stats_get_max(stats); /* FIXME: Need traceback of failure */ cpl_test(line > 0); cpl_test_leq( fabs(mymin), tol ); cpl_test_leq( fabs(mymax), tol ); cpl_stats_delete(stats); }