/* $Id: cpl_vector-test.c,v 1.57 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.57 $
 * $Name:  $
 */

/*-----------------------------------------------------------------------------
                                   Includes
 -----------------------------------------------------------------------------*/

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <float.h>
#include <assert.h>
#include <math.h>

#include "cpl_tools.h"
#include "cpl_math_const.h"
#include "cpl_memory.h"

#define VECTOR_SIZE 1024
#define VECTOR_CUT 24

/* alphaev56 SIGFPE's with more than 20 */
#define POLY_SIZE 20

/*-----------------------------------------------------------------------------
                        Private function prototypes
 -----------------------------------------------------------------------------*/

static double cpl_vector_get_diff(const cpl_vector *, const cpl_vector *);

static void cpl_vector_save_bench(int);

/*-----------------------------------------------------------------------------
                                  Main
 -----------------------------------------------------------------------------*/
int main(void)
{
    double        xc;
    double        emax = 0; /* Maximum observed xc-error */
    cpl_vector  * sinus;
    cpl_vector  * cosinus;
    cpl_vector  * tmp_vec;
    cpl_vector  * taylor;
    cpl_vector  * tmp_vec2;
    cpl_vector  * vxc;
    cpl_vector  * vxc1;
    cpl_vector  * vxc3;
    double      * data;
    FILE        * f_out;
    char          filename[1024];
    const int     vdif = VECTOR_SIZE - VECTOR_CUT > VECTOR_CUT
                       ? VECTOR_CUT : VECTOR_SIZE - VECTOR_CUT;
    const int     vdif2 = (VECTOR_SIZE - VECTOR_CUT)/2 > VECTOR_CUT/2
                        ? VECTOR_CUT/2 : (VECTOR_SIZE - VECTOR_CUT)/2;
    int           delta;
    int           half_search;
    int           i,k;
    cpl_boolean   do_bench;

    cpl_test_init(PACKAGE_BUGREPORT, CPL_MSG_WARNING);

    /* Insert tests below */

    do_bench = cpl_msg_get_level() <= CPL_MSG_INFO ? CPL_TRUE : CPL_FALSE;

    /* Create the vector sinus */
    cpl_test_nonnull( sinus = cpl_vector_new(VECTOR_SIZE) );
    
    /* Test cpl_vector_get_size() */
    cpl_test( cpl_vector_get_size(sinus) == VECTOR_SIZE );

    /* Fill the vector sinus */
    /* Test cpl_vector_get_data(), cpl_vector_set(), cpl_vector_get() */
    data = cpl_vector_get_data(sinus);
    cpl_test_nonnull( data );
    for (i=0 ; i < VECTOR_SIZE ; i++) {
        const double value = sin(i*CPL_MATH_2PI/VECTOR_SIZE);
        cpl_test_zero( cpl_vector_set(sinus, i, value) );
        cpl_test(  cpl_vector_get(sinus, i) == data[i] );
    }

    /* Create a Taylor-expansion of exp() */
    cpl_test_nonnull( taylor = cpl_vector_new(POLY_SIZE) );
    i = 0;
    cpl_vector_set(taylor, i, 1);
    for (i=1 ; i<POLY_SIZE ; i++)
        cpl_vector_set(taylor, i, cpl_vector_get(taylor, i-1)/i);

    /* Evaluate exp(sinus) using Horners scheme on the Taylor expansion */
    cpl_test( tmp_vec2 = cpl_vector_new(VECTOR_SIZE) );
    cpl_test_zero( cpl_vector_fill(tmp_vec2, cpl_vector_get(taylor, POLY_SIZE-1)) );
    
    for (k=POLY_SIZE-1; k > 0 ; k--) {
        cpl_test_zero( cpl_vector_multiply(tmp_vec2, sinus) );
        if (k&1) {
          cpl_test_zero( cpl_vector_add_scalar(tmp_vec2,
                                         cpl_vector_get(taylor, k-1)) );
        } else {
          cpl_test_zero( cpl_vector_subtract_scalar(tmp_vec2,
                                              -cpl_vector_get(taylor, k-1)) );
        }
    }        

    /* Verify the result (against cpl_vector_exponential() ) */
    cpl_test( tmp_vec = cpl_vector_duplicate(sinus) );

    cpl_test_zero( cpl_vector_exponential(tmp_vec, CPL_MATH_E) );

    cpl_test_zero( cpl_vector_subtract(tmp_vec2, tmp_vec) );
    cpl_test_zero( cpl_vector_divide(tmp_vec2, tmp_vec) );
    cpl_test_zero( cpl_vector_divide_scalar(tmp_vec2, DBL_EPSILON) );

    cpl_test_leq( fabs(cpl_vector_get_max(tmp_vec2)), 2.60831 );
    cpl_test_leq( fabs(cpl_vector_get_min(tmp_vec2)), 2.03626 );

    /* Evaluate exp() using cpl_vector_pow() on the Taylor expansion */
    cpl_test_zero( cpl_vector_fill(tmp_vec2, cpl_vector_get(taylor, 0)) );

    /* POLY_SIZE > 20 on alphaev56:
        Program received signal SIGFPE, Arithmetic exception.
          0x200000a3ff0 in cpl_vector_multiply_scalar ()
    */
    for (k=1; k < POLY_SIZE ; k++) {
        cpl_vector * vtmp = cpl_vector_duplicate(sinus);
        cpl_test_zero( cpl_vector_power(vtmp, k) );
        cpl_test_zero( cpl_vector_multiply_scalar(vtmp, cpl_vector_get(taylor, k)) );
        cpl_test_zero( cpl_vector_add(tmp_vec2, vtmp) );
        cpl_vector_delete(vtmp);
    }        

    /* Much less precise than Horner ... */
    cpl_test_leq(cpl_vector_get_diff(tmp_vec, tmp_vec2), 8 * DBL_EPSILON);

    cpl_vector_delete(taylor);

    /* Verify cpl_vector_logarithm() ) */
    cpl_test_zero( cpl_vector_logarithm(tmp_vec, CPL_MATH_E) );
    for (i=0 ; i < VECTOR_SIZE ; i++) {
        const double value = cpl_vector_get(sinus, i);
        double error = value - cpl_vector_get(tmp_vec,  i);
        if (2*i == VECTOR_SIZE) {
            /* value should really be zero */
          cpl_test_leq( fabs(value), 0.552 * DBL_EPSILON );
        } else {
          if (value != 0) error /= value;
        }
        cpl_test_leq( fabs(error), 43 * DBL_EPSILON );
    }

    /* Verify cpl_vector_power() */
    cpl_test_zero(  cpl_vector_copy(tmp_vec, sinus) );
    /* Just be positive */
    cpl_test_zero( cpl_vector_exponential(tmp_vec, CPL_MATH_E) );

    cpl_test_zero( cpl_vector_copy(tmp_vec2, tmp_vec) );

    cpl_test_zero( cpl_vector_sqrt(tmp_vec2) );
    cpl_test_zero( cpl_vector_power(tmp_vec, 0.5) );

    /* Necessary on AMD 64 (x86_64) Linux */
    cpl_test_leq(cpl_vector_get_diff(tmp_vec, tmp_vec2), 1.1*DBL_EPSILON);

    cpl_test_zero(  cpl_vector_copy(tmp_vec, sinus) );
    cpl_test_zero( cpl_vector_exponential(tmp_vec, CPL_MATH_E) );

    cpl_test_zero(  cpl_vector_multiply(tmp_vec2, tmp_vec) );

    cpl_test_zero( cpl_vector_power(tmp_vec, 1.5) );

    cpl_test_leq(cpl_vector_get_diff(tmp_vec, tmp_vec2), 8 * DBL_EPSILON);

    cpl_test_zero(  cpl_vector_copy(tmp_vec2, tmp_vec) );

    cpl_test_zero( cpl_vector_power(tmp_vec, 2) );

    cpl_test_zero(  cpl_vector_divide(tmp_vec2, tmp_vec) );

    cpl_test_zero( cpl_vector_power(tmp_vec, -0.5) );

    cpl_test_leq(cpl_vector_get_diff(tmp_vec, tmp_vec2), 8 * DBL_EPSILON);

    cpl_test_zero( cpl_vector_fill(tmp_vec, 1) );

    cpl_test_zero( cpl_vector_power(tmp_vec2, 0) );

    cpl_test_leq(cpl_vector_get_diff(tmp_vec, tmp_vec2), 0);

    cpl_vector_delete(tmp_vec);
    cpl_vector_delete(tmp_vec2);

    /* Test cpl_vector_dump() */
    sprintf(filename, "cpl_vector_dump.txt");
    cpl_test_nonnull( filename);
    cpl_test_nonnull( f_out = fopen(filename, "w") );

    /* Will not print any values */
    cpl_vector_dump(NULL, f_out);

    cpl_vector_dump(sinus, f_out);
    fclose(f_out);

    tmp_vec = cpl_vector_read(NULL);
    cpl_test_error(CPL_ERROR_NULL_INPUT);
    cpl_test_null( tmp_vec );

    tmp_vec = cpl_vector_read("/dev/null");
    cpl_test_error(CPL_ERROR_BAD_FILE_FORMAT);
    cpl_test_null( tmp_vec );

    /* Test cpl_vector_read() */
    tmp_vec = cpl_vector_read("cpl_vector_dump.txt");
    cpl_test( remove("cpl_vector_dump.txt") == 0 );
    cpl_test_nonnull( tmp_vec );

    cpl_test( cpl_vector_get_size(tmp_vec) == cpl_vector_get_size(sinus) );

    /* Test cpl_vector_save() / cpl_vector_load() */
    cpl_vector_save(tmp_vec, "cpl_vector_save.fits", CPL_BPP_IEEE_DOUBLE, NULL,
            CPL_IO_DEFAULT); 
    cpl_test_error( CPL_ERROR_NONE );
    tmp_vec2 = cpl_vector_load("cpl_vector_save.fits", 0) ;
    cpl_test( remove("cpl_vector_save.fits") == 0 );
    cpl_test_nonnull( tmp_vec2 );

    /* Verify that the save/load did not change the vector */
    cpl_test_leq(cpl_vector_get_diff(tmp_vec, tmp_vec2), 0.0);

    cpl_vector_delete(tmp_vec2);

    /* Loss of precision in cpl_vector_dump() */
    cpl_test_leq(cpl_vector_get_diff(tmp_vec, sinus), 10*FLT_EPSILON);

    cpl_vector_subtract(tmp_vec, sinus);

    /* Same loss for positive as for negative numbers */
    cpl_test_abs( cpl_vector_get_max(tmp_vec)+cpl_vector_get_min(tmp_vec), 0.0,
        2.5 * DBL_EPSILON);

    cpl_vector_delete(tmp_vec);

    /* Test cpl_vector_duplicate */
    tmp_vec = cpl_vector_duplicate(sinus);
    cpl_test_leq(cpl_vector_get_diff(tmp_vec, sinus), 0);

    /* Test fill function */
    cpl_test( cpl_vector_fill(tmp_vec, 1.0) == CPL_ERROR_NONE );
    cpl_test_leq( cpl_vector_get_mean(tmp_vec) - 1, DBL_EPSILON );

    /* Test extract function */
    tmp_vec2 = cpl_vector_extract(tmp_vec, 0, VECTOR_SIZE/2, 1);
    cpl_test_nonnull( tmp_vec2 );

    cpl_vector_delete(tmp_vec2);
    
    /* Create the vector cosinus */
    cosinus = cpl_vector_new(VECTOR_SIZE);
    cpl_test( cpl_vector_get_size(sinus) == VECTOR_SIZE );
    
    /* Fill the vector cosinus */
    data = cpl_vector_get_data(cosinus);
    cpl_test_nonnull( data );
    for (i=0 ; i<VECTOR_SIZE ; i++) data[i] = cos(i*CPL_MATH_2PI/VECTOR_SIZE);

    /* Test mean function */
    cpl_test_abs( cpl_vector_get_mean(cosinus), 0.0, 0.68*DBL_EPSILON );
    /* Test stdev function  (NB: the mean-value of cosinus-squared is 1/2) */
    cpl_test_leq( fabs(cpl_vector_get_stdev(cosinus) - sqrt(VECTOR_SIZE*0.5
                  /(VECTOR_SIZE - 1))), 0.036*DBL_EPSILON * sqrt(VECTOR_SIZE));

    /* Test copy function */
    cpl_test( cpl_vector_copy(tmp_vec, cosinus) == CPL_ERROR_NONE );
    cpl_test_leq(cpl_vector_get_diff(tmp_vec, cosinus), 0);

    cpl_vector_delete(tmp_vec);

    /* Test add & sub functions */
    tmp_vec = cpl_vector_duplicate(sinus);
    cpl_test_leq(cpl_vector_get_diff(tmp_vec, sinus), 0);
    cpl_vector_add(tmp_vec, cosinus);
    cpl_vector_subtract(tmp_vec, sinus);

    cpl_test_leq(cpl_vector_get_diff(tmp_vec, cosinus), DBL_EPSILON);

    /* Test cpl_vector_subtract_scalar() function */
    cpl_test( cpl_vector_subtract_scalar(tmp_vec, 2) == CPL_ERROR_NONE );

   /* Test div function */
    cpl_test( cpl_vector_divide(tmp_vec, tmp_vec) == CPL_ERROR_NONE );
    cpl_test_leq( cpl_vector_get_mean(tmp_vec) - 1, DBL_EPSILON );

    cpl_vector_delete(tmp_vec);

    /* Test dot-product - using orthogonal vectors and pythagoras */
    cpl_test_leq( cpl_vector_product(sinus, cosinus),
                       DBL_EPSILON*VECTOR_SIZE);

    cpl_test_leq( fabs(cpl_vector_product(  sinus,   sinus) - VECTOR_SIZE
                          + cpl_vector_product(cosinus, cosinus)),
                       DBL_EPSILON*VECTOR_SIZE );

    /* Test filtering */
    tmp_vec = cpl_vector_filter_lowpass_create(sinus, CPL_LOWPASS_LINEAR, 2);
    cpl_test( cpl_vector_get_size(tmp_vec) == cpl_vector_get_size(sinus) );
    cpl_vector_delete(tmp_vec);

    tmp_vec = cpl_vector_filter_median_create(sinus, 2);
    cpl_test( cpl_vector_get_size(tmp_vec) == cpl_vector_get_size(sinus) );
    cpl_vector_delete(tmp_vec);

    /* Test existence of cpl_vector_fit_gaussian() */
    cpl_test( cpl_vector_fit_gaussian(NULL, NULL, NULL, NULL, CPL_FIT_ALL,
					       NULL, NULL, NULL, NULL, NULL, NULL,
					       NULL) == CPL_ERROR_NULL_INPUT );

    /* sinus <- sinus*sinus */
    cpl_test( cpl_vector_multiply(sinus, sinus) == CPL_ERROR_NONE );

    /* Multiply by -1 */
    cpl_test( cpl_vector_multiply_scalar(sinus, -1) == CPL_ERROR_NONE );
    
    /* Add 1 */
    cpl_test( cpl_vector_add_scalar(sinus, 1) == CPL_ERROR_NONE );
   
    /* sinus <- sqrt(1-sinus^2) */
    cpl_test( cpl_vector_sqrt(sinus) == CPL_ERROR_NONE );

    /* Compute the absolute value of cosinus */
    data = cpl_vector_get_data(cosinus);
    cpl_test_nonnull( data );
    for (i=0 ; i<VECTOR_SIZE ; i++) data[i] = fabs(data[i]);

    /* Compare fabs(cosinus) with sqrt(1-sinus^2) */
    cpl_test_leq(cpl_vector_get_diff(sinus, cosinus), 1000*DBL_EPSILON);

    cpl_vector_copy(sinus, cosinus);

    cpl_vector_sort(cosinus, 1);
    for (i=1 ; i<VECTOR_SIZE ; i++) cpl_test( data[i] >= data[i-1] );

    cpl_vector_sort(sinus, -1);
    data = cpl_vector_get_data(sinus);
    cpl_test_nonnull( data );
    for (i=1 ; i<VECTOR_SIZE ; i++) cpl_test( data[i] <= data[i-1] );

    cpl_test_leq( fabs(cpl_vector_get_mean(cosinus)
                             - cpl_vector_get_mean(sinus)),
                             5.5*DBL_EPSILON );

    cpl_vector_delete(cosinus);
    cpl_vector_delete(sinus);

    /* Create the double-length vector sinus */
    cpl_test_nonnull( sinus = cpl_vector_new(2*VECTOR_SIZE) );
    
    /* Fill the vector sinus */
    data = cpl_vector_get_data(sinus);
    cpl_test_nonnull( data );
    for (i=0 ; i<2*VECTOR_SIZE ; i++) data[i] = sin(i*CPL_MATH_2PI/VECTOR_SIZE);


    /* Create the vector cosinus */
    cpl_test_nonnull( cosinus = cpl_vector_new(VECTOR_SIZE) );
    
    /* Fill the vector cosinus */
    data = cpl_vector_get_data(cosinus);
    cpl_test_nonnull( data );
    for (i=0 ; i<VECTOR_SIZE ; i++) data[i] = cos(i*CPL_MATH_2PI/VECTOR_SIZE);

    /* Create the vector tmp_vec */
    tmp_vec = cpl_vector_new(VECTOR_SIZE-1);
    cpl_test_nonnull( tmp_vec );
    cpl_test( cpl_vector_fill(tmp_vec, 1.0) == CPL_ERROR_NONE );

    vxc1 = cpl_vector_new(1);
    vxc3 = cpl_vector_new(3);
   /*
    Various error conditions
    */

    delta = cpl_vector_correlate(NULL, sinus, sinus);
    cpl_test_error( CPL_ERROR_NULL_INPUT );
    cpl_test( delta < 0 );

    delta = cpl_vector_correlate(vxc1, NULL, sinus);
    cpl_test_error( CPL_ERROR_NULL_INPUT );
    cpl_test( delta < 0 );

    delta = cpl_vector_correlate(vxc1, sinus, NULL);
    cpl_test_error( CPL_ERROR_NULL_INPUT );
    cpl_test( delta < 0 );

    delta = cpl_vector_correlate(cosinus, sinus, sinus);
    cpl_test_error( CPL_ERROR_ILLEGAL_INPUT );
    cpl_test( delta < 0 );

    delta = cpl_vector_correlate(vxc1, cosinus, sinus);
    cpl_test_error( CPL_ERROR_ILLEGAL_INPUT );
    cpl_test( delta < 0 );

    delta = cpl_vector_correlate(vxc3, cosinus, tmp_vec);
    cpl_test_error( CPL_ERROR_ILLEGAL_INPUT );
    cpl_test( delta < 0 );

    /* Normal conditions */
    cpl_tools_reset_flops();

    cpl_test( cpl_vector_correlate(vxc1, cosinus, tmp_vec) == 0);
    cpl_test( cpl_vector_get(vxc1, 0) == 0 );

    cpl_vector_delete(tmp_vec);

    cpl_test_zero( cpl_vector_multiply_scalar(sinus,   CPL_MATH_SQRT2) );
    cpl_test_zero( cpl_vector_add_scalar(     sinus,   CPL_MATH_PI   ) );
    cpl_test_zero( cpl_vector_multiply_scalar(cosinus, CPL_MATH_SQRT2) );

    cpl_test( cpl_vector_correlate(vxc1, sinus, sinus) == 0);
    /* without -O3 a zero-tolereance would be OK */
    cpl_test_leq( 1.0 - cpl_vector_get(vxc1, 0), 144.0*DBL_EPSILON );

    cpl_test( cpl_vector_correlate(vxc1, cosinus, cosinus) == 0 );
    xc = cpl_vector_get(vxc1, 0);
    cpl_test_leq( fabs(1-xc), 5 * DBL_EPSILON );
    if (fabs(1-xc) > emax) emax = fabs(1-xc);

    if (VECTOR_SIZE % 2 == 0) {
        /* Sinus and cosinus have zero cross-correlation with zero shift */
        cpl_test( cpl_vector_correlate(vxc1, sinus, cosinus) == 0 );
        xc = cpl_vector_get(vxc1, 0);
        cpl_test_leq( fabs(xc), 0.82*DBL_EPSILON );
        if (fabs(xc) > emax) emax = fabs(xc);
    }

    /*  cosinus and -cosinus have cross-correlation -1 with zero shift */
    tmp_vec = cpl_vector_duplicate(cosinus);
    cpl_test_leq(cpl_vector_get_diff(tmp_vec, cosinus), 0);
    cpl_test_zero( cpl_vector_divide_scalar(tmp_vec, -1) );

    cpl_test( cpl_vector_correlate(vxc1, tmp_vec, cosinus) == 0 );
    xc = cpl_vector_get(vxc1, 0);
    cpl_vector_delete(tmp_vec);
    cpl_test_leq( fabs(1+xc), 5*DBL_EPSILON );
    if (fabs(1+xc) > emax) emax = fabs(1+xc);

    vxc = cpl_vector_new( 1 );
    if (VECTOR_SIZE % 2 == 0) {
        /* Cross-correlation between sinus and cosinus grows to maximum
           at shift of pi/2 */
        for (i=0; i<VECTOR_SIZE/4; i++) {
            const double xcp = xc;
            half_search = i+1;

            cpl_test_zero( cpl_vector_set_size(vxc, 2*half_search + 1) );

            delta = cpl_vector_correlate(vxc, sinus, cosinus);
            xc = cpl_vector_get(vxc, delta);
            cpl_test( xc > xcp );
            cpl_test( abs(delta-(i+1)) == i+1 );
        }
        cpl_test_leq( fabs(1-xc), 16*DBL_EPSILON);
        half_search = VECTOR_SIZE/3;
        cpl_test_zero( cpl_vector_set_size(vxc, 2*half_search + 1) );
        delta = cpl_vector_correlate(vxc, sinus, cosinus);
        xc = cpl_vector_get(vxc, delta );
        cpl_test( abs(delta-VECTOR_SIZE/3) == VECTOR_SIZE/4 );
        if (fabs(1-xc) > emax) emax = fabs(1-xc);
    }

    cpl_vector_delete(sinus);

    /* Vectors of almost the same length - no commutativity */

    /* Create the vector sinus */
    cpl_test_nonnull( sinus = cpl_vector_new(VECTOR_SIZE-VECTOR_CUT) );
    
    /* Fill the vector sinus */
    data = cpl_vector_get_data(sinus);
    cpl_test_nonnull( data );
    for (i=0 ; i<cpl_vector_get_size(sinus) ; i++)
        data[i] = cos(i*CPL_MATH_2PI/VECTOR_SIZE);

    /* Compare with no shift - other than half the length difference */
    half_search = VECTOR_SIZE;
    cpl_test_zero( cpl_vector_set_size(vxc, 2*half_search + 1) );

    delta = cpl_vector_correlate(vxc, cosinus, sinus);
    xc = cpl_vector_get(vxc, delta);
    delta -= VECTOR_SIZE;
    cpl_test( delta + VECTOR_CUT/2 == 0);
    cpl_test_leq( fabs(1-xc), 6.5 * DBL_EPSILON );
    if (fabs(1-xc) > emax) emax = fabs(1-xc);

    /* Compare with increasing shift and increasing drop of elements
       - only up to the length-difference */
    for (k = 1; k < vdif; k++) {

        for (i=0 ; i<cpl_vector_get_size(sinus) ; i++)
            data[i] = cos((i+k)*CPL_MATH_2PI/VECTOR_SIZE);

        delta = cpl_vector_correlate(vxc, cosinus, sinus);
        xc = cpl_vector_get(vxc, delta);
        delta -= VECTOR_SIZE;
        cpl_test( delta + VECTOR_CUT/2 == k );
        cpl_test_leq( fabs(1-xc), 8.5 * DBL_EPSILON );
        if (fabs(1-xc) > emax) emax = fabs(1-xc);

    }

    /* Continue - maximum xc found with drop */
    for (; k < vdif; k++) {
        half_search = k-VECTOR_CUT/2;
        cpl_test_zero( cpl_vector_set_size(vxc, 2*half_search + 1) );

        for (i=0 ; i<cpl_vector_get_size(sinus) ; i++)
            data[i] = cos((i+k)*CPL_MATH_2PI/VECTOR_SIZE);

        delta = cpl_vector_correlate(vxc, cosinus, sinus);
        xc = cpl_vector_get(vxc, delta);
        delta -= half_search;
        cpl_test_leq( fabs(1-xc), 25 * DBL_EPSILON );
        cpl_test( delta + VECTOR_CUT/2 == k );
        if (fabs(1-xc) > emax) emax = fabs(1-xc);
    }

    /* Compare with increasing negative shift and increasing drop of elements
       - only up to half the length-difference */
    half_search = VECTOR_CUT;
    cpl_test_zero( cpl_vector_set_size(vxc, 2*half_search + 1) );
    xc = 1;
    for (k = 1; k < vdif2; k++) {
        const double xcp = xc;

        for (i=0 ; i<cpl_vector_get_size(sinus) ; i++)
            data[i] = cos((i-k)*CPL_MATH_2PI/VECTOR_SIZE);

        delta = cpl_vector_correlate(vxc, cosinus, sinus);
        xc = cpl_vector_get(vxc, delta);
        delta -= half_search;
        cpl_test( xc <= xcp );
        cpl_test( delta + k + VECTOR_CUT/2 >= 0 );
    }

    cpl_vector_delete(sinus);

    /* Vectors of the same length - commutativity */

    sinus = cpl_vector_duplicate(cosinus);
    cpl_test_leq(cpl_vector_get_diff(sinus, cosinus), 0);

    half_search = VECTOR_CUT;
    cpl_test_zero( cpl_vector_set_size(vxc, 2*half_search + 1) );
    delta = cpl_vector_correlate(vxc, sinus, cosinus);
    xc = cpl_vector_get(vxc, delta);
    delta -= half_search;
    cpl_test( delta == 0 );
    cpl_test_leq( fabs(1-xc), 5 * DBL_EPSILON );
    if (fabs(1-xc) > emax) emax = fabs(1-xc);

    /* Verify commutativity */
    cpl_test( delta+half_search == cpl_vector_correlate(vxc, cosinus, sinus) );
    cpl_test( xc == cpl_vector_get(vxc, delta+half_search) );

    data = cpl_vector_get_data(sinus);
    cpl_test_nonnull( data );

    half_search = VECTOR_SIZE/2;
    cpl_test_zero( cpl_vector_set_size(vxc, 2*half_search + 1) );
    /* Compare with increasing shift and increasing drop of elements
       -  delta tests will not hold for large shifts */
    xc = 1;
    for (k = 1; k < VECTOR_SIZE/50; k+=7) {
        const double xcp = xc;
        double xcn;

        for (i=0 ; i<VECTOR_SIZE ; i++) data[i] = cos((i+k)*CPL_MATH_2PI/VECTOR_SIZE);

        delta = cpl_vector_correlate(vxc, cosinus, sinus);
        xc = cpl_vector_get(vxc, delta);
        delta -= half_search;
        cpl_test( k == delta );
        cpl_test( xc < xcp );

        /* Commutativity */
        delta = cpl_vector_correlate(vxc, sinus, cosinus);
        xcn = cpl_vector_get(vxc, delta);
        delta -= half_search;
        cpl_test( k + delta == 0);
        cpl_test_leq( fabs(xcn - xc), 7*DBL_EPSILON); /* SUSE 9.0 */

        /* Shift in opposite direction, i.e. reverse sign on k */
        for (i=0 ; i<VECTOR_SIZE ; i++) data[i] = cos((i-k)*CPL_MATH_2PI/VECTOR_SIZE);

        delta = cpl_vector_correlate(vxc, cosinus, sinus);
        xc = cpl_vector_get(vxc, delta);
        delta -= half_search;
        cpl_test( k + delta == 0);
        cpl_test( xc < xcp );
        
    }

    half_search = VECTOR_SIZE;
    cpl_test_zero( cpl_vector_set_size(vxc, 2*half_search + 1) );

    /* Check with pseudo-random data */
    srand(1);
    for (i=0 ; i<VECTOR_SIZE ; i++) data[i] = 2.0*cpl_drand() - 1.0;

    cpl_vector_copy(cosinus, sinus);

    cpl_test( cpl_vector_correlate(vxc, cosinus, sinus) == half_search );
    /* without -O3 a zero-tolereance would be OK */
    cpl_test_leq( 1.0 - cpl_vector_get(vxc, half_search), 3.5*DBL_EPSILON );

    half_search = VECTOR_SIZE/2;
    cpl_test_zero( cpl_vector_set_size(vxc, 2*half_search + 1) );

#if 0
    /* Execute this to only measure time and FLOPS for the following code */
    cpl_tools_get_cputime(CPL_CLOCK_RESET);
    cpl_tools_reset_flops();
    cpl_tools_get_cputime(CPL_CLOCK_START);
#endif
    for (k = 2; k < VECTOR_SIZE-2; k+=2) {
        double * pcosinus;

        cpl_vector_delete(cosinus);
        cosinus  = cpl_vector_new(VECTOR_SIZE-k);
        pcosinus = cpl_vector_get_data(cosinus);
        cpl_test_nonnull( pcosinus );

        for (i=0 ; i<VECTOR_SIZE-k ; i++) pcosinus[i] = data[i];

        delta = cpl_vector_correlate(vxc, sinus, cosinus);
        xc = cpl_vector_get(vxc, delta);
        delta -= half_search;
        cpl_test( delta <= 0 );
        cpl_test_abs( xc, 1.0, 13.5 * DBL_EPSILON );
        if (fabs(1-xc) > emax) emax = fabs(1-xc);

    }

    cpl_msg_info("","Largest cross-correlation rounding error [DBL_EPSILON]: "
                 "%g", emax/DBL_EPSILON);

    if (do_bench) {
        /* cpl_msg_set_component_on(); */
        cpl_vector_save_bench(200);
        /* cpl_msg_set_component_off(); */
    } else {
        cpl_vector_save_bench(1);
    }

    /* Free and return */
    cpl_vector_delete(cosinus);
    cpl_vector_delete(sinus);

    cpl_vector_delete(vxc);
    cpl_vector_delete(vxc1);
    cpl_vector_delete(vxc3);

    /* End of tests */
    return cpl_test_end(0);
}

/**@}*/

/*----------------------------------------------------------------------------*/
/**
  @brief    Compute the inf-norm of a vector difference
  @param    v1          First cpl_vector
  @param    v2          Second cpl_vector of the same size as v1
  @return   The (non-negative) inf-norm of the vector difference
  @note     The function will fail assert() on illegal input.

  The inf-norm of a vector is the maximum absolute value of among its elements.
  
 */
/*----------------------------------------------------------------------------*/
static double cpl_vector_get_diff(
        const cpl_vector * v1,
        const cpl_vector * v2)
{

    const double *p1;
    const double *p2;

    double norm = 0;
    int n, i;

    /* Test inputs */
    assert(v1 != NULL);
    assert(v2 != NULL);

    n = cpl_vector_get_size(v1);

    assert(cpl_vector_get_size(v2) == n);

    p1 = cpl_vector_get_data_const(v1);
    p2 = cpl_vector_get_data_const(v2);

    assert( p1 );
    assert( p2 );

    for (i = 0; i < n; i++) {
        const double diff = fabs(p1[i] - p2[i]);
        if (diff > norm) norm = diff;
    }

    return norm;

}

/*----------------------------------------------------------------------------*/
/**
  @brief    Benchmark the CPL function
  @param n  The number of repeats
  @return   void
 */
/*----------------------------------------------------------------------------*/
static void cpl_vector_save_bench(int n)
{

    const int          nprops = 100;
    int                i;
    double             secs;
    const char       * filename = "cpl_vector_save.fits";
    const double       vval[] = {0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0};
    const int          nvals = (int)(sizeof(vval)/sizeof(double));
    const cpl_vector * vec = cpl_vector_wrap(nvals, (double*)vval);
    cpl_propertylist * qclist = cpl_propertylist_new();
    char               key[81];


    cpl_msg_info(cpl_func, "Benchmarking with %d %d-length vectors", n, nvals);

    for (i = 0; i < nprops; i++) {

        const int nlen = snprintf(key, 81, "ESO QC CARD%04d", i);

        cpl_test( nlen > 0 && nlen < 81);

        cpl_test_zero( cpl_propertylist_append_int(qclist, key, i));

    }

    cpl_test( cpl_propertylist_get_size(qclist) == nprops);


    cpl_tools_get_cputime(CPL_CLOCK_STOP);
    secs  = cpl_tools_get_cputime(CPL_CLOCK_TIME);
    cpl_tools_get_cputime(CPL_CLOCK_START);

    for (i = 0; i < n; i++) {

       cpl_test_zero( cpl_vector_save(vec, filename, CPL_BPP_IEEE_DOUBLE,
                                           qclist, CPL_IO_DEFAULT));
    }

    cpl_tools_get_cputime(CPL_CLOCK_STOP);
    secs = cpl_tools_get_cputime(CPL_CLOCK_TIME) - secs;
    cpl_tools_get_cputime(CPL_CLOCK_START);

    cpl_msg_info(cpl_func,"Time spent saving %d %d-sized vectors [s]: %g",
                 n, nvals, secs);

    cpl_test( remove(filename) == 0 );

    cpl_vector_unwrap((cpl_vector*)vec);
    cpl_propertylist_delete(qclist);

    return;

}