/* $Id: cpl_vector.c,v 1.133 2007/12/06 14:30:47 yjung 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: yjung $ * $Date: 2007/12/06 14:30:47 $ * $Revision: 1.133 $ * $Name: $ */ #ifdef HAVE_CONFIG_H #include #endif /*----------------------------------------------------------------------------- Includes -----------------------------------------------------------------------------*/ #include #include #include #include /* Needed by memcpy() */ #include #include #include "cpl_memory.h" #include "cpl_matrix.h" #include "cpl_vector.h" #include "cpl_vector_fit_impl.h" #include "cpl_tools.h" #include "cpl_math_const.h" /*----------------------------------------------------------------------------*/ /** * @defgroup cpl_vector Vector * * This module provides functions to handle @em cpl_vector. * * A @em cpl_vector is an object containing a list of values (as doubles) * and the (always positive) number of values. The functionalities provided * here are simple ones like sorting, statistics, or simple operations. * The @em cpl_bivector object is composed of two of these vectors. * * @par Synopsis: * @code * #include "cpl_vector.h" * @endcode */ /*----------------------------------------------------------------------------*/ /**@{*/ /*----------------------------------------------------------------------------- Type definition -----------------------------------------------------------------------------*/ struct _cpl_vector_ { int size; double * data; }; /*----------------------------------------------------------------------------- Private function prototypes -----------------------------------------------------------------------------*/ static cpl_error_code cpl_vector_append(const cpl_vector *, const char *, cpl_type_bpp, cpl_propertylist *) ; static cpl_vector * cpl_vector_gen_lowpass_kernel(cpl_lowpass, int); static double cpl_tools_sinc(double); static void cpl_vector_fill_tanh_kernel(cpl_vector *, double); static void cpl_vector_fill_alpha_kernel(cpl_vector *, double, double); static void cpl_vector_reverse_tanh_kernel(double *, int); static int cpl_fit_gaussian_1d_compare(const void *left, const void *right); static int gauss(const double x[], const double a[], double *result); static int gauss_derivative(const double x[], const double a[], double result[]); /*----------------------------------------------------------------------------- Private types -----------------------------------------------------------------------------*/ typedef struct { double x; double y; } cpl_vector_fit_gaussian_input; /*----------------------------------------------------------------------------- Function codes -----------------------------------------------------------------------------*/ /*----------------------------------------------------------------------------*/ /** @brief Create a new cpl_vector @param n Number of element of the cpl_vector @return 1 newly allocated cpl_vector or NULL in case of an error The returned object must be deallocated using cpl_vector_delete(). There is no default values assigned to the created object, they are undefined intil they are set. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_ILLEGAL_INPUT if n is negative or zero */ /*----------------------------------------------------------------------------*/ cpl_vector * cpl_vector_new(int n) { cpl_vector * v; /* Test input */ cpl_ensure(n > 0, CPL_ERROR_ILLEGAL_INPUT, NULL); /* Allocate memory */ v = cpl_malloc(sizeof(cpl_vector)); v->size = n; v->data = cpl_malloc(n * sizeof(double)); /* Return */ return v; } /*----------------------------------------------------------------------------*/ /** @brief Create a cpl_vector from existing data @param n Number of elements in the vector @param data Pointer to array of n doubles @return 1 newly allocated cpl_vector or NULL in case of an error The returned object must be deallocated using cpl_vector_unwrap(). Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if an input pointer is NULL - CPL_ERROR_ILLEGAL_INPUT if n is negative or zero */ /*----------------------------------------------------------------------------*/ cpl_vector * cpl_vector_wrap(int n, double * data) { cpl_vector * v; /* Test input */ cpl_ensure(data, CPL_ERROR_NULL_INPUT, NULL); cpl_ensure(n > 0, CPL_ERROR_ILLEGAL_INPUT, NULL); /* Allocate memory */ v = cpl_malloc(sizeof(cpl_vector)); v->size = n; v->data = data; /* Return */ return v; } /*----------------------------------------------------------------------------*/ /** @brief Delete a cpl_vector @param v cpl_vector to delete @return void */ /*----------------------------------------------------------------------------*/ void cpl_vector_delete(cpl_vector * v) { if (v == NULL) return; if (v->data != NULL) cpl_free(v->data); cpl_free(v); return; } /*----------------------------------------------------------------------------*/ /** @brief Delete a cpl_vector except the data array @param v cpl_vector to delete @return A pointer to the data array or NULL if the input is NULL. @note The data array must subsequently be deallocated. Failure to do so will result in a memory leak. */ /*----------------------------------------------------------------------------*/ void * cpl_vector_unwrap(cpl_vector * v) { void * data; if (v == NULL) return NULL; data = (void *) v->data; cpl_free(v); return data; } /*----------------------------------------------------------------------------*/ /** @brief Read a list of values from an ASCII file and create a cpl_vector @param filename Name of the input ASCII file @return 1 newly allocated cpl_vector or NULL in case of an error @see cpl_vector_dump Parse an input ASCII file values and create a cpl_vector from it Lines beginning with a hash are ignored, blank lines also. In valid lines the value is preceeded by an integer, which is ignored. The returned object must be deallocated using cpl_vector_delete() In addition to normal files, FIFO (see man mknod) are also supported. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if an input pointer is NULL - CPL_ERROR_FILE_IO if the file cannot be read - CPL_ERROR_BAD_FILE_FORMAT if the file contains no valid lines */ /*----------------------------------------------------------------------------*/ cpl_vector * cpl_vector_read(const char * filename) { FILE * in; cpl_vector * v; int np = 0; int size = 1000; /* Default size */ char line[1024]; double x; cpl_ensure(filename != NULL, CPL_ERROR_NULL_INPUT, NULL); /* Open the file */ in = fopen(filename, "r"); cpl_ensure(in != NULL, CPL_ERROR_FILE_IO, NULL); /* Create and fill the vector */ v = cpl_vector_new(size); while (fgets(line, 1024, in) != NULL) { if (line[0] != '#' && sscanf(line, "%*d %lg", &x) == 1) { /* Found new element - increase vector size if necessary, - insert element at end and - increment size counter */ if (np == size) cpl_vector_set_size(v, size *= 2); cpl_vector_set(v, np++, x); } } /* Check that the loop ended due to eof and not an error */ if (ferror(in)) { fclose(in); cpl_vector_delete(v); cpl_ensure(0, CPL_ERROR_FILE_IO, NULL); } fclose(in); /* Check that the file was not empty and set the size to its true value */ if (np == 0 || cpl_vector_set_size(v, np)) { cpl_vector_delete(v); cpl_ensure(0, CPL_ERROR_BAD_FILE_FORMAT, NULL); } return v; } /*----------------------------------------------------------------------------*/ /** @brief Dump a cpl_vector as ASCII to a stream @param v Input cpl_vector to dump @param stream Output stream, accepts @c stdout or @c stderr @return void Each element is preceded by its index number (starting with 1!) and written on a single line. Comment lines start with the hash character. stream may be NULL in which case @c stdout is used. @note In principle a cpl_vector can be saved using cpl_vector_dump() and re-read using cpl_vector_read(). This will however introduce significant precision loss due to the limited accuracy of the ASCII representation. */ /*----------------------------------------------------------------------------*/ void cpl_vector_dump( const cpl_vector * v, FILE * stream) { int i; if (stream == NULL) stream = stdout; fprintf(stream, "#----- vector -----\n"); if (v == NULL) { fprintf(stream, "#NULL vector\n"); return; } fprintf(stream, "#Index\t\tX\n"); for (i=0; isize; i++) fprintf(stream, "%d\t\t%g\n", i+1, v->data[i]); fprintf(stream, "#------------------\n"); return; } /*----------------------------------------------------------------------------*/ /** @brief Load a list of values from a FITS file @param filename Name of the input file @param xtnum Extension number in the file. @return 1 newly allocated cpl_vector or NULL in case of an error @see cpl_vector_save This function loads a vector from a FITS file (NAXIS=1), using cfitsio. The returned image has to be deallocated with cpl_vector_delete(). 'xtnum' specifies from which extension the image should be loaded. This could be 0 for the main data section (files without extension), or any number between 1 and N, where N is the number of extensions present in the file. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if an input pointer is NULL - CPL_ERROR_ILLEGAL_INPUT if the extension is not valid - CPL_ERROR_FILE_IO if the file cannot be read */ /*----------------------------------------------------------------------------*/ cpl_vector * cpl_vector_load( const char * filename, int xtnum) { fitsfile * fptr ; int fio_status=0 ; cpl_vector * vec ; char * extname ; long nx ; long fpixel[1] ; int naxis ; /* Test entries */ cpl_ensure(filename != NULL, CPL_ERROR_NULL_INPUT, NULL); cpl_ensure(xtnum>=0, CPL_ERROR_ILLEGAL_INPUT, NULL) ; /* Open the file extension */ if (xtnum == 0) extname = cpl_sprintf(filename) ; else extname = cpl_sprintf("%s[%d]", filename, xtnum) ; fits_open_file(&fptr, extname, READONLY, &fio_status) ; cpl_free(extname) ; cpl_ensure(fio_status==0, CPL_ERROR_FILE_IO, NULL) ; /* Get the file dimension */ fits_get_img_dim(fptr, &naxis, &fio_status) ; if (naxis != 1) { fits_close_file(fptr, &fio_status) ; cpl_ensure(0, CPL_ERROR_ILLEGAL_INPUT, NULL) ; } fits_get_img_size(fptr, 1, &nx, &fio_status) ; /* Create the vector */ vec = cpl_vector_new((int)nx) ; /* Read */ fpixel[0] = 1 ; fits_read_pix(fptr,TDOUBLE,fpixel, nx, NULL, vec->data, NULL, &fio_status) ; /* Close */ fits_close_file(fptr, &fio_status) ; if (fio_status != 0) { cpl_free(vec) ; cpl_ensure(0, CPL_ERROR_FILE_IO, NULL) ; } /* Return */ return vec ; } /*----------------------------------------------------------------------------*/ /** @brief Save a vector to a FITS file @param to_save Vector to write to disk or NULL @param filename Name of the file to write @param bpp Bits per pixel @param pl Property list for the output header or NULL @param mode The desired output options (combined with bitwise or) @return CPL_ERROR_NONE or the relevant #_cpl_error_code_ on error This function saves a vector to a FITS file (NAXIS=1), using cfitsio. If a property list is provided, it is written to the named file before the pixels are written. If the image is not provided, the created file will only contain the primary header. This can be useful to create multiple extension files. The requested pixel depth (bpp) follows the FITS convention. Possible values are CPL_BPP_8_UNSIGNED (8), CPL_BPP_16_SIGNED (16), CPL_BPP_32_SIGNED (32), CPL_BPP_IEEE_FLOAT (-32), CPL_BPP_IEEE_DOUBLE (-64) Supported output modes are CPL_IO_DEFAULT (create a new file) and CPL_IO_EXTEND (append to an existing file) If you are in append mode, make sure that the file has writing permissions. You may have problems if you create a file in your application and append something to it with the umask set to 222. In this case, the file created by your application would not be writable, and the append would fail. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if (one of) the input pointer(s) is NULL - CPL_ERROR_ILLEGAL_INPUT if the bpp value or the mode is not supported - CPL_ERROR_FILE_NOT_CREATED if the output file cannot be created - CPL_ERROR_FILE_IO if the data cannot be written in the file */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_vector_save( const cpl_vector * to_save, const char * filename, cpl_type_bpp bpp, const cpl_propertylist * pl, unsigned mode) { fitsfile * fptr ; char sval[1024] ; int fio_status=0 ; long naxes[1] ; int ival ; double dval ; /* Test entries */ cpl_ensure_code(filename, CPL_ERROR_NULL_INPUT) ; cpl_ensure_code(bpp==CPL_BPP_8_UNSIGNED || bpp==CPL_BPP_16_SIGNED || bpp==CPL_BPP_32_SIGNED || bpp==CPL_BPP_IEEE_FLOAT || bpp==CPL_BPP_IEEE_DOUBLE, CPL_ERROR_ILLEGAL_INPUT); cpl_ensure_code(mode <= CPL_IO_MAX, CPL_ERROR_ILLEGAL_INPUT) ; /* In append mode, call the dedicated function */ if (mode & CPL_IO_EXTEND) return cpl_vector_append(to_save, filename, bpp, (cpl_propertylist*)pl); /* Create the file */ sprintf(sval, "!%s", filename) ; fits_create_file(&fptr, sval, &fio_status); cpl_ensure_code(fio_status==0, CPL_ERROR_FILE_IO) ; /* Create the vector in the primary unit */ if (to_save != NULL) { naxes[0] = to_save->size ; fits_create_img(fptr, bpp, 1, naxes, &fio_status); /* Write the data */ fits_write_img(fptr, TDOUBLE, 1, naxes[0], to_save->data, &fio_status); } else { /* Create minimal header */ fits_create_hdu(fptr, &fio_status); ival = 1 ; /* FIXME: Cast string literal to avoid compiler warnings - and pray that CFITSIO does not try to modify them :-(((( */ fits_write_key(fptr, TLOGICAL, (char*)"SIMPLE", &ival, (char*)"Fits format", &fio_status); fits_write_key(fptr, TINT, (char*)"BITPIX", &bpp, (char*)"Bits per pixel", &fio_status); ival = 0 ; fits_write_key(fptr, TINT, (char*)"NAXIS", &ival, (char*)"Empty unit", &fio_status); ival = 1 ; fits_write_key(fptr, TLOGICAL, (char*)"EXTEND", &ival, (char*)"There may be FITS extensions", &fio_status); } /* Add BSCALE BZERO DATE */ dval = 1.0 ; fits_write_key(fptr, TDOUBLE, (char*)"BSCALE", &dval, (char*)"Scale factor", &fio_status); dval = 0.0 ; fits_write_key(fptr, TDOUBLE, (char*)"BZERO", &dval, (char*)"Offset factor", &fio_status); fits_write_date(fptr, &fio_status); /* Check */ if (fio_status != 0) { fits_close_file(fptr, &fio_status) ; cpl_ensure_code(0, CPL_ERROR_FILE_IO) ; } /* Add the property list */ if (cpl_fits_add_properties(fptr, pl, CPL_FITS_BADKEYS_PRIM "|" CPL_FITS_COMPRKEYS)!=CPL_ERROR_NONE) { fits_close_file(fptr, &fio_status) ; cpl_ensure_code(0, CPL_ERROR_ILLEGAL_INPUT) ; } /* Close and write on disk */ fits_close_file(fptr, &fio_status) ; return CPL_ERROR_NONE ; } /*----------------------------------------------------------------------------*/ /** @brief This function duplicates an existing vector and allocates memory @param v the input cpl_vector @return a newly allocated cpl_vector or NULL in case of an error The returned object must be deallocated using cpl_vector_delete() Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if an input pointer is NULL */ /*----------------------------------------------------------------------------*/ cpl_vector * cpl_vector_duplicate(const cpl_vector * v) { cpl_vector * out; /* Test inputs */ cpl_ensure(v != NULL, CPL_ERROR_NULL_INPUT, NULL); /* Create a new cpl_vector */ out = cpl_vector_new(v->size); /* Copy the contents */ cpl_vector_copy(out, v); return out; } /*----------------------------------------------------------------------------*/ /** @brief This function copies contents of a vector into another vector @param destination destination cpl_vector @param source source cpl_vector @return CPL_ERROR_NONE or the relevant #_cpl_error_code_ on error @see cpl_vector_set_size() if source and destination have different sizes. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if an input pointer is NULL */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_vector_copy( cpl_vector * destination, const cpl_vector * source) { cpl_ensure_code(source != NULL, CPL_ERROR_NULL_INPUT); cpl_ensure_code(!cpl_vector_set_size(destination, source->size), cpl_error_get_code()); memcpy(destination->data, source->data, source->size * sizeof(double)); return CPL_ERROR_NONE; } /*----------------------------------------------------------------------------*/ /** @brief Get the size of the vector @param in the input vector @return The size or -1 in case of an error Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if an input pointer is NULL */ /*----------------------------------------------------------------------------*/ int cpl_vector_get_size(const cpl_vector * in) { /* Test entries */ cpl_ensure(in!=NULL, CPL_ERROR_NULL_INPUT, -1); return in->size; } /*----------------------------------------------------------------------------*/ /** @brief Resize the vector @param in The vector to be resized @param newsize The new (positive) number of elements in the vector @return CPL_ERROR_NONE or the relevant #_cpl_error_code_ on error @note On succesful return the value of the elements of the vector is unchanged to the minimum of the old and new sizes; any newly allocated elements are undefined. The pointer to the vector data buffer may change, therefore pointers previously retrieved by calling @c cpl_vector_get_data() should be discarded. If the vector was created with cpl_vector_wrap() the argument pointer to that call should be discarded as well. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if an input pointer is NULL - CPL_ERROR_ILLEGAL_INPUT if newsize is negative or zero */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_vector_set_size(cpl_vector * in, int newsize) { /* Test entries */ cpl_ensure_code(in != NULL, CPL_ERROR_NULL_INPUT); cpl_ensure_code(newsize > 0, CPL_ERROR_ILLEGAL_INPUT); if (in->size == newsize) return CPL_ERROR_NONE; assert( in->data ); in->data = cpl_realloc(in->data, newsize * sizeof(double)); in->size = newsize; return CPL_ERROR_NONE; } /*----------------------------------------------------------------------------*/ /** @brief Get a pointer to the data part of the vector @param in the input vector @return Pointer to the data or NULL in case of an error The returned pointer refers to already allocated data. @note Use at your own risk: direct manipulation of vector data rules out any check performed by the vector object interface, and may introduce inconsistencies between the information maintained internally, and the actual vector data and structure. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if an input pointer is NULL */ /*----------------------------------------------------------------------------*/ double * cpl_vector_get_data(cpl_vector * in) { /* Test entries */ cpl_ensure(in!=NULL, CPL_ERROR_NULL_INPUT, NULL); assert( in->data ); return in->data; } /*----------------------------------------------------------------------------*/ /** @brief Get a pointer to the data part of the vector @param in the input vector @return Pointer to the data or NULL in case of an error @see cpl_vector_get_data */ /*----------------------------------------------------------------------------*/ const double * cpl_vector_get_data_const(const cpl_vector * in) { return cpl_vector_get_data((cpl_vector *)in); } /*----------------------------------------------------------------------------*/ /** @brief Get an element of the vector @param in the input vector @param index the index of the element (0 to nelem-1) @return The element value In case of error, the #_cpl_error_code_ code is set, and the returned double is undefined. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if an input pointer is NULL - CPL_ERROR_ILLEGAL_INPUT if index is negative - CPL_ERROR_ACCESS_OUT_OF_RANGE if index is out of the vector bounds */ /*----------------------------------------------------------------------------*/ double cpl_vector_get(const cpl_vector * in, int index) { /* Test entries */ cpl_ensure(in, CPL_ERROR_NULL_INPUT, -1.0); cpl_ensure(index >= 0, CPL_ERROR_ILLEGAL_INPUT, -2.0); cpl_ensure(index < in->size, CPL_ERROR_ACCESS_OUT_OF_RANGE, -3.0); return in->data[index]; } /*----------------------------------------------------------------------------*/ /** @brief Set an element of the vector @param in the input vector @param index the index of the element (0 to nelem-1) @param value the value to set in the vector @return CPL_ERROR_NONE or the relevant #_cpl_error_code_ on error Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if an input pointer is NULL - CPL_ERROR_ILLEGAL_INPUT if index is negative - CPL_ERROR_ACCESS_OUT_OF_RANGE if index is out of the vector bounds */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_vector_set( cpl_vector * in, int index, double value) { /* Test entries */ cpl_ensure_code(in != NULL, CPL_ERROR_NULL_INPUT); cpl_ensure_code(index >= 0, CPL_ERROR_ILLEGAL_INPUT); cpl_ensure_code(index < in->size, CPL_ERROR_ACCESS_OUT_OF_RANGE); /* Assign the value */ in->data[index] = value; return CPL_ERROR_NONE; } /*----------------------------------------------------------------------------*/ /** @brief Add a cpl_vector to another @param v1 First cpl_vector (modified) @param v2 Second cpl_vector @return CPL_ERROR_NONE or the relevant #_cpl_error_code_ on error The second vector is added to the first one. The input first vector is modified. The input vectors must have the same size. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if an input pointer is NULL - CPL_ERROR_INCOMPATIBLE_INPUT if v1 and v2 have different sizes */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_vector_add( cpl_vector * v1, const cpl_vector * v2) { int i; /* Test inputs */ cpl_ensure_code(v1 != NULL, CPL_ERROR_NULL_INPUT); cpl_ensure_code(v2 != NULL, CPL_ERROR_NULL_INPUT); cpl_ensure_code(v1->size == v2->size, CPL_ERROR_INCOMPATIBLE_INPUT); /* Add v2 to v1 */ for (i=0; i < v1->size; i++) v1->data[i] += v2->data[i]; cpl_tools_add_flops(v1->size); return CPL_ERROR_NONE; } /*----------------------------------------------------------------------------*/ /** @brief Subtract a cpl_vector from another @param v1 First cpl_vector (modified) @param v2 Second cpl_vector @return CPL_ERROR_NONE or the relevant #_cpl_error_code_ on error @see cpl_vector_add() */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_vector_subtract( cpl_vector * v1, const cpl_vector * v2) { int i; /* Test inputs */ cpl_ensure_code(v1 != NULL, CPL_ERROR_NULL_INPUT); cpl_ensure_code(v2 != NULL, CPL_ERROR_NULL_INPUT); cpl_ensure_code(v1->size == v2->size, CPL_ERROR_INCOMPATIBLE_INPUT); /* subtract v2 from v1 */ for (i=0; i < v1->size; i++) v1->data[i] -= v2->data[i]; cpl_tools_add_flops(v1->size); return CPL_ERROR_NONE; } /*----------------------------------------------------------------------------*/ /** @brief Multiply two vectors component-wise @param v1 First cpl_vector @param v2 Second cpl_vector @return CPL_ERROR_NONE or the relevant #_cpl_error_code_ on error @see cpl_vector_add() */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_vector_multiply( cpl_vector * v1, const cpl_vector * v2) { int i; /* Test inputs */ cpl_ensure_code(v1 != NULL, CPL_ERROR_NULL_INPUT); cpl_ensure_code(v2 != NULL, CPL_ERROR_NULL_INPUT); cpl_ensure_code(v1->size == v2->size, CPL_ERROR_INCOMPATIBLE_INPUT); /* Multiply the cpl_vector fields */ for (i=0; i < v1->size; i++) v1->data[i] *= v2->data[i]; cpl_tools_add_flops(v1->size); return CPL_ERROR_NONE; } /*----------------------------------------------------------------------------*/ /** @brief Divide two vectors element-wise @param v1 First cpl_vector @param v2 Second cpl_vector @return CPL_ERROR_NONE or the relevant #_cpl_error_code_ on error @see cpl_vector_add() If an element in v2 is zero v1 is not modified and CPL_ERROR_DIVISION_BY_ZERO is returned. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if an input pointer is NULL - CPL_ERROR_INCOMPATIBLE_INPUT if v1 and v2 have different sizes - CPL_ERROR_DIVISION_BY_ZERO if a division by 0 would occur */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_vector_divide( cpl_vector * v1, const cpl_vector * v2) { int i; /* Test inputs */ cpl_ensure_code(v1 != NULL, CPL_ERROR_NULL_INPUT); cpl_ensure_code(v2 != NULL, CPL_ERROR_NULL_INPUT); cpl_ensure_code(v1->size == v2->size, CPL_ERROR_INCOMPATIBLE_INPUT); for (i=0; i < v2->size; i++) cpl_ensure_code(v2->data[i] != 0, CPL_ERROR_DIVISION_BY_ZERO); /* Divide the cpl_vector fields */ for (i=0; i < v1->size; i++) v1->data[i] /= v2->data[i]; cpl_tools_add_flops(v1->size); return CPL_ERROR_NONE; } /*----------------------------------------------------------------------------*/ /** @brief Compute the vector dot product. @param v1 One vector @param v2 Another vector of the same size @return The (non-negative) product or negative on error. The same vector may be passed twice, in which case the square of its 2-norm is computed. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if an input pointer is NULL - CPL_ERROR_INCOMPATIBLE_INPUT if v1 and v2 have different sizes */ /*----------------------------------------------------------------------------*/ double cpl_vector_product(const cpl_vector * v1, const cpl_vector * v2) { const double * x = cpl_vector_get_data_const(v1); const double * y = cpl_vector_get_data_const(v2); const int n = cpl_vector_get_size(v1); double sum = 0; int i; cpl_ensure(x, CPL_ERROR_NULL_INPUT, -1); cpl_ensure(y, CPL_ERROR_NULL_INPUT, -2); cpl_ensure(cpl_vector_get_size(v2) == n, CPL_ERROR_INCOMPATIBLE_INPUT, -3); for (i = 0; i < n; i++) sum += x[i] * y[i]; cpl_tools_add_flops(2*n); return sum; } /*----------------------------------------------------------------------------*/ /** @brief Sort a cpl_vector by increasing/decreasing data @param v cpl_vector to sort @param c Sorting criterion to use @return CPL_ERROR_NONE or the relevant #_cpl_error_code_ on error The input cpl_vector is modified to sort out all values following the sorting criterion. Possible sorting criteria are: - +1 to sort by increasing data - -1 to sort by decreasing data Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if an input pointer is NULL - CPL_ERROR_ILLEGAL_INPUT if c is neither 1 nor -1 */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_vector_sort( cpl_vector * v, int c) { int i; cpl_ensure_code( v != NULL, CPL_ERROR_NULL_INPUT); if (v->size == 1) return CPL_ERROR_NONE; /* Sort by increasing data */ cpl_tools_sort_double(v->data, v->size); switch (c) { case 1 : /* Nothing to do */ break; case -1 : /* Invert the order */ for (i=0; i < v->size/2; i++) { const double tmp = v->data[i]; v->data[i] = v->data[v->size - i - 1]; v->data[v->size - i - 1] = tmp; } break; default: cpl_ensure_code(0, CPL_ERROR_ILLEGAL_INPUT); break; } return CPL_ERROR_NONE; } /*----------------------------------------------------------------------------*/ /** @brief Fill a cpl_vector @param v cpl_vector to be filled with the value val @param val Value used to fill the cpl_vector @return CPL_ERROR_NONE or the relevant #_cpl_error_code_ on error Input vector is modified Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if an input pointer is NULL */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_vector_fill(cpl_vector * v, double val) { int i; cpl_ensure_code(v, CPL_ERROR_NULL_INPUT); for (i=0; isize; i++) v->data[i] = val; return CPL_ERROR_NONE; } /*----------------------------------------------------------------------------*/ /** @brief Compute the sqrt of a cpl_vector @param v cpl_vector @return CPL_ERROR_NONE or the relevant #_cpl_error_code_ on error The sqrt of the data is computed. The input cpl_vector is modified If an element in v is negative v is not modified and CPL_ERROR_ILLEGAL_INPUT is returned. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if an input pointer is NULL - CPL_ERROR_ILLEGAL_INPUT if one of the vector values is negative */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_vector_sqrt(cpl_vector * v) { int i; /* Test inputs */ cpl_ensure_code(v!=NULL, CPL_ERROR_NULL_INPUT); for (i=0; isize; i++) cpl_ensure_code(v->data[i] >= 0, CPL_ERROR_ILLEGAL_INPUT); /* Compute the sqrt */ for (i=0; i < v->size; i++) v->data[i] = sqrt(v->data[i]); cpl_tools_add_flops(v->size); return CPL_ERROR_NONE; } /*----------------------------------------------------------------------------*/ /** @brief In a sorted vector find the element closest to the given value @param sorted Sorted cpl_vector @param key Value to find @return The index that minimizes fabs(sorted[index] - key) or -1 on error Bisection is used to find the element. If two (neighboring) elements with different values both minimize fabs(sorted[index] - key) the index of the larger element is returned. If the vector contains identical elements that minimize fabs(sorted[index] - key) then it is undefined which element has its index returned. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if an input pointer is NULL */ /*----------------------------------------------------------------------------*/ int cpl_vector_find( const cpl_vector * sorted, double key) { int low, high; /* Test inputs */ cpl_ensure(sorted != NULL, CPL_ERROR_NULL_INPUT, -1); /* Initialize */ low = 0; high = sorted->size - 1; if (key <= sorted->data[ low]) return low; if (key >= sorted->data[high]) return high; /* key is in the range of the vector values */ while (high - low > 1) { const int middle = (high + low) >> 1; if (sorted->data[middle] > key) high = middle; else if (sorted->data[middle] < key) low = middle; else low = high = middle; } /* assert( low == high || low == high - 1); */ return key - sorted->data[low] < sorted->data[high] - key ? low : high; } /*----------------------------------------------------------------------------*/ /** @brief Extract a sub_vector from a vector @param v Input cpl_vector @param istart Start index (from 0 to number of elements - 1) @param istop Stop index (from 0 to number of elements - 1) @param istep Extract every step element @return A newly allocated cpl_vector or NULL in case of an error The returned object must be deallocated using cpl_vector_delete() FIXME: Currently istop must be greater than istart. FIXME: Currently istep must equal 1. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if an input pointer is NULL - CPL_ERROR_ILLEGAL_INPUT if istart, istop, istep are not as requested */ /*----------------------------------------------------------------------------*/ cpl_vector * cpl_vector_extract( const cpl_vector * v, int istart, int istop, int istep) { cpl_vector * out; /* Test inputs */ cpl_ensure(v, CPL_ERROR_NULL_INPUT, NULL); cpl_ensure(istart >= 0, CPL_ERROR_ACCESS_OUT_OF_RANGE, NULL); cpl_ensure(istart <= istop, CPL_ERROR_ILLEGAL_INPUT, NULL); cpl_ensure(istep == 1, CPL_ERROR_ILLEGAL_INPUT, NULL); cpl_ensure(istop < v->size, CPL_ERROR_ACCESS_OUT_OF_RANGE, NULL); /* Allocate and fill the output vector */ out = cpl_vector_new(istop - istart + 1); assert( out ); memcpy(out->data, &(v->data[istart]), out->size * sizeof(double)); return out; } /*----------------------------------------------------------------------------*/ /** @brief Get the minimum of the cpl_vector @param v const cpl_vector @return The minimum value of the vector or undefined on error In case of error, the #_cpl_error_code_ code is set, and the returned double is undefined. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if an input pointer is NULL */ /*----------------------------------------------------------------------------*/ double cpl_vector_get_min(const cpl_vector * v) { double min; int i; /* Test entries */ cpl_ensure(v != NULL, CPL_ERROR_NULL_INPUT, -1.0); min = v->data[0]; for (i=1; isize; i++) if (v->data[i] < min) min = v->data[i]; return min; } /*----------------------------------------------------------------------------*/ /** @brief Get the maximum of the cpl_vector @param v const cpl_vector @return the maximum value of the vector or undefined on error @see cpl_vector_get_min() */ /*----------------------------------------------------------------------------*/ double cpl_vector_get_max(const cpl_vector * v) { double max; int i; /* Test entries */ cpl_ensure(v != NULL, CPL_ERROR_NULL_INPUT, -1.0); max = v->data[0]; for (i=1; i < v->size; i++) if (v->data[i] > max) max = v->data[i]; return max; } /*----------------------------------------------------------------------------*/ /** @brief Compute the mean value of vector elements. @param v Input const cpl_vector @return Mean value of vector elements or undefined on error. @see cpl_vector_get_min() */ /*----------------------------------------------------------------------------*/ double cpl_vector_get_mean(const cpl_vector * v) { double mean = 0; int i; cpl_ensure(v != NULL, CPL_ERROR_NULL_INPUT, -1.0); /* Compute the arithmetic mean of a dataset using the recurrence relation mean_(n) = mean(n-1) + (v[n] - mean(n-1))/(n+1) - this has a measurable impact on the output of cpl_polynomial_fit_{1,2}d_create() */ for (i = 0; i < v->size; i++) mean += (v->data[i] - mean) / (i + 1); cpl_tools_add_flops(3 * v->size); return mean; } /*----------------------------------------------------------------------------*/ /** @brief Compute the median of the elements of a vector. @param v Input cpl_vector @return Median value of the vector elements or undefined on error. @see cpl_vector_get_median_const() @note For efficiency reasons, this function modifies the order of the elements of the input vector. */ /*----------------------------------------------------------------------------*/ double cpl_vector_get_median(cpl_vector * v) { cpl_ensure(v != NULL, CPL_ERROR_NULL_INPUT, 0.0); return cpl_tools_get_median_double(v->data, v->size); } /*----------------------------------------------------------------------------*/ /** @brief Compute the median of the elements of a vector. @param v Input const cpl_vector @return Median value of the vector elements or undefined on error. @see cpl_vector_get_min() For a finite population or sample, the median is the middle value of an odd number of values (arranged in ascending order) or any value between the two middle values of an even number of values. The computed median should actually be a value of the input array, so, in case of an even number of values in the input array, the median would be the lowest of the two central values. */ /*----------------------------------------------------------------------------*/ double cpl_vector_get_median_const(const cpl_vector * v) { double * darr; double med; cpl_ensure(v != NULL, CPL_ERROR_NULL_INPUT, 0.0); /* Create a separate data buffer because */ /* cpl_tools_get_median_double() modifies its input */ darr = cpl_malloc(v->size * sizeof(double)); memcpy(darr, v->data, v->size * sizeof(double)); med = cpl_tools_get_median_double(darr, v->size); cpl_free(darr); return med; } /*----------------------------------------------------------------------------*/ /** @brief Compute the bias-corrected standard deviation of a vectors elements @param v Input const cpl_vector @return standard deviation of the elements or a negative number on error. @see cpl_vector_get_min() S(n-1) = sqrt((1/n-1) sum((xi-mean)^2) (i=1 -> n) The length of v must be at least 2. */ /*----------------------------------------------------------------------------*/ double cpl_vector_get_stdev(const cpl_vector * v) { const double mean = cpl_vector_get_mean(v); double variance; cpl_ensure(v != NULL, CPL_ERROR_NULL_INPUT,-1); cpl_ensure(v->size > 1, CPL_ERROR_ILLEGAL_INPUT,-2); variance = cpl_tools_get_variance_double(v->data, v->size, mean); /* Compute the bias-corrected standard deviation. - With the recurrence relation rounding can likely not cause the variance to become negative, but check just to be safe */ return variance > 0 ? sqrt(variance * v->size / (double) (v->size-1)) : 0; } /*----------------------------------------------------------------------------*/ /** @brief Cross-correlation of two vectors @param vxc Odd-sized vector with the computed cross-correlations @param v1 1st vector to correlate @param v2 2nd vector to correlate @return Index of maximum cross-correlation, or negative on error vxc must have an odd number of elements, 2*half_search+1, where half_search is the half-size of the search domain. The length of v2 may not exceed that of v1. If the difference in length between v1 and v2 is less than half_search then this difference must be even (if the difference is odd resampling of v2 may be useful). The cross-correlation is computed with shifts ranging from -half_search to half_search. On succesful return element i (starting with 0) of vxc contains the cross- correlation at offset i-half_search. On error vxc is unmodified. The cross-correlation is in fact the dot-product of two unit-vectors and ranges therefore from -1 to 1. The cross-correlation is, in absence of rounding errors, commutative only for equal-sized vectors, i.e. changing the order of v1 and v2 will move element j in vxc to 2*half_search - j and thus change the return value from i to 2*half_search - i. If, in absence of rounding errors, more than one shift would give the maximum cross-correlation, rounding errors may cause any one of those shifts to be returned. If rounding errors have no effect the index corresponding to the shift with the smallest absolute value is returned (with preference given to the smaller of two indices that correspond to the same absolute shift). If v1 is longer than v2, the first element in v1 used for the resulting cross-correlation is max(0,shift + (cpl_vector_get_size(v1)-cpl_vector_get_size(v2))/2). Cross-correlation with half_search == 0 requires about 8n FLOPs, where n = cpl_vector_get_size(v2). Each increase of half_search by 1 requires about 4n FLOPs more, when all of v2's elements can be cross-correlated, otherwise the extra cost is about 4m, where m is the number of elements in v2 that can be cross-correlated, n - half_search <= m < n. In case of error, the #_cpl_error_code_ code is set, and the returned delta and cross-correlation is undefined. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if an input pointer is NULL - CPL_ERROR_ILLEGAL_INPUT if v1 and v2 have different sizes or if vxc is not as requested */ /*----------------------------------------------------------------------------*/ int cpl_vector_correlate(cpl_vector * vxc, const cpl_vector * v1, const cpl_vector * v2) { double xc = 0; double mean1 = 0; double mean2 = 0; double var1 = 0; double var2 = 0; int half_search; int delta; int i; int i1; int size; cpl_ensure(vxc != NULL, CPL_ERROR_NULL_INPUT, -10); cpl_ensure(v1 != NULL, CPL_ERROR_NULL_INPUT, -20); cpl_ensure(v2 != NULL, CPL_ERROR_NULL_INPUT, -30); cpl_ensure(v1->size >= v2->size, CPL_ERROR_ILLEGAL_INPUT, -50); half_search = cpl_vector_get_size(vxc)-1; cpl_ensure(!(half_search & 1), CPL_ERROR_ILLEGAL_INPUT, -60); half_search >>=1; size = v2->size; /* Number of elements used in cross-correlation */ /* 1st v1 index used - non-zero iff v1 is longer than v2 */ i1 = (v1->size - v2->size)/2; cpl_ensure(half_search <= i1 || 2*i1 == v1->size - v2->size, CPL_ERROR_ILLEGAL_INPUT, -70); /* The current implementation requires O(size * half_search) FLOPs. If size and half_search are large enough it should be preferable to compute the 1+2*half_search cross-correlations as a Toeplitz matrix vector product in terms of FFTs - since this requires O(size * log(size)) FLOPs */ /* 1st v1 index not used : i1 + size == (v1->size + size)/2 */ /* Compute mean of vectors, normalization factors and cross-correlation - over those elements used in the no-shift cross-correlation */ #if 0 /* This approach requires 10n FLOPS */ for (i=0; idata[i+i1]; mean2 += v2->data[i ]; } mean1 /= size; mean2 /= size; for (i=0; idata[i+i1] - mean1) * (v1->data[i+i1] - mean1); var2 += (v2->data[i ] - mean2) * (v2->data[i ] - mean2); xc += (v1->data[i+i1] - mean1) * (v2->data[i ] - mean2); } cpl_tools_add_flops(10 * size + 2); #else /* This approach requires only 8n FLOPS. It is only faster if v1 and v2 are already in cache - or if they do not fit in the cache. The round-off can be 10 times bigger */ for (i=0; idata[i+i1]; mean2 += v2->data[i ]; var1 += v1->data[i+i1] * v1->data[i+i1]; var2 += v2->data[i ] * v2->data[i ]; xc += v1->data[i+i1] * v2->data[i ]; } mean1 /= size; mean2 /= size; var1 -= mean1 * mean1 * size; var2 -= mean2 * mean2 * size; xc -= mean1 * mean2 * size; cpl_tools_add_flops(8 * size + 11); #endif /* var can only be zero with a constant vector - in which case xc is zero */ if ( var1 > 0 && var2 > 0) { xc /= sqrt(var1 * var2); /* xc can only be outside [-1;1] due to rounding errors */ if (xc < -1) { assert( -1-xc < size*DBL_EPSILON); xc = -1; } else if (xc > 1) { assert( xc-1 < size*DBL_EPSILON); xc = 1; } } else { /* Remove some rounding errors */ if (var1 < 0) var1 = 0; if (var2 < 0) var2 = 0; xc = 0; } delta = 0; (void)cpl_vector_set(vxc, half_search, xc); if (half_search > 0) { /* less than maximal precision OK here */ const double rsize = (double) 1 / size; const double dsize = 1 + rsize; double mean1_p = mean1; double mean1_n = mean1; double mean2_p = mean2; double mean2_n = mean2; /* The normalization factors: size * variance They are updated by use of the relation variance = sum of squares - square of sums / size */ double var1_p = var1; double var1_n = var1; double var2_p = var2; double var2_n = var2; /* No use to search further than i1 + size - 1 */ const int max_hs = i1 + size - 2; const int hs = half_search > max_hs ? max_hs : half_search; const int hs1 = hs > i1 ? i1 : hs; /* Skip non-drop part if var2 == 0 */ int step = var2 > 0 ? 1 : hs1 + 1; /* Number of elements in shortened cross-correlation */ int istop = i1 + size - step; for (; step<=hs1; step++, istop--) { /* Can cross-correlate over all of v2 - mean2 and var2 are unchanged */ double xc_p = 0; double xc_n = 0; /* Subtract term from normalization factors */ var1_p -= (v1->data[i1+step-1] * dsize - 2 * mean1_p) * v1->data[i1+step-1]; var1_n -= (v1->data[istop ] * dsize - 2 * mean1_n) * v1->data[istop ]; /* Update mean1 */ mean1_p += (v1->data[i1+size+step-1] - v1->data[i1+step-1 ]) * rsize; mean1_n += (v1->data[i1-step ] - v1->data[i1+size-step]) * rsize; /* Add term to normalization factors */ var1_n += (v1->data[i1-step ] * dsize - 2 * mean1_n) * v1->data[i1-step ]; var1_p += (v1->data[i1+size+step-1] * dsize - 2 * mean1_p) * v1->data[i1+size+step-1]; /* Compute xc - exploiting that the dot-product is distributive over subtraction */ for (i=0; idata[i] * v1->data[i+i1-step]; xc_p += v2->data[i] * v1->data[i+i1+step]; } if (var1_n > 0) { /* Subtract the mean-term */ xc_n -= mean2 * mean1_n * size; /* - and divide by the norm of the mean-corrected vectors */ xc_n /= sqrt(var1_n * var2); /* xc_n can only be outside [-1;1] due to rounding errors */ if (xc_n < -1) { assert( -1-xc_n < 4.0 * size * DBL_EPSILON); xc_n = -1; } else if (xc_n > 1) { assert( xc_n-1 < 4.0 * size * DBL_EPSILON); xc_n = 1; } } else xc_n = 0; if (xc_n > xc) { xc = xc_n; delta = -step; } (void)cpl_vector_set(vxc, half_search - step, xc_n); if (var1_p > 0) { /* Subtract the mean-term */ xc_p -= mean2 * mean1_p * size; /* - and divide by the norm of the mean-corrected vectors */ xc_p /= sqrt(var1_p * var2); /* xc_p can only be outside [-1;1] due to rounding errors */ if (xc_p < -1) { assert( -1-xc_p < 4.0 * size * DBL_EPSILON); xc_p = -1; } else if (xc_p > 1) { assert( xc_p-1 < 4.0 * size * DBL_EPSILON); xc_p = 1; } } else xc_p = 0; if (xc_p > xc) { xc = xc_p; delta = step; } (void)cpl_vector_set(vxc, half_search + step, xc_p); } assert( step == hs1 + 1 ); if (var2 > 0) cpl_tools_add_flops( hs1 * (4*size + 26) ); for (; step<=hs; step++, istop--) { /* v1 is too short - Cross-correlate on reduced number of elements - elements out of range are defined to be zero */ double xc_p = 0; double xc_n = 0; /* Subtract dropped term from normalization factors */ var1_p -= (v1->data[step+i1-1] * dsize - 2 * mean1_p) * v1->data[step+i1-1]; var1_n -= (v1->data[istop ] * dsize - 2 * mean1_n) * v1->data[istop ]; var2_p -= (v2->data[istop ] * dsize - 2 * mean2_p) * v2->data[istop ]; var2_n -= (v2->data[step-i1-1] * dsize - 2 * mean2_n) * v2->data[step-i1-1]; /* Subtract dropped elements from mean */ mean1_p -= v1->data[step+i1-1] * rsize; mean1_n -= v1->data[istop ] * rsize; mean2_p -= v2->data[istop ] * rsize; mean2_n -= v2->data[step-i1-1] * rsize; /* Compute xc */ for (i=0; idata[i] * v2->data[i-i1+step]; xc_p += v2->data[i] * v1->data[i+i1+step]; } if (var1_n * var2_n > 0) { /* Subtract the mean-term */ xc_n -= mean1_n * mean2_n * (2*size - istop); /* - and divide by the norm of the mean-corrected vectors */ xc_n /= sqrt(var1_n * var2_n); /* xc_n can only be outside [-1;1] due to rounding errors */ if (xc_n < -1) { assert( -1-xc_n < 4.0 * size * DBL_EPSILON); xc_n = -1; } else if (xc_n > 1) { assert( xc_n-1 < 4.0 * size * DBL_EPSILON); xc_n = 1; } } else xc_n = 0; if (xc_n > xc) { xc = xc_n; delta = -step; } (void)cpl_vector_set(vxc, half_search - step, xc_n); if (var1_p * var2_p > 0) { /* Subtract the mean-term */ xc_p -= mean1_p * mean2_p * (2*size - istop); /* - and divide by the norm of the mean-corrected vectors */ xc_p /= sqrt(var1_p * var2_p); /* xc_p can only be outside [-1;1] due to rounding errors */ if (xc_p < -1) { assert( -1-xc_p < 4.0 * size * DBL_EPSILON); xc_p = -1; } else if (xc_p > 1) { assert( xc_p-1 < 4.0 * size * DBL_EPSILON); xc_p = 1; } } else xc_p = 0; if (xc_p > xc) { xc = xc_p; delta = step; } (void)cpl_vector_set(vxc, half_search + step, xc_p); } if (hs > hs1) cpl_tools_add_flops( (hs-hs1) * ( 4*istop + 28 ) ); } return half_search + delta; } /*----------------------------------------------------------------------------*/ /** @brief Apply a low-pass filter to a cpl_vector @param v cpl_vector @param filter_type Type of filter to use @param hw Filter half-width @return Pointer to newly allocated cpl_vector or NULL in case of an error This type of low-pass filtering consists in a convolution with a given kernel. The chosen filter type determines the kind of kernel to apply for convolution. Supported kernels are CPL_LOWPASS_LINEAR and CPL_LOWPASS_GAUSSIAN. In the case of CPL_LOWPASS_GAUSSIAN, the gaussian sigma used is 1/sqrt(2). As this function is not meant to be general and cover all possible cases, this sigma is hardcoded and cannot be changed. The returned smooth cpl_vector must be deallocated using cpl_vector_delete(). The returned signal has exactly as many samples as the input signal. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if an input pointer is NULL - CPL_ERROR_ILLEGAL_INPUT if filter_type is not supported or if hw is negative or bigger than half the vector v size */ /*----------------------------------------------------------------------------*/ cpl_vector * cpl_vector_filter_lowpass_create( const cpl_vector * v, cpl_lowpass filter_type, int hw) { cpl_vector * filtered; cpl_vector * kernel; double replace; int i, j; cpl_ensure(v != NULL, CPL_ERROR_NULL_INPUT, NULL); cpl_ensure(hw >= 0, CPL_ERROR_ILLEGAL_INPUT, NULL); cpl_ensure(2*hw <= v->size, CPL_ERROR_ILLEGAL_INPUT, NULL); /* allocate output cpl_vector */ filtered = cpl_vector_new(v->size); /* generate low-pass filter kernel */ kernel = cpl_vector_gen_lowpass_kernel(filter_type, hw); /* This will catch illegal values for filter_type */ cpl_ensure(kernel != NULL, CPL_ERROR_ILLEGAL_INPUT, NULL); /* compute edge effects for the first hw elements */ for (i=0; idata[hw+j] * v->data[0]; } else { replace += kernel->data[hw+j] * v->data[i+j]; } } filtered->data[i] = replace; } /* compute edge effects for the last hw elements */ for (i=v->size-hw; isize; i++) { replace = 0.0; for (j=-hw; j<=hw; j++) { if (i+j>v->size-1) { replace += kernel->data[hw+j] * v->data[v->size-1]; } else { replace += kernel->data[hw+j] * v->data[i+j]; } } filtered->data[i] = replace; } /* compute all other elements */ for (i=hw; isize-hw; i++) { replace = 0.0; for (j=-hw; j<=hw; j++) { replace += kernel->data[hw+j] * v->data[i+j]; } filtered->data[i] = replace; } cpl_vector_delete(kernel); cpl_tools_add_flops(4*hw*v->size); return filtered; } /*----------------------------------------------------------------------------*/ /** @brief Apply a 1d median filter of given half-width to a cpl_vector @param v cpl_vector @param hw Filter half-width @return Pointer to newly allocated cpl_vector or NULL in case of an error This function applies a median smoothing to a cpl_vector and returns a newly allocated cpl_vector containing a median-smoothed version of the input. The returned cpl_vector has exactly as many samples as the input one. It must be deallocated using cpl_vector_delete(). For half-widths of 1,2,3,4, the filtering is optimised for speed. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if an input pointer is NULL - CPL_ERROR_ILLEGAL_INPUT - CPL_ERROR_ILLEGAL_INPUT if hw is negative or bigger than half the vector v size */ /*----------------------------------------------------------------------------*/ cpl_vector * cpl_vector_filter_median_create( const cpl_vector * v, int hw) { cpl_vector * filtered; double * row; int i, j; cpl_ensure(v != NULL, CPL_ERROR_NULL_INPUT, NULL); cpl_ensure(hw >= 0, CPL_ERROR_ILLEGAL_INPUT, NULL); cpl_ensure(2*hw <= v->size, CPL_ERROR_ILLEGAL_INPUT, NULL); /* Create the filtered vector */ filtered = cpl_vector_new(v->size); /* simply copy first hw and last hw items */ for (i=0; idata[i] = v->data[i]; for (i=v->size-hw; isize; i++) filtered->data[i] = v->data[i]; /* median filter on all central items */ row = cpl_malloc((2*hw+1) * sizeof(double)); for (i=hw; isize-hw; i++) { for (j=-hw; j<=hw; j++) row[j+hw] = v->data[i+j]; filtered->data[i] = cpl_tools_get_median_double(row, 2*hw+1); } cpl_free(row); return filtered; } /*----------------------------------------------------------------------------*/ /** @brief Elementwise addition of a scalar to a vector @param v cpl_vector to modify @param addend Number to add @return CPL_ERROR_NONE or the relevant #_cpl_error_code_ on error Add a number to each element of the cpl_vector. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if an input pointer is NULL */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_vector_add_scalar( cpl_vector * v, double addend) { int i; /* Test inputs */ cpl_ensure_code(v != NULL, CPL_ERROR_NULL_INPUT); /* Perform the addition */ for (i=0; isize; i++) v->data[i] += addend; cpl_tools_add_flops(v->size); return CPL_ERROR_NONE; } /*----------------------------------------------------------------------------*/ /** @brief Elementwise subtraction of a scalar from a vector @param v cpl_vector to modify @param subtrahend Number to subtract @return CPL_ERROR_NONE or the relevant #_cpl_error_code_ on error Subtract a number from each element of the cpl_vector. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if an input pointer is NULL */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_vector_subtract_scalar( cpl_vector * v, double subtrahend) { int i; /* Test inputs */ cpl_ensure_code(v != NULL, CPL_ERROR_NULL_INPUT); /* Perform the subtraction */ for (i=0; isize; i++) v->data[i] -= subtrahend; cpl_tools_add_flops(v->size); return CPL_ERROR_NONE; } /*----------------------------------------------------------------------------*/ /** @brief Elementwise multiplication of a vector with a scalar @param v cpl_vector to modify @param factor Number to multiply with @return CPL_ERROR_NONE or the relevant #_cpl_error_code_ on error Multiply each element of the cpl_vector with a number. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if an input pointer is NULL */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_vector_multiply_scalar( cpl_vector * v, double factor) { int i; /* Test inputs */ cpl_ensure_code(v != NULL, CPL_ERROR_NULL_INPUT); /* Perform the subtraction */ for (i=0; isize; i++) v->data[i] *= factor; cpl_tools_add_flops(v->size); return CPL_ERROR_NONE; } /*----------------------------------------------------------------------------*/ /** @brief Elementwise division of a vector with a scalar @param v cpl_vector to modify @param divisor Non-zero number to divide with @return CPL_ERROR_NONE or the relevant #_cpl_error_code_ on error Divide each element of the cpl_vector with a number. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if an input pointer is NULL - CPL_ERROR_DIVISION_BY_ZERO if divisor is 0.0 */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_vector_divide_scalar( cpl_vector * v, double divisor) { int i; /* Test inputs */ cpl_ensure_code(v != NULL, CPL_ERROR_NULL_INPUT); cpl_ensure_code(divisor != 0, CPL_ERROR_DIVISION_BY_ZERO); /* Perform the subtraction */ for (i=0; isize; i++) v->data[i] /= divisor; cpl_tools_add_flops(v->size); return CPL_ERROR_NONE; } /*----------------------------------------------------------------------------*/ /** @brief Compute the element-wise logarithm. @param v cpl_vector to modify. @param base Logarithm base. @return CPL_ERROR_NONE or the relevant #_cpl_error_code_ on error The base and all the vector elements must be positive and the base must be different from 1, or a cpl_error_code will be returned and the vector will be left unmodified. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if an input pointer is NULL - CPL_ERROR_ILLEGAL_INPUT if base is negative or zero or if one of the vector values is negative or zero - CPL_ERROR_DIVISION_BY_ZERO if a division by zero occurs */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_vector_logarithm( cpl_vector * v, double base) { int i; cpl_ensure_code(v != NULL, CPL_ERROR_NULL_INPUT); cpl_ensure_code(base > 0, CPL_ERROR_ILLEGAL_INPUT); cpl_ensure_code(base != 1, CPL_ERROR_DIVISION_BY_ZERO); for (i = 0; i < v->size; i++) cpl_ensure_code(v->data[i] > 0, CPL_ERROR_ILLEGAL_INPUT); for (i = 0; i < v->size; i++) v->data[i] = log(v->data[i]) / log(base); cpl_tools_add_flops(2*v->size); return CPL_ERROR_NONE; } /*----------------------------------------------------------------------------*/ /** @brief Compute the exponential of all vector elements. @param v Target cpl_vector. @param base Exponential base. @return CPL_ERROR_NONE or the relevant #_cpl_error_code_ on error If the base is zero all vector elements must be positive and if the base is negative all vector elements must be integer, otherwise a cpl_error_code is returned and the vector is unmodified. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if an input pointer is NULL - CPL_ERROR_ILLEGAL_INPUT base and v are not as requested - CPL_ERROR_DIVISION_BY_ZERO if one of the v values is negative or 0 */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_vector_exponential( cpl_vector * v, double base) { int i; cpl_ensure_code(v != NULL, CPL_ERROR_NULL_INPUT); if (base == 0) { for (i = 0; i < v->size; i++) cpl_ensure_code(v->data[i] > 0, CPL_ERROR_DIVISION_BY_ZERO); } else if (base < 0) { for (i = 0; i < v->size; i++) cpl_ensure_code(v->data[i] == ceil(v->data[i]), CPL_ERROR_ILLEGAL_INPUT); } for (i = 0; i < v->size; i++) v->data[i] = pow(base, v->data[i]); cpl_tools_add_flops(v->size); return CPL_ERROR_NONE; } /*----------------------------------------------------------------------------*/ /** @brief Compute the power of all vector elements. @param v Target cpl_vector. @param exponent Constant exponent. @return CPL_ERROR_NONE or the relevant #_cpl_error_code_ on error If the exponent is non-positive all vector elements must be non-zero and if the exponent is non-integer all vector elements must be non-negative, otherwise a cpl_error_code is returned and the vector is unmodified. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if an input pointer is NULL - CPL_ERROR_ILLEGAL_INPUT if v and exponent are not as requested - CPL_ERROR_DIVISION_BY_ZERO if one of the v values is 0 */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_vector_power( cpl_vector * v, double exponent) { int i; cpl_ensure_code(v != NULL, CPL_ERROR_NULL_INPUT); if (exponent <= 0) { for (i = 0; i < v->size; i++) cpl_ensure_code(v->data[i] != 0, CPL_ERROR_DIVISION_BY_ZERO); } else if (exponent != ceil(exponent)) { for (i = 0; i < v->size; i++) cpl_ensure_code(v->data[i] >= 0, CPL_ERROR_ILLEGAL_INPUT); } for (i = 0; i < v->size; i++) v->data[i] = pow(v->data[i], exponent); cpl_tools_add_flops(v->size); return CPL_ERROR_NONE; } /*----------------------------------------------------------------------------*/ /** @brief Fill a vector with a kernel profile @param profile cpl_vector to be filled @param type Type of kernel profile @param radius Radius of the profile in pixels @return CPL_ERROR_NONE or the relevant #_cpl_error_code_ on error @see cpl_image_get_interpolated A number of predefined kernel profiles are available: - CPL_KERNEL_DEFAULT: default kernel, currently CPL_KERNEL_TANH - CPL_KERNEL_TANH: Hyperbolic tangent - CPL_KERNEL_SINC: Sinus cardinal - CPL_KERNEL_SINC2: Square sinus cardinal - CPL_KERNEL_LANCZOS: Lanczos2 kernel - CPL_KERNEL_HAMMING: Hamming kernel - CPL_KERNEL_HANN: Hann kernel Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if an input pointer is NULL - CPL_ERROR_ILLEGAL_INPUT if radius is non-positive */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_vector_fill_kernel_profile(cpl_vector * profile, cpl_kernel type, double radius) { const int size = cpl_vector_get_size(profile); double dx; int i; cpl_ensure_code(size > 0, cpl_error_get_code()); cpl_ensure_code(radius > 0, CPL_ERROR_ILLEGAL_INPUT); dx = radius / (size-1); /* TABS per pixel == 1 / dx */ /* Switch on the requested kernel type */ switch (type) { case CPL_KERNEL_TANH: cpl_vector_fill_tanh_kernel(profile, 0.5 / dx); break; case CPL_KERNEL_SINC: for (i=0; i < size; i++) { const double x = i * dx; profile->data[i] = cpl_tools_sinc(x); } cpl_tools_add_flops( 4 * size ); break; case CPL_KERNEL_SINC2: for (i=0; i < size; i++) { const double x = i * dx; profile->data[i] = cpl_tools_sinc(x); profile->data[i] *= profile->data[i]; } cpl_tools_add_flops( 5 * size ); break; case CPL_KERNEL_LANCZOS: for (i=0; i < size; i++) { const double x = i * dx; if (x >= 2) break; profile->data[i] = cpl_tools_sinc(x) * cpl_tools_sinc(x/2); } for (; i < size; i++) profile->data[i] = 0; cpl_tools_add_flops( 8 * i ); break; case CPL_KERNEL_HAMMING: cpl_vector_fill_alpha_kernel(profile, dx, 0.54); break; case CPL_KERNEL_HANN: cpl_vector_fill_alpha_kernel(profile, dx, 0.50); break; default: /* Only reached if cpl_kernel is extended in error */ assert( 0 ); } return CPL_ERROR_NONE; } /*----------------------------------------------------------------------------*/ /** @brief Apply a 1d gaussian fit @param x Positions to fit @param sigma_x Uncertainty (one sigma, gaussian errors assumed) assosiated with @em x. Taking into account the uncertainty of the independent variable is currently unsupported, and this parameter must therefore be set to NULL. @param y Values to fit @param sigma_y Uncertainty (one sigma, gaussian errors assumed) associated with @em y. If NULL, constant uncertainties are assumed. @param fit_pars Specifies which parameters participate in the fit (any other parameters will be held constant). Possible values are CPL_FIT_CENTROID, CPL_FIT_STDEV, CPL_FIT_AREA, CPL_FIT_OFFSET and any bitwise combination of these. As a shorthand for including all four parameters in the fit, use CPL_FIT_ALL. @param x0 (output) Center of best fit gaussian. If CPL_FIT_CENTROID is not set, this is also an input parameter. @param sigma (output) Width of best fit gaussian. A positive number on success. If CPL_FIT_STDEV is not set, this is also an input parameter. @param area (output) Area of gaussian. A positive number on succes. If CPL_FIT_AREA is not set, this is also an input parameter. @param offset (output) Fitted background level. If CPL_FIT_OFFSET is not set, this is also an input parameter. @param mse (output) If non-NULL, the mean squared error of the best fit is returned. @param red_chisq (output) If non-NULL, the reduced chi square of the best fit is returned. This requires the noise vector to be specified. @param covariance (output) If non-NULL, the formal covariance matrix of the best fit is returned. This requires @em sigma_y to be specified. The order of fit parameters in the covariance matrix is defined as (@em x0, @em sigma, @em area, @em offset), for example the (3,3) element of the matrix (counting from zero) is the variance of the fitted @em offset. The matrix must be deallocated by calling @c cpl_matrix_delete() . On error, NULL is returned. @return CPL_ERROR_NONE iff okay This function fits to the input vectors a 1d gaussian function of the form f(x) = @em area / sqrt(2 pi @em sigma^2) * exp( -(@em x - @em x0)^2/(2 @em sigma^2)) + @em offset (@em area > 0) by minimizing chi^2 using a Levenberg-Marquardt algorithm. The values to fit are read from the input vector @em x. Optionally, a vector @em sigma_x (of same size as @em x) may be specified. Optionally, the mean squared error, the reduced chi square and the covariance matrix of the best fit are computed . Set corresponding parameter to NULL to ignore. If the covariance matrix is requested and successfully computed, the diagonal elements (the variances) are guaranteed to be positive. Occasionally, the Levenberg-Marquardt algorithm fails to converge to a set of sensible parameters. In this case (and only in this case), a CPL_ERROR_CONTINUE is set. To allow the caller to recover from this particular error, the parameters @em x0, @em sigma, @em area and @em offset will on output contain estimates of the best fit parameters, specifically estimated as the median position, the median of the absolute residuals multiplied by 1.4828, the minimum flux value and the maximum flux difference multiplied by sqrt(2 pi sigma^2), respectively. In this case, @em mse, @em red_chisq and @em covariance are not computed. Note that the variance of @em x0 (the (0,0) element of the covariance matrix) in this case can be estimated by @em sigma^2 / @em area . A CPL_ERROR_SINGULAR_MATRIX occurs if the covariance matrix cannot be computed. In that case all other output parameters are valid. Current limitations - Taking into account the uncertainties of the independent variable is not supported. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if @em x, @em y, @em x0, @em sigma, @em area or @em offset is NULL. - CPL_ERROR_INVALID_TYPE if the specified @em fit_pars is not a bitwise combination of the allowed values (e.g. 0 or 1). - CPL_ERROR_UNSUPPORTED_MODE @em sigma_x is non-NULL. - CPL_ERROR_INCOMPATIBLE_INPUT if the sizes of any input vectors are different, or if the computation of reduced chi square or covariance is requested, but @em sigma_y is not provided. - CPL_ERROR_ILLEGAL_INPUT if any input noise values, @em sigma or @em area is non-positive, or if chi square computation is requested and there are less than 5 data points to fit. - CPL_ERROR_ILLEGAL_OUTPUT if memory allocation failed. - CPL_ERROR_CONTINUE if the fitting algorithm failed. - CPL_ERROR_SINGULAR_MATRIX if the covariance matrix could not be calculated. */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_vector_fit_gaussian(const cpl_vector *x, const cpl_vector *sigma_x, const cpl_vector *y, const cpl_vector *sigma_y, cpl_fit_mode fit_pars, double *x0, double *sigma, double *area, double *offset, double *mse, double *red_chisq, cpl_matrix **covariance) { cpl_matrix *x_matrix = NULL; /* LM algorithm needs a matrix, not a vector */ int N; /* Number of data points */ double xlo, xhi; /* Min/max x */ /* Initial parameter values */ double x0_guess = 0; /* Avoid warnings about uninitialized variables */ double sigma_guess = 0; double area_guess; double offset_guess; /* Validate input */ cpl_ensure_code( x != NULL, CPL_ERROR_NULL_INPUT); cpl_ensure_code( sigma_x == NULL, CPL_ERROR_UNSUPPORTED_MODE); cpl_ensure_code( y != NULL, CPL_ERROR_NULL_INPUT); /* sigma_y may be NULL or non-NULL */ /* Bits in 'fit_pars' must be non-empty subset of bits in 'CPL_FIT_ALL' */ cpl_ensure_code( fit_pars != 0 && (CPL_FIT_ALL | fit_pars) == CPL_FIT_ALL, CPL_ERROR_INVALID_TYPE); N = cpl_vector_get_size(x); cpl_ensure_code( N == cpl_vector_get_size(y), CPL_ERROR_INCOMPATIBLE_INPUT); if (sigma_x != NULL) { cpl_ensure_code( N == cpl_vector_get_size(sigma_x), CPL_ERROR_INCOMPATIBLE_INPUT); } if (sigma_y != NULL) { cpl_ensure_code( N == cpl_vector_get_size(sigma_y), CPL_ERROR_INCOMPATIBLE_INPUT); } cpl_ensure_code( x0 != NULL && sigma != NULL && area != NULL && offset != NULL, CPL_ERROR_NULL_INPUT); if (! (fit_pars & CPL_FIT_STDEV)) { cpl_ensure_code( *sigma > 0, CPL_ERROR_ILLEGAL_INPUT); } if (! (fit_pars & CPL_FIT_AREA)) { cpl_ensure_code( *area > 0, CPL_ERROR_ILLEGAL_INPUT); } /* mse, red_chisq may be NULL */ /* Need more than number_of_parameters points to calculate chi^2. * There are less than 5 parameters. */ cpl_ensure_code( red_chisq == NULL || N >= 5, CPL_ERROR_ILLEGAL_INPUT); if (covariance != NULL) *covariance = NULL; /* If covariance computation is requested, then * return either the covariance matrix or NULL * (don't return undefined pointer). */ /* Cannot compute chi square & covariance without sigma_y */ cpl_ensure_code( (red_chisq == NULL && covariance == NULL) || sigma_y != NULL, CPL_ERROR_INCOMPATIBLE_INPUT); /* Create (non-modified) matrix from x-data */ x_matrix = cpl_matrix_wrap(N, 1, (double*)cpl_vector_get_data_const(x)); if (x_matrix == NULL) { cpl_ensure_code( CPL_FALSE, CPL_ERROR_ILLEGAL_OUTPUT); } /* Check that any provided sigmas are positive. */ if (sigma_x != NULL && cpl_vector_get_min(sigma_x) <= 0) { cpl_matrix_unwrap(x_matrix); cpl_ensure_code( CPL_FALSE, CPL_ERROR_ILLEGAL_INPUT); } if (sigma_y != NULL && cpl_vector_get_min(sigma_y) <= 0) { cpl_matrix_unwrap(x_matrix); cpl_ensure_code( CPL_FALSE, CPL_ERROR_ILLEGAL_INPUT); } /* Compute guess parameters using robust estimation * (This step might be improved by taking into account the sigmas, * but for simplicity do not) */ { double sum = 0.0; double quartile[3]; double fraction[3] = {0.25 , 0.50 , 0.75}; const double *y_data = cpl_vector_get_data_const(y); if (fit_pars & CPL_FIT_OFFSET) { /* Estimate offset as 25% percentile of y-values. (The minimum value may be too low for noisy input, the median is too high if there is not much background in the supplied data, so use something inbetween). */ cpl_vector *y_dup = cpl_vector_duplicate(y); if (y_dup == NULL) { cpl_matrix_unwrap(x_matrix); cpl_ensure_code( CPL_FALSE, CPL_ERROR_ILLEGAL_OUTPUT); } offset_guess = cpl_tools_get_kth_double( cpl_vector_get_data(y_dup), N, N/4); cpl_vector_delete(y_dup); } else { offset_guess = *offset; } /* Get quartiles of distribution (only bother if it's needed for estimation of x0 or sigma) */ if ( (fit_pars & CPL_FIT_CENTROID) || (fit_pars & CPL_FIT_STDEV ) ) { /* The algorithm requires the input to be sorted as function of x, so do that (using qsort), and work on the sorted copy. Of course, the y-vector must be re-ordered along with x. sigma_x and sigma_y are not used, so don't copy those. */ cpl_vector_fit_gaussian_input *sorted_input = cpl_malloc(N * sizeof(*sorted_input)); double *x_data = cpl_matrix_get_data(x_matrix); cpl_boolean is_sorted = CPL_TRUE; int i; if (sorted_input == NULL) { cpl_matrix_unwrap(x_matrix); cpl_ensure_code( CPL_FALSE, CPL_ERROR_ILLEGAL_INPUT); } for (i = 0; i < N; i++) { sorted_input[i].x = x_data[i]; sorted_input[i].y = y_data[i]; is_sorted = is_sorted && (i==0 || (x_data[i-1] < x_data[i])); } /* For efficiency, sort only if necessary. This is to * keep the asymptotic time complexity at O(n) if the * input was already sorted. (Otherwise the time * complexity is dominated by the following qsort call.) */ if (!is_sorted) { qsort(sorted_input, N, sizeof(*sorted_input), &cpl_fit_gaussian_1d_compare); } for(i = 0; i < N; i++) { double flux = sorted_input[i].y; sum += (flux - offset_guess); } /* Note that 'sum' must be calculated the same way as 'running_sum' below, Otherwise (due to round-off error) 'running_sum' might end up being different from 'sum'(!). Specifically, it will not work to calculate 'sum' as (flux1 + ... + fluxN) - N*offset_guess */ if (sum > 0.0) { double flux, x1, x2; double running_sum = 0.0; int j; i = 0; flux = sorted_input[i].y - offset_guess; for (j = 0; j < 3; j++) { double limit = fraction[j] * sum; double k; while (running_sum + flux < limit && i < N-1) { running_sum += flux; i++; flux = sorted_input[i].y - offset_guess; } /* Fraction [0;1] of current flux needed to reach the quartile */ k = (limit - running_sum)/flux; if (k <= 0.5) { /* Interpolate linearly between current and previous position */ if (0 < i) { x1 = sorted_input[i-1].x; x2 = sorted_input[i].x; quartile[j] = x1*(0.5-k) + x2*(0.5+k); /* k=0 => quartile = midpoint, k=0.5 => quartile = x2 */ } else { quartile[j] = sorted_input[i].x; } } else { /* Interpolate linearly between current and next position */ if (i < N-1) { x1 = sorted_input[i].x; x2 = sorted_input[i+1].x; quartile[j] = x1*( 1.5-k) + x2*(-0.5+k); /* k=0.5 => quartile = x1, k=1.0 => quartile = midpoint */ } else { quartile[j] = sorted_input[i].x; } } } } else { /* If there's no flux (sum = 0) then set quartiles to something that's not completely insensible. */ quartile[1] = cpl_matrix_get_mean(x_matrix); quartile[2] = quartile[1]; quartile[0] = quartile[2] - 1.0; /* Then sigma_guess = 1.0 */ } cpl_free(sorted_input); } /* x0_guess = median of distribution */ x0_guess = (fit_pars & CPL_FIT_CENTROID) ? quartile[1] : *x0; /* sigma_guess = median of absolute residuals * * (68% is within 1 sigma, and 50% is within 0.6744 * sigma, => quartile3-quartile1 = 2*0.6744 sigma) */ sigma_guess = (fit_pars & CPL_FIT_STDEV) ? (quartile[2] - quartile[0]) / (2*0.6744) : *sigma; area_guess = (fit_pars & CPL_FIT_AREA) ? (cpl_vector_get_max(y) - offset_guess) * CPL_MATH_SQRT2PI * sigma_guess : *area; /* Make sure that the area/sigma are positive number */ if (area_guess <= 0) area_guess = 1.0; if (sigma_guess <= 0) sigma_guess = 1.0; } /* Wrap parameters, fit, unwrap */ { cpl_vector *a = cpl_vector_new(4); /* There are four parameters */ int ia[4]; cpl_error_code ec; if (a == NULL) { cpl_matrix_unwrap(x_matrix); cpl_ensure_code( CPL_FALSE, CPL_ERROR_ILLEGAL_OUTPUT); } cpl_vector_set(a, 0, x0_guess); cpl_vector_set(a, 1, sigma_guess); cpl_vector_set(a, 2, area_guess); cpl_vector_set(a, 3, offset_guess); ia[0] = fit_pars & CPL_FIT_CENTROID; ia[1] = fit_pars & CPL_FIT_STDEV; ia[2] = fit_pars & CPL_FIT_AREA; ia[3] = fit_pars & CPL_FIT_OFFSET; ec = _cpl_fit_lvmq(x_matrix, NULL, y, sigma_y, a, ia, gauss, gauss_derivative, CPL_FIT_LVMQ_TOLERANCE, CPL_FIT_LVMQ_COUNT, CPL_FIT_LVMQ_MAXITER, mse, red_chisq, covariance); cpl_matrix_unwrap(x_matrix); /* Check return status of fitting routine */ if (ec == CPL_ERROR_NONE || ec == CPL_ERROR_SINGULAR_MATRIX) { /* The LM algorithm converged. The computation * of the covariance matrix might have failed. */ /* In principle, the LM algorithm might have converged * to a negative sigma (even if the guess value was * positive). Make sure that the returned sigma is positive * (by convention). */ if (CPL_FIT_CENTROID) *x0 = cpl_vector_get(a, 0); if (CPL_FIT_STDEV ) *sigma = fabs(cpl_vector_get(a, 1)); if (CPL_FIT_AREA ) *area = cpl_vector_get(a, 2); if (CPL_FIT_OFFSET ) *offset = cpl_vector_get(a, 3); } cpl_vector_delete(a); xlo = cpl_vector_get_min(x); xhi = cpl_vector_get_max(x); if (ec == CPL_ERROR_CONTINUE || !( xlo <= *x0 && *x0 <= xhi && 0 < *sigma && *sigma < (xhi - xlo + 1) && 0 < *area /* This extra check on the background level makes sense iff the input flux is assumed to be positive && *offset > - *area */ ) ) { /* The LM algorithm did not converge, or it converged to * a non-sensical result. Return the guess parameter values * in order to enable the caller to recover. */ *x0 = x0_guess; *sigma = sigma_guess; *area = area_guess; *offset = offset_guess; /* In this case the covariance matrix will not make sense (because the LM algorithm failed), so delete it */ if (covariance != NULL && *covariance != NULL) { cpl_matrix_delete(*covariance); *covariance = NULL; } /* Return CPL_ERROR_CONTINUE if the fitting routine failed */ cpl_ensure_code( CPL_FALSE, CPL_ERROR_CONTINUE); } } return CPL_ERROR_NONE; } /**@}*/ /*----------------------------------------------------------------------------*/ /** @brief Append a vector to a FITS file @param to_save Vector to write in extension @param filename Name of the file to write @param bpp Bits per pixel @param pl Property list for the output header @return the #_cpl_error_code_ or CPL_ERROR_NONE @see cpl_vector_save() This function appends a vector to a FITS file, using cfitsio. If a property list is provided, it is used as the extension header. If the vector is not provided, the extension will only contain the header. The requested pixel depth (bpp) follows the FITS convention. Possible values are CPL_BPP_8_UNSIGNED (8), CPL_BPP_16_SIGNED (16), CPL_BPP_32_SIGNED (32), CPL_BPP_IEEE_FLOAT (-32), CPL_BPP_IEEE_DOUBLE (-64) Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if (one of) the input pointer(s) is NULL - CPL_ERROR_ILLEGAL_INPUT if the bpp value or the mode is not supported - CPL_ERROR_FILE_NOT_CREATED if the output file cannot be created - CPL_ERROR_FILE_IO if the data cannot be written in the file */ /*----------------------------------------------------------------------------*/ static cpl_error_code cpl_vector_append( const cpl_vector * to_save, const char * filename, cpl_type_bpp bpp, cpl_propertylist * pl) { fitsfile * fptr ; int fio_status=0 ; long naxes[1] ; int ival ; double dval ; /* Test entries */ cpl_ensure_code(filename, CPL_ERROR_NULL_INPUT) ; cpl_ensure_code(bpp==CPL_BPP_8_UNSIGNED || bpp==CPL_BPP_16_SIGNED || bpp==CPL_BPP_32_SIGNED || bpp==CPL_BPP_IEEE_FLOAT || bpp==CPL_BPP_IEEE_DOUBLE, CPL_ERROR_ILLEGAL_INPUT) ; /* Open the file */ fits_open_file(&fptr, filename, READWRITE, &fio_status); cpl_ensure_code(fio_status==0, CPL_ERROR_FILE_IO) ; /* Create the vector and append it */ if (to_save != NULL) { naxes[0] = to_save->size ; fits_create_img(fptr, bpp, 1, naxes, &fio_status); /* Write the data */ fits_write_img(fptr, TDOUBLE, 1, naxes[0], to_save->data, &fio_status); } else { /* Create minimal header */ fits_create_hdu(fptr, &fio_status); /* FIXME: Cast string literal to avoid compiler warnings - and pray that CFITSIO does not try to modify them :-(((( */ fits_write_key(fptr, TSTRING, (char*)"XTENSION", (char*)"VEXTOR", (char*)"FITS Vector Extension", &fio_status); fits_write_key(fptr, TINT, (char*)"BITPIX", &bpp, (char*)"Bits per pixel", &fio_status); ival = 0 ; fits_write_key(fptr, TINT, (char*)"NAXIS", &ival, (char*)"Empty unit", &fio_status); } /* Add BSCALE BZERO */ dval = 1.0 ; fits_write_key(fptr, TDOUBLE, (char*)"BSCALE", &dval, (char*)"Scale factor", &fio_status); dval = 0.0 ; fits_write_key(fptr, TDOUBLE, (char*)"BZERO", &dval, (char*)"Offset factor", &fio_status); /* Check */ if (fio_status != 0) { fits_close_file(fptr, &fio_status) ; cpl_ensure_code(0, CPL_ERROR_FILE_IO) ; } /* Add the property list */ if (cpl_fits_add_properties(fptr, pl, CPL_FITS_BADKEYS_EXT "|" CPL_FITS_COMPRKEYS)!=CPL_ERROR_NONE) { fits_close_file(fptr, &fio_status) ; cpl_ensure_code(0, CPL_ERROR_ILLEGAL_INPUT) ; } /* Close and write on disk */ fits_close_file(fptr, &fio_status) ; return CPL_ERROR_NONE ; } /*----------------------------------------------------------------------------*/ /** @brief Define order of input data @param left Left operand @param right Right operand This function uses void pointers because it is used as input for ANSI's qsort(). */ /*----------------------------------------------------------------------------*/ static int cpl_fit_gaussian_1d_compare(const void *left, const void *right) { return (((cpl_vector_fit_gaussian_input *)left )->x < ((cpl_vector_fit_gaussian_input *)right)->x) ? -1 : (((cpl_vector_fit_gaussian_input *)left )->x == ((cpl_vector_fit_gaussian_input *)right)->x) ? 0 : 1; } /*----------------------------------------------------------------------------*/ /** @brief Generate a kernel for smoothing filters (low-pass) @param filt_type Type of kernel to generate @param hw Kernel half-width. @return Pointer to newly allocated cpl_vector or NULL in case of an error @see cpl_vector_filter_lowpass_create() The kernel is returned as a cpl_vector (values of the Y field) The returned cpl_vector must be deallocated using cpl_vector_delete(). The returned cpl_vector's size is 2hw+1. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_ILLEGAL_INPUT if hw is negative */ /*----------------------------------------------------------------------------*/ static cpl_vector * cpl_vector_gen_lowpass_kernel( cpl_lowpass filt_type, int hw) { cpl_vector * kernel; double norm; int i; cpl_ensure(hw >= 0, CPL_ERROR_ILLEGAL_INPUT, NULL); /* Create the kernel */ kernel = cpl_vector_new(2*hw+1); switch(filt_type) { case CPL_LOWPASS_LINEAR: /* flat kernel */ for (i=-hw; i<=hw; i++) kernel->data[hw+i] = 1.0/(double)(2*hw+1); break; case CPL_LOWPASS_GAUSSIAN: norm = 0.00; for (i=-hw; i<=hw; i++) { /* gaussian kernel */ kernel->data[hw+i] = exp(-(double)(i*i)); norm += kernel->data[hw+i]; } for (i=0; i<2*hw+1; i++) kernel->data[i] /= norm; cpl_tools_add_flops(2*hw*4); break; default: /* Only reached if cpl_lowpass is extended in error */ assert( 0 ); break; } return kernel; } /*----------------------------------------------------------------------------*/ /** @brief Cardinal sine, sin(pi*x)/(pi*x). @param x double value. @return The cardinal sine of x. */ /*----------------------------------------------------------------------------*/ static double cpl_tools_sinc(double x) { const double xpi = x * CPL_MATH_PI; /* sinc(0) == 1, sinx(x*pi) = 0, for integral x */ return x == 0 ? 1 : (ceil(x) == x ? 0 : sin(xpi) / xpi); } /*----------------------------------------------------------------------------*/ /** @brief Fill a vector with an alpha kernel. @param dx x increment @param alpha Alpha @return void The function will SIGABRT on invalid input. */ /*----------------------------------------------------------------------------*/ static void cpl_vector_fill_alpha_kernel(cpl_vector * profile, double dx, double alpha) { int i; assert( profile ); for (i=0; i < profile->size; i++) { const double x = i * dx; if (x >= 1) break; profile->data[i] = alpha + (1-alpha) * cos (x * CPL_MATH_PI); } cpl_tools_add_flops( 5 * i ); for (; i < profile->size; i++) profile->data[i] = 0; return; } /*----------------------------------------------------------------------------*/ /** @brief Fill a vector with a hyperbolic tangent kernel. @param width Half of the number of profile values per pixel @return void The following function builds up a good approximation of a box filter. It is built from a product of hyperbolic tangents. It has the following properties: \begin{itemize} \item It converges very quickly towards +/- 1. \item The converging transition is very sharp. \item It is infinitely differentiable everywhere (i.e. smooth). \item The transition sharpness is scalable. \end{itemize} The function will SIGABRT on invalid input. */ /*----------------------------------------------------------------------------*/ /* Kernel definitions */ #ifndef CPL_KERNEL_TANH_STEEPNESS #define CPL_KERNEL_TANH_STEEPNESS 5 #endif #define CPL_hk_gen(x,s) (((tanh(s*(x+0.5))+1)/2)*((tanh(s*(-x+0.5))+1)/2)) static void cpl_vector_fill_tanh_kernel(cpl_vector * profile, double width) { double * x; const int np = 32768; /* Hardcoded: should never be changed */ const double inv_np = 1 / (double)np; const double steep = CPL_KERNEL_TANH_STEEPNESS; double ind; const int samples = cpl_vector_get_size(profile); int i; assert( samples > 0); /* * Generate the kernel expression in Fourier space * with a correct frequency ordering to allow standard FT */ x = cpl_malloc((2*np+1)*sizeof(double)); for (i=0; i < np/2; i++) { ind = i * 2 * width * inv_np; x[2*i] = CPL_hk_gen(ind, steep); x[2*i+1] = 0; } for (i=np/2; idata[i] = 2.0 * width * x[2*i] * inv_np; } cpl_free(x); cpl_tools_add_flops( (13 * np + samples * 4)/2 ); return; } /*----------------------------------------------------------------------------*/ /** @brief Bring a hyperbolic tangent kernel from Fourier to normal space. @param data Kernel samples in Fourier space. @param nn Number of samples in the input kernel. @return void Bring back a hyperbolic tangent kernel from Fourier to normal space. Do not try to understand the implementation and DO NOT MODIFY THIS FUNCTION. */ /*----------------------------------------------------------------------------*/ #define CPL_KERNEL_SW(a,b) tempr=(a);(a)=(b);(b)=tempr static void cpl_vector_reverse_tanh_kernel(double * data, int nn) { unsigned long n, mmax, m, i, j, istep; double wtemp, wr, wpr, wpi, wi, theta; double tempr, tempi; n = (unsigned long)nn << 1; j = 1; for (i=1; i i) { CPL_KERNEL_SW(data[j-1],data[i-1]); CPL_KERNEL_SW(data[j],data[i]); } m = n >> 1; while (m>=2 && j>m) { j -= m; m >>= 1; } j += m; } mmax = 2; while (n > mmax) { istep = mmax << 1; theta = CPL_MATH_2PI / mmax; wtemp = sin(0.5 * theta); wpr = -2.0 * wtemp * wtemp; wpi = sin(theta); wr = 1.0; wi = 0.0; for (m=1; m