using namespace std;
-template <typename obstype>
+template <typename meastype, typename restype>
class obstat
{
public:
- void addMeas(const obstype& val);
- void addMeas(obstype val[], int valsize);
+ void addMeas(const meastype& val);
+ void addMeas(meastype val[], int valsize);
- void computeEasy(const string& compid="one", const int& ival=0) { mean(compid, &measurements, ival); }
+ int computeEasy(const int& ival=0) { return mean(&measurements, ival); }
- void computeJack(const string& compid, double (*func)(vector< vector <obstype> > *vals, void *para)=NULL, void *para=NULL);
- void computeJack(double (*func)(vector< vector <obstype> > *vals, void *para)=NULL, void *para=NULL) { computeJack("one",func,para); }
+ int computeJack(restype (*func)(vector< vector <meastype> > *vals, void *para), void *para=NULL);
+ 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);
- double getMean(const string& compid="one") { return computations[compid].val; }
- double getErr(const string& compid="one") { return computations[compid].err; }
+ restype getMean(int compid) { return computations[compid].val; }
+ restype getErr(int compid) { return computations[compid].err; }
void reset();
private:
struct result{
- double val;
- double err;
+ restype val;
+ restype err;
};
- vector< vector<obstype> > measurements;
- map<string,result> computations;
+ vector< vector<meastype> > measurements;
+ vector<result> computations;
- void mean(const string& compid, vector< vector <double> > *meas, const int& ival);
- void mean(const string& compid, vector< vector <int> > *meas, const int& ival);
+ int mean(vector< vector <meastype> > *meas, const int& ival);
};
-
-template <typename obstype>
-void obstat<obstype>::reset()
+template <typename meastype, typename restype>
+void obstat<meastype,restype>::reset()
{
computations.clear();
measurements.clear();
}
-template <typename obstype>
-void obstat<obstype>::addMeas(const obstype& val)
+template <typename meastype, typename restype>
+void obstat<meastype,restype>::addMeas(const meastype& val)
{
- measurements.push_back( vector<obstype>(1,val) );
+ measurements.push_back( vector<meastype>(1,val) );
}
-template <typename obstype>
-void obstat<obstype>::addMeas(obstype val[], int valsize)
+template <typename meastype, typename restype>
+void obstat<meastype,restype>::addMeas(meastype val[], int valsize)
{
- vector<obstype> tmpvec;
- for(int i=0; i<valsize; i++) tmpvec.push_back(val[i]);
+ vector<meastype> tmpvec;
+ for (int i=0; i<valsize; i++)
+ tmpvec.push_back(val[i]);
measurements.push_back(tmpvec);
}
-template <typename obstype>
-void obstat<obstype>::mean(const string& compid, vector< vector<double> > *meas, const int& ival)
+template <typename meastype, typename restype>
+int obstat<meastype,restype>::mean(vector< vector<meastype> > *meas, const int& ival)
{
- computations[compid].val = 0;
- computations[compid].err = 0;
+ result meanres;
+
+ meanres.val = 0;
+ meanres.err = 0;
int nmeas = meas->size();
- for(vector< vector<double> >::iterator measIt = meas->begin(); measIt != meas->end(); ++measIt)
- computations[compid].val += (*measIt)[ival];
- computations[compid].val /= nmeas;
+ for (typename vector< vector<meastype> >::iterator measIt = meas->begin(); measIt != meas->end(); ++measIt)
+ meanres.val += (*measIt)[ival];
+ meanres.val /= nmeas;
- for(vector< vector<double> >::iterator measIt = meas->begin(); measIt != meas->end(); ++measIt)
- computations[compid].err += pow( (*measIt)[ival] - computations[compid].val, 2 );
- computations[compid].err = sqrt( computations[compid].err ) / 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;
-template <typename obstype>
-void obstat<obstype>::mean(const string& compid, vector< vector<int> > *meas, const int& ival)
-{
- computations[compid].val = 0;
- computations[compid].err = 0;
-
- int nmeas = meas->size();
+ computations.push_back(meanres);
- for(vector< vector<int> >::iterator measIt = meas->begin(); measIt != meas->end(); ++measIt)
- computations[compid].val += (*measIt)[ival];
- computations[compid].val /= nmeas;
-
- for(vector< vector<int> >::iterator measIt = meas->begin(); measIt != meas->end(); ++measIt)
- computations[compid].err += pow( (*measIt)[ival] - computations[compid].val, 2 );
- computations[compid].err = sqrt( computations[compid].err ) / nmeas;
-}
-
-template <typename obstype>
-void obstat<obstype>::computeJack(const string& compid, double (*func)(vector< vector<obstype> > *vals, void *para), void *para)
+ return computations.size()-1;
+}
+
+template <typename meastype, typename restype>
+int obstat<meastype,restype>::computeJack(restype (*func)(vector< vector<meastype> > *vals, void *para), void *para)
{
int nmeas=measurements.size();
- double manymeans[nmeas];
+ restype *manymeans = new restype[nmeas];
+ result jackres;
- computations[compid].val = 0;
- computations[compid].err = 0;
+ jackres.val = 0;
+ jackres.err = 0;
int imeas=0;
- for(typename vector< vector<obstype> >::iterator removedIt = measurements.begin(); removedIt != measurements.end(); ++removedIt, imeas++)
+ for(typename vector< vector<meastype> >::iterator removedIt = measurements.begin(); removedIt != measurements.end(); ++removedIt, imeas++)
{
- vector<obstype> removed = *removedIt;
+ vector<meastype> removed = *removedIt;
*removedIt = measurements.back();
measurements.pop_back();
manymeans[imeas] = func(&measurements, para);
- computations[compid].val += manymeans[imeas];
+ jackres.val += manymeans[imeas];
measurements.push_back( *removedIt );
*removedIt = removed;
}
- computations[compid].val /= nmeas;
+ 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;
+}
+
+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++)
- computations[compid].err += pow( manymeans[imean] - computations[compid].val, 2 );
- computations[compid].err *= (double)(nmeas-1)/nmeas;
- computations[compid].err = sqrt(computations[compid].err);
+ 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
projects / phys / latlib.git / blobdiff