/* $Id: cpl_bivector-test.c,v 1.25 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.25 $
 * $Name:  $
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <float.h>

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

#define FUNCTION_SIZE  1024

/*-----------------------------------------------------------------------------
                                  Main
 -----------------------------------------------------------------------------*/
int main(void)
{
    cpl_bivector    *   sinus;
    cpl_bivector    *   cosinus;
    cpl_bivector    *   tmp_fun;
    double          *   data_x,
                    *   data_y;
    cpl_vector      *   vec1,
                    *   vec2;
    double              scale;
    double              error;
    FILE            *   f_out;
    char                filename[1024];
    int                 i;
    

    cpl_test_init(PACKAGE_BUGREPORT, CPL_MSG_WARNING);

    /* Insert actual test code here */

    scale = CPL_MATH_2PI / FUNCTION_SIZE;
    
    /* Create the first 1d function : sinus function */
    cpl_test( sinus = cpl_bivector_new(FUNCTION_SIZE) );

    /* Fill the sinus function */
    data_x = cpl_bivector_get_x_data(sinus);
    data_y = cpl_bivector_get_y_data(sinus);
    for (i=0 ; i<FUNCTION_SIZE ; i++) {
        data_x[i] = i * scale;
        data_y[i] = sin(data_x[i]);
    }
    
    /* Create the second 1d function : cosinus function */
    cpl_test (cosinus = cpl_bivector_new(FUNCTION_SIZE));

    /* Fill the cosinus function */
    data_x = cpl_bivector_get_x_data(cosinus);
    data_y = cpl_bivector_get_y_data(cosinus);
    for (i=0 ; i<FUNCTION_SIZE ; i++) {
        data_x[i] = i * scale - CPL_MATH_PI;
        data_y[i] = cos(data_x[i]);
    }
  
    /* Test cpl_bivector_get_y() */
    cpl_test( vec1 = cpl_bivector_get_y(sinus) );


    /* Test cpl_bivector_get_size() */
    cpl_test( cpl_bivector_get_size(sinus) == FUNCTION_SIZE );
   
    /* Test cpl_bivector_dump() */
    cpl_test( sprintf(filename, "cpl_bivector_dump.txt") );
    cpl_test( f_out = fopen(filename, "w") );

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

    cpl_bivector_dump(sinus, f_out);
    cpl_test_zero( fclose(f_out) );

    /* Test cpl_bivector_read() */

    tmp_fun = cpl_bivector_read(NULL);
    cpl_test_error( CPL_ERROR_NULL_INPUT );
    cpl_test_null( tmp_fun );

    tmp_fun = cpl_bivector_read("/dev/null");
    cpl_test_error(CPL_ERROR_BAD_FILE_FORMAT);
    cpl_test_null( tmp_fun );

    cpl_test( tmp_fun = cpl_bivector_duplicate(sinus) );
    cpl_bivector_delete(sinus);

    cpl_test( sinus = cpl_bivector_read("cpl_bivector_dump.txt") );

    remove("cpl_bivector_dump.txt");

    /* Measure the error incurred through the dump/read */
    vec1 = cpl_bivector_get_x(tmp_fun);
    vec2 = cpl_bivector_get_y(tmp_fun);

    cpl_test_zero( cpl_vector_subtract(vec1, cpl_bivector_get_x(sinus)) );
    cpl_test_zero( cpl_vector_subtract(vec2, cpl_bivector_get_y(sinus)) );

    /* The 2-norm divided by sqrt(size) */
    error = sqrt(cpl_tools_get_variance_double(cpl_vector_get_data(vec1),
                                          FUNCTION_SIZE, 0));
    cpl_test_leq( error,  23*FLT_EPSILON);

    error = sqrt(cpl_tools_get_variance_double(cpl_vector_get_data(vec2),
                                          FUNCTION_SIZE, 0));
    cpl_test_leq( error,  2.4*FLT_EPSILON);

    cpl_bivector_delete(tmp_fun);
    
    /* Test cpl_bivector_interpolate_linear */
    tmp_fun = cpl_bivector_new(FUNCTION_SIZE);
    vec1 = cpl_bivector_get_x(tmp_fun);
    vec2 = cpl_bivector_get_x(sinus);
    cpl_vector_copy(vec1, vec2);
    cpl_vector_add_scalar(vec1, 0.5);
    cpl_test_eq( cpl_bivector_interpolate_linear(tmp_fun, sinus),
                 CPL_ERROR_ILLEGAL_INPUT);
    cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
    cpl_bivector_delete(tmp_fun);

    /* Free */
    cpl_bivector_delete(sinus);
    cpl_bivector_delete(cosinus);

    sinus = cpl_bivector_new(FUNCTION_SIZE+1);

    /* Create a function with sine values, with abscissa offset by 0.5 */
    data_x = cpl_bivector_get_x_data(sinus);
    data_y = cpl_bivector_get_y_data(sinus);
    for (i=0 ; i<FUNCTION_SIZE+1 ; i++) {
        data_x[i] = (i-0.5) * scale;
        data_y[i] = sin(data_x[i]);
    }
    
    tmp_fun = cpl_bivector_new(FUNCTION_SIZE);
    vec1 = cpl_bivector_get_x(tmp_fun);
    data_x = cpl_vector_get_data(vec1);
    for (i=0 ; i<FUNCTION_SIZE ; i++) data_x[i] = i * scale;

    cpl_test_zero( cpl_bivector_interpolate_linear(tmp_fun, sinus) );

    data_y = cpl_bivector_get_y_data(tmp_fun);

    /* Check interpolation error at 0, pi/2, pi, 3*pi/2 */
    i = 0;
    cpl_test_leq( fabs(data_y[i]),  DBL_EPSILON);

    i+= FUNCTION_SIZE/2;
    cpl_test_leq( fabs(data_y[i]), 10 * DBL_EPSILON);

    i = FUNCTION_SIZE/4;
    cpl_test_leq( fabs(1 - data_y[i]) * FUNCTION_SIZE, 1 );
    i+= FUNCTION_SIZE/2;
    cpl_test_leq( fabs(1 + data_y[i]) * FUNCTION_SIZE, 1 );

    cpl_bivector_delete(sinus);
    cpl_bivector_delete(tmp_fun);

    /* End of actual test code */
    return cpl_test_end(0);
}