#ifndef OBSTAT_HPP
#define OBSTAT_HPP

#include <vector>
#include <string>
#include <map>
#include <math.h>

using namespace std;

template <typename meastype, typename restype>
class obstat
{ 
public:
  void addMeas(const meastype& val);
  void addMeas(meastype val[], int valsize);

  int computeEasy(const int& ival=0) { return mean(&measurements, ival); }

  int computeJack(void (*preMeasFunc)(vector< vector <meastype> > *allVals, vector <meastype> *preCalculated, void *para), 
		  restype (*measFunc)(vector <meastype> *preCalculated, vector <meastype> *excludedmeas, int nmeas, void *para), void *para=NULL);

  restype getMean(int compid) { return computations[compid].val; }
  restype getErr(int compid) { return computations[compid].err; }

  void reset();

private:
  struct result{
    restype val;
    restype err;
  };

  vector< vector<meastype> > measurements;
  vector<result> computations;
  
  int mean(vector< vector <meastype> > *meas, const int& ival);
};

template <typename meastype, typename restype>
void obstat<meastype,restype>::reset()
{
  computations.clear();
  measurements.clear();
}

template <typename meastype, typename restype>
void obstat<meastype,restype>::addMeas(const meastype& val)
{
  measurements.push_back( vector<meastype>(1,val) );
}

template <typename meastype, typename restype>
void obstat<meastype,restype>::addMeas(meastype val[], int valsize)
{
  vector<meastype> tmpvec;
  for (int i=0; i<valsize; i++) 
    tmpvec.push_back(val[i]);
  measurements.push_back(tmpvec);
}

template <typename meastype, typename restype>
int obstat<meastype,restype>::mean(vector< vector<meastype> > *meas, const int& ival)
{
  result meanres;

  meanres.val = 0;
  meanres.err = 0;
  
  int nmeas = meas->size();

  for (typename vector< vector<meastype> >::iterator measIt = meas->begin(); measIt != meas->end(); ++measIt)
    meanres.val += (*measIt)[ival];
  meanres.val /= nmeas;
  
  for (typename vector< vector<meastype> >::iterator measIt = meas->begin(); measIt != meas->end(); ++measIt)
    meanres.err += pow( (*measIt)[ival] - meanres.val, 2 );
  meanres.err = sqrt( meanres.err ) / (restype)nmeas;	 

  computations.push_back(meanres);

  return computations.size()-1;
}

template <typename meastype, typename restype>
int obstat<meastype,restype>::computeJack(void (*preMeasFunc)(vector< vector <meastype> > *allVals, vector<meastype> *preCalculated, void *para), 
					  restype (*measFunc)(vector<meastype> *preCalculated, vector<meastype> *excludedmeas, int nmeas, void *para),
					  void *para) {
  int nmeas=measurements.size();
  restype *manymeans = new restype[nmeas];
  result jackres;
  vector<meastype> preCalculated;

  jackres.val = 0;
  jackres.err = 0;

  preMeasFunc(&measurements, &preCalculated, para);

  int imeas=0;
  for(typename vector< vector<meastype> >::iterator removedIt = measurements.begin(); removedIt != measurements.end(); ++removedIt, imeas++) {
    manymeans[imeas] = measFunc(&preCalculated, &(*removedIt), measurements.size(), para);
    jackres.val += manymeans[imeas]; 
  }
  jackres.val /= nmeas;
  
  for(int imean=0; imean<nmeas; imean++)
    jackres.err += pow( manymeans[imean] - jackres.val, 2 );
  jackres.err *= (double)(nmeas-1)/nmeas;
  jackres.err = sqrt(jackres.err);

  computations.push_back(jackres);

  delete [] manymeans;

  return computations.size()-1;
}

#endif