Index: trunk/tests/test_HistoND.cc
===================================================================
--- trunk/tests/test_HistoND.cc	(revision 615)
+++ trunk/tests/test_HistoND.cc	(revision 616)
@@ -1,791 +1,816 @@
 #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.inPlaceMul(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.inPlaceMul(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(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.7*0.4*1 + 0.7*0.6*2 + 0.3*0.4*4 + 0.3*0.6*8;
 
         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);
+    }
 }