/** @file psFunctions.h
 *  @brief Standard Mathematical Functions.
 *  @ingroup Stats
 *
 *  This file will hold the prototypes for procedures which allocate, free,
 *  and evaluate various polynomials.  Those polynomial structures are also
 *  defined here.
 *
 *  @ingroup Stats
 *
 *  @author Someone at IfA
 *  @author GLG, MHPCC
 *
 *  @version $Revision: 1.46 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2005-06-09 19:26:48 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
 */

#ifndef PS_FUNCTIONS_H
#define PS_FUNCTIONS_H

#include <stdio.h>
#include <stdbool.h>
#include <float.h>
#include <math.h>

#include "psVector.h"
#include "psScalar.h"

/** \addtogroup Stats
 *  \{
 */

/** Evaluate a non-normalized Gaussian with the given mean and sigma at the
 *  given coordianate.
 *
 *  Note that this is not a Gaussian deviate.  The evaluated Gaussian is:
 *        \f[ exp(-\frac{(x-mean)^2}{2\sigma^2}) \f]
 *
 *  @return psF32      value on the gaussian curve given the input parameters
 */
psF32 psGaussian(
    psF32 x,                           ///< Value at which to evaluate
    psF32 mean,                        ///< Mean for the Gaussian
    psF32 stddev,                      ///< Standard deviation for the Gaussian
    psBool normal                      ///< Indicates whether result should be normalized
);

/** Produce a vector of random numbers from a Gaussian distribution with
 *  the specified mean and sigma
 *
 *  @return psVector*    vector of random numbers
 *
 */
psVector* p_psGaussianDev(
    psF32 mean,                        ///< The mean of the Gaussian
    psF32 sigma,                       ///< The sigma of the Gaussian
    psS32 Npts                         ///< The size of the vector
);

/** Polynomial Type.
 *
 *  Enumeration for Polynomial types.
 */
typedef enum {
    PS_POLYNOMIAL_ORD,                 ///< Ordinary Polynomial
    PS_POLYNOMIAL_CHEB                 ///< Chebyshev Polynomial
} psPolynomialType;

/** One-dimensional polynomial */
typedef struct
{
    psPolynomialType type;             ///< Polynomial type
    psS32 n;                             ///< Number of terms
    psF32 *coeff;                      ///< Coefficients
    psF32 *coeffErr;                   ///< Error in coefficients
    psU8 *mask;                        ///< Coefficient mask
}
psPolynomial1D;

/** Two-dimensional polynomial */
typedef struct
{
    psPolynomialType type;             ///< Polynomial type
    psS32 nX;                            ///< Number of terms in x
    psS32 nY;                            ///< Number of terms in y
    psF32 **coeff;                     ///< Coefficients
    psF32 **coeffErr;                  ///< Error in coefficients
    psU8 **mask;                       ///< Coefficients mask
}
psPolynomial2D;

/** Three-dimensional polynomial */
typedef struct
{
    psPolynomialType type;             ///< Polynomial type
    psS32 nX;                            ///< Number of terms in x
    psS32 nY;                            ///< Number of terms in y
    psS32 nZ;                            ///< Number of terms in z
    psF32 ***coeff;                    ///< Coefficients
    psF32 ***coeffErr;                 ///< Error in coefficients
    psU8 ***mask;                      ///< Coefficients mask
}
psPolynomial3D;

/** Four-dimensional polynomial */
typedef struct
{
    psPolynomialType type;             ///< Polynomial type
    psS32 nW;                            ///< Number of terms in w
    psS32 nX;                            ///< Number of terms in x
    psS32 nY;                            ///< Number of terms in y
    psS32 nZ;                            ///< Number of terms in z
    psF32 ****coeff;                   ///< Coefficients
    psF32 ****coeffErr;                ///< Error in coefficients
    psU8 ****mask;                     ///< Coefficients mask
}
psPolynomial4D;


/** Allocates a psPolynomial1D structure with n terms
 *
 *  @return  psPolynomial1D*    new 1-D polynomial struct
 */
psPolynomial1D* psPolynomial1DAlloc(
    psS32 n,                            ///< Number of terms
    psPolynomialType type              ///< Polynomial Type
);

/** Allocates a 2-D polynomial structure
 *
 *  @return  psPolynomial2D*    new 2-D polynomial struct
 */
psPolynomial2D* psPolynomial2DAlloc(
    psS32 nX,                           ///< Number of terms in x
    psS32 nY,                           ///< Number of terms in y
    psPolynomialType type              ///< Polynomial Type
);

/** Allocates a 3-D polynomial structure
 *
 *  @return  psPolynomial3D*    new 3-D polynomial struct
 */
psPolynomial3D* psPolynomial3DAlloc(
    psS32 nX,                           ///< Number of terms in x
    psS32 nY,                           ///< Number of terms in y
    psS32 nZ,                           ///< Number of terms in z
    psPolynomialType type              ///< Polynomial Type
);

/** Allocates a 4-D polynomial structure
 *
 *  @return  psPolynomial4D*    new 4-D polynomial struct
 */
psPolynomial4D* psPolynomial4DAlloc(
    psS32 nW,                           ///< Number of terms in w
    psS32 nX,                           ///< Number of terms in x
    psS32 nY,                           ///< Number of terms in y
    psS32 nZ,                           ///< Number of terms in z
    psPolynomialType type              ///< Polynomial Type
);

/** Evaluates a 1-D polynomial at specific coordinates.
 *
 *  @return psF32    result of polynomial at given location
 */
psF32 psPolynomial1DEval(
    const psPolynomial1D* myPoly,      ///< Coefficients for the polynomial
    psF32 x                            ///< location at which to evaluate
);

/** Evaluates a 2-D polynomial at specific coordinates.
 *
 *  @return psF32    result of polynomial at given location
 */
psF32 psPolynomial2DEval(
    const psPolynomial2D* myPoly,      ///< Coefficients for the polynomial
    psF32 x,                            ///< x location at which to evaluate
    psF32 y                             ///< y location at which to evaluate
);

/** Evaluates a 3-D polynomial at specific coordinates.
 *
 *  @return psF32    result of polynomial at given location
 */
psF32 psPolynomial3DEval(
    const psPolynomial3D* myPoly,      ///< Coefficients for the polynomial
    psF32 x,                            ///< x location at which to evaluate
    psF32 y,                            ///< y location at which to evaluate
    psF32 z                             ///< z location at which to evaluate
);

/** Evaluates a 4-D polynomial at specific coordinates.
 *
 *  @return psF32    result of polynomial at given location
 */
psF32 psPolynomial4DEval(
    const psPolynomial4D* myPoly,      ///< Coefficients for the polynomial
    psF32 w,                            ///< w location at which to evaluate
    psF32 x,                            ///< x location at which to evaluate
    psF32 y,                            ///< y location at which to evaluate
    psF32 z                             ///< z location at which to evaluate
);

/** Evaluates a 1-D polynomial at specific sets of coordinates
 *  
 *  @return psVector*    results of polynomials at given locations
 */
psVector *psPolynomial1DEvalVector(
    const psPolynomial1D *myPoly,      ///< Coefficients for the polynomial
    const psVector *x                   ///< x locations at which to evaluate
);

/** Evaluates a 2-D polynomial at specific sets of coordinates
 *  
 *  @return psVector*    results of polynomial at given locations
 */
psVector *psPolynomial2DEvalVector(
    const psPolynomial2D *myPoly,      ///< Coefficients for the polynomial
    const psVector *x,                  ///< x locations at which to evaluate
    const psVector *y                   ///< y locations at which to evaluate
);

/** Evaluates a 3-D polynomial at specific sets of coordinates
 *  
 *  @return psVector*    results of polynomial at given locations
 */
psVector *psPolynomial3DEvalVector(
    const psPolynomial3D *myPoly,      ///< Coefficients for the polynomial
    const psVector *x,                  ///< x locations at which to evaluate
    const psVector *y,                  ///< y locations at which to evaluate
    const psVector *z                   ///< z locations at which to evaluate
);

/** Evaluates a 4-D polynomial at specific sets of coordinates
 *  
 *  @return psVector*    results of polynomial at given locations
 */
psVector *psPolynomial4DEvalVector(
    const psPolynomial4D *myPoly,      ///< Coefficients for the polynomial
    const psVector *w,                  ///< w locations at which to evaluate
    const psVector *x,                  ///< x locations at which to evaluate
    const psVector *y,                  ///< y locations at which to evaluate
    const psVector *z                   ///< z locations at which to evaluate
);

/*****************************************************************************/

/* Double-precision polynomials, mainly for use in astrometry */

/** Double-precision one-dimensional polynomial */
typedef struct
{
    psPolynomialType type;             ///< Polynomial type
    psS32 n;                            ///< Number of terms
    psF64 *coeff;                      ///< Coefficients
    psF64 *coeffErr;                   ///< Error in coefficients
    psU8 *mask;                        ///< Coefficient mask
}
psDPolynomial1D;

/** Double-precision two-dimensional polynomial */
typedef struct
{
    psPolynomialType type;             ///< Polynomial type
    psS32 nX;                           ///< Number of terms in x
    psS32 nY;                           ///< Number of terms in y
    psF64 **coeff;                     ///< Coefficients
    psF64 **coeffErr;                  ///< Error in coefficients
    psU8 **mask;                       ///< Coefficients mask
}
psDPolynomial2D;

/** Double-precision three-dimensional polynomial */
typedef struct
{
    psPolynomialType type;             ///< Polynomial type
    psS32 nX;                           ///< Number of terms in x
    psS32 nY;                           ///< Number of terms in y
    psS32 nZ;                           ///< Number of terms in z
    psF64 ***coeff;                    ///< Coefficients
    psF64 ***coeffErr;                 ///< Error in coefficients
    psU8 ***mask;                      ///< Coefficient mask
}
psDPolynomial3D;

/** Double-precision four-dimensional polynomial */
typedef struct
{
    psPolynomialType type;             ///< Polynomial type
    psS32 nW;                           ///< Number of terms in w
    psS32 nX;                           ///< Number of terms in x
    psS32 nY;                           ///< Number of terms in y
    psS32 nZ;                           ///< Number of terms in z
    psF64 ****coeff;                   ///< Coefficients
    psF64 ****coeffErr;                ///< Error in coefficients
    psU8 ****mask;                     ///< Coefficients mask
}
psDPolynomial4D;

/** Allocates a double-precision 1-D polynomial structure with n terms
 *
 *  @return  psPolynomial1D*    new double-precision 1-D polynomial struct
 */
psDPolynomial1D* psDPolynomial1DAlloc(
    psS32 n,                              ///< Number of terms
    psPolynomialType type                ///< Polynomial Type
);

/** Allocates a double-precision 2-D polynomial structure
 *
 *  @return  psPolynomial2D*    new double-precision 2-D polynomial struct
 */
psDPolynomial2D* psDPolynomial2DAlloc(
    psS32 nX,                           ///< Number of terms in x
    psS32 nY,                           ///< Number of terms in y
    psPolynomialType type              ///< Polynomial Type
);

/** Allocates a double-precision 3-D polynomial structure
 *
 *  @return  psPolynomial3D*    new double-precision 3-D polynomial struct
 */
psDPolynomial3D* psDPolynomial3DAlloc(
    psS32 nX,                           ///< Number of terms in x
    psS32 nY,                           ///< Number of terms in y
    psS32 nZ,                           ///< Number of terms in z
    psPolynomialType type              ///< Polynomial Type
);

/** Allocates a double-precision 4-D polynomial structure
 *
 *  @return  psPolynomial4D*    new double-precision 4-D polynomial struct
 */
psDPolynomial4D* psDPolynomial4DAlloc(
    psS32 nW,                           ///< Number of terms in w
    psS32 nX,                           ///< Number of terms in x
    psS32 nY,                           ///< Number of terms in y
    psS32 nZ,                           ///< Number of terms in z
    psPolynomialType type              ///< Polynomial Type
);

/** Evaluates a double-precision 1-D polynomial at specific coordinates.
 *
 *  @return psF32    result of polynomial at given location
 */
psF64 psDPolynomial1DEval(
    const psDPolynomial1D* myPoly,     ///< Coefficients for the polynomial
    psF64 x                             ///< Value at which to evaluate
);

/** Evaluates a double-precision 2-D polynomial at specific coordinates.
 *
 *  @return psF32    result of polynomial at given location
 */
psF64 psDPolynomial2DEval(
    const psDPolynomial2D* myPoly,      ///< Coefficients for the polynomial
    psF64 x,                             ///< Value x at which to evaluate
    psF64 y                              ///< Value y at which to evaluate
);

/** Evaluates a double-precision 3-D polynomial at specific coordinates.
 *
 *  @return psF32    result of polynomial at given location
 */
psF64 psDPolynomial3DEval(
    const psDPolynomial3D* myPoly,     ///< Coefficients for the polynomial
    psF64 x,                            ///< Value x at which to evaluate
    psF64 y,                            ///< Value y at which to evaluate
    psF64 z                             ///< Value z at which to evaluate
);

/** Evaluates a double-precision 4-D polynomial at specific coordinates.
 *
 *  @return psF32    result of polynomial at given location
 */
psF64 psDPolynomial4DEval(
    const psDPolynomial4D* myPoly,     ///< Coefficients for the polynomial
    psF64 w,                            ///< Value w at which to evaluate
    psF64 x,                            ///< Value x at which to evaluate
    psF64 y,                            ///< Value y at which to evaluate
    psF64 z                             ///< Value z at which to evaluate
);

/** Evaluates a double-precision 1-D polynomial at specific sets of coordinates.
 *  
 *  @return psVector*    results of polynomial at given locations
 */
psVector *psDPolynomial1DEvalVector(
    const psDPolynomial1D *myPoly,     ///< Coefficients for the polynomial
    const psVector *x                   ///< x locations at which to evaluate
);

/** Evaluates a double-precision 2-D polynomial at specific sets of coordinates.
 *  
 *  @return psVector*    results of polynomial at given locations
 */
psVector *psDPolynomial2DEvalVector(
    const psDPolynomial2D *myPoly,     ///< Coefficients for the polynomial
    const psVector *x,                  ///< x locations at which to evaluate
    const psVector *y                   ///< y locations at which to evaluate
);

/** Evaluates a double-precision 3-D polynomial at specific sets of coordinates.
 *  
 *  @return psVector*    results of polynomial at given locations
 */
psVector *psDPolynomial3DEvalVector(
    const psDPolynomial3D *myPoly,     ///< Coefficients for the polynomial
    const psVector *x,                  ///< x locations at which to evaluate
    const psVector *y,                  ///< y locations at which to evaluate
    const psVector *z                   ///< z locations at which to evaluate
);

/** Evaluates a double-precision 4-D polynomial at specific sets of coordinates.
 *  
 *  @return psVector*    results of polynomial at given locations
 */
psVector *psDPolynomial4DEvalVector(
    const psDPolynomial4D *myPoly,     ///< Coefficients for the polynomial
    const psVector *w,                  ///< w locations at which to evaluate
    const psVector *x,                  ///< x locations at which to evaluate
    const psVector *y,                  ///< y locations at which to evaluate
    const psVector *z                   ///< z locations at which to evaluate
);

/** One-Dimensional Spline */
typedef struct
{
    psS32 n;                            ///< The number of spline polynomials
    psPolynomial1D **spline;           ///< An array of n pointers to the spline polynomials
    psF32 *p_psDeriv2;                 ///< For cubic splines, the second derivative at each domain point.  Size is n+1.
    psF32 *domains;                    ///< The boundaries between each spline piece.  Size is n+1.
    psVector *knots;                   ///< The boundaries between each spline piece.  Size is n+1.
}
psSpline1D;

/** Allocates a psSpline1D structure
 *  
 *  Allocator for psSpline1D where the bounds are implicitly specified through specifying
 *  min and max values along with the number of splines.
 *
 *  @return psSpline1D*    new 1-D spline struct
 */
psSpline1D *psSpline1DAlloc(psS32 n,             ///< Number of spline polynomials
                            psS32 order,         ///< Order of spline polynomials
                            psF32 min,           ///< Lower boundary value of spline polynomials
                            psF32 max);          ///< Upper boundary value of spline polynomials

/** Allocates a psSpline1D structure
 *  
 *  Allocator for psSpline1D where the bounds are explicitly specified. 
 *
 *  @return psSpline1D*    new 1-D spline struct
 */
psSpline1D *psSpline1DAllocGeneric(const psVector *bounds,   ///< Bounds for spline polynomials
                                   psS32 order);             ///< Order of spline polynomials

/** Evaluates 1-D spline polynomials at a specific coordinate.
 *  
 *  @return psF32    result of spline polynomials evaluated at given location
 */
psF32 psSpline1DEval(
    const psSpline1D *spline,          ///< Coefficients for spline polynomials
    psF32 x                            ///< location at which to evaluate
);

/** Evaluates 1-D spline polynomials at a set of specific coordinates.
 *  
 *  @return psVector*    results of spline polynomials evaluated at given locations
 */
psVector *psSpline1DEvalVector(
    const psSpline1D *spline,          ///< Coefficients of spline polynomials
    const psVector *x                  ///< locations at which to evaluate
);

/** Performs a binary disection on a given vector.
 *  Searches through an array of data for a specified value.
 *  
 *  @return psS32    corresponding index number of specified value
 */
psS32 p_psVectorBinDisect(psVector *bins,        ///< Array of non-decreasing values
                          psScalar *x);          ///< Target value to find

/** Interpolates a series of data points for evaluation at a specific coordinate.  Uses a
 *  Lagrange interpolation method.
 *
 *  @return psScalar*    Lagrange interpolation value at given location
 */
psScalar *p_psVectorInterpolate(psVector *domain,     ///< Domain (x coords) for interpolation
                                psVector *range,      ///< Range (y coords) for interpolation
                                psS32 order,          ///< Order of interpolation function
                                psScalar *x);         ///< Location at which to evaluate

#if 0
psF32 p_psNRSpline1DEval(psSpline1D *spline,
                         const psVector* x,
                         const psVector* y,
                         psF32 X);
#endif // #if 0

/** \} */ // End of MathGroup Functions

#endif // #ifndef PS_FUNCTIONS_H

