#if !defined(PS_STATS_H)
#define PS_STATS_H

/** \file psStats.h
 *  \brief basic statistical operations
 *  \ingroup MathGroup
 */

/** statistics which may be calculated */
typedef enum {
    PS_STAT_SAMPLE_MEAN           = 0x000001,
    PS_STAT_SAMPLE_MEDIAN         = 0x000002,
    PS_STAT_SAMPLE_STDEV          = 0x000004,
    PS_STAT_SAMPLE_QUARTILE       = 0x000008, 
    PS_STAT_ROBUST_MEAN           = 0x000010,
    PS_STAT_ROBUST_MEDIAN         = 0x000020,
    PS_STAT_ROBUST_MODE           = 0x000040,
    PS_STAT_ROBUST_STDEV          = 0x000080,
    PS_STAT_ROBUST_QUARTILE       = 0x000100, 
    PS_STAT_CLIPPED_MEAN          = 0x000200,
    PS_STAT_CLIPPED_STDEV         = 0x000400,
    PS_STAT_MAX                   = 0x000800,     
    PS_STAT_MIN                   = 0x001000,
    PS_STAT_USE_RANGE             = 0x002000,
    PS_STAT_USE_BINSIZE           = 0x004000,
} psStatsOptions;

/** generic statistics structure */
typedef struct {
    double sampleMean;                  ///< formal mean of sample
    double sampleMedian;                ///< formal median of sample
    double sampleStdev;                 ///< standard deviation of sample
    double sampleUQ;                    ///< upper quartile of sample
    double sampleLQ;                    ///< lower quartile of sample
    double robustMean;                  ///< robust mean of data
    double robustMedian;                ///< robust median of data
    double robustMode;                  ///< Robust mode of data
    double robustStdev;                 ///< robust standard deviation of data
    double robustUQ;                    ///< robust upper quartile
    double robustLQ;                    ///< robust lower quartile
    int    robustN50;                   ///< Number of points UQ-LQ
    int    robustNfit;                  ///< Number of points in Gauss. fit
    double clippedMean;                 ///< Nsigma clipped mean
    double clippedStdev;                ///< standard deviation after clipping
    int    clippedNvalues;              ///< number of data points used for clipped mean
    double clipSigma;                   ///< Nsigma used for clipping; user input
    int    clipIter;                    ///< Number of clipping iterations; user input
    double min;                         ///< minimum data value in data; input range
    double max;                         ///< maximum data value in data; input range
    double binsize;		        ///< binsize for robust fit (input/output)
    psStatsOptions options;             ///< bitmask of calculated values
} psStats;


/** Do Statistics on a vector.  Returns a status value. \ingroup MathGroup */
psStats *
psVectorStats(psStats *stats,		///< stats structure defines stats to be calculated, how and result
	      const psVector *in,	///< Array to be analysed
	      const psVector *mask,	///< Ignore elements where (maskArray & maskVal) != 0.  May be NULL
	      unsigned int maskVal)	///< Only mask elements with one of these bits set in maskArray
;

/** Constructor */
psStats *
psStatsAlloc(psStatsOptions options)	///< Statistics to measure
;

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

/** Histograms  */
typedef struct {
    const psVector *bounds;		///< Bounds for the bins
    psVector *nums;			///< Number in each of the bins
    int minNum, maxNum;			///< Number below the minimum and above the maximum
    bool uniform;			///< Is it a uniform distribution?
} psHistogram;

/** Constructor \ingroup MathGroup */
psHistogram *
psHistogramAlloc(float lower,		///< Lower limit for the bins
		 float upper,		///< Upper limit for the bins
		 int n)			///< Number of the bins
;

/** Generic constructor \ingroup MathGroup */
psHistogram *
psHistogramAllocGeneric(const psVector *bounds) ///< Bounds for the bins
;

/** Calculate a histogram \ingroup MathGroup **/
psHistogram *
psVectorHistogram(psHistogram *out, 	///< Input and output histogram
		  const psVector *in	///< Array to analyse
		  const psVector *mask,	///< Ignore elements where (maskArray & maskVal) != 0.  May be NULL
		  unsigned int maskVal	///< Only mask elements with one of these bits set in maskArray
    );

#endif

