/* $Id: cpl_polynomial-test.c,v 1.80 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.80 $ * $Name: $ */ /*----------------------------------------------------------------------------- Includes -----------------------------------------------------------------------------*/ #ifdef HAVE_CONFIG_H #include #endif #include #include #include #include #include #include "cpl_polynomial.h" #include "cpl_error.h" #include "cpl_tools.h" #include "cpl_memory.h" #include "cpl_math_const.h" /*----------------------------------------------------------------------------- Defines -----------------------------------------------------------------------------*/ #define POLY_COEF_NB 5 #define VECTOR_SIZE 1024 #define POLY_SIZE 24 #define POLY_DIM 11 /*----------------------------------------------------------------------------- Main -----------------------------------------------------------------------------*/ int main(void) { cpl_polynomial * poly1; cpl_polynomial * poly2; cpl_polynomial * poly3; double x, y, z1, z2; double z, y_1, y2, x1, x2; cpl_bivector * xy_pos; cpl_vector * x_pos; cpl_vector * y_pos; cpl_vector * vec; cpl_vector * taylor; cpl_vector * xpoint; double * vec_data; const double xy_offset = 100; /* Some binomial coefficients */ double p15[8] = {1,15,105,455,1365,3003,5005,6435}; double xmax = 0; /* Maximum rounding error on x*/ double ymax = 0; /* Maximum rounding error on y*/ double zmax = 0; double eps; int expo[POLY_DIM]; int degree; int i, j; int k; cpl_test_init(PACKAGE_BUGREPORT, CPL_MSG_WARNING); /* FIXME: Move test of non-exported cpl_tools_ipow() to cpl_tools-test.c */ cpl_test_abs( cpl_tools_ipow(1.0, 0), 1.0, 0.0); cpl_test_abs( cpl_tools_ipow(0.0, 0), 1.0, 0.0); cpl_test_abs( cpl_tools_ipow(0.0, 1), 0.0, 0.0); cpl_test_abs( cpl_tools_ipow(0.0, 2), 0.0, 0.0); x = 3.14; y = 2.72; cpl_test_null( cpl_polynomial_extract(NULL, 0, NULL) ); cpl_test_error(CPL_ERROR_NULL_INPUT); cpl_test( cpl_polynomial_derivative(NULL, 0) ); cpl_test_error(CPL_ERROR_NULL_INPUT); /* Create a polynomial */ poly1 = cpl_polynomial_new(1); if (cpl_msg_get_level() <= CPL_MSG_DEBUG) cpl_polynomial_dump(poly1, stdout); cpl_test_eq( cpl_polynomial_get_degree(poly1), 0); /* Various tests regarding 0-degree polynomials */ cpl_test( cpl_polynomial_extract(poly1, 0, poly1) == NULL); cpl_test_error(CPL_ERROR_INVALID_TYPE); cpl_test( cpl_polynomial_derivative(poly1, 1) ); cpl_test_error(CPL_ERROR_ACCESS_OUT_OF_RANGE); cpl_test_zero( cpl_polynomial_derivative(poly1, 0) ); cpl_test_eq( cpl_polynomial_eval_1d(poly1, 5, &x), 0); cpl_test( x == 0 ); cpl_test_zero( cpl_polynomial_solve_1d(poly1, 5, &x, 1) ); cpl_test( x == 0 ); cpl_test_eq( cpl_polynomial_eval_1d(poly1, 5, &x), 0); cpl_test( x == 0 ); /* Properly react to a 0-degree polynomial with no roots */ i = 0; cpl_test_zero( cpl_polynomial_set_coeff(poly1, &i, 1) ); cpl_test_eq( cpl_polynomial_get_coeff(poly1, &i), 1 ); cpl_test_eq( cpl_polynomial_solve_1d(poly1, 1, &x, 1), CPL_ERROR_DIVISION_BY_ZERO ); cpl_test_error(CPL_ERROR_DIVISION_BY_ZERO); /* Equivalent to setting all coefficients to zero */ cpl_test_zero( cpl_polynomial_derivative(poly1, 0) ); cpl_test_eq( cpl_polynomial_eval_1d(poly1, 5, &x), 0); cpl_test( x == 0 ); i = VECTOR_SIZE; cpl_test_zero( cpl_polynomial_set_coeff(poly1, &i, 1) ); cpl_test_eq( cpl_polynomial_get_coeff(poly1, &i), 1 ); cpl_test_zero( cpl_polynomial_set_coeff(poly1, &i, 0) ); cpl_test_eq( cpl_polynomial_get_dimension(poly1), 1); cpl_test_eq( cpl_polynomial_get_degree(poly1), 0); cpl_test( poly2 = cpl_polynomial_duplicate(poly1) ); cpl_test_eq( cpl_polynomial_get_degree(poly2), 0); cpl_polynomial_delete(poly2); cpl_test( poly2 = cpl_polynomial_new(VECTOR_SIZE) ); cpl_test_eq( cpl_polynomial_get_dimension(poly2), VECTOR_SIZE ); cpl_test_eq( cpl_polynomial_get_degree(poly2), 0); cpl_test_zero( cpl_polynomial_copy(poly2, poly1)); cpl_test_eq( cpl_polynomial_get_dimension(poly2), 1 ); cpl_test_eq( cpl_polynomial_get_degree(poly2), 0); cpl_polynomial_delete(poly2); /* Set the coefficients : 1 + 2.x + 3.x^2 + 4.x^3 + 5.x^4 */ for (i=0; i < POLY_COEF_NB ; i++) { if (i % (POLY_COEF_NB-2) != 0) continue; cpl_test_zero( cpl_polynomial_set_coeff(poly1, &i, i+1) ); if (cpl_msg_get_level() <= CPL_MSG_DEBUG) cpl_polynomial_dump(poly1, stdout); } for (i=0; i < POLY_COEF_NB ; i++) { if (i % (POLY_COEF_NB-2) == 0) continue; if (cpl_msg_get_level() <= CPL_MSG_DEBUG) cpl_polynomial_dump(poly1, stdout); cpl_test_zero( cpl_polynomial_set_coeff(poly1, &i, i+1) ); } if (cpl_msg_get_level() <= CPL_MSG_DEBUG) cpl_polynomial_dump(poly1, stdout); /* Test cpl_polynomial_get_degree() */ cpl_test_eq( cpl_polynomial_get_degree(poly1), POLY_COEF_NB-1); /* Test cpl_polynomial_get_coeff() */ for (i=0; i < POLY_COEF_NB ; i++) cpl_test_eq( cpl_polynomial_get_coeff(poly1, &i), i+1 ); /* Test cpl_polynomial_eval() */ vec = cpl_vector_new(1); vec_data = cpl_vector_get_data(vec); vec_data[0] = x; z1 = cpl_polynomial_eval(poly1, vec); cpl_vector_delete(vec); /* Dump polynomials */ if (cpl_msg_get_level() <= CPL_MSG_DEBUG) cpl_polynomial_dump(poly1, stdout); /* Test cpl_polynomial_duplicate() and cpl_polynomial_copy() */ poly2 = cpl_polynomial_duplicate(poly1); cpl_test_zero( cpl_polynomial_compare(poly1, poly2, 0) ); cpl_polynomial_copy(poly1, poly2); cpl_test_zero( cpl_polynomial_compare(poly1, poly2, 0) ); /* Free */ cpl_polynomial_delete(poly1); cpl_polynomial_delete(poly2); /* Properly react to a polynomial with one real root at x = 1 and a proper (positive) 1st guess */ poly1 = cpl_polynomial_new(1); i = 0; cpl_test_zero( cpl_polynomial_set_coeff(poly1, &i, -1.0) ); i++; cpl_test_zero( cpl_polynomial_set_coeff(poly1, &i, 1.0) ); i = 2 * POLY_SIZE; cpl_test_zero( cpl_polynomial_set_coeff(poly1, &i, -1.0) ); i++; cpl_test_zero( cpl_polynomial_set_coeff(poly1, &i, 1.0) ); cpl_test_zero( cpl_polynomial_solve_1d(poly1, 5, &x, 1) ); cpl_test_abs( x, 1.0, 0.0 ); cpl_polynomial_delete(poly1); poly1 = cpl_polynomial_new(1); /* Create a tmporary vector - with a small length */ cpl_test_nonnull( vec = cpl_vector_new(POLY_SIZE) ); /* Fill the vector with a exp() Taylor series */ vec_data = cpl_vector_get_data(vec); cpl_test_nonnull( vec_data ); i = 0; vec_data[i] = 1; for (i=1 ; i < POLY_SIZE ; i++) { vec_data[i] = vec_data[i-1]/i; } for (i=POLY_SIZE-1 ; i >= 0; i--) { cpl_test_zero( cpl_polynomial_set_coeff(poly1, &i, vec_data[i]) ); } cpl_test_eq( cpl_polynomial_get_degree(poly1), POLY_SIZE-1); cpl_test_abs(cpl_polynomial_eval_1d(poly1, 0, &y), 1.0, DBL_EPSILON); cpl_test_abs( y, 1.0, DBL_EPSILON); /* See how far away from zero the approximation holds */ x = DBL_EPSILON; while ( fabs(cpl_polynomial_eval_1d(poly1, x, NULL)-exp( x)) < DBL_EPSILON * exp( x) && fabs(cpl_polynomial_eval_1d(poly1,-x, &y)-exp(-x)) < DBL_EPSILON * exp(-x) ) { /* Differentation of exp() does not change anything - but in the case of a Taylor approximation one term is lost and with it a bit of precision */ cpl_test_rel( exp(-x), y, 7.5 * DBL_EPSILON); x *= 2; } x /= 2; /* FIXME: Verify the correctness of this test */ cpl_test_leq( -x, -FLT_EPSILON); /* OK for POLY_SIZE >= 4 */ z2 = 2 * x / VECTOR_SIZE; /* Evaluate a sequence of exp()-approximations */ taylor = cpl_vector_new(VECTOR_SIZE); cpl_test_zero( cpl_vector_fill_polynomial(taylor, poly1, -x, z2)); vec_data = cpl_vector_get_data(taylor); cpl_test_nonnull( vec_data ); for (i=0 ; i < VECTOR_SIZE ; i++) { const double xapp = -x + i * z2; const double yapp = exp( xapp ); /* cpl_vector_fill_polynomial() is just a wrapper */ cpl_test_abs( vec_data[i], cpl_polynomial_eval_1d(poly1, xapp, &y), 2*DBL_EPSILON); if ( fabs(y - yapp ) > ymax * yapp ) ymax = fabs(y - yapp ) / yapp; if ( fabs(vec_data[i] - yapp ) > xmax * yapp ) xmax = fabs(vec_data[i] - yapp ) / yapp; } cpl_msg_info("","Rounding on %d-term Taylor-exp() in range [%g; %g] " "[DBL_EPSILON]: %7.5f %7.5f", POLY_SIZE, -x, x, xmax/DBL_EPSILON, ymax/DBL_EPSILON); cpl_test_leq( xmax, 2.0*DBL_EPSILON); cpl_test_leq( xmax, ymax); cpl_test_leq( ymax, 7.39*DBL_EPSILON); vec_data = cpl_vector_get_data(vec); eps = vec_data[0]; /* Solve p(y) = exp(x/2) - i.e. compute a logarithm */ i = 0; vec_data[i] -= exp(x/2); cpl_polynomial_set_coeff(poly1, &i, vec_data[i]); cpl_test_zero( cpl_polynomial_solve_1d(poly1, -x * x, &y, 1) ); vec_data[0] = eps; cpl_polynomial_set_coeff(poly1, &i, vec_data[i]); /* Check solution - allow up to 2 * meps rounding */ cpl_test_rel( y, x/2, 2.0*DBL_EPSILON); /* Check Residual - allow up to 2 * meps rounding */ cpl_test_rel( cpl_polynomial_eval_1d(poly1, y, NULL), exp(x/2), 2.0 * DBL_EPSILON ); /* Free */ cpl_vector_delete(vec); cpl_vector_delete(taylor); cpl_polynomial_delete(poly1); poly1 = cpl_polynomial_new(1); i = 0; cpl_polynomial_set_coeff(poly1, &i, 2); i++; cpl_polynomial_set_coeff(poly1, &i, 2); i++; cpl_polynomial_set_coeff(poly1, &i, 1); cpl_test_eq( cpl_polynomial_get_degree(poly1), 2); /* Properly react on a polynomial with no real roots */ cpl_test_eq( cpl_polynomial_solve_1d(poly1, 0, &x, 1), CPL_ERROR_DIVISION_BY_ZERO ); i = 0; cpl_polynomial_set_coeff(poly1, &i, 4); /* Properly react on a polynomial with no real roots */ cpl_test_eq( cpl_polynomial_solve_1d(poly1, 0, &x, 1), CPL_ERROR_CONTINUE ); cpl_test_error(CPL_ERROR_CONTINUE); cpl_polynomial_delete(poly1); poly1 = cpl_polynomial_new(1); /* The simplest 15-degree polynomial */ i =15; cpl_polynomial_set_coeff(poly1, &i, 1); cpl_test_eq( cpl_polynomial_get_degree(poly1), i); cpl_test_zero( cpl_polynomial_solve_1d(poly1, 10, &x, i) ); cpl_test_abs( x, 0.0, DBL_EPSILON); cpl_test_abs( cpl_polynomial_eval_1d(poly1, x, &y), 0.0, DBL_EPSILON); cpl_test_abs( y, 0.0, DBL_EPSILON); /* -1 is a root with multiplicity 15 */ x1 = -1; poly2 = cpl_polynomial_duplicate(poly1); cpl_test_zero( cpl_polynomial_shift_1d(poly2, -x1) ); /* poly2 now holds the binomial coefficients for n = 15 */ i = 0; for (i = 0; i < 8; i++) { cpl_polynomial_set_coeff(poly1, &i, p15[i]); j = 15 - i; cpl_polynomial_set_coeff(poly1, &j, p15[i]); } i = 15; cpl_test_eq( cpl_polynomial_get_degree(poly1), i); for (i = 0; i < 8; i++) { j = 15 - i; z1 = cpl_polynomial_get_coeff(poly1, &i); z2 = cpl_polynomial_get_coeff(poly2, &i); z = cpl_polynomial_get_coeff(poly2, &j); cpl_test_rel( z1, z2, 1.0 * DBL_EPSILON); cpl_test_rel( z, z2, 1.0 * DBL_EPSILON); } cpl_test_zero( cpl_polynomial_compare(poly1,poly2,0) ); i = 15; cpl_polynomial_solve_1d(poly1, 10*x1, &x, i); z = cpl_polynomial_eval_1d(poly1, x, &y); cpl_msg_info("", "(X+1)^15 (%d): %g %g %g", i, (x-x1)/DBL_EPSILON, z/DBL_EPSILON, y/DBL_EPSILON); cpl_test_rel( x, x1, DBL_EPSILON ); cpl_test_abs( z, 0.0, DBL_EPSILON ); cpl_test_abs( y, 0.0, DBL_EPSILON); /* Lots of round-off here, which depends on the long double precision */ i = 5; cpl_test_zero( cpl_polynomial_solve_1d(poly1, -10*x1, &x, i) ); z = cpl_polynomial_eval_1d(poly1, x, &y); cpl_msg_info("", "(X+1)^15 (%d) " CPL_LDBL_EPSILON_STR ": %g %g %g", i, x-x1, (double)(z/CPL_LDBL_EPSILON), (double)(y/CPL_LDBL_EPSILON)); cpl_test_rel( x, x1, 0.35 ); /* alphaev56 */ cpl_test_abs( z, 0.0, 18202.0 * DBL_EPSILON); cpl_test_abs( y, 0.0, 1554616.0 * DBL_EPSILON); /* alphaev56 */ i = 15; eps = 2 * DBL_EPSILON; cpl_polynomial_set_coeff(poly1, &i, 1 + eps); cpl_test( cpl_polynomial_compare(poly1,poly2,-eps/2) < 0); cpl_test_error(CPL_ERROR_ILLEGAL_INPUT); cpl_test_zero( cpl_polynomial_compare(poly1,poly2,eps) ); cpl_test( cpl_polynomial_compare(poly1,poly2,eps/2) > 0); cpl_polynomial_delete(poly2); poly2 = cpl_polynomial_new(POLY_DIM); if (cpl_msg_get_level() <= CPL_MSG_DEBUG) cpl_polynomial_dump(poly2, stdout); cpl_test_eq( cpl_polynomial_get_dimension(poly2), POLY_DIM); cpl_test_eq( cpl_polynomial_get_degree(poly2), 0); /* Set and reset some (outragous) coefficient */ for (j = 0; j =0;i--) cpl_test( cpl_polynomial_get_coeff(poly1, &i) == 0 ); cpl_test_zero( cpl_polynomial_copy(poly1, poly2) ); cpl_test_eq( cpl_polynomial_get_dimension(poly2), POLY_DIM); cpl_test_zero( cpl_polynomial_compare(poly1, poly2, 0) ); i = 12; for (j = 0; j 0 ); cpl_test_zero( cpl_polynomial_set_coeff(poly2, expo, cpl_polynomial_get_coeff(poly1, expo)) ); cpl_test_zero( cpl_polynomial_compare(poly1, poly2, 0) ); x = cpl_polynomial_get_coeff(poly1, expo); for (j = 0; j 0 ); cpl_polynomial_delete(poly2); poly3 = cpl_polynomial_new(POLY_DIM-1); for (j = 0; j 0) cpl_test_leq( eps/xmax, 9016); xmax = eps; i = 0; eps = cpl_polynomial_get_coeff(poly1, &i); if (fabs(eps) > fabs(zmax)) { cpl_msg_info("FIXME","1D-Polynomial with degree %d fit of a %d-point " "dataset has a greater error on P0 than a %d-degree " "fit to a %d-point dataset: fabs(%g) > fabs(%g)", j, j-1, j-2, j-1, eps, zmax); } /* Should loose at most one decimal digit per degree fabs(eps) < DBL_EPSILON * 10 ** (degree-2) */ cpl_test_leq( fabs(eps), DBL_EPSILON * pow(10, j-3) ); zmax = eps; /* Compute approximation to sqrt(0.25) - this will systematically be too low, but less and less so with increasing degree - until the approximation goes bad */ eps = cpl_polynomial_eval_1d(poly1, 0.25, NULL); cpl_polynomial_delete(poly1); if (j < 18) cpl_test_leq( fabs(0.5-eps), fabs(0.5 - ymax)); if (eps <= ymax) { /* Should be able to fit at least an 18-degree polynomial with increased accuracy - and without error-margin */ cpl_msg_info("FIXME","1D-Polynomial with degree %d fit of a %d-" "point dataset has a greater residual than a %d-" "degree fit to a %d-point dataset: %g > %g", j, j-1, j-2, j-1, 0.5-eps, 0.5-ymax); break; } ymax = eps; } /* And the mse itself must be bound */ cpl_test_leq( eps, 0.411456 * cpl_error_margin); /* Try to fit increasing degrees to a sqrt() */ cpl_test_zero( cpl_vector_set_size(vec, VECTOR_SIZE) ); cpl_test_zero( cpl_vector_set_size(taylor, VECTOR_SIZE) ); for (i=0; i < VECTOR_SIZE; i++) { cpl_test_zero( cpl_vector_set(vec, i, i + 1e4) ); cpl_test_zero( cpl_vector_set(taylor, i, sqrt((double)i)) ); } eps = FLT_MAX; for (i = 0; i < VECTOR_SIZE-1; i++) { xmax = eps; poly1 = cpl_polynomial_fit_1d_create(vec, taylor, i, &eps); cpl_test_nonnull( poly1 ); cpl_test_error( CPL_ERROR_NONE ); cpl_polynomial_delete(poly1); if (eps > xmax) { /* Should be able to fit at least a 8-degree polynomial with no degradation - and without error-margin (i686 can manage one degree more) */ if (i < 9) cpl_test_leq( eps, xmax); cpl_msg_info("FIXME","1D-Polynomial with degree %d fit of a %d-" "point dataset has a higher mean square error than a " "%d-degree fit: %g > %g", VECTOR_SIZE, i, i-1, eps, xmax); break; } } cpl_vector_delete(vec); cpl_vector_delete(taylor); /* Try to fit increasing degrees to f(x,y) = sqrt(x)*log(1+y) - with exactly as many points as there are coefficient to determine thus the residual should be exactly zero. */ poly1 = NULL; xpoint = cpl_vector_new(2); /* f(1/4,sqrt(e)-1) = 1/4 */ cpl_vector_set(xpoint, 0, 0.25+xy_offset); cpl_vector_set(xpoint, 1, sqrt(CPL_MATH_E)-1+xy_offset); xmax = 0; ymax = 0; for (degree = 0; degree < POLY_SIZE; degree++) { const int nc = (degree+1)*(degree+2)/2; k = 0; vec = cpl_vector_new(nc); xy_pos = cpl_bivector_new(nc); x_pos = cpl_bivector_get_x(xy_pos); y_pos = cpl_bivector_get_y(xy_pos); for (i=0; i <= degree; i++) { for (j=0; j <= i; j++, k++) { cpl_test_zero( cpl_vector_set(x_pos, k, i+xy_offset) ); cpl_test_zero( cpl_vector_set(y_pos, k, j+xy_offset) ); cpl_test_zero( cpl_vector_set(vec, k, sqrt((double)i)*log((double)(j+1)))); } } cpl_test( k == nc ); poly2 = cpl_polynomial_fit_2d_create(xy_pos, vec, degree, &eps); cpl_test_nonnull( poly2 ); cpl_test_error( CPL_ERROR_NONE ); if (cpl_msg_get_level() <= CPL_MSG_DEBUG) cpl_polynomial_dump(poly2, stdout); /* The increase in mse must be bound */ if (xmax > 0) cpl_test_leq( eps/xmax, 5.39953e+09); cpl_vector_delete(vec); cpl_bivector_delete(xy_pos); cpl_msg_info("FIXME","2D-Polynomial with degree %d fit of a %d-point " "dataset has a mean square error ratio to a %d-degree " "fit: %g (%g > %g)", degree, nc, degree-1, xmax > 0 ? eps/xmax : 0, eps, xmax); xmax = eps; eps = cpl_polynomial_eval(poly2, xpoint); cpl_polynomial_delete(poly1); poly1 = poly2; if (degree < 8) cpl_test_leq( fabs(0.25-eps), fabs(0.25 - ymax)); if (fabs(0.25-eps) >= fabs(0.25 - ymax) && degree > 0) { /* Should be able to fit at least a 5-degree polynomial with increased accuracy - and without error-margin */ cpl_msg_info("FIXME","2D-Polynomial with degree %d fit of a %d-" "point dataset has a greater residual than a %d-" "degree fit to a %d-point dataset: fabs(%g) > " "fabs(%g)", degree, nc, degree-1, degree*(degree+1)/2, eps-0.25, ymax-0.25); break; } ymax = eps; } /* Try to fit increasing degrees to f(x,y) = sqrt(x)*log(1+y) - with a constant, overdetermined number of points - The residual should decrease with increasing degree until the system becomes near singular */ /* f(1/4,sqrt(e)-1) = 1/4 */ cpl_vector_set(xpoint, 0, 0.25+xy_offset); cpl_vector_set(xpoint, 1, sqrt(CPL_MATH_E)-1-xy_offset); k = 0; vec = cpl_vector_new(2 * POLY_SIZE * 2 * POLY_SIZE); xy_pos = cpl_bivector_new(2 * POLY_SIZE * 2 * POLY_SIZE); x_pos = cpl_bivector_get_x(xy_pos); y_pos = cpl_bivector_get_y(xy_pos); for (i=0, k = 0; i < 2 * POLY_SIZE; i++) { for (j=0; j < 2 * POLY_SIZE; j++, k++) { x = (double) i * 0.5; y = (double) j * 2.0; cpl_test_zero( cpl_vector_set(x_pos, k, i+xy_offset) ); cpl_test_zero( cpl_vector_set(y_pos, k, j-xy_offset) ); cpl_test_zero( cpl_vector_set(vec, k, sqrt(x)*log(y+1.0))); } } cpl_test( k == 2 * POLY_SIZE * 2 * POLY_SIZE ); ymax = 0; for (degree = 0; degree < POLY_SIZE; degree++) { const int nc = 2 * POLY_SIZE * 2 * POLY_SIZE; double mse; poly2 = cpl_polynomial_fit_2d_create(xy_pos, vec, degree, &mse); cpl_test_nonnull( poly2 ); cpl_test_error( CPL_ERROR_NONE ); if (cpl_msg_get_level() <= CPL_MSG_DEBUG) cpl_polynomial_dump(poly2, stdout); eps = cpl_polynomial_eval(poly2, xpoint); cpl_msg_debug("FIXME","2D-Polynomial with degree %d fit of a %d-point " "dataset has a mean square error: %g (P0-error=%g)", degree, nc, mse, 0.25-eps); /* mse must decrease */ if (degree > 0) { if (fabs(eps-0.25) > fabs(ymax - 0.25)) cpl_msg_info("FIXME","2D-Polynomial with degree %d fit of a %d-" "point dataset has a larger error than a %d-degree " "fit: fabs(%g-0.25) > fabs(%g-0.25)", degree, nc, degree - 1, eps, ymax); if (mse > xmax) { cpl_msg_info("FIXME","2D-Polynomial with degree %d fit of a %d-" "point dataset has a larger mean square error than" " a %d-degree fit: %g > %g", degree, nc, degree - 1, mse, xmax); break; } } xmax = mse; ymax = eps; cpl_polynomial_delete(poly1); poly1 = poly2; } cpl_vector_delete(vec); cpl_bivector_delete(xy_pos); cpl_vector_delete(xpoint); cpl_polynomial_delete(poly2); /* A small performance test - twice to test DFS02166 */ for (k=0; k < 1 ; k++) { vec = cpl_vector_new(POLY_DIM); cpl_test_zero( cpl_vector_fill(vec, 1)); poly2 = cpl_polynomial_new(POLY_DIM); for (i = 0; i < POLY_SIZE; i++) { for (j = 0; j < POLY_DIM; j++) expo[j] = j + i; cpl_test_zero( cpl_polynomial_set_coeff(poly2, expo, 1)); } for (i = 0; i < POLY_SIZE*VECTOR_SIZE; i++) { eps = cpl_polynomial_eval(poly2, vec); } cpl_polynomial_delete(poly2); cpl_vector_delete(vec); } for (k=0; k < 1 ; k++) { vec = cpl_vector_new(POLY_DIM); cpl_test_zero( cpl_vector_fill(vec, 1)); poly2 = cpl_polynomial_new(POLY_DIM); for (i = 0; i < POLY_SIZE; i++) { for (j = 0; j < POLY_DIM; j++) expo[j] = j + i; cpl_test_zero( cpl_polynomial_set_coeff(poly2, expo, 1)); } for (i = 0; i < POLY_SIZE*VECTOR_SIZE; i++) { eps = cpl_polynomial_eval(poly2, vec); } cpl_polynomial_delete(poly2); cpl_vector_delete(vec); } #if 0 if (1) { /* A performance test */ cpl_flops flops; double secs; const int nc = 5 * POLY_SIZE; k = 0; vec = cpl_vector_new(nc*nc); xy_pos = cpl_bivector_new(nc*nc); x_pos = cpl_bivector_get_x(xy_pos); y_pos = cpl_bivector_get_y(xy_pos); for (i=0, k = 0; i < nc; i++) { for (j=0; j < nc; j++, k++) { const double x = (double) i * 0.5; const double y = (double) j * 2.0; cpl_test_zero( cpl_vector_set(x_pos, k, i+xy_offset) ); cpl_test_zero( cpl_vector_set(y_pos, k, j-xy_offset) ); cpl_test_zero( cpl_vector_set(vec, k, sqrt(x)*log(y+1.0))); } } cpl_test( k == nc * nc ); cpl_tools_get_cputime(CPL_CLOCK_STOP); flops = cpl_tools_get_flops(); secs = cpl_tools_get_cputime(CPL_CLOCK_TIME); for (degree = 0; degree < POLY_SIZE; degree++) { double mse; cpl_tools_get_cputime(CPL_CLOCK_START); poly2 = cpl_polynomial_fit_2d_create(xy_pos, vec, degree, &mse); cpl_tools_get_cputime(CPL_CLOCK_STOP); if ( poly2 == NULL ) break; if (cpl_msg_get_level() <= CPL_MSG_DEBUG) cpl_polynomial_dump(poly2, stdout); cpl_polynomial_delete(poly2); if (degree > 0 && mse > xmax) break; xmax = mse; } flops = cpl_tools_get_flops() - flops; secs = cpl_tools_get_cputime(CPL_CLOCK_TIME) - secs; cpl_tools_get_cputime(CPL_CLOCK_START); cpl_vector_delete(vec); cpl_bivector_delete(xy_pos); cpl_msg_info("","Speed while fitting %d points with up to degree %d in " "%g secs [Mflop/s]: %g", nc*nc, degree, secs, flops/secs/1e6); } #endif poly3 = cpl_polynomial_new(1); cpl_test_null( cpl_polynomial_extract(poly1, -1, poly3) ); cpl_test_error(CPL_ERROR_ILLEGAL_INPUT); cpl_test_null( cpl_polynomial_extract(poly1, 2, poly3) ); cpl_test_error(CPL_ERROR_ACCESS_OUT_OF_RANGE); if (cpl_msg_get_level() <= CPL_MSG_DEBUG) cpl_polynomial_dump(poly1, stdout); poly2 = cpl_polynomial_extract(poly1, 0, poly3); cpl_test( cpl_polynomial_get_dimension(poly2) == cpl_polynomial_get_dimension(poly1) - 1); /* poly1 comes from cpl_polynomial_fit_2d_create() */ cpl_test( cpl_polynomial_get_degree(poly2) == cpl_polynomial_get_degree(poly1)); if (cpl_msg_get_level() <= CPL_MSG_DEBUG) cpl_polynomial_dump(poly2, stdout); i = 0; eps = cpl_polynomial_get_coeff(poly2, &i); cpl_polynomial_delete(poly2); poly2 = cpl_polynomial_extract(poly1, 1, poly3); cpl_test( cpl_polynomial_get_dimension(poly2) == cpl_polynomial_get_dimension(poly1) - 1); /* poly1 comes from cpl_polynomial_fit_2d_create() */ cpl_test( cpl_polynomial_get_degree(poly2) == cpl_polynomial_get_degree(poly1)); if (cpl_msg_get_level() <= CPL_MSG_DEBUG) cpl_polynomial_dump(poly2, stdout); i = 0; xmax = eps - cpl_polynomial_get_coeff(poly2, &i); cpl_polynomial_delete(poly2); cpl_test_leq( fabs(xmax), 0.0); /* Use this block to benchmark 2D cpl_polynomial_eval() */ vec = cpl_vector_new(2); cpl_test_zero( cpl_vector_fill(vec, 1)); for (k=0; k < VECTOR_SIZE*VECTOR_SIZE; k++) eps = cpl_polynomial_eval(poly1, vec); cpl_vector_delete(vec); cpl_polynomial_delete(poly1); /* Create and differentiate a 1d-polynomial with uniform coefficients */ poly1 = cpl_polynomial_new(1); x = CPL_MATH_PI_4; for (i=0; i < POLY_SIZE; i++) { cpl_test_zero( cpl_polynomial_set_coeff(poly1, &i, x)); } x1 = x; for (i=POLY_SIZE-1; i > 0; i--) { x1 *= i; cpl_test_zero( cpl_polynomial_derivative(poly1, 0) ); } cpl_test( cpl_polynomial_get_degree(poly1) == 0); /* Verify the value of the constant term */ cpl_test_leq( fabs(x1 - cpl_polynomial_get_coeff(poly1, &i)), 0.0); cpl_polynomial_delete(poly1); /* Create and collapse and differentiate a 2d-polynomial with uniform coefficients */ poly1 = cpl_polynomial_new(2); x = CPL_MATH_PI_4; y = CPL_MATH_E; x1 = 1.0; for (i=0; i < POLY_SIZE; i++) { for (j=0; j < POLY_SIZE; j++) { expo[0] = i; expo[1] = j; cpl_test_zero( cpl_polynomial_set_coeff(poly1, expo, x)); } if (i > 0) x1 *= i; } i = 0; cpl_test_zero( cpl_polynomial_set_coeff(poly3, &i, y)); for (j = 0; j < 2; j++) { poly2 = cpl_polynomial_extract(poly1, j, poly3); z1 = 0.0; for (i=0; i < POLY_SIZE; i++) { z1 += pow(y, i); } z1 *= x; for (i=0; i < POLY_SIZE; i++) { z = cpl_polynomial_get_coeff(poly2, &i); cpl_test_leq( fabs(z - z1 ), 2*POLY_SIZE*FLT_EPSILON); } cpl_test_zero( cpl_polynomial_derivative(poly2, 0) ); i = 0; z = cpl_polynomial_get_coeff(poly2, &i); for (i=1; i < POLY_SIZE-1; i++) { z1 = cpl_polynomial_get_coeff(poly2, &i); cpl_test_leq( fabs(z * (i+1) - z1 ), 11*POLY_SIZE*FLT_EPSILON); } cpl_polynomial_delete(poly2); } cpl_polynomial_delete(poly3); for (i=POLY_SIZE-1; i > 0; i--) { cpl_test_zero( cpl_polynomial_derivative(poly1, 0) ); cpl_test_zero( cpl_polynomial_derivative(poly1, 1) ); } cpl_test( cpl_polynomial_get_degree(poly1) == 0); expo[0] = expo[1] = 0; x2 = cpl_polynomial_get_coeff(poly1, expo); /* The constant term is huge at this point - reduce by a factor eps which is inaccurate */ eps = x2 * DBL_EPSILON > 1 ? x2 * DBL_EPSILON : 1; xmax = x2/eps - x * x1 * x1/eps; cpl_test_leq( fabs(xmax), 1.0); cpl_polynomial_delete(poly1); /* End of tests */ return cpl_test_end(0); }