test_HistoND.cc
No OneTemporary
Actions

Size

29 KB

Subscribers

None

test_HistoND.cc
View Options

	#include <string>
	#include <sstream>
	#include <functional>

	#include "UnitTest++.h"

	#include "test_utils.hh"

	#include "npstat/stat/HistoNDCdf.hh"
	#include "npstat/stat/arrayStats.hh"
	#include "npstat/stat/histoStats.hh"
	#include "npstat/stat/StatAccumulator.hh"
	#include "npstat/stat/interpolateHistoND.hh"
	#include "npstat/nm/SimpleFunctors.hh"
	#include "npstat/stat/ArrayProjectors.hh"
	#include "npstat/stat/NUHistoAxis.hh"
	#include "npstat/stat/DualHistoAxis.hh"
	#include "npstat/nm/EquidistantSequence.hh"

	#include "rk/geom3.hh"

	#include "geners/CPP11_array.hh"

	using namespace npstat;
	using namespace std;

	namespace {
	template <typename T1, typename T2>
	struct pluseq_left_ptr
	{
	inline T1& operator()(T1*& left, const T2& right) const
	{return *left += right;}
	};

	TEST(HistoAxis)
	{
	const unsigned nsteps = 331;
	const HistoAxis axis(10, 1.0, 5.0);
	const double step = 4.0/nsteps;
	const LinearMapper1d mapper(axis.binNumberMapper());
	const LinearMapper1d mapper2(axis.binNumberMapper(false));
	const double halfbin = axis.binWidth()/2.0;

	for (unsigned i=0; i<nsteps; ++i)
	{
	const double x = 1.0 + step*(i + 0.5);
	const int ibin = axis.binNumber(x);
	CHECK_CLOSE(mapper(x), axis.fltBinNumber(x), 1.0e-14);
	CHECK_CLOSE(mapper2(x), axis.fltBinNumber(x, false), 1.0e-14);
	CHECK_EQUAL(ibin, static_cast<int>(mapper(x)));
	CHECK_EQUAL(ibin, static_cast<int>(mapper2(x + halfbin)));
	}
	}

	TEST(HistoND)
	{
	HistoND<double> h1d(std::vector<HistoAxis>(1, HistoAxis(10,0.0,5.0)));
	HistoND<double> h1d_fcn(h1d);
	HistoND<double*> h1d_ptr(h1d.axes());
	HistoND<double> h1d_base(h1d);
	HistoND<double> h1dcop(HistoAxis(10,0.0,5.0));

	CHECK(h1dcop == h1d);

	for (unsigned i=0; i<h1d_base.nBins(); ++i)
	h1d_ptr.setBin(i, const_cast<double*>(&h1d_base.binContents()(i)));

	double one = 1.0;
	pluseq_left<double,double> fcn;
	pluseq_left_ptr<double,double> fcn2;
	for (unsigned i=0; i<1000; ++i)
	{
	const double x = test_rng()*5.0;
	h1d.fill(x, 1.0);
	h1d_fcn.dispatch(x, one, fcn);
	h1d_ptr.dispatch(x, one, fcn2);
	}
	h1d_fcn.recalculateNFillsFromData();
	h1d_base.recalculateNFillsFromData();
	CHECK(h1d == h1d_fcn);
	CHECK(h1d == h1d_base);
	CPP11_auto_ptr<BinnedDensity1D> dens1 = histoDensity1D(h1d);

	const double xinter = 1.12;
	CHECK_EQUAL(interpolateHistoND(h1d, xinter, 1U),
	interpolateHistoND(h1d, &xinter, 1U, 1U));

	std::vector<HistoAxis> axes;
	axes.push_back(HistoAxis(10, 0.0, 5.0));
	axes.push_back(HistoAxis(15, 0.0, 3.0));

	const double xmin = axes[0].min();
	const double bwx = axes[0].binWidth();
	const unsigned nx = axes[0].nBins();
	const double ymin = axes[1].min();
	const double bwy = axes[1].binWidth();
	const unsigned ny = axes[1].nBins();

	HistoND<unsigned> hu1(axes);
	HistoND<unsigned> hu1_cop(HistoAxis(10, 0.0, 5.0),
	HistoAxis(15, 0.0, 3.0));
	CHECK(hu1 == hu1_cop);

	std::vector<HistoAxis> moreAxes(axes);
	moreAxes.push_back(HistoAxis(20, -1.0, 4.0, "haha"));
	HistoND<unsigned> hu33(moreAxes);
	HistoND<unsigned> hu33_3(HistoAxis(10, 0.0, 5.0),
	HistoAxis(15, 0.0, 3.0),
	HistoAxis(20, -1.0, 4.0, "haha"));
	CHECK(hu33 == hu33_3);

	moreAxes.push_back(HistoAxis(2, -10.0, -1.0, "hoho"));
	HistoND<unsigned> hu44(moreAxes);
	HistoND<unsigned> hu44_4(HistoAxis(10, 0.0, 5.0),
	HistoAxis(15, 0.0, 3.0),
	HistoAxis(20, -1.0, 4.0, "haha"),
	HistoAxis(2, -10.0, -1.0, "hoho"));
	CHECK(hu44 == hu44_4);

	hu1.fill(2.0, 1.0, 1);
	hu1.fill(2.0, 1.0, 1);
	CHECK_EQUAL(2U, hu1.examine(2.0, 1.0));
	CHECK_EQUAL(2U, hu1.closestBin(2.0, 1.0));
	const double dindex[2] = {2.0, 1.0};
	CHECK_EQUAL(2U, hu1.examine(dindex, 2U));
	CHECK_EQUAL(2U, hu1.closestBin(dindex, 2U));
	CHECK_CLOSE(5.0*3.0/150, hu1.binVolume(), 1.e-15);
	CHECK_EQUAL(5.0*3.0, hu1.volume());

	HistoND<unsigned> hu2(hu1);
	CHECK(hu1 == hu2);

	hu2 += hu1;
	hu2 -= hu1;
	CHECK(hu1.isSameData(hu2));

	hu2 *= 7;
	CHECK(!hu1.isSameData(hu2));
	hu2 /= 7;
	CHECK(hu1.isSameData(hu2));

	HistoND<unsigned> hu3(axes);
	double coords[2] = {2.0, 1.0};
	hu3.fill(coords, 2, 1);
	hu3.fill(coords, 2, 1);
	CHECK(hu1 == hu3);

	// HistoND<double> hu4(axes);
	HistoND<float> hu4(axes);
	hu4.fillC(2.0, 1.0, 5.0);

	HistoND<unsigned> hu5(axes);
	for (unsigned ix=0; ix<hu1.axis(0).nBins(); ++ix)
	for (unsigned iy=0; iy<hu1.axis(1).nBins(); ++iy)
	hu5.setBinAt(ix, iy, hu1.binContents().at(ix, iy));
	hu5.recalculateNFillsFromData();
	CHECK(hu1 == hu5);

	double mean[2], hmean[2];
	arrayCoordMean(hu4.binContents(), hu4.boundingBox(), mean, 2U);
	CHECK_CLOSE(2.0, mean[0], 1.0e-15);
	CHECK_CLOSE(1.0, mean[1], 1.0e-15);

	histoMean(hu4, hmean, 2U);
	CHECK_CLOSE(2.0, hmean[0], 1.0e-15);
	CHECK_CLOSE(1.0, hmean[1], 1.0e-15);

	for (unsigned i=0; i<10000; ++i)
	hu4.fill(test_rng()5.0, test_rng()3.0, 1.0);
	CHECK(hu4.nFillsTotal() == 10001);

	CHECK_EQUAL(interpolateHistoND(hu4, dindex, 2U, 1),
	interpolateHistoND(hu4, dindex[0], dindex[1], 1));

	Matrix<double> cov1(2, 2);
	Matrix<double> cov2(2, 2);
	arrayCoordCovariance(hu4.binContents(), hu4.boundingBox(), &cov1);
	histoCovariance(hu4, &cov2);
	for (unsigned i=0; i<2; ++i)
	for (unsigned j=0; j<2; ++j)
	CHECK_CLOSE(cov1[i][j], cov2[i][j], 1.0e-15);
	CPP11_auto_ptr<BinnedDensityND> dens = histoDensityND(hu4);

	ostringstream os;
	CHECK(hu4.classId().write(os));
	CHECK(hu4.write(os));

	HistoND<geom3::Vector3> hvec(axes);
	hvec.fill(1.0, 1.0, geom3::Vector3(1.0, 2.0, 3.0));
	hvec.fill(2.0, 1.5, geom3::Vector3(1.1, 2.3, 4.0));
	const geom3::Vector3& bin = hvec.closestBin(1.0, 1.5);
	bin.x();

	HistoND<string> hs1(axes);
	hs1.fill(2.0, 2.0, "abc");
	hs1.fill(1.0, 1.0, "hook");
	hs1.fill(1.0, 1.0, "cat");
	hs1.fill(100.0, 199.0, "baba");
	gs::write_item(os, hs1, false);

	HistoND<string> hs2(hs1);
	hs2.addToBinContents("gqlosc");

	HistoND<string> hs3(hs1);
	for (unsigned ix=0; ix<nx; ++ix)
	{
	const double x = xmin + bwx*(ix + 0.5);
	for (unsigned iy=0; iy<ny; ++iy)
	{
	const double y = ymin + bwy*(iy + 0.5);
	hs3.fill(x, y, "gqlosc");
	}
	}
	CHECK(hs2 == hs3);

	istringstream is(os.str());
	gs::ClassId hid(is, 1);
	CHECK(hid == hu4.classId());
	HistoND<float>* hu4_restored = HistoND<float>::read(hid, is);
	CHECK(hu4_restored);
	CHECK(hu4 == *hu4_restored);
	delete hu4_restored;

	HistoND<string>* hs1_restored = HistoND<string>::read(
	gs::ClassId::makeId<HistoND<string> >(), is);
	CHECK(hs1_restored);
	CHECK(hs1 == *hs1_restored);
	delete hs1_restored;

	// Check the code which gets the bin centers
	std::vector<CPP11_array<float,2> > binCenters;
	hu1.allBinCenters(&binCenters);
	unsigned ibin = 0;
	for (unsigned ix=0; ix<hu1.axis(0).nBins(); ++ix)
	{
	double cnt[2];
	cnt[0] = hu1.axis(0).binCenter(ix);
	for (unsigned iy=0; iy<hu1.axis(1).nBins(); ++iy, ++ibin)
	{
	cnt[1] = hu1.axis(1).binCenter(iy);
	for (unsigned idim=0; idim<2; ++idim)
	CHECK_CLOSE(cnt[idim], binCenters[ibin][idim], 1.0e-7);
	}
	}

	HistoND<StatAccumulator> hacc(axes);
	hacc.fill(2.0, 1.0, 1);

	HistoND<double> hd(axes);
	for (unsigned i=0; i<1000; ++i)
	hd.fill(test_rng()5.0, test_rng()3.0, 1.0);
	CHECK(hd.nFillsTotal() == 1000);

	ArrayND<double> arrhd(hd.binContents());
	hd.scaleBinContents(hu4.binContents().data(), hu4.nBins());
	arrhd *= hu4.binContents();
	const unsigned long lenb = hu4.nBins();
	for (unsigned long i=0; i<lenb; ++i)
	CHECK_EQUAL(arrhd.data()[i], hd.binContents().data()[i]);

	HistoND<StatAccumulator> hacc2(hd.axes());
	hacc2 += hd;
	hacc2 += hd;
	HistoND<double> hmeanh(hacc2, mem_fun_ref(&StatAccumulator::mean));
	HistoND<float> fmean(hmeanh, SameRef<double>());

	hd.setBinsToConst(1.0f);
	hd.setOverflowsToConst(1U);

	hacc2 /= -2;
	hacc2 *= -2;
	}

	TEST(HistoND_dual)
	{
	HistoND<double,DualHistoAxis> h1d(DualHistoAxis(10,0.0,5.0));
	HistoND<double,DualHistoAxis> h1d_fcn(h1d);
	HistoND<double*,DualHistoAxis> h1d_ptr(h1d.axes());
	HistoND<double,DualHistoAxis> h1d_base(h1d);
	HistoND<double,DualHistoAxis> h1dcop(DualHistoAxis(10,0.0,5.0));

	CHECK(h1dcop == h1d);

	for (unsigned i=0; i<h1d_base.nBins(); ++i)
	h1d_ptr.setBin(i, const_cast<double*>(&h1d_base.binContents()(i)));

	double one = 1.0;
	pluseq_left<double,double> fcn;
	pluseq_left_ptr<double,double> fcn2;
	for (unsigned i=0; i<1000; ++i)
	{
	const double x = test_rng()*5.0;
	h1d.fill(x, 1.0);
	h1d_fcn.dispatch(x, one, fcn);
	h1d_ptr.dispatch(x, one, fcn2);
	}
	h1d_fcn.recalculateNFillsFromData();
	h1d_base.recalculateNFillsFromData();
	CHECK(h1d == h1d_fcn);
	CHECK(h1d == h1d_base);
	CPP11_auto_ptr<BinnedDensity1D> dens1 = histoDensity1D(h1d);

	const unsigned nCdf = 1000;
	std::vector<double> cdf(nCdf);
	std::vector<double> quantiles(nCdf);
	for (unsigned i=0; i<nCdf; ++i)
	{
	const double r = test_rng();
	cdf[i] = r;
	quantiles[i] = dens1->quantile(r);
	}
	histoQuantiles(h1d, 0, &cdf[0], &cdf[0], nCdf);

	for (unsigned i=0; i<nCdf; ++i)
	CHECK_CLOSE(quantiles[i], cdf[i], 1.0e-12);

	const double xinter = 1.12;
	CHECK_EQUAL(interpolateHistoND(h1d, xinter, 1U),
	interpolateHistoND(h1d, &xinter, 1U, 1U));

	std::vector<DualHistoAxis> axes;
	axes.push_back(DualHistoAxis(10, 0.0, 5.0));
	axes.push_back(DualHistoAxis(15, 0.0, 3.0));

	const double xmin = axes[0].min();
	const double bwx = axes[0].binWidth(0);
	const unsigned nx = axes[0].nBins();
	const double ymin = axes[1].min();
	const double bwy = axes[1].binWidth(0);
	const unsigned ny = axes[1].nBins();

	HistoND<unsigned,DualHistoAxis> hu1(axes);
	HistoND<unsigned,DualHistoAxis> hu1_cop(DualHistoAxis(10, 0.0, 5.0),
	DualHistoAxis(15, 0.0, 3.0));
	CHECK(hu1 == hu1_cop);

	std::vector<DualHistoAxis> moreAxes(axes);
	moreAxes.push_back(DualHistoAxis(20, -1.0, 4.0, "haha"));
	HistoND<unsigned,DualHistoAxis> hu33(moreAxes);
	HistoND<unsigned,DualHistoAxis> hu33_3(DualHistoAxis(10, 0.0, 5.0),
	DualHistoAxis(15, 0.0, 3.0),
	DualHistoAxis(20, -1.0, 4.0, "haha"));
	CHECK(hu33 == hu33_3);

	moreAxes.push_back(DualHistoAxis(2, -10.0, -1.0, "hoho"));
	HistoND<unsigned,DualHistoAxis> hu44(moreAxes);
	HistoND<unsigned,DualHistoAxis> hu44_4(DualHistoAxis(10, 0.0, 5.0),
	DualHistoAxis(15, 0.0, 3.0),
	DualHistoAxis(20, -1.0, 4.0, "haha"),
	DualHistoAxis(2, -10.0, -1.0, "hoho"));
	CHECK(hu44 == hu44_4);

	hu1.fill(2.0, 1.0, 1);
	hu1.fill(2.0, 1.0, 1);
	CHECK_EQUAL(2U, hu1.examine(2.0, 1.0));
	CHECK_EQUAL(2U, hu1.closestBin(2.0, 1.0));
	const double dindex[2] = {2.0, 1.0};
	CHECK_EQUAL(2U, hu1.examine(dindex, 2U));
	CHECK_EQUAL(2U, hu1.closestBin(dindex, 2U));
	CHECK_CLOSE(5.0*3.0/150, hu1.binVolume(), 1.e-15);
	CHECK_EQUAL(5.0*3.0, hu1.volume());

	HistoND<unsigned,DualHistoAxis> hu2(hu1);
	CHECK(hu1 == hu2);

	hu2 += hu1;
	hu2 -= hu1;
	CHECK(hu1.isSameData(hu2));

	hu2 *= 7;
	CHECK(!hu1.isSameData(hu2));
	hu2 /= 7;
	CHECK(hu1.isSameData(hu2));

	HistoND<unsigned,DualHistoAxis> hu3(axes);
	double coords[2] = {2.0, 1.0};
	hu3.fill(coords, 2, 1);
	hu3.fill(coords, 2, 1);
	CHECK(hu1 == hu3);

	// HistoND<double> hu4(axes);
	HistoND<float,DualHistoAxis> hu4(axes);
	hu4.fill(0.75, 0.3, 5.0);

	HistoND<unsigned,DualHistoAxis> hu5(axes);
	for (unsigned ix=0; ix<hu1.axis(0).nBins(); ++ix)
	for (unsigned iy=0; iy<hu1.axis(1).nBins(); ++iy)
	hu5.setBinAt(ix, iy, hu1.binContents().at(ix, iy));
	hu5.recalculateNFillsFromData();
	CHECK(hu1 == hu5);

	double mean[2], hmean[2];
	arrayCoordMean(hu4.binContents(), hu4.boundingBox(), mean, 2U);
	CHECK_CLOSE(0.75, mean[0], 1.0e-15);
	CHECK_CLOSE(0.3, mean[1], 1.0e-15);

	histoMean(hu4, hmean, 2U);
	CHECK_CLOSE(0.75, hmean[0], 1.0e-15);
	CHECK_CLOSE(0.3, hmean[1], 1.0e-15);

	for (unsigned i=0; i<10000; ++i)
	hu4.fill(test_rng()5.0, test_rng()3.0, 1.0);
	CHECK(hu4.nFillsTotal() == 10001);

	CHECK_EQUAL(interpolateHistoND(hu4, dindex, 2U, 1),
	interpolateHistoND(hu4, dindex[0], dindex[1], 1));

	Matrix<double> cov1(2, 2);
	Matrix<double> cov2(2, 2);
	arrayCoordCovariance(hu4.binContents(), hu4.boundingBox(), &cov1);
	histoCovariance(hu4, &cov2);
	for (unsigned i=0; i<2; ++i)
	for (unsigned j=0; j<2; ++j)
	CHECK_CLOSE(cov1[i][j], cov2[i][j], 1.0e-15);
	CPP11_auto_ptr<BinnedDensityND> dens = histoDensityND(hu4);

	ostringstream os;
	CHECK(hu4.classId().write(os));
	CHECK(hu4.write(os));

	HistoND<geom3::Vector3,DualHistoAxis> hvec(axes);
	hvec.fill(1.0, 1.0, geom3::Vector3(1.0, 2.0, 3.0));
	hvec.fill(2.0, 1.5, geom3::Vector3(1.1, 2.3, 4.0));
	const geom3::Vector3& bin = hvec.closestBin(1.0, 1.5);
	bin.x();

	HistoND<string,DualHistoAxis> hs1(axes);
	hs1.fill(2.0, 2.0, "abc");
	hs1.fill(1.0, 1.0, "hook");
	hs1.fill(1.0, 1.0, "cat");
	hs1.fill(100.0, 199.0, "baba");
	gs::write_item(os, hs1, false);

	HistoND<string,DualHistoAxis> hs2(hs1);
	hs2.addToBinContents("gqlosc");

	HistoND<string,DualHistoAxis> hs3(hs1);
	for (unsigned ix=0; ix<nx; ++ix)
	{
	const double x = xmin + bwx*(ix + 0.5);
	for (unsigned iy=0; iy<ny; ++iy)
	{
	const double y = ymin + bwy*(iy + 0.5);
	hs3.fill(x, y, "gqlosc");
	}
	}
	CHECK(hs2 == hs3);

	istringstream is(os.str());
	gs::ClassId hid(is, 1);
	CHECK(hid == hu4.classId());
	HistoND<float,DualHistoAxis>* hu4_restored = HistoND<float,DualHistoAxis>::read(hid, is);
	CHECK(hu4_restored);
	CHECK(hu4 == *hu4_restored);
	delete hu4_restored;

	HistoND<string,DualHistoAxis>* hs1_restored = HistoND<string,DualHistoAxis>::read(
	gs::ClassId::makeId<HistoND<string,DualHistoAxis> >(), is);
	CHECK(hs1_restored);
	CHECK(hs1 == *hs1_restored);
	delete hs1_restored;

	// Check the code which gets the bin centers
	std::vector<CPP11_array<float,2> > binCenters;
	hu1.allBinCenters(&binCenters);
	unsigned ibin = 0;
	for (unsigned ix=0; ix<hu1.axis(0).nBins(); ++ix)
	{
	double cnt[2];
	cnt[0] = hu1.axis(0).binCenter(ix);
	for (unsigned iy=0; iy<hu1.axis(1).nBins(); ++iy, ++ibin)
	{
	cnt[1] = hu1.axis(1).binCenter(iy);
	for (unsigned idim=0; idim<2; ++idim)
	CHECK_CLOSE(cnt[idim], binCenters[ibin][idim], 1.0e-7);
	}
	}

	HistoND<StatAccumulator,DualHistoAxis> hacc(axes);
	hacc.fill(2.0, 1.0, 1);

	HistoND<double,DualHistoAxis> hd(axes);
	for (unsigned i=0; i<1000; ++i)
	hd.fill(test_rng()5.0, test_rng()3.0, 1.0);
	CHECK(hd.nFillsTotal() == 1000);

	ArrayND<double> arrhd(hd.binContents());
	hd.scaleBinContents(hu4.binContents().data(), hu4.nBins());
	arrhd *= hu4.binContents();
	const unsigned long lenb = hu4.nBins();
	for (unsigned long i=0; i<lenb; ++i)
	CHECK_EQUAL(arrhd.data()[i], hd.binContents().data()[i]);

	HistoND<StatAccumulator,DualHistoAxis> hacc2(hd.axes());
	hacc2 += hd;
	hacc2 += hd;
	HistoND<double,DualHistoAxis> hmeanh(hacc2, mem_fun_ref(&StatAccumulator::mean));
	HistoND<float,DualHistoAxis> fmean(hmeanh, SameRef<double>());

	hd.setBinsToConst(1.0f);
	hd.setOverflowsToConst(1U);

	hacc2 /= -2;
	hacc2 *= -2;
	}

	TEST(HistoND_2)
	{
	HistoND<double> h1(std::vector<HistoAxis>(1, HistoAxis(10,0.0,5.0)));
	EquidistantInLinearSpace bins(0.0,5.0,11);
	HistoND<double,NUHistoAxis> h2(std::vector<NUHistoAxis>(1, bins));

	for (unsigned i=0; i<10000; ++i)
	{
	const double x = test_rng()*12 - 1.0;
	h1.fill(x, 1.0);
	h2.fill(x, 1.0);
	}

	CHECK_CLOSE(h1.integral(), 0.5*h1.nFillsInRange(), 1.0e-12);
	CHECK_CLOSE(h2.integral(), 0.5*h2.nFillsInRange(), 1.0e-12);
	CHECK_EQUAL(h1.nBins(), h2.nBins());
	const ArrayND<double>& a1 = h1.binContents();
	const ArrayND<double>& a2 = h2.binContents();
	for (unsigned i=0; i<h2.nBins(); ++i)
	{
	CHECK(a1(i) == a2(i));
	CHECK_CLOSE(h1.axis(0).leftBinEdge(i),
	h2.axis(0).leftBinEdge(i), 1.0e-14);
	CHECK_CLOSE(h1.axis(0).rightBinEdge(i),
	h2.axis(0).rightBinEdge(i), 1.0e-14);
	CHECK_CLOSE(h1.axis(0).binWidth(i),
	h2.axis(0).binWidth(i), 1.0e-14);
	CHECK_CLOSE(h1.axis(0).binCenter(i),
	h2.axis(0).binCenter(i), 1.0e-14);
	}
	const ArrayND<double>& o1 = h1.overflows();
	const ArrayND<double>& o2 = h2.overflows();
	for (unsigned i=0; i<o1.length(); ++i)
	CHECK(o1(i) == o2(i));

	ostringstream os;
	CHECK(h2.classId().write(os));
	CHECK(h2.write(os));

	istringstream is(os.str());
	gs::ClassId hid(is, 1);
	CHECK(hid == h2.classId());
	HistoND<double,NUHistoAxis>* h2_restored =
	HistoND<double,NUHistoAxis>::read(hid, is);
	CHECK(h2_restored);
	CHECK(h2 == *h2_restored);
	delete h2_restored;
	}

	TEST(HistoND_2_dual)
	{
	HistoND<double> h1(std::vector<HistoAxis>(1, HistoAxis(10,0.0,5.0)));
	EquidistantInLinearSpace bins(0.0,5.0,11);
	HistoND<double,DualHistoAxis> h2(std::vector<DualHistoAxis>(1, bins));

	for (unsigned i=0; i<10000; ++i)
	{
	const double x = test_rng()*12 - 1.0;
	h1.fill(x, 1.0);
	h2.fill(x, 1.0);
	}

	CHECK_CLOSE(h1.integral(), 0.5*h1.nFillsInRange(), 1.0e-12);
	CHECK_CLOSE(h2.integral(), 0.5*h2.nFillsInRange(), 1.0e-12);
	CHECK_EQUAL(h1.nBins(), h2.nBins());
	const ArrayND<double>& a1 = h1.binContents();
	const ArrayND<double>& a2 = h2.binContents();
	for (unsigned i=0; i<h2.nBins(); ++i)
	{
	CHECK(a1(i) == a2(i));
	CHECK_CLOSE(h1.axis(0).leftBinEdge(i),
	h2.axis(0).leftBinEdge(i), 1.0e-14);
	CHECK_CLOSE(h1.axis(0).rightBinEdge(i),
	h2.axis(0).rightBinEdge(i), 1.0e-14);
	CHECK_CLOSE(h1.axis(0).binWidth(i),
	h2.axis(0).binWidth(i), 1.0e-14);
	CHECK_CLOSE(h1.axis(0).binCenter(i),
	h2.axis(0).binCenter(i), 1.0e-14);
	}
	const ArrayND<double>& o1 = h1.overflows();
	const ArrayND<double>& o2 = h2.overflows();
	for (unsigned i=0; i<o1.length(); ++i)
	CHECK(o1(i) == o2(i));

	ostringstream os;
	CHECK(h2.classId().write(os));
	CHECK(h2.write(os));

	istringstream is(os.str());
	gs::ClassId hid(is, 1);
	CHECK(hid == h2.classId());
	HistoND<double,DualHistoAxis>* h2_restored =
	HistoND<double,DualHistoAxis>::read(hid, is);
	CHECK(h2_restored);
	CHECK(h2 == *h2_restored);
	delete h2_restored;
	}

	TEST(HistoND_overflows)
	{
	const double eps = 1.0e-10;

	HistoND<double,HistoAxis> h2(HistoAxis(10, 0.25, 0.75),
	HistoAxis(7, 0.1, 0.5));
	double cnt0[3] = {0.,};
	double cnt1[3] = {0.,};
	for (unsigned i=0; i<1000; ++i)
	{
	const double x1 = test_rng();
	const double x2 = test_rng();
	const double w = test_rng();
	h2.fill(x1, x2, w);
	if (x1 < 0.25 \|\| x1 >= 0.75 \|\| x2 < 0.1 \|\| x2 >= 0.5)
	{
	if (x1 < 0.25)
	cnt0[0] += w;
	else if (x1 >= 0.75)
	cnt0[2] += w;
	else
	cnt0[1] += w;
	if (x2 < 0.1)
	cnt1[0] += w;
	else if (x2 >= 0.5)
	cnt1[2] += w;
	else
	cnt1[1] += w;
	}
	}
	CHECK_CLOSE(cnt0[0], h2.underflowWeight(0), eps);
	CHECK_CLOSE(cnt0[1], h2.inRangeOverWeight(0), eps);
	CHECK_CLOSE(cnt0[2], h2.overflowWeight(0), eps);
	CHECK_CLOSE(cnt1[0], h2.underflowWeight(1), eps);
	CHECK_CLOSE(cnt1[1], h2.inRangeOverWeight(1), eps);
	CHECK_CLOSE(cnt1[2], h2.overflowWeight(1), eps);
	}

	TEST(accumulateBinsInBox)
	{
	const double eps = 1.0e-14;
	const double eps2 = 1.0e-10;

	std::vector<HistoAxis> axes;
	axes.push_back(HistoAxis(2, 0.0, 2.0));
	axes.push_back(HistoAxis(2, 0.0, 2.0));

	HistoND<double> hd(axes);
	hd.setBin(0U, 0U, 1.0);
	hd.setBin(0U, 1U, 2.0);
	hd.setBin(1U, 0U, 4.0);
	hd.setBin(1U, 1U, 8.0);

	BoxND<double> box;
	box.push_back(Interval<double>(0.3, 1.3));
	box.push_back(Interval<double>(0.6, 1.6));

	const double expect = 0.70.41 + 0.70.62 + 0.30.44 + 0.30.68;

	double sum0 = 0.0;
	for (unsigned i=0; i<2; ++i)
	{
	for (unsigned j=0; j<2; ++j)
	{
	BoxND<double> bin;
	bin.push_back(hd.axis(0).binInterval(i));
	bin.push_back(hd.axis(1).binInterval(j));
	const double over = bin.overlapFraction(box);
	sum0 += over*hd.binContents()(i, j);
	}
	}
	CHECK_CLOSE(expect, sum0, eps);

	long double sum = 0.0;
	hd.accumulateBinsInBox(box, &sum);
	CHECK_CLOSE(expect, static_cast<double>(sum), eps);

	double center[2];
	center[0] = hd.axis(0).fltBinNumber(box[0].midpoint(), false);
	center[1] = hd.axis(1).fltBinNumber(box[1].midpoint(), false);
	const double v1 = hd.binContents().interpolate1(center, 2);
	CHECK_CLOSE(expect, v1, eps);

	double scales[2] = {10.0, 10.0};
	HistoNDCdf hcdf(hd, scales, 2);
	const double nall = hd.binContents().sum<long double>();
	const double coverage = hcdf.boxCoverage(box);
	CHECK_CLOSE(expect, nall*coverage, eps);

	center[0] = box[0].midpoint();
	center[1] = box[1].midpoint();
	BoxND<double> cobox;
	hcdf.coveringBox(coverage, center, 2, &cobox);
	for (unsigned i=0; i<2; ++i)
	{
	CHECK_CLOSE(box[i].min(), cobox[i].min(), eps2);
	CHECK_CLOSE(box[i].max(), cobox[i].max(), eps2);
	}
	}

	TEST(HistoNDCdf)
	{
	const unsigned npoints = 1000;
	const double eps = 1.0e-11;

	for (unsigned dim = 1; dim < 4; ++dim)
	{
	std::vector<HistoAxis> axes;
	axes.push_back(HistoAxis(5, -2.0, 5.0));
	if (dim > 1)
	axes.push_back(HistoAxis(10, -4.0, 8.0));
	if (dim > 2)
	axes.push_back(HistoAxis(7, -1.0, 2.0));

	HistoND<double> hd(axes);

	for (unsigned icycle=0; icycle<100; ++icycle)
	{
	hd.clear();
	for (unsigned i=0; i<npoints; ++i)
	{
	double values[3];
	values[0] = test_rng()*7.0-2.0;
	values[1] = test_rng()*12.0-4.0;
	values[2] = test_rng()*3.0-1.0;
	hd.fill(values, dim, 1.0);
	}
	const double nInside = hd.binContents().sum<long double>();
	CHECK(nInside > 0.0);

	double scales[3];
	scales[0] = 3.5*test_rng() + 0.1;
	scales[1] = 6.0*test_rng() + 0.1;
	scales[2] = 1.5*test_rng() + 0.1;

	HistoNDCdf hcdf(hd, scales, dim);

	BoxND<double> box(BoxND<double>::sizeTwoBox(dim));
	double shift[3];
	shift[0] = test_rng()*7.0-2.0;
	shift[1] = test_rng()*12.0-4.0;
	shift[2] = test_rng()*3.0-1.0;
	box.shift(shift, dim);
	box.expand(scales, dim);

	long double sum = 0.0L;
	hd.accumulateBinsInBox(box, &sum);
	const double expectedCover = sum/nInside;

	const double cover1 = hcdf.boxCoverage(box);

	CHECK_CLOSE(expectedCover, cover1, eps);

	BoxND<double> cobox;
	hcdf.coveringBox(cover1, shift, dim, &cobox);
	for (unsigned i=0; i<dim; ++i)
	{
	CHECK_CLOSE(box[i].midpoint(),cobox[i].midpoint(),1.0e-15);
	CHECK_CLOSE(box[i].min(), cobox[i].min(), eps);
	CHECK_CLOSE(box[i].max(), cobox[i].max(), eps);
	}
	}
	}
	}

	TEST(HistoND_addToProjection)
	{
	const double eps = 1.0e-12;
	const unsigned npoints = 1000;

	std::vector<HistoAxis> axes;
	axes.push_back(HistoAxis(5, 0, 10, "X"));
	axes.push_back(HistoAxis(10, -5, 5, "Y"));

	HistoND<unsigned> expected(axes, "Some Title", "V");

	axes.push_back(HistoAxis(7, 0, 7, "Z"));
	HistoND<unsigned> hd(axes, "Dummy", "V");

	HistoND<unsigned> marg0(axes[0]);
	HistoND<unsigned> marg1(axes[1]);
	HistoND<unsigned> marg2(axes[2]);

	for (unsigned i=0; i<npoints; ++i)
	{
	double values[3];
	values[0] = test_rng()*10.0;
	values[1] = test_rng()*10.0-5.0;
	values[2] = test_rng()*7.0;
	hd.fill(values, 3, 1U);
	expected.fill(values, 2, 1U);
	marg0.fill(values[0], 1U);
	marg1.fill(values[1], 1U);
	marg2.fill(values[2], 1U);
	}
	CHECK(hd.nFillsInRange() == npoints);
	CHECK(expected.nFillsInRange() == npoints);

	unsigned exclude = 2;
	HistoND<unsigned> proj(hd, &exclude, 1U, "Some Title");
	ArraySumProjector<unsigned, unsigned, unsigned> summer;
	hd.addToProjection(&proj, summer, &exclude, 1U);
	proj.recalculateNFillsFromData();

	CHECK(expected == proj);

	const double qvalues[3] = {0.25, 0.5, 0.75};
	double quantiles[3], quantComp[3];

	histoQuantiles(hd, 0, qvalues, quantiles, 3);
	histoQuantiles(marg0, 0, qvalues, quantComp, 3);
	for (unsigned i=0; i<3; ++i)
	CHECK_CLOSE(quantComp[i], quantiles[i], eps);

	histoQuantiles(hd, 1, qvalues, quantiles, 3);
	histoQuantiles(marg1, 0, qvalues, quantComp, 3);
	for (unsigned i=0; i<3; ++i)
	CHECK_CLOSE(quantComp[i], quantiles[i], eps);

	histoQuantiles(hd, 2, qvalues, quantiles, 3);
	histoQuantiles(marg2, 0, qvalues, quantComp, 3);
	for (unsigned i=0; i<3; ++i)
	CHECK_CLOSE(quantComp[i], quantiles[i], eps);
	}

	TEST(histoQuantiles)
	{
	HistoND<int> h1(HistoAxis(100,0.0,5.0));
	for (unsigned i=0; i<10; ++i)
	{
	const double x = test_rng()*5.0;
	h1.fill(x, 1.0);
	}
	CPP11_auto_ptr<BinnedDensity1D> dens1 = histoDensity1D(h1);

	const unsigned nCdf = 1000;
	std::vector<double> cdf(nCdf);
	std::vector<double> quantiles(nCdf);
	for (unsigned i=0; i<nCdf; ++i)
	{
	const double r = test_rng();
	cdf[i] = r;
	quantiles[i] = dens1->quantile(r);
	}
	histoQuantiles(h1, 0, &cdf[0], &cdf[0], nCdf);

	for (unsigned i=0; i<nCdf; ++i)
	CHECK_CLOSE(quantiles[i], cdf[i], 1.0e-12);
	}
	}

File Metadata

Mime Type: text/x-c
Expires: Wed, May 14, 11:27 AM (16 h, 37 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 5085035
Default Alt Text: test_HistoND.cc (29 KB)

test_HistoND.ccNo OneTemporaryActions

test_HistoND.ccView Options

File Metadata

Event Timeline

test_HistoND.cc
No OneTemporary
Actions

test_HistoND.cc
View Options