diff --git a/build/CMakeLists.txt b/build/CMakeLists.txt index 45db2eb1..f61635fd 100644 --- a/build/CMakeLists.txt +++ b/build/CMakeLists.txt @@ -98,9 +98,6 @@ elseif(CMAKE_BUILD_TYPE MATCHES release) endif() endif() -# sync examples in documentation with example folder -add_custom_target(sync_examples ALL COMMAND python ${PROJECT_SOURCE_DIR}/../doc/sync_code.py) - # checks if(BUILD_CHECKS) add_executable(axis_size ../test/axis_size.cpp) @@ -143,7 +140,6 @@ foreach(SRC IN ITEMS ${TEST_SOURCES}) endif() else() add_executable(${BASENAME} ${SRC}) - add_dependencies(${BASENAME} sync_examples) if (SANITIZE) target_compile_options(${BASENAME} PRIVATE -fsanitize=address,undefined) target_link_libraries(${BASENAME} -lasan -lubsan ${LIBRARIES}) diff --git a/doc/getting_started.qbk b/doc/getting_started.qbk index f0ac53f1..38e0a49a 100644 --- a/doc/getting_started.qbk +++ b/doc/getting_started.qbk @@ -6,87 +6,8 @@ To get you started quickly, here are some heavily commented examples to copy pas If possible, use the static histogram. It is faster and user errors are caught at compile time. -[c++]`` -#include -#include - -int main(int, char**) { - namespace bh = boost::histogram; - using namespace bh::literals; // enables _c suffix - - /* - create a static 1d-histogram with an axis that has 10 equidistant - bins on the real line from -1.0 to 2.0, and label it as "x" - */ - auto h = bh::make_static_histogram( - bh::axis::regular<>(10, -1.0, 2.0, "x") - ); - - // fill histogram with data, typically this would happen in a loop - h.fill(-1.5); // put in underflow bin - h.fill(-1.0); // included in first bin, bin interval is semi-open - h.fill(-0.5); - h.fill(1.1); - h.fill(0.3); - h.fill(1.7); - h.fill(2.0); // put in overflow bin, bin interval is semi-open - h.fill(20.0); // put in overflow bin - - /* - do a weighted fill using bh::weight, which accepts a double - */ - h.fill(0.1, bh::weight(2.5)); - - /* - iterate over bins, loop excludes under- and overflow bins - - index 0_c is a compile-time number, the only way in C++ to make - axis(...) to return a different type for each index - - for-loop yields instances of `bin_type`, usually is a semi-open - interval representing the bin, whose edges can be accessed with - methods `lower()` and `upper()`, but the choice depends on the - axis type, please look it up in the reference - - `operator()` returns the bin counter at index, you can then - access its `value() and `variance()` methods; the first returns the - actual count, the second returns a variance estimate of the count; - a bin_type is convertible into an index - (see Rationale section for what this means) - */ - for (auto ai : h.axis(0_c)) { - std::cout << "bin " << ai.idx() << " x in [" - << ai.lower() << ", " << ai.upper() << "): " - << h(ai).value() << " +/- " - << std::sqrt(h(ai).variance()) - << std::endl; - } - - // accessing under- and overflow bins is easy, use indices -1 and 10 - std::cout << "underflow bin [" << h.axis(0_c)[-1].lower() - << ", " << h.axis(0_c)[-1].upper() << "): " - << h(-1).value() << " +/- " << std::sqrt(h(-1).variance()) - << std::endl; - std::cout << "overflow bin [" << h.axis(0_c)[10].lower() - << ", " << h.axis(0_c)[10].upper() << "): " - << h(10).value() << " +/- " << std::sqrt(h(10).variance()) - << std::endl; - - /* program output: - - bin 0 x in [-1, -0.7): 1 +/- 1 - bin 1 x in [-0.7, -0.4): 1 +/- 1 - bin 2 x in [-0.4, -0.1): 0 +/- 0 - bin 3 x in [-0.1, 0.2): 2.5 +/- 2.5 - bin 4 x in [0.2, 0.5): 1 +/- 1 - bin 5 x in [0.5, 0.8): 0 +/- 0 - bin 6 x in [0.8, 1.1): 4 +/- 2 - bin 7 x in [1.1, 1.4): 1 +/- 1 - bin 8 x in [1.4, 1.7): 0 +/- 0 - bin 9 x in [1.7, 2): 1 +/- 1 - underflow bin [-inf, -1): 1 +/- 1 - overflow bin [2, inf): 2 +/- 1.41421 - - */ -} -`` +[import ../examples/getting_started_listing_01.cpp] +[getting_started_listing_01] [endsect] @@ -94,55 +15,9 @@ int main(int, char**) { Dynamic histograms are a bit slower than static histograms, but still faster than other libraries. Use a dynamic histogram when you only know at runtime which and how many axis are going to be used, for example, because you wrote a graphical user interface that uses Boost.Histogram underneath. -[c++]`` -#include -#include -#include -#include -#include +[import ../examples/getting_started_listing_02.cpp] +[getting_started_listing_02] -namespace br = boost::random; -namespace bh = boost::histogram; - -int main() { - /* - create a dynamic histogram with the factory `make_dynamic_histogram` - - axis can be passed directly just like for `make_static_histogram` - - in addition, the factory also accepts iterators over a sequence of - axis::any, the polymorphic type that can hold concrete axis types - */ - std::vector> axes; - axes.emplace_back(bh::axis::category({"red", "blue"})); - axes.emplace_back(bh::axis::regular<>(5, -5, 5, "x")); - axes.emplace_back(bh::axis::regular<>(5, -5, 5, "y")); - auto h = bh::make_dynamic_histogram(axes.begin(), axes.end()); - - // fill histogram with random numbers - br::mt19937 gen; - br::normal_distribution<> norm; - for (int i = 0; i < 1000; ++i) - h.fill(i % 2 ? "red" : "blue", norm(gen), norm(gen)); - - /* - print dynamic histogram by iterating over bins - - for most axis types, the for loop looks just like for a static - histogram, except that we can pass runtime numbers, too - - if the [bin type] of the axis is not convertible to a - double interval, one needs to cast axis::any before looping; - this is here the case for the category axis - */ - using cas = bh::axis::category; - for (auto cbin : bh::axis::cast(h.axis(0))) { - std::printf("%s\n", cbin.value().c_str()); - for (auto ybin : h.axis(2)) { // rows - for (auto xbin : h.axis(1)) { // columns - std::printf("%3.0f ", h(cbin, xbin, ybin).value()); - } - std::printf("\n"); - } - } -} -`` [note If you care about maximum performance: In this example, `axis::category` is used with two string labels "red" and "blue". It is faster to use an enum, `enum { red, blue };` and a `axis::category<>` axis. ] @@ -153,60 +28,8 @@ If you care about maximum performance: In this example, `axis::category - -namespace bh = boost::histogram; - -int main() { - /* - create a 1d-histogram in default configuration which - covers the real line from -1 to 1 in 100 bins, the same - call with `make_dynamic_histogram` would also work - */ - auto h = bh::make_static_histogram(bh::axis::regular<>(100, -1, 1)); - - // do something with h -} -`` +[import ../examples/guide_make_static_histogram.cpp] +[guide_make_static_histogram] An axis object defines how input values are mapped to bins, which means that it defines the number of bins along that axis and a mapping function from input values to bins. If you provide one axis, the histogram is one-dimensional. If you provide two, it is two-dimensional, and so on. When you work with dynamic histograms, you can also create a sequence of axes at run-time and pass them to the factory: -[c++]`` -#include -#include - -namespace bh = boost::histogram; - -int main() { - // create vector of axes, axis::any is a polymorphic axis type - auto v = std::vector>(); - v.push_back(bh::axis::regular<>(100, -1, 1)); - v.push_back(bh::axis::integer<>(1, 7)); - - // create dynamic histogram (make_static_histogram be used with iterators) - auto h = bh::make_dynamic_histogram(v.begin(), v.end()); - - // do something with h -} -`` +[import ../examples/guide_make_dynamic_histogram.cpp] +[guide_make_dynamic_histogram] [funcref boost::histogram::make_static_histogram make_static_histogram] cannot handle this case because a static histogram can only be constructed when the number of types of all axes are known already at compile time. While strictly speaking that is also true in this example, you could have filled the vector also at run-time, based on run-time user input. @@ -88,38 +58,15 @@ The [classref boost::histogram::axis::regular regular axis] also accepts a secon In addition to the required parameters for an axis, you can assign an optional label to any axis, which helps to remember what the axis is about. Example: you have census data and you want to investigate how yearly income correlates with age, you could do: -[c++]`` -#include - -namespace bh = boost::histogram; - -int main() { - // create a 2d-histogram with an "age" and an "income" axis - auto h = bh::make_static_histogram( - bh::axis::regular<>(20, 0, 100, "age in years"), - bh::axis::regular<>(20, 0, 100, "yearly income in $1000") - ); - - // do something with h -} -`` +[import ../examples/guide_axis_with_labels.cpp] +[guide_axis_with_labels] Without the labels it would be difficult to remember which axis was covering which quantity, because they look the same otherwise. Labels are the only axis property that can be changed later. Axes objects with different label do not compare equal with `operator==`. By default, additional under- and overflow bins are added automatically for each axis where that makes sense. If you create an axis with 20 bins, the histogram will actually have 22 bins along that axis. The two extra bins are generally very good to have, as explained in [link histogram.rationale.uoflow the rationale]. If you are certain that the input cannot exceed the axis range, you can disable the extra bins to save memory. This is done by passing the enum value [enumref boost::histogram::axis::uoflow uoflow::off] to the axis constructor: -[c++]`` -#include - -namespace bh = boost::histogram; - -int main() { - // create a 1d-histogram for dice throws with integer values from 1 to 6 - auto h = bh::make_static_histogram(bh::axis::integer<>(1, 7, "eyes", bh::axis::uoflow::off)); - - // do something with h -} -`` +[import ../examples/guide_axis_with_uoflow_off.cpp] +[guide_axis_with_uoflow_off] We use an [classref boost::histogram::axis::integer integer axis] here, because the input values are integers and we want one bin for each eye value. @@ -131,60 +78,20 @@ We use an [classref boost::histogram::axis::integer integer axis] here, because [section Fill histogram] -After you created a histogram, you want to insert (possibly multi-dimensional) data. This is done with the flexible `histogram::fill(...)` method, which you typically call in a loop. Some extra parameters can be passed to the method as shown in the next example. +After you created a histogram, you want to insert (possibly multi-dimensional) data. This is done with the flexible `histogram(...)` call, which you typically do in a loop. Some extra parameters can be passed to the method as shown in the next example. -[c++]`` -#include +[import ../examples/guide_fill_histogram.cpp] +[guide_fill_histogram] -namespace bh = boost::histogram; +`histogram(...)` either takes `N` arguments or one container with `N` elements, where `N` is equal to the number of axes of the histogram, finds the corresponding bin, and increments the bin counter by one. -int main() { - auto h = bh::make_dynamic_histogram(bh::axis::integer<>(0, 4), - bh::axis::regular<>(10, 0, 5)); +This library was designed to feel familiar to users of [@boost:/libs/accumulators/index.html Boost.Accumulators] and to make it compatible with `std::for_each`, that's why values are filled with `operator()`. If you have the input values in a collection, `std::for_each` fills the histogram in one line. - // fill histogram, number of arguments must be equal to number of axes - h.fill(0, 4.1); // increases bin counter by one - h.fill(3, 3.4, bh::weight(2)); // increase bin counter by 2 - h.fill(bh::weight(2), 3, 3.4); // the same as the previous call +`histogram(weight(x), ...)` does the same as the first call, but increments the bin counter by the value `x`. The type of `x` is not restricted, usually it is a real number. The `weight(x)` helper class must be first argument. You can freely mix calls with and without a `weight`. Calls without a `weight` act like the weight is `1`. - // a dynamic histogram also supports fills from an interator range while - // a static histogram does not allow it; using a range of wrong length - // is an error - std::vector values = {4, 3.1}; - h.fill(values.begin(), values.end()); -} -`` +Why weighted increments are sometimes useful, especially in a scientific context, is explained [link histogram.rationale.weights in the rationale]. If you don't see the point, you can just ignore this type of call. This feature does not affect the performance of the histogram if you don't use it. -Here is an explanation for these calls. - -* `histogram::fill(...)` takes `N` arguments, where `N` is equal to the number of axes of the histogram, finds the corresponding bin, and increments the bin counter by one. - -* `histogram::fill(..., weight(x))` does the same as the first call, but increments the bin counter by a weight `x`. The type of `x` is not restricted, but usually it is a real number. The position of the `weight` helper class is not restricted, it can appear at any place. - -Having more than one `weight` instance in the same call is an error. You can freely mix calls with and without a `weight`. Calls without a `weight` act like the weight is `1`. - -Why weighted increments are sometimes useful, mostly in a scientific context. How weights are handled is explained [link histogram.rationale.weights in the rationale]. If you don't see the point, you can just ignore this type of call, this feature does not affect the preformance of the histogram if you don't use it. - -[note The first call to a weighted fill internally switches the default storage from integral counters to another type, which holds two real numbers per bin, one for the sum of weights (the weighted count), and another for the sum of weights squared (the variance of the weighted count). This is not necessary for unweighted fills, because the two sums are identical is all weights are `1`. The default storage automatically optimizes this case by using only one integral number per bin.] - -In contrast to [@boost:/libs/accumulators/index.html Boost.Accumulators], this library asks you to write a loop yourself to pass data to the histogram, because that is the most flexible solution. If you prefer a functional programming style, you can use a lambda, as shown in this example. - -[c++]`` -#include -#include -#include - -namespace bh = boost::histogram; - -int main() { - auto h = bh::make_static_histogram(bh::axis::integer<>(0, 9)); - std::vector v{0, 1, 2, 3, 4, 5, 6, 7, 8, 9}; - std::for_each(v.begin(), v.end(), - [&h](int x) { h.fill(x, bh::weight(2.0)); } - ); - // h is now filled -} -`` +[note The first call to a weighted fill internally switches the default storage from integral counters to another type, which holds two real numbers per bin, one for the sum of weights (the weighted count), and another for the sum of weights squared (the variance of the weighted count). This is not necessary for unweighted fills, because the two sums are identical is all weights are `1`. The default storage automatically optimizes this case by using only one integral number per bin as long as no weights are encountered.] [endsect] @@ -196,60 +103,10 @@ To access each bin, you use a multi-dimensional index, which consists of a seque The calls are demonstrated in the next example. -[c++]`` -#include -#include +[import ../examples/guide_access_bin_counts.cpp] +[guide_access_bin_counts] -namespace bh = boost::histogram; - -int main() { - // make histogram with 2 x 2 = 4 bins (not counting under-/overflow bins) - auto h = bh::make_static_histogram( - bh::axis::regular<>(2, -1, 1), - bh::axis::regular<>(2, 2, 4) - ); - - h.fill(-0.5, 2.5, bh::weight(1)); // bin index 0, 0 - h.fill(-0.5, 3.5, bh::weight(2)); // bin index 0, 1 - h.fill( 0.5, 2.5, bh::weight(3)); // bin index 1, 0 - h.fill( 0.5, 3.5, bh::weight(4)); // bin index 1, 1 - - // access count value, number of indices must match number of axes - std::cout << h(0, 0).value() << " " - << h(0, 1).value() << " " - << h(1, 0).value() << " " - << h(1, 1).value() - << std::endl; - - // prints: 1 2 3 4 - - // access count variance, number of indices must match number of axes - std::cout << h(0, 0).variance() << " " - << h(0, 1).variance() << " " - << h(1, 0).variance() << " " - << h(1, 1).variance() - << std::endl; - // prints: 1 4 9 16 - - // you can also make a copy of the type that holds the bin count - auto c11 = h(1, 1); - std::cout << c11.value() << " " << c11.variance() << std::endl; - // prints: 4 16 - - // a dynamic histogram also supports access via an interator range, - // using a range of wrong length is an error - auto h2 = bh::make_dynamic_histogram( - bh::axis::regular<>(2, -1, 1), - bh::axis::regular<>(2, 2, 4) - ); - - std::vector idx(2); - idx = {0, 1}; - std::cout << h2(idx.begin(), idx.end()).value() << std::endl; -} -`` - -[note The numbers returned by `value()` and `variance()` are equal, if weighted fills are not used. The internal structure, which handles the bin counters, has been optimised for this common case. Internally only a single integral number per bin is used until a weighted fill, then the counters internally switch to storing two real numbers per bin. If the very first call to `histogram::fill` is already a weighted increment, the two real numbers per bin are allocated directly without any superfluous conversion from integral counters to double counters. This special case is efficiently handled.] +[note The numbers returned by `value()` and `variance()` are always equal, if weighted fills are not used. The internal structure, which handles the bin counters, has been optimised for this common case. Internally only a single integral number per bin is used until a weighted fill, then the counters internally switch to storing two real numbers per bin. If the very first call to `histogram(...)` is already a weighted increment, the two real numbers per bin are allocated directly without any superfluous conversion from integral counters to double counters. This special case is efficiently handled.] [endsect] @@ -272,68 +129,14 @@ Adding histograms is useful, if you want to parallelize the filling of a histogr Multiplying by a number is useful to re-weight histograms before adding them, for those who need to work with weights. Multiplying by a factor `x` has a different effect on value and variance of each bin counter. The value is multiplied by `x`, but the variance is multiplied by `x*x`. This follows from the properties of the variance, as explained in [link histogram.rationale.variance the rationale]. -[warning Because of behavior of the variance, adding a histogram to itself is not identical to multiplying the original histogram by two, as far as the variance is concerned.] +[warning Because of special behavior of the variance, adding a histogram to itself is not identical to multiplying the original histogram by two, as far as the variance is concerned.] [note Scaling a histogram automatically converts the bin counters from an integral number per bin to two real numbers per bin, if that has not happened already, because value and variance are different after the multiplication.] Here is an example which demonstrates the supported operators. -[c++]`` -#include -#include - -namespace bh = boost::histogram; - -int main() { - // make two histograms - auto h1 = bh::make_static_histogram(bh::axis::regular<>(2, -1, 1)); - auto h2 = bh::make_static_histogram(bh::axis::regular<>(2, -1, 1)); - - h1.fill(-0.5); // counts are: 1 0 - h2.fill(0.5); // counts are: 0 1 - - // add them - auto h3 = h1; - h3 += h2; // counts are: 1 1 - - // adding multiple histograms at once is efficient and does not create - // superfluous temporaries since operator+ functions are overloaded to - // accept and return rvalue references where possible - auto h4 = h1 + h2 + h3; // counts are: 2 2 - - std::cout << h4(0).value() << " " << h4(1).value() << std::endl; - // prints: 2 2 - - // multiply by number - h4 *= 2; // counts are: 4 4 - - // divide by number - auto h5 = h4 / 4; // counts are: 1 1 - - std::cout << h5(0).value() << " " << h5(1).value() << std::endl; - // prints: 1 1 - - // compare histograms - std::cout << (h4 == 4 * h5) << " " << (h4 != h5) << std::endl; - // prints: 1 1 - - // note: special effect of multiplication on counter variance - auto h = bh::make_static_histogram(bh::axis::regular<>(2, -1, 1)); - h.fill(-0.5); // counts are: 1 0 - std::cout << "value " << (2 * h)(0).value() - << " " << (h + h)(0).value() << "\n" - << "variance " << (2 * h)(0).variance() - << " " << (h + h)(0).variance() << std::endl; - // equality operator also checks variances, so the statement is false - std::cout << (h + h == 2 * h) << std::endl; - - /* prints: - value 2 2 - variance 4 2 - 0 - */ -} -`` +[import ../examples/guide_histogram_operators.cpp] +[guide_histogram_operators] [endsect] @@ -345,51 +148,8 @@ For this purpose use the `histogram::reduce_to(...)` method, which returns a new Here is an example to illustrates this. -[c++]`` -#include -#include - -namespace bh = boost::histogram; - -int main() { - using namespace bh::literals; // enables _c suffix - - // make a 2d histogram - auto h = bh::make_static_histogram(bh::axis::regular<>(4, -1, 1), - bh::axis::integer<>(0, 4)); - - h.fill(-0.9, 0); - h.fill(0.9, 3); - h.fill(0.1, 2); - - auto hr0 = h.reduce_to(0_c); // keep only first axis - auto hr1 = h.reduce_to(1_c); // keep only second axis - - /* - reduce does not remove counts; returned histograms are summed over - the removed axes, so h, hr0, and hr1 have same number of total counts - */ - std::cout << h.sum().value() << " " - << hr0.sum().value() << " " - << hr1.sum().value() << std::endl; - // prints: 3 3 3 - - for (auto yi : h.axis(1_c)) { - for (auto xi : h.axis(0_c)) { - std::cout << h(xi, yi).value() << " "; - } - std::cout << std::endl; - } - - for (auto xi : hr0.axis()) - std::cout << hr0(xi).value() << " "; - std::cout << std::endl; - - for (auto yi : hr1.axis()) - std::cout << hr1(yi).value() << " "; - std::cout << std::endl; -} -`` +[import ../examples/guide_histogram_reduction.cpp] +[guide_histogram_reduction] [endsect] @@ -397,43 +157,8 @@ int main() { Simple ostream operators are shipped with the library, which are internally used by the Python interface bindings. These give text representations of axis and histogram configurations, but do not show the histogram content. For users, the builtin streaming operators may be useful for debugging. They are not included by the super header `#include `, so that users can use their own implementations. The following example shows the effect of output streaming. -[c++]`` -#include -#include -#include - -namespace bh = boost::histogram; - -int main() { - namespace axis = boost::histogram::axis; - - enum { A, B, C }; - - auto h = bh::make_static_histogram( - axis::regular<>(2, -1, 1, "regular1", axis::uoflow::off), - axis::regular(2, 1, 10, "regular2"), - axis::circular<>(4, 0.1, 1.0, "polar"), - axis::variable<>({-1, 0, 1}, "variable", axis::uoflow::off), - axis::category<>({A, B, C}, "category"), - axis::integer<>(-1, 1, "integer", axis::uoflow::off) - ); - - std::cout << h << std::endl; - - /* prints: - - histogram( - regular(2, -1, 1, label='regular1', uoflow=False), - regular_log(2, 1, 10, label='regular2'), - circular(4, phase=0.1, perimeter=1, label='polar'), - variable(-1, 0, 1, label='variable', uoflow=False), - category(0, 1, 2, label='category'), - integer(-1, 1, label='integer', uoflow=False), - ) - - */ -} -`` +[import ../examples/guide_histogram_streaming.cpp] +[guide_histogram_streaming] [endsect] @@ -441,46 +166,8 @@ int main() { The library supports serialization via [@boost:/libs/serialization/index.html Boost.Serialization]. The serialization machinery is used in the Python module to enable pickling of histograms. The streaming code is not included by the super header `#include `, so that the library can be used without having Boost.Serialization installed. -[c++]`` -#include -#include // includes serialization code -#include -#include -#include - -namespace bh = boost::histogram; - -int main() { - auto a = bh::make_static_histogram(bh::axis::regular<>(3, -1, 1, "r"), - bh::axis::integer<>(0, 2, "i")); - a.fill(0.5, 1); - - std::string buf; // holds persistent representation - - // store histogram - { - std::ostringstream os; - boost::archive::text_oarchive oa(os); - oa << a; - buf = os.str(); - } - - auto b = decltype(a)(); // create a default-constructed second histogram - - std::cout << "before restore " << (a == b) << std::endl; - // prints: before restore 0 - - // load histogram - { - std::istringstream is(buf); - boost::archive::text_iarchive ia(is); - ia >> b; - } - - std::cout << "after restore " << (a == b) << std::endl; - // prints: after restore 1 -} -`` +[import ../examples/guide_histogram_serialization.cpp] +[guide_histogram_serialization] [endsect] @@ -504,36 +191,8 @@ The documentation of the Python interface is best explored from within the Pytho Here is an example of the Python module in action. -[python]`` -import histogram as hg - -# make 1-d histogram with 5 logarithmic bins from 1e0 to 1e5 -h = hg.histogram(hg.axis.regular_log(5, 1e0, 1e5, "x")) - -# fill histogram with numbers -for x in (2e0, 2e1, 2e2, 2e3, 2e4): - h.fill(x, weight=4) # increment bin counter by 4 - -# iterate over bins and access bin counter -for idx, (lower, upper) in enumerate(h.axis(0)): - print "bin {0} x in [{1}, {2}): {3} +/- {4}".format( - idx, lower, upper, h(idx).value, h(idx).variance ** 0.5) - -# under- and overflow bins are accessed like in C++ -lo, up = h.axis(0)[-1] -print "underflow [{0}, {1}): {2} +/- {3}".format(lo, up, h(-1).value, h(-1).variance) -lo, up = h.axis(0)[5] -print "overflow [{0}, {1}): {2} +/- {3}".format(lo, up, h(5).value, h(5).variance) - -# prints: -# bin 0 x in [1.0, 10.0): 4.0 +/- 4.0 -# bin 1 x in [10.0, 100.0): 4.0 +/- 4.0 -# bin 2 x in [100.0, 1000.0): 4.0 +/- 4.0 -# bin 3 x in [1000.0, 10000.0): 4.0 +/- 4.0 -# bin 4 x in [10000.0, 100000.0): 4.0 +/- 4.0 -# underflow [0.0, 1.0): 0.0 +/- 0.0 -# overflow [100000.0, inf): 0.0 +/- 0.0 -`` +[import ../examples/guide_python_histogram.py] +[guide_python_histogram] [endsect] @@ -541,91 +200,15 @@ print "overflow [{0}, {1}): {2} +/- {3}".format(lo, up, h(5).value, h(5).varian Looping over a large collection of numbers in Python to fill a histogram is very slow. If the Python module was build with Numpy support, the input values can be passed as Numpy arrays. This is very fast and convenient. The performance is almost as good as using C++ directly, and usually better than Numpy's histogram functions, see the [link histogram.benchmarks benchmark]. Here is an example to illustrate: -[python]`` -import histogram as hg -import numpy as np +[import ../examples/guide_numpy_support.py] +[guide_numpy_support] -# create 2d-histogram with two axes with 10 equidistant bins from -3 to 3 -h = hg.histogram(hg.axis.regular(10, -3, 3, "x"), - hg.axis.regular(10, -3, 3, "y")) - -# generate some numpy arrays with data to fill into histogram, -# in this case normal distributed random numbers in x and y -x = np.random.randn(1000) -y = 0.5 * np.random.randn(1000) - -# fill histogram with numpy arrays, this is very fast -h.fill(x, y) # call looks the same as if x, y were values - -# get representations of the bin edges as Numpy arrays; this representation -# differs from `list(h.axis(0))` as explained in the next example -x = np.array(h.axis(0)) -y = np.array(h.axis(1)) - -# creates a view of the counts (no copy involved) -count_matrix = np.asarray(h) - -# cut off the under- and overflow bins to not confuse matplotib (no copy) -reduced_count_matrix = count_matrix[:-2,:-2] - -try: - # draw the count matrix - import matplotlib.pyplot as plt - plt.pcolor(x, y, reduced_count_matrix.T) - plt.xlabel(h.axis(0).label) - plt.ylabel(h.axis(1).label) - plt.savefig("example_2d_python.png") -except ImportError: - # ok, no matplotlib, then just print the full count matrix - print count_matrix - - # output of the print looks something like this, the two right-most rows - # and two down-most columns represent under-/overflow bins - # [[ 0 0 0 1 5 0 0 1 0 0 0 0] - # [ 0 0 0 1 17 11 6 0 0 0 0 0] - # [ 0 0 0 5 31 26 4 1 0 0 0 0] - # [ 0 0 3 20 59 62 26 4 0 0 0 0] - # [ 0 0 1 26 96 89 16 1 0 0 0 0] - # [ 0 0 4 21 86 84 20 1 0 0 0 0] - # [ 0 0 1 24 71 50 15 2 0 0 0 0] - # [ 0 0 0 6 26 37 7 0 0 0 0 0] - # [ 0 0 0 0 11 10 2 0 0 0 0 0] - # [ 0 0 0 1 2 3 1 0 0 0 0 0] - # [ 0 0 0 0 0 2 0 0 0 0 0 0] - # [ 0 0 0 0 0 1 0 0 0 0 0 0]] -`` - -[note The `fill(...)` method accepts not only Numpy arrays as arguments, but also any sequence that can be converted into a Numpy array with `dtype=float`. The best performance is achieved, if numpy arrays are used with dtype `float` directly, so that no conversion happens.] +[note `histogram(...)` accepts any sequence that can be converted into a Numpy array with `dtype=float` for convenience. To achieve the best performance, use numpy arrays with dtype `float` as input to avoid a conversion.] The conversion of an axis to a Numpy array has been optimised for interoperability with other Numpy functions and matplotlib. For an axis with `N` bins, a sequence of bin edges is generated of length `N+1`. The Numpy array is a unique ordered sequence of all bin edges, including the upper edge of the last bin. The following code illustrates how the sequence is constructed. -[python]`` -import histogram as hg -import numpy as np - -ax = hg.axis.regular(5, 0, 1) -xedge1 = np.array(ax) # this is equivalent to... -xedge2 = [] -for idx, (lower, upper) in enumerate(ax): - xedge2.append(lower) - if idx == len(ax)-1: - xedge2.append(upper) - -print xedge1 -print xedge2 - -# prints: -# [ 0. 0.2 0.4 0.6 0.8 1. ] -# [0.0, 0.2, 0.4, 0.6, 0.8, 1.0] - -# sequences constructed from an axis use its iterator, the result differs -xedge3 = list(ax) - -print xedge3 - -# prints: -# [(0., 0.2), (0.2, 0.4), (0.4, 0.6), (0.8, 1.0)] -`` +[import ../examples/guide_python_axis_representation.py] +[guide_python_axis_representation] Using a Numpy array to view the count matrix reveals an implementation detail of how under- and overflow bins are stored internally. For a 1-d histogram with 5 bins, the counter structure looks like this: @@ -641,31 +224,8 @@ It is very easy to crop the extra counters by slicing each dimension with the sl Histograms in Python support the same operators as in C++. In addition, histograms support the Pickle protocol. An example: -[python]`` -import histogram as hg -import cPickle - -h1 = hg.histogram(hg.axis.regular(2, -1, 1)) -h2 = hg.histogram(h1) # creates copy -h1.fill(-0.5) -h2.fill(0.5) - -# arithmetic operators (see performance note below) -h3 = h1 + h2 -h4 = h3 * 2 - -print h4(0).value, h4(1).value -# prints: 2.0 2.0 - -# now save the histogram -with open("h4_saved.pkl", "w") as f: - cPickle.dump(h4, f) -with open("h4_saved.pkl", "r") as f: - h5 = cPickle.load(f) - -print h4 == h5 -# prints: True -`` +[import ../examples/guide_histogram_pickle.py] +[guide_histogram_pickle] [note Python has no concept of rvalue references and therefore cannot avoid creating temporary objects in arithmetic operations like C++ can. A statement `h3 = h1 + h2` is equivalent to `tmp = hg.histogram(h1); tmp += h2; h3 = hg.histogram(tmp)`, which amounts to two allocations, one in the first and one in the last copy-assignment, while only one allocation is really necessary. @@ -685,48 +245,13 @@ Here is an example, consisting of a C++ part, which needs be compiled as a share C++ code that is compiled into a shared library. -[c++]`` -#include -#include - -namespace bh = boost::histogram; -namespace bp = boost::python; - -// function that runs in C++ and accepts reference to dynamic histogram -void process(bh::dynamic_histogram<>& h) { - // fill histogram, in reality this would be arbitrarily complex code - for (int i = 0; i < 4; ++i) - h.fill(0.25 * i, i); -} - -// a minimal Python module, which exposes the process function to Python -BOOST_PYTHON_MODULE(cpp_filler) { - bp::def("process", process); -} -`` +[import ../examples/guide_mixed_cpp_python.cpp] +[guide_mixed_cpp_python_part_cpp] Python script which uses the shared library. -[python]`` -import histogram as bh -import cpp_filler - -h = bh.histogram(bh.axis.regular(4, 0, 1), - bh.axis.integer(0, 4)) - -cpp_filler.process(h) # histogram is filled with input values in C++ - -for iy, y in enumerate(h.axis(1)): - for ix, x in enumerate(h.axis(0)): - print h(ix, iy).value, - print - -# prints: -# 1.0 0.0 0.0 0.0 -# 0.0 1.0 0.0 0.0 -# 0.0 0.0 1.0 0.0 -# 0.0 0.0 0.0 1.0 -`` +[import ../examples/guide_mixed_cpp_python.py] +[guide_mixed_cpp_python_part_py] [endsect] @@ -742,46 +267,8 @@ In C++, users can implement their own axis class without touching any library co The simplest way to make a custom axis is to derive from a builtin class. Here is a contrieved example of a custom axis that inherits from the [classref boost::histogram::axis::integer integer axis] and accepts c-strings representing numbers. -[c++]`` -#include -#include - -namespace bh = boost::histogram; - -// custom axis, which adapts builtin integer axis -struct custom_axis : public bh::axis::integer<> { - using value_type = const char*; // type that is fed to the axis - - using integer::integer; // inherit ctors of base - - // the customization point - // - accept const char* and convert to int - // - then call index method of base class - int index(value_type s) const { - return integer::index(std::atoi(s)); - } -}; - -int main() { - auto h = bh::make_static_histogram(custom_axis(0, 3)); - h.fill("-10"); - h.fill("0"); - h.fill("1"); - h.fill("9"); - - for (auto xi : h.axis()) { - std::cout << "bin " << xi.idx() << " [" << xi.lower() << ", " - << xi.upper() << ") " << h(xi).value() - << std::endl; - } - - /* prints: - bin 0 [0, 1) 1 - bin 1 [1, 2] 1 - bin 2 [2, 3] 0 - */ -} -`` +[import ../examples/guide_custom_axis.cpp] +[guide_custom_axis] [note The axis types available in Python cannot be changed without changing library code and recompiling the library.] @@ -797,21 +284,8 @@ If you work exclusively with weighted fills, [classref boost::histogram::array_s Here is an example of a histogram construced with an alternative storage policy. -[c++]`` -#include -#include - -namespace bh = boost::histogram; - -int main() { - // create static histogram with array_storage, using int as counter type - auto h = bh::make_static_histogram_with>( - bh::axis::regular<>(10, 0, 1) - ); - - // do something with h -} -`` +[import ../examples/guide_custom_storage.cpp] +[guide_custom_storage] [warning The guarantees regarding protection against overflow and count capping are only valid, if the default [classref boost::histogram::adaptive_storage adaptive_storage] is used. If you change the storage policy, you need know what you are doing.] diff --git a/examples/getting_started_listing_01.cpp b/examples/getting_started_listing_01.cpp new file mode 100644 index 00000000..2eb8040f --- /dev/null +++ b/examples/getting_started_listing_01.cpp @@ -0,0 +1,84 @@ +//[ getting_started_listing_01 + +#include +#include + +int main(int, char**) { + namespace bh = boost::histogram; + using namespace bh::literals; // enables _c suffix + + /* + create a static 1d-histogram with an axis that has 10 equidistant + bins on the real line from -1.0 to 2.0, and label it as "x" + */ + auto h = bh::make_static_histogram( + bh::axis::regular<>(10, -1.0, 2.0, "x") + ); + + // fill histogram with data, typically this happens in a loop + h(-1.5); // put in underflow bin + h(-1.0); // included in first bin, bin interval is semi-open + h(2.0); // put in overflow bin, bin interval is semi-open + h(20.0); // put in overflow bin + + // STL algorithms are supported + auto data = { -0.5, 1.1, 0.3, 1.7 }; + std::for_each(data.begin(), data.end(), h); + + /* + do a weighted fill using bh::weight, a wrapper for any type, + which may appear at the beginning of the argument list + */ + h(bh::weight(1.0), 0.1); + + /* + iterate over bins, loop excludes under- and overflow bins + - index 0_c is a compile-time number, the only way in C++ to make + axis(...) to return a different type for each index + - for-loop yields instances of `bin_type`, usually is a semi-open + interval representing the bin, whose edges can be accessed with + methods `lower()` and `upper()`, but the choice depends on the + axis type, please look it up in the reference + - `operator()` returns the bin counter at index, you can then + access its `value() and `variance()` methods; the first returns the + actual count, the second returns a variance estimate of the count; + a bin_type is convertible into an index + (see Rationale section for what this means) + */ + for (auto bin : h.axis(0_c)) { + std::cout << "bin " << bin.idx() << " x in [" + << bin.lower() << ", " << bin.upper() << "): " + << h.bin(bin).value() << " +/- " + << std::sqrt(h.bin(bin).variance()) + << std::endl; + } + + // accessing under- and overflow bins is easy, use indices -1 and 10 + std::cout << "underflow bin [" << h.axis(0_c)[-1].lower() + << ", " << h.axis(0_c)[-1].upper() << "): " + << h.bin(-1).value() << " +/- " << std::sqrt(h.bin(-1).variance()) + << std::endl; + std::cout << "overflow bin [" << h.axis(0_c)[10].lower() + << ", " << h.axis(0_c)[10].upper() << "): " + << h.bin(10).value() << " +/- " << std::sqrt(h.bin(10).variance()) + << std::endl; + + /* program output: + + bin 0 x in [-1, -0.7): 1 +/- 1 + bin 1 x in [-0.7, -0.4): 1 +/- 1 + bin 2 x in [-0.4, -0.1): 0 +/- 0 + bin 3 x in [-0.1, 0.2): 2.5 +/- 2.5 + bin 4 x in [0.2, 0.5): 1 +/- 1 + bin 5 x in [0.5, 0.8): 0 +/- 0 + bin 6 x in [0.8, 1.1): 4 +/- 2 + bin 7 x in [1.1, 1.4): 1 +/- 1 + bin 8 x in [1.4, 1.7): 0 +/- 0 + bin 9 x in [1.7, 2): 1 +/- 1 + underflow bin [-inf, -1): 1 +/- 1 + overflow bin [2, inf): 2 +/- 1.41421 + + */ +} + +//] diff --git a/examples/getting_started_listing_02.cpp b/examples/getting_started_listing_02.cpp new file mode 100644 index 00000000..475bbd5c --- /dev/null +++ b/examples/getting_started_listing_02.cpp @@ -0,0 +1,51 @@ +//[ getting_started_listing_02 + +#include +#include +#include +#include +#include + +namespace br = boost::random; +namespace bh = boost::histogram; + +int main() { + /* + create a dynamic histogram with the factory `make_dynamic_histogram` + - axis can be passed directly just like for `make_static_histogram` + - in addition, the factory also accepts iterators over a sequence of + axis::any, the polymorphic type that can hold concrete axis types + */ + std::vector> axes; + axes.emplace_back(bh::axis::category({"red", "blue"})); + axes.emplace_back(bh::axis::regular<>(5, -5, 5, "x")); + axes.emplace_back(bh::axis::regular<>(5, -5, 5, "y")); + auto h = bh::make_dynamic_histogram(axes.begin(), axes.end()); + + // fill histogram with random numbers + br::mt19937 gen; + br::normal_distribution<> norm; + for (int i = 0; i < 1000; ++i) + h(i % 2 ? "red" : "blue", norm(gen), norm(gen)); + + /* + print dynamic histogram by iterating over bins + - for most axis types, the for loop looks just like for a static + histogram, except that we can pass runtime numbers, too + - if the [bin type] of the axis is not convertible to a + double interval, one needs to cast axis::any before looping; + this is here the case for the category axis + */ + using cas = bh::axis::category; + for (auto cbin : bh::axis::cast(h.axis(0))) { + std::printf("%s\n", cbin.value().c_str()); + for (auto ybin : h.axis(2)) { // rows + for (auto xbin : h.axis(1)) { // columns + std::printf("%3.0f ", h.bin(cbin, xbin, ybin).value()); + } + std::printf("\n"); + } + } +} + +//] diff --git a/examples/getting_started_listing_03.py b/examples/getting_started_listing_03.py new file mode 100644 index 00000000..733db994 --- /dev/null +++ b/examples/getting_started_listing_03.py @@ -0,0 +1,56 @@ +#[ getting_started_listing_03 + +import histogram as hg +import numpy as np + +# create 2d-histogram with two axes with 10 equidistant bins from -3 to 3 +h = hg.histogram(hg.axis.regular(10, -3, 3, "x"), + hg.axis.regular(10, -3, 3, "y")) + +# generate some numpy arrays with data to fill into histogram, +# in this case normal distributed random numbers in x and y +x = np.random.randn(1000) +y = 0.5 * np.random.randn(1000) + +# fill histogram with numpy arrays, this is very fast +h(x, y) + +# get representations of the bin edges as Numpy arrays, this representation +# differs from `list(h.axis(0))`, because it is optimised for compatibility +# with existing Numpy code, i.e. to replace numpy.histogram +x = np.array(h.axis(0)) +y = np.array(h.axis(1)) + +# creates a view of the counts (no copy involved) +count_matrix = np.asarray(h) + +# cut off the under- and overflow bins to not confuse matplotib (no copy) +reduced_count_matrix = count_matrix[:-2,:-2] + +try: + # draw the count matrix + import matplotlib.pyplot as plt + plt.pcolor(x, y, reduced_count_matrix.T) + plt.xlabel(h.axis(0).label) + plt.ylabel(h.axis(1).label) + plt.savefig("example_2d_python.png") +except ImportError: + # ok, no matplotlib, then just print the full count matrix + print count_matrix + + # output of the print looks something like this, the two right-most rows + # and two down-most columns represent under-/overflow bins + # [[ 0 0 0 1 5 0 0 1 0 0 0 0] + # [ 0 0 0 1 17 11 6 0 0 0 0 0] + # [ 0 0 0 5 31 26 4 1 0 0 0 0] + # [ 0 0 3 20 59 62 26 4 0 0 0 0] + # [ 0 0 1 26 96 89 16 1 0 0 0 0] + # [ 0 0 4 21 86 84 20 1 0 0 0 0] + # [ 0 0 1 24 71 50 15 2 0 0 0 0] + # [ 0 0 0 6 26 37 7 0 0 0 0 0] + # [ 0 0 0 0 11 10 2 0 0 0 0 0] + # [ 0 0 0 1 2 3 1 0 0 0 0 0] + # [ 0 0 0 0 0 2 0 0 0 0 0 0] + # [ 0 0 0 0 0 1 0 0 0 0 0 0]] + +#] \ No newline at end of file diff --git a/examples/getting_started_listing_04.py b/examples/getting_started_listing_04.py new file mode 100644 index 00000000..3471b661 --- /dev/null +++ b/examples/getting_started_listing_04.py @@ -0,0 +1,18 @@ +#[ getting_started_listing_04 + +from __future__ import print_function +import histogram as hg + +# make 1-d histogram with 5 logarithmic bins from 1e0 to 1e5 +h = hg.histogram(hg.axis.regular_log(5, 1e0, 1e5, "x")) + +# fill histogram with numbers +for x in (2e0, 2e1, 2e2, 2e3, 2e4): + h(x, weight=2) + +# iterate over bins and access bin counter +for idx, (lower, upper) in enumerate(h.axis(0)): + print("bin {0} x in [{1}, {2}): {3} +/- {4}".format( + idx, lower, upper, h.bin(idx).value, h.bin(idx).variance ** 0.5)) + +#] diff --git a/examples/guide_access_bin_counts.cpp b/examples/guide_access_bin_counts.cpp new file mode 100644 index 00000000..c34d0cfd --- /dev/null +++ b/examples/guide_access_bin_counts.cpp @@ -0,0 +1,56 @@ +//[ guide_access_bin_counts + +#include +#include +#include + +namespace bh = boost::histogram; + +int main() { + // make histogram with 2 x 2 = 4 bins (not counting under-/overflow bins) + auto h = bh::make_static_histogram( + bh::axis::regular<>(2, -1, 1), + bh::axis::regular<>(2, 2, 4) + ); + + h(bh::weight(1), -0.5, 2.5); // bin index 0, 0 + h(bh::weight(2), -0.5, 3.5); // bin index 0, 1 + h(bh::weight(3), 0.5, 2.5); // bin index 1, 0 + h(bh::weight(4), 0.5, 3.5); // bin index 1, 1 + + // access count value, number of indices must match number of axes + std::cout << h.bin(0, 0).value() << " " + << h.bin(0, 1).value() << " " + << h.bin(1, 0).value() << " " + << h.bin(1, 1).value() + << std::endl; + + // prints: 1 2 3 4 + + // access count variance, number of indices must match number of axes + std::cout << h.bin(0, 0).variance() << " " + << h.bin(0, 1).variance() << " " + << h.bin(1, 0).variance() << " " + << h.bin(1, 1).variance() + << std::endl; + // prints: 1 4 9 16 + + // you can also make a copy of the type that holds the bin count + auto c11 = h.bin(1, 1); + std::cout << c11.value() << " " << c11.variance() << std::endl; + // prints: 4 16 + + // histogram also supports access via container; using a container of + // wrong size trips an assertion in debug mode + auto idx = {0, 1}; + std::cout << h.bin(idx).value() << std::endl; + // prints: 2 + + // histogram also provides extended iterators + auto sum = std::accumulate(h.begin(), h.end(), + typename decltype(h)::element_type(0)); + std::cout << sum.value() << " " << sum.variance() << std::endl; + // prints: 10 30 +} + +//] diff --git a/examples/guide_axis_with_labels.cpp b/examples/guide_axis_with_labels.cpp new file mode 100644 index 00000000..e2fbfcda --- /dev/null +++ b/examples/guide_axis_with_labels.cpp @@ -0,0 +1,17 @@ +//[ guide_axis_with_labels + +#include + +namespace bh = boost::histogram; + +int main() { + // create a 2d-histogram with an "age" and an "income" axis + auto h = bh::make_static_histogram( + bh::axis::regular<>(20, 0, 100, "age in years"), + bh::axis::regular<>(20, 0, 100, "yearly income in $1000") + ); + + // do something with h +} + +//] diff --git a/examples/guide_axis_with_uoflow_off.cpp b/examples/guide_axis_with_uoflow_off.cpp new file mode 100644 index 00000000..786c61ea --- /dev/null +++ b/examples/guide_axis_with_uoflow_off.cpp @@ -0,0 +1,14 @@ +//[ guide_axis_with_uoflow_off + +#include + +namespace bh = boost::histogram; + +int main() { + // create a 1d-histogram for dice throws with integer values from 1 to 6 + auto h = bh::make_static_histogram(bh::axis::integer<>(1, 7, "eyes", bh::axis::uoflow::off)); + + // do something with h +} + +//] diff --git a/examples/guide_custom_axis.cpp b/examples/guide_custom_axis.cpp new file mode 100644 index 00000000..07edec06 --- /dev/null +++ b/examples/guide_custom_axis.cpp @@ -0,0 +1,42 @@ +//[ guide_custom_axis + +#include +#include + +namespace bh = boost::histogram; + +// custom axis, which adapts builtin integer axis +struct custom_axis : public bh::axis::integer<> { + using value_type = const char*; // type that is fed to the axis + + using integer::integer; // inherit ctors of base + + // the customization point + // - accept const char* and convert to int + // - then call index method of base class + int index(value_type s) const { + return integer::index(std::atoi(s)); + } +}; + +int main() { + auto h = bh::make_static_histogram(custom_axis(0, 3)); + h("-10"); + h("0"); + h("1"); + h("9"); + + for (auto xi : h.axis()) { + std::cout << "bin " << xi.idx() << " [" << xi.lower() << ", " + << xi.upper() << ") " << h.bin(xi).value() + << std::endl; + } + + /* prints: + bin 0 [0, 1) 1 + bin 1 [1, 2] 1 + bin 2 [2, 3] 0 + */ +} + +//] diff --git a/examples/guide_custom_storage.cpp b/examples/guide_custom_storage.cpp new file mode 100644 index 00000000..aaa4a000 --- /dev/null +++ b/examples/guide_custom_storage.cpp @@ -0,0 +1,17 @@ +//[ guide_custom_storage + +#include +#include + +namespace bh = boost::histogram; + +int main() { + // create static histogram with array_storage, using int as counter type + auto h = bh::make_static_histogram_with>( + bh::axis::regular<>(10, 0, 1) + ); + + // do something with h +} + +//] diff --git a/examples/guide_fill_histogram.cpp b/examples/guide_fill_histogram.cpp new file mode 100644 index 00000000..0d4e247e --- /dev/null +++ b/examples/guide_fill_histogram.cpp @@ -0,0 +1,38 @@ +//[ guide_fill_histogram + +#include +#include +#include + +namespace bh = boost::histogram; + +int main() { + auto h = bh::make_static_histogram(bh::axis::integer<>(0, 4), + bh::axis::regular<>(10, 0, 5)); + + // fill histogram, number of arguments must be equal to number of axes + h(0, 1.1); // increases bin counter by one + h(bh::weight(2), 3, 3.4); // increase bin counter by 2 instead of 1 + + // histogram also supports fills from a container of values; a container + // of wrong size trips an assertion in debug mode + auto xy1 = std::make_pair(4, 3.1); + h(xy1); + auto xy2 = std::vector({3.0, 4.9}); + h(xy2); + + // functional-style processing is also supported + std::vector> input_data{ + {0, 1.2}, {2, 3.4}, {4, 5.6} + }; + // Note that std::for_each takes the functor by value, thus it makes a + // potentially expensive copy of your histogram. Passing freshly created + // histograms is ok, though, because of return-value-optimization + auto h2 = std::for_each(input_data.begin(), input_data.end(), + bh::make_static_histogram( + bh::axis::integer<>(0, 4), + bh::axis::regular<>(10, 0, 5))); + // h is filled +} + +//] diff --git a/examples/guide_histogram_operators.cpp b/examples/guide_histogram_operators.cpp new file mode 100644 index 00000000..f9fd53cb --- /dev/null +++ b/examples/guide_histogram_operators.cpp @@ -0,0 +1,57 @@ +//[ guide_histogram_operators + +#include +#include + +namespace bh = boost::histogram; + +int main() { + // make two histograms + auto h1 = bh::make_static_histogram(bh::axis::regular<>(2, -1, 1)); + auto h2 = bh::make_static_histogram(bh::axis::regular<>(2, -1, 1)); + + h1(-0.5); // counts are: 1 0 + h2(0.5); // counts are: 0 1 + + // add them + auto h3 = h1; + h3 += h2; // counts are: 1 1 + + // adding multiple histograms at once is efficient and does not create + // superfluous temporaries since operator+ functions are overloaded to + // accept and return rvalue references where possible + auto h4 = h1 + h2 + h3; // counts are: 2 2 + + std::cout << h4.bin(0).value() << " " << h4.bin(1).value() << std::endl; + // prints: 2 2 + + // multiply by number + h4 *= 2; // counts are: 4 4 + + // divide by number + auto h5 = h4 / 4; // counts are: 1 1 + + std::cout << h5.bin(0).value() << " " << h5.bin(1).value() << std::endl; + // prints: 1 1 + + // compare histograms + std::cout << (h4 == 4 * h5) << " " << (h4 != h5) << std::endl; + // prints: 1 1 + + // note: special effect of multiplication on counter variance + auto h = bh::make_static_histogram(bh::axis::regular<>(2, -1, 1)); + h(-0.5); // counts are: 1 0 + std::cout << "value " << (2 * h).bin(0).value() + << " " << (h + h).bin(0).value() << "\n" + << "variance " << (2 * h).bin(0).variance() + << " " << (h + h).bin(0).variance() << std::endl; + // equality operator also checks variances, so the statement is false + std::cout << (h + h == 2 * h) << std::endl; + /* prints: + value 2 2 + variance 4 2 + 0 + */ +} + +//] diff --git a/examples/guide_histogram_pickle.py b/examples/guide_histogram_pickle.py new file mode 100644 index 00000000..686f9089 --- /dev/null +++ b/examples/guide_histogram_pickle.py @@ -0,0 +1,28 @@ +#[ guide_histogram_pickle + +from __future__ import print_function +import histogram as hg +import pickle + +h1 = hg.histogram(hg.axis.regular(2, -1, 1)) +h2 = hg.histogram(h1) # creates copy +h1(-0.5) +h2(0.5) + +# arithmetic operators (see performance note below) +h3 = h1 + h2 +h4 = h3 * 2 + +print(h4.bin(0).value, h4.bin(1).value) +# prints: 2.0 2.0 + +# now save the histogram +with open("h4_saved.pkl", "w") as f: + pickle.dump(h4, f) +with open("h4_saved.pkl", "r") as f: + h5 = pickle.load(f) + +print(h4 == h5) +# prints: True + +#] diff --git a/examples/guide_histogram_reduction.cpp b/examples/guide_histogram_reduction.cpp new file mode 100644 index 00000000..1e91c8cf --- /dev/null +++ b/examples/guide_histogram_reduction.cpp @@ -0,0 +1,62 @@ +//[ guide_histogram_reduction + +#include +#include + +namespace bh = boost::histogram; + +// example of a generic function for histograms, this one sums all entries +template +typename bh::histogram::element_type sum(const bh::histogram& h) { + auto result = typename bh::histogram::element_type(0); + for (auto x : h) + result += x; + return result; +} + +int main() { + using namespace bh::literals; // enables _c suffix + + // make a 2d histogram + auto h = bh::make_static_histogram(bh::axis::regular<>(3, -1, 1), + bh::axis::integer<>(0, 4)); + + h(-0.9, 0); + h(0.9, 3); + h(0.1, 2); + + auto hr0 = h.reduce_to(0_c); // keep only first axis + auto hr1 = h.reduce_to(1_c); // keep only second axis + + /* + reduce does not remove counts; returned histograms are summed over + the removed axes, so h, hr0, and hr1 have same number of total counts + */ + std::cout << sum(h).value() << " " + << sum(hr0).value() << " " + << sum(hr1).value() << std::endl; + // prints: 3 3 3 + + for (auto yi : h.axis(1_c)) { + for (auto xi : h.axis(0_c)) { + std::cout << h.bin(xi, yi).value() << " "; + } + std::cout << std::endl; + } + // prints: 1 0 0 + // 0 0 0 + // 0 1 0 + // 0 0 1 + + for (auto xi : hr0.axis()) + std::cout << hr0.bin(xi).value() << " "; + std::cout << std::endl; + // prints: 1 1 1 + + for (auto yi : hr1.axis()) + std::cout << hr1.bin(yi).value() << " "; + std::cout << std::endl; + // prints: 1 0 1 1 +} + +//] diff --git a/examples/guide_histogram_serialization.cpp b/examples/guide_histogram_serialization.cpp new file mode 100644 index 00000000..754297bc --- /dev/null +++ b/examples/guide_histogram_serialization.cpp @@ -0,0 +1,42 @@ +//[ guide_histogram_serialization + +#include +#include // includes serialization code +#include +#include +#include + +namespace bh = boost::histogram; + +int main() { + auto a = bh::make_static_histogram(bh::axis::regular<>(3, -1, 1, "r"), + bh::axis::integer<>(0, 2, "i")); + a(0.5, 1); + + std::string buf; // holds persistent representation + + // store histogram + { + std::ostringstream os; + boost::archive::text_oarchive oa(os); + oa << a; + buf = os.str(); + } + + auto b = decltype(a)(); // create a default-constructed second histogram + + std::cout << "before restore " << (a == b) << std::endl; + // prints: before restore 0 + + // load histogram + { + std::istringstream is(buf); + boost::archive::text_iarchive ia(is); + ia >> b; + } + + std::cout << "after restore " << (a == b) << std::endl; + // prints: after restore 1 +} + +//] diff --git a/examples/guide_histogram_streaming.cpp b/examples/guide_histogram_streaming.cpp new file mode 100644 index 00000000..fb357499 --- /dev/null +++ b/examples/guide_histogram_streaming.cpp @@ -0,0 +1,39 @@ +//[ guide_histogram_streaming + +#include +#include +#include + +namespace bh = boost::histogram; + +int main() { + namespace axis = boost::histogram::axis; + + enum { A, B, C }; + + auto h = bh::make_static_histogram( + axis::regular<>(2, -1, 1, "regular1", axis::uoflow::off), + axis::regular(2, 1, 10, "regular2"), + axis::circular<>(4, 0.1, 1.0, "polar"), + axis::variable<>({-1, 0, 1}, "variable", axis::uoflow::off), + axis::category<>({A, B, C}, "category"), + axis::integer<>(-1, 1, "integer", axis::uoflow::off) + ); + + std::cout << h << std::endl; + + /* prints: + + histogram( + regular(2, -1, 1, label='regular1', uoflow=False), + regular_log(2, 1, 10, label='regular2'), + circular(4, phase=0.1, perimeter=1, label='polar'), + variable(-1, 0, 1, label='variable', uoflow=False), + category(0, 1, 2, label='category'), + integer(-1, 1, label='integer', uoflow=False), + ) + + */ +} + +//] diff --git a/examples/guide_make_dynamic_histogram.cpp b/examples/guide_make_dynamic_histogram.cpp new file mode 100644 index 00000000..c8f37d67 --- /dev/null +++ b/examples/guide_make_dynamic_histogram.cpp @@ -0,0 +1,20 @@ +//[ guide_make_dynamic_histogram + +#include +#include + +namespace bh = boost::histogram; + +int main() { + // create vector of axes, axis::any is a polymorphic axis type + auto v = std::vector>(); + v.push_back(bh::axis::regular<>(100, -1, 1)); + v.push_back(bh::axis::integer<>(1, 7)); + + // create dynamic histogram (make_static_histogram be used with iterators) + auto h = bh::make_dynamic_histogram(v.begin(), v.end()); + + // do something with h +} + +//] diff --git a/examples/guide_make_static_histogram.cpp b/examples/guide_make_static_histogram.cpp new file mode 100644 index 00000000..6e353813 --- /dev/null +++ b/examples/guide_make_static_histogram.cpp @@ -0,0 +1,18 @@ +//[ guide_make_static_histogram + +#include + +namespace bh = boost::histogram; + +int main() { + /* + create a 1d-histogram in default configuration which + covers the real line from -1 to 1 in 100 bins, the same + call with `make_dynamic_histogram` would also work + */ + auto h = bh::make_static_histogram(bh::axis::regular<>(100, -1, 1)); + + // do something with h +} + +//] diff --git a/examples/guide_mixed_cpp_python.cpp b/examples/guide_mixed_cpp_python.cpp new file mode 100644 index 00000000..9662b63a --- /dev/null +++ b/examples/guide_mixed_cpp_python.cpp @@ -0,0 +1,21 @@ +//[ guide_mixed_cpp_python_part_cpp + +#include +#include + +namespace bh = boost::histogram; +namespace bp = boost::python; + +// function that runs in C++ and accepts reference to dynamic histogram +void process(bh::dynamic_histogram<>& h) { + // fill histogram, in reality this would be arbitrarily complex code + for (int i = 0; i < 4; ++i) + h(0.25 * i, i); +} + +// a minimal Python module, which exposes the process function to Python +BOOST_PYTHON_MODULE(cpp_filler) { + bp::def("process", process); +} + +//] diff --git a/examples/guide_mixed_cpp_python.py b/examples/guide_mixed_cpp_python.py new file mode 100644 index 00000000..42286208 --- /dev/null +++ b/examples/guide_mixed_cpp_python.py @@ -0,0 +1,23 @@ +#[ guide_mixed_cpp_python_part_py + +from __future__ import print_function +import histogram as bh +import cpp_filler + +h = bh.histogram(bh.axis.regular(4, 0, 1), + bh.axis.integer(0, 4)) + +cpp_filler.process(h) # histogram is filled with input values in C++ + +for iy, y in enumerate(h.axis(1)): + for ix, x in enumerate(h.axis(0)): + print(h(ix, iy).value, end=' ') + print() + +# prints: +# 1.0 0.0 0.0 0.0 +# 0.0 1.0 0.0 0.0 +# 0.0 0.0 1.0 0.0 +# 0.0 0.0 0.0 1.0 + +#] diff --git a/examples/guide_numpy_support.py b/examples/guide_numpy_support.py new file mode 100644 index 00000000..acb3c402 --- /dev/null +++ b/examples/guide_numpy_support.py @@ -0,0 +1,55 @@ +#[ guide_numpy_support + +import histogram as hg +import numpy as np + +# create 2d-histogram with two axes with 10 equidistant bins from -3 to 3 +h = hg.histogram(hg.axis.regular(10, -3, 3, "x"), + hg.axis.regular(10, -3, 3, "y")) + +# generate some numpy arrays with data to fill into histogram, +# in this case normal distributed random numbers in x and y +x = np.random.randn(1000) +y = 0.5 * np.random.randn(1000) + +# fill histogram with numpy arrays, this is very fast +h(x, y) # call looks the same as if x, y were values + +# get representations of the bin edges as Numpy arrays; this representation +# differs from `list(h.axis(0))` as explained in the next example +x = np.array(h.axis(0)) +y = np.array(h.axis(1)) + +# creates a view of the counts (no copy involved) +count_matrix = np.asarray(h) + +# cut off the under- and overflow bins to not confuse matplotib (no copy) +reduced_count_matrix = count_matrix[:-2,:-2] + +try: + # draw the count matrix + import matplotlib.pyplot as plt + plt.pcolor(x, y, reduced_count_matrix.T) + plt.xlabel(h.axis(0).label) + plt.ylabel(h.axis(1).label) + plt.savefig("example_2d_python.png") +except ImportError: + # ok, no matplotlib, then just print the full count matrix + print count_matrix + + # output of the print looks something like this, the two right-most rows + # and two down-most columns represent under-/overflow bins + # [[ 0 0 0 1 5 0 0 1 0 0 0 0] + # [ 0 0 0 1 17 11 6 0 0 0 0 0] + # [ 0 0 0 5 31 26 4 1 0 0 0 0] + # [ 0 0 3 20 59 62 26 4 0 0 0 0] + # [ 0 0 1 26 96 89 16 1 0 0 0 0] + # [ 0 0 4 21 86 84 20 1 0 0 0 0] + # [ 0 0 1 24 71 50 15 2 0 0 0 0] + # [ 0 0 0 6 26 37 7 0 0 0 0 0] + # [ 0 0 0 0 11 10 2 0 0 0 0 0] + # [ 0 0 0 1 2 3 1 0 0 0 0 0] + # [ 0 0 0 0 0 2 0 0 0 0 0 0] + # [ 0 0 0 0 0 1 0 0 0 0 0 0]] + +#] diff --git a/examples/guide_python_axis_representation.py b/examples/guide_python_axis_representation.py new file mode 100644 index 00000000..882e6bb8 --- /dev/null +++ b/examples/guide_python_axis_representation.py @@ -0,0 +1,29 @@ +#[ guide_python_axis_representation + +import histogram as hg +import numpy as np + +ax = hg.axis.regular(5, 0, 1) +xedge1 = np.array(ax) # this is equivalent to... +xedge2 = [] +for idx, (lower, upper) in enumerate(ax): + xedge2.append(lower) + if idx == len(ax)-1: + xedge2.append(upper) + +print xedge1 +print xedge2 + +# prints: +# [ 0. 0.2 0.4 0.6 0.8 1. ] +# [0.0, 0.2, 0.4, 0.6, 0.8, 1.0] + +# sequences constructed from an axis use its iterator, the result differs +xedge3 = list(ax) + +print xedge3 + +# prints: +# [(0., 0.2), (0.2, 0.4), (0.4, 0.6), (0.8, 1.0)] + +#] diff --git a/examples/guide_python_histogram.py b/examples/guide_python_histogram.py new file mode 100644 index 00000000..43ee48d8 --- /dev/null +++ b/examples/guide_python_histogram.py @@ -0,0 +1,32 @@ +#[ guide_python_histogram + +import histogram as hg + +# make 1-d histogram with 5 logarithmic bins from 1e0 to 1e5 +h = hg.histogram(hg.axis.regular_log(5, 1e0, 1e5, "x")) + +# fill histogram with numbers +for x in (2e0, 2e1, 2e2, 2e3, 2e4): + h(x, weight=4) # increment bin counter by 4 + +# iterate over bins and access bin counter +for idx, (lower, upper) in enumerate(h.axis(0)): + print "bin {0} x in [{1}, {2}): {3} +/- {4}".format( + idx, lower, upper, h.bin(idx).value, h.bin(idx).variance ** 0.5) + +# under- and overflow bins are accessed like in C++ +lo, up = h.axis(0)[-1] +print "underflow [{0}, {1}): {2} +/- {3}".format(lo, up, h(-1).value, h(-1).variance) +lo, up = h.axis(0)[5] +print "overflow [{0}, {1}): {2} +/- {3}".format(lo, up, h(5).value, h(5).variance) + +# prints: +# bin 0 x in [1.0, 10.0): 4.0 +/- 4.0 +# bin 1 x in [10.0, 100.0): 4.0 +/- 4.0 +# bin 2 x in [100.0, 1000.0): 4.0 +/- 4.0 +# bin 3 x in [1000.0, 10000.0): 4.0 +/- 4.0 +# bin 4 x in [10000.0, 100000.0): 4.0 +/- 4.0 +# underflow [0.0, 1.0): 0.0 +/- 0.0 +# overflow [100000.0, inf): 0.0 +/- 0.0 + +#] diff --git a/include/boost/histogram/axis/axis.hpp b/include/boost/histogram/axis/axis.hpp index 90064906..02e78d16 100644 --- a/include/boost/histogram/axis/axis.hpp +++ b/include/boost/histogram/axis/axis.hpp @@ -534,7 +534,7 @@ public: } template > + typename = ::boost::histogram::detail::requires_iterator> category(Iterator begin, Iterator end, string_view label = {}) : base_type(std::distance(begin, end), label), map_(new map_type()) { int index = 0; diff --git a/include/boost/histogram/detail/meta.hpp b/include/boost/histogram/detail/meta.hpp index be75a7d9..2324dc07 100644 --- a/include/boost/histogram/detail/meta.hpp +++ b/include/boost/histogram/detail/meta.hpp @@ -25,46 +25,59 @@ #include #include #include +#include namespace boost { namespace histogram { namespace detail { +#define BOOST_HISTOGRAM_MAKE_SFINAE(name, cond) \ +template struct name { \ + template \ + struct SFINAE {}; \ + template static std::true_type Test(SFINAE *); \ + template static std::false_type Test(...); \ + using type = decltype(Test(nullptr)); \ +}; \ +template \ +using name##_t = typename name::type + +BOOST_HISTOGRAM_MAKE_SFINAE(has_variance_support, + (std::declval().value(), std::declval().variance())); + +BOOST_HISTOGRAM_MAKE_SFINAE(has_method_lower, + (std::declval().lower(0))); + +BOOST_HISTOGRAM_MAKE_SFINAE(is_dynamic_container, + (std::begin(std::declval()))); + +BOOST_HISTOGRAM_MAKE_SFINAE(is_static_container, + (std::get<0>(std::declval()))); + +struct static_container_tag {}; +struct dynamic_container_tag {}; +struct no_container_tag {}; + +template +using classify_container_t = + typename std::conditional< + is_static_container_t::value, + static_container_tag, + typename std::conditional< + is_dynamic_container_t::value, + dynamic_container_tag, + no_container_tag + >::type + >::type; + template ().size(), std::declval().increase(0), std::declval()[0])> struct requires_storage {}; -template struct has_variance_support { - template ().value(), - std::declval().variance())> - struct SFINAE {}; - template static std::true_type Test(SFINAE *); - template static std::false_type Test(...); - using type = decltype(Test(nullptr)); -}; - -template -using has_variance_support_t = typename has_variance_support::type; - -template struct has_method_lower { - template ().lower(0))> - struct SFINAE {}; - template static std::true_type Test(SFINAE *); - template static std::false_type Test(...); - using type = decltype(Test(nullptr)); -}; - -template -using has_method_lower_t = typename has_method_lower::type; - template (), ++std::declval())> -struct is_iterator {}; - -template ()), - std::end(std::declval()))> -struct is_sequence {}; +struct requires_iterator {}; template struct union_ @@ -110,6 +123,15 @@ template using axes_select_t = typename mpl::transform>::type; +template +using size_of = std::tuple_size::type>; + +template +using type_of = typename std::tuple_element::type>::type; + +template +using if_else = typename std::conditional::type; + } // namespace detail } // namespace histogram } // namespace boost diff --git a/include/boost/histogram/detail/utility.hpp b/include/boost/histogram/detail/utility.hpp index c2076949..1ef7e0af 100644 --- a/include/boost/histogram/detail/utility.hpp +++ b/include/boost/histogram/detail/utility.hpp @@ -8,6 +8,7 @@ #define _BOOST_HISTOGRAM_DETAIL_UTILITY_HPP_ #include +#include #include #include #include @@ -34,25 +35,8 @@ inline void lin(std::size_t &out, std::size_t &stride, const Axis &a, int j) noexcept { // the following is highly optimized code that runs in a hot loop; // please measure the performance impact of changes - const int uoflow = a.uoflow(); - // set stride to zero if 'j' is not in range, - // this communicates the out-of-range condition to the caller - stride *= (j >= -uoflow) & (j < (a.size() + uoflow)); - j += (j < 0) * (a.size() + 2); // wrap around if in < 0 - out += j * stride; -#ifndef _MSC_VER -#pragma GCC diagnostic ignored "-Wstrict-overflow" -#endif - stride *= a.shape(); -} - -template -inline void xlin(std::size_t &out, std::size_t &stride, const Axis &a, - const typename Axis::value_type &x) noexcept { - // the following is highly optimized code that runs in a hot loop; - // please measure the performance impact of changes - int j = a.index(x); - // j is guaranteed to be in range [-1, bins] + BOOST_ASSERT_MSG(-1 <= j && j <= a.size(), + "index not in range [-1, size], check your logic"); j += (j < 0) * (a.size() + 2); // wrap around if j < 0 out += j * stride; #ifndef _MSC_VER diff --git a/include/boost/histogram/dynamic_histogram.hpp b/include/boost/histogram/dynamic_histogram.hpp index 914250ed..ee612376 100644 --- a/include/boost/histogram/dynamic_histogram.hpp +++ b/include/boost/histogram/dynamic_histogram.hpp @@ -9,6 +9,7 @@ #include #include +#include #include #include #include @@ -34,6 +35,7 @@ #include #include #include +#include // forward declaration for serialization namespace boost { @@ -75,7 +77,7 @@ public: storage_ = Storage(bincount_from_axes()); } - template > + template > histogram(Iterator begin, Iterator end) : axes_(std::distance(begin, end)) { std::copy(begin, end, axes_.begin()); storage_ = Storage(bincount_from_axes()); @@ -136,69 +138,67 @@ public: return *this; } - template void fill(const Args &... args) { - using n_weight = - typename mpl::count_if, detail::is_weight>; - using n_sample = - typename mpl::count_if, detail::is_sample>; - static_assert(n_weight::value <= 1, - "more than one weight argument is not allowed"); - static_assert(n_sample::value <= 1, - "more than one sample argument is not allowed"); - if (dim() != sizeof...(args) - n_weight::value - n_sample::value) - throw std::invalid_argument( - "fill arguments does not match histogram dimension"); - fill_impl(mpl::bool_(), mpl::bool_(), - args...); - } - - template > - void fill(Iterator begin, Iterator end) { - if (dim() != std::distance(begin, end)) - throw std::invalid_argument( - "fill iterator range does not match histogram dimension"); + template void operator()(Ts&&... ts) { + // case with one argument is ambiguous, is specialized below + BOOST_ASSERT_MSG(dim() == sizeof...(Ts), + "fill arguments does not match histogram dimension " + "(did you use weight() in the wrong place?)"); std::size_t idx = 0, stride = 1; - xlin_iter(idx, stride, begin); + xlin<0>(idx, stride, std::forward(ts)...); if (stride) { - storage_.increase(idx); + fill_storage_impl(idx); } } - template > - void fill(Iterator begin, Iterator end, const detail::weight_t &w) { - if (dim() != std::distance(begin, end)) - throw std::invalid_argument( - "fill iterator range does not match histogram dimension"); - std::size_t idx = 0, stride = 1; - xlin_iter(idx, stride, begin); - if (stride) { - storage_.add(idx, w); + template void operator()(T && t) { + // check whether T is unpackable + if (dim() == 1) { + fill_impl(detail::no_container_tag(), std::forward(t)); + } else { + fill_impl(detail::classify_container_t(), std::forward(t)); } } - template - const_reference operator()(Indices &&... indices) const { - if (dim() != sizeof...(indices)) - throw std::invalid_argument( - "value arguments does not match histogram dimension"); + template void operator()(detail::weight&& w, + Ts&&... ts) { + // case with one argument is ambiguous, is specialized below + BOOST_ASSERT_MSG(dim() == sizeof...(Ts), + "fill arguments does not match histogram dimension"); std::size_t idx = 0, stride = 1; - lin<0>(idx, stride, indices...); - if (stride == 0) - throw std::out_of_range("invalid index"); + xlin<0>(idx, stride, std::forward(ts)...); + if (stride) { + fill_storage_impl(idx, std::move(w)); + } + } + + template + void operator()(detail::weight&& w, T&& t) { + // check whether T is unpackable + if (dim() == 1) { + fill_impl(detail::no_container_tag(), std::forward(t), std::move(w)); + } else { + fill_impl(detail::classify_container_t(), std::forward(t), std::move(w)); + } + } + + template + const_reference bin(Ts &&... ts) const { + // case with one argument is ambiguous, is specialized below + BOOST_ASSERT_MSG(dim() == sizeof...(Ts), + "bin arguments does not match histogram dimension"); + std::size_t idx = 0, stride = 1; + lin<0>(idx, stride, static_cast(ts)...); + BOOST_ASSERT_MSG(stride > 0, "invalid index"); return storage_[idx]; } - template > - const_reference operator()(Iterator begin, Iterator end) const { - if (dim() != std::distance(begin, end)) - throw std::invalid_argument( - "iterator range in operator() does not match histogram dimension"); - std::size_t idx = 0, stride = 1; - lin_iter(idx, stride, begin); - if (stride == 0) - throw std::out_of_range("invalid index"); - return storage_[idx]; + template + const_reference bin(T&&t) const { + // check whether T is unpackable + if (dim() == 1) + return bin_impl(detail::no_container_tag(), std::forward(t)); + else + return bin_impl(detail::classify_container_t(), std::forward(t)); } /// Number of axes (dimensions) of histogram @@ -207,16 +207,6 @@ public: /// Total number of bins in the histogram (including underflow/overflow) std::size_t bincount() const noexcept { return storage_.size(); } - /// Sum of all counts in the histogram - element_type sum() const noexcept { - element_type result(0); - // don't use bincount() here, so sum() still works in a moved-from object - for (std::size_t i = 0, n = storage_.size(); i < n; ++i) { - result += storage_[i]; - } - return result; - } - /// Reset bin counters to zero void reset() { storage_ = Storage(bincount_from_axes()); } @@ -258,7 +248,7 @@ public: } /// Return a lower dimensional histogram - template > + template > histogram reduce_to(Iterator begin, Iterator end) const { std::vector b(dim(), false); for (; begin != end; ++begin) @@ -267,11 +257,11 @@ public: } const_iterator begin() const noexcept { - return const_iterator(*this, storage_, true); + return const_iterator(*this, storage_, 0); } const_iterator end() const noexcept { - return const_iterator(*this, storage_); + return const_iterator(*this, storage_, storage_.size()); } private: @@ -284,34 +274,87 @@ private: return v.value; } - template - inline void fill_impl(mpl::false_, mpl::false_, const Args &... args) { + template + void fill_storage_impl(std::size_t idx, detail::weight && w) { + storage_.add(idx, w); + } + + void fill_storage_impl(std::size_t idx) { + storage_.increase(idx); + } + + template + void fill_impl(detail::dynamic_container_tag, T && t, Ts&&... ts) { + BOOST_ASSERT_MSG(dim() == std::distance(std::begin(t), std::end(t)), + "fill container does not match histogram dimension"); std::size_t idx = 0, stride = 1; - xlin<0>(idx, stride, args...); + xlin_iter(idx, stride, std::begin(t)); if (stride) { - storage_.increase(idx); + fill_storage_impl(idx, std::forward(ts)...); } } - template - inline void fill_impl(mpl::true_, mpl::false_, const Args &... args) { + template + void fill_impl(detail::static_container_tag, T && t, Ts&&... ts) { + BOOST_ASSERT_MSG(dim() == detail::size_of::value, + "fill container does not match histogram dimension"); std::size_t idx = 0, stride = 1; - typename mpl::deref, detail::is_weight>::type>::type w; - wxlin<0>(idx, stride, w, args...); + xlin_get(mpl::int_::value>(), idx, stride, std::forward(t)); if (stride) { - storage_.add(idx, w); + fill_storage_impl(idx, std::forward(ts)...); } } - template - inline void fill_impl(mpl::false_, mpl::true_, const Args &... args) { - // not implemented + template + void fill_impl(detail::no_container_tag, T && t, Ts&&... ts) { + BOOST_ASSERT_MSG(dim() == 1, + "fill argument does not match histogram dimension"); + std::size_t idx = 0, stride = 1; + xlin<0>(idx, stride, t); + if (stride) { + fill_storage_impl(idx, std::forward(ts)...); + } } - template - inline void fill_impl(mpl::true_, mpl::true_, const Args &... args) { - // not implemented + template + const_reference bin_impl(detail::dynamic_container_tag, T && t) const { + BOOST_ASSERT_MSG(dim() == std::distance(std::begin(t), std::end(t)), + "bin container does not match histogram dimension"); + std::size_t idx = 0, stride = 1; + lin_iter(idx, stride, std::begin(t)); + BOOST_ASSERT_MSG(stride > 0, "invalid index"); + return storage_[idx]; + } + + template + const_reference bin_impl(detail::static_container_tag, T && t) const { + BOOST_ASSERT_MSG(dim() == detail::size_of::value, + "bin container does not match histogram dimension"); + std::size_t idx = 0, stride = 1; + lin_get(mpl::int_::value>(), idx, stride, std::forward(t)); + BOOST_ASSERT_MSG(stride > 0, "invalid index"); + return storage_[idx]; + } + + template + typename std::enable_if::value, const_reference>::type + bin_impl(detail::no_container_tag, T && t) const { + BOOST_ASSERT_MSG(dim() == 1, + "bin argument does not match histogram dimension"); + std::size_t idx = 0, stride = 1; + lin<0>(idx, stride, static_cast(t)); + BOOST_ASSERT_MSG(stride > 0, "invalid index"); + return storage_[idx]; + } + + template + typename std::enable_if::value), const_reference>::type + bin_impl(detail::no_container_tag, T && t) const { + BOOST_ASSERT_MSG(dim() == 1, + "bin argument does not match histogram dimension"); + BOOST_ASSERT_MSG(false, + "bin argument not convertible to int"); + return storage_[0]; } struct lin_visitor : public static_visitor { @@ -328,11 +371,11 @@ private: template inline void lin(std::size_t &, std::size_t &) const noexcept {} - template - inline void lin(std::size_t &idx, std::size_t &stride, const First &x, - const Rest &... rest) const noexcept { - apply_visitor(lin_visitor{idx, stride, static_cast(x)}, axes_[D]); - lin(idx, stride, rest...); + template + inline void lin(std::size_t &idx, std::size_t &stride, int x, + Ts... ts) const noexcept { + apply_visitor(lin_visitor{idx, stride, x}, axes_[D]); + lin(idx, stride, ts...); } template struct xlin_visitor : public static_visitor { @@ -346,7 +389,7 @@ private: } template void impl(std::true_type, const Axis &a) const { - detail::xlin(idx, stride, a, val); + detail::lin(idx, stride, a, a.index(val)); } template void impl(std::false_type, const Axis &) const { @@ -359,31 +402,11 @@ private: template inline void xlin(std::size_t &, std::size_t &) const {} - template - inline void xlin(std::size_t &idx, std::size_t &stride, const First &first, - const Rest &... rest) const { - apply_visitor(xlin_visitor{idx, stride, first}, axes_[D]); - xlin(idx, stride, rest...); - } - - template - inline void wxlin(std::size_t &, std::size_t &, Weight &) const {} - - // enable_if needed, because gcc thinks the overloads are ambiguous - template - inline typename std::enable_if::value)>::type - wxlin(std::size_t &idx, std::size_t &stride, Weight &w, const First &first, - const Rest &... rest) const { - apply_visitor(xlin_visitor{idx, stride, first}, axes_[D]); - wxlin(idx, stride, w, rest...); - } - - template - inline void wxlin(std::size_t &idx, std::size_t &stride, Weight &w, - const detail::weight_t &first, - const Rest &... rest) const { - w = first; - wxlin(idx, stride, w, rest...); + template + inline void xlin(std::size_t &idx, std::size_t &stride, T &&t, + Ts &&... ts) const { + apply_visitor(xlin_visitor{idx, stride, t}, axes_[D]); + xlin<(D+1)>(idx, stride, std::forward(ts)...); } template @@ -403,6 +426,25 @@ private: } } + template void xlin_get(mpl::int_<0>, std::size_t&, std::size_t &, T&&t) const {} + + template void xlin_get(mpl::int_, std::size_t& idx, + std::size_t & stride, T&&t) const { + constexpr unsigned D = detail::size_of::value - N; + apply_visitor(xlin_visitor>{idx, stride, std::get(t)}, axes_[D]); + xlin_get(mpl::int_<(N-1)>(), idx, stride, std::forward(t)); + } + + template void lin_get(mpl::int_<0>, std::size_t& , + std::size_t & , T&&) const {} + + template void lin_get(mpl::int_, std::size_t& idx, + std::size_t & stride, T&&t) const { + constexpr unsigned D = detail::size_of::value - N; + apply_visitor(lin_visitor{idx, stride, static_cast(std::get(t))}, axes_[D]); + lin_get(mpl::int_<(N-1)>(), idx, stride, std::forward(t)); + } + histogram reduce_impl(const std::vector &b) const { axes_type axes; std::vector n(b.size()); @@ -456,7 +498,7 @@ make_dynamic_histogram_with(Axes &&... axes) { Storage>(std::forward(axes)...); } -template > +template > histogram> make_dynamic_histogram(Iterator begin, Iterator end) { @@ -466,7 +508,7 @@ make_dynamic_histogram(Iterator begin, Iterator end) { begin, end); } -template > +template > histogram, array_storage>> @@ -478,7 +520,7 @@ make_dynamic_weighted_histogram(Iterator begin, Iterator end) { } template > + typename = detail::requires_iterator> histogram, Storage> diff --git a/include/boost/histogram/histogram_fwd.hpp b/include/boost/histogram/histogram_fwd.hpp index d82b6b4d..4b4f0cf8 100644 --- a/include/boost/histogram/histogram_fwd.hpp +++ b/include/boost/histogram/histogram_fwd.hpp @@ -61,21 +61,18 @@ template using static_histogram = histogram; namespace detail { -template struct weight_t { - T value; - operator const T &() const { return value; } -}; -template struct is_weight : std::false_type {}; -template struct is_weight> : std::true_type {}; +template struct weight { T value; }; +// template struct is_weight : std::false_type {}; +// template struct is_weight> : std::true_type {}; -template struct sample_t { T value; }; -template struct is_sample : std::false_type {}; -template struct is_sample> : std::true_type {}; +template struct sample { T value; }; +// template struct is_sample : std::false_type {}; +// template struct is_sample> : std::true_type {}; } // namespace detail -template detail::weight_t weight(T &&t) { return {t}; } +template detail::weight weight(T &&t) { return {t}; } -template detail::sample_t sample(T &&t) { return {t}; } +template detail::sample sample(T &&t) { return {t}; } } // namespace histogram } // namespace boost diff --git a/include/boost/histogram/iterator.hpp b/include/boost/histogram/iterator.hpp index c73f1a68..ca2b8211 100644 --- a/include/boost/histogram/iterator.hpp +++ b/include/boost/histogram/iterator.hpp @@ -25,23 +25,6 @@ struct dim_t { std::size_t stride; }; -std::size_t inc(dim_t *iter, dim_t *end) noexcept { - for (; iter != end; ++iter) { - ++(iter->idx); - if (iter->idx < iter->size) - break; - iter->idx = 0; - } - - if (iter == end) - return std::numeric_limits::max(); - - std::size_t idx = 0; - for (; iter != end; ++iter) - idx += iter->idx * iter->stride; - return idx; -} - struct dim_visitor { mutable std::size_t stride; mutable dim_t *dims; @@ -51,25 +34,72 @@ struct dim_visitor { } }; +void decode(std::size_t idx, unsigned dim, dim_t* dims) { + dims += dim; + while ((--dims, --dim)) { + dims->idx = idx / dims->stride; + idx -= dims->idx * dims->stride; + dims->idx -= (dims->idx > dims->size) * (dims->size + 2); + } + dims->idx = idx; + dims->idx -= (dims->idx > dims->size) * (dims->size + 2); +} + class multi_index { public: - int idx(unsigned dim = 0) const noexcept { return dims_[dim].idx; } - unsigned dim() const noexcept { return dims_.size(); } + unsigned dim() const noexcept { return dim_; } + + int idx(unsigned dim=0) const noexcept { + if (idx_ != last_) { + const_cast(last_) = idx_; + decode(idx_, dim_, const_cast(dims_.get())); + } + return dims_[dim].idx; + } protected: multi_index() = default; template - explicit multi_index(const Histogram &h) : idx_(0), dims_(h.dim()) { - h.for_each_axis(dim_visitor{1, dims_.data()}); + multi_index(const Histogram &h, std::size_t idx) : + dim_(h.dim()), idx_(idx), last_(idx), + dims_(new dim_t[h.dim()]) { + h.for_each_axis(dim_visitor{1, dims_.get()}); + if (0 < idx && idx < h.bincount()) { + decode(idx, dim_, dims_.get()); + } } - void increment() noexcept { - idx_ = inc(dims_.data(), dims_.data() + dims_.size()); + multi_index(const multi_index& o) : + dim_(o.dim_), idx_(o.idx_), last_(o.last_), + dims_(new dim_t[o.dim_]) + { + std::copy(o.dims_.get(), o.dims_.get() + o.dim_, dims_.get()); } + multi_index& operator=(const multi_index& o) { + if (this != &o) { + if (dim_ != o.dim_) { + dims_.reset(new dim_t[o.dim_]); + } + dim_ = o.dim_; + idx_ = o.idx_; + last_ = o.last_; + std::copy(o.dims_.get(), o.dims_.get() + o.dim_, dims_.get()); + } + return *this; + } + + multi_index(multi_index&&) = default; + multi_index& operator=(multi_index&&) = default; + + void increment() noexcept { ++idx_; } + void decrement() noexcept { --idx_; } + + unsigned dim_; std::size_t idx_; - std::vector dims_; + std::size_t last_; + std::unique_ptr dims_; }; } // namespace detail @@ -77,22 +107,17 @@ template class iterator_over : public iterator_facade< iterator_over, typename Storage::element_type, - forward_traversal_tag, typename Storage::const_reference>, + bidirectional_traversal_tag, typename Storage::const_reference>, public detail::multi_index { public: - /// begin iterator - iterator_over(const Histogram &h, const Storage &s, bool) - : detail::multi_index(h), histogram_(h), storage_(s) {} - - /// end iterator - iterator_over(const Histogram &h, const Storage &s) - : histogram_(h), storage_(s) { - idx_ = std::numeric_limits::max(); - } + iterator_over(const Histogram &h, const Storage &s, std::size_t idx) + : detail::multi_index(h, idx), histogram_(h), storage_(s) {} iterator_over(const iterator_over &) = default; iterator_over &operator=(const iterator_over &) = default; + iterator_over(iterator_over &&) = default; + iterator_over &operator=(iterator_over &&) = default; auto bin() const -> decltype(std::declval().axis(mpl::int_<0>())[0]) { @@ -105,7 +130,7 @@ public: return histogram_.axis(mpl::int_())[dims_[N].idx]; } - template + template // TODO: solve this with a std::enable_if auto bin(Int dim) const -> decltype(std::declval().axis(dim)[0]) { return histogram_.axis(dim)[dims_[dim].idx]; @@ -115,7 +140,7 @@ private: bool equal(const iterator_over &rhs) const noexcept { return &storage_ == &rhs.storage_ && idx_ == rhs.idx_; } - typename Storage::const_reference dereference() const { + typename Storage::const_reference dereference() const noexcept { return storage_[idx_]; } diff --git a/include/boost/histogram/static_histogram.hpp b/include/boost/histogram/static_histogram.hpp index 130e8274..b7bcbaae 100644 --- a/include/boost/histogram/static_histogram.hpp +++ b/include/boost/histogram/static_histogram.hpp @@ -7,6 +7,7 @@ #ifndef _BOOST_HISTOGRAM_HISTOGRAM_IMPL_STATIC_HPP_ #define _BOOST_HISTOGRAM_HISTOGRAM_IMPL_STATIC_HPP_ +#include #include #include #include @@ -36,6 +37,7 @@ #include #include #include +#include // forward declaration for serialization namespace boost { @@ -147,48 +149,65 @@ public: return *this; } - template void fill(const Args &... args) { - using n_weight = - typename mpl::count_if, detail::is_weight>; - using n_sample = - typename mpl::count_if, detail::is_sample>; - static_assert(n_weight::value <= 1, - "more than one weight argument is not allowed"); - static_assert(n_sample::value <= 1, - "more than one sample argument is not allowed"); - static_assert(sizeof...(args) == - (axes_size::value + n_weight::value + n_sample::value), - "number of arguments does not match histogram dimension"); - fill_impl(mpl::bool_(), mpl::bool_(), - args...); + template void operator()(Ts &&... ts) { + // case with one argument is ambiguous, is specialized below + static_assert(sizeof...(Ts) == axes_size::value, + "fill arguments do not match histogram dimension"); + std::size_t idx = 0, stride = 1; + xlin<0>(idx, stride, ts...); + if (stride) + fill_storage_impl(idx); + } + + template void operator()(T && t) { + // check whether we need to unpack argument + fill_impl(detail::if_else<(axes_size::value == 1), + detail::no_container_tag, + detail::classify_container_t>(), std::forward(t)); + } + + template void operator()(detail::weight&& w, + Ts &&... ts) { + // case with one argument is ambiguous, is specialized below + std::size_t idx = 0, stride = 1; + xlin<0>(idx, stride, ts...); + if (stride) + fill_storage_impl(idx, std::move(w)); + } + + template void operator()(detail::weight&& w, + T&&t) { + // check whether we need to unpack argument + fill_impl(detail::if_else<(axes_size::value == 1), + detail::no_container_tag, + detail::classify_container_t>(), std::forward(t), std::move(w)); } template - const_reference operator()(const Indices &... indices) const { + const_reference bin(Indices &&... indices) const { + // case with one argument is ambiguous, is specialized below static_assert(sizeof...(indices) == axes_size::value, - "number of arguments does not match histogram dimension"); + "bin arguments do not match histogram dimension"); std::size_t idx = 0, stride = 1; - lin<0>(idx, stride, indices...); - if (stride == 0) { - throw std::out_of_range("invalid index"); - } + lin<0>(idx, stride, static_cast(indices)...); + BOOST_ASSERT_MSG(stride > 0, "invalid index"); return storage_[idx]; } + template + const_reference bin(T&&t) const { + // check whether we need to unpack argument + return bin_impl(detail::if_else<(axes_size::value == 1), + detail::no_container_tag, + detail::classify_container_t>(), std::forward(t)); + } + /// Number of axes (dimensions) of histogram constexpr unsigned dim() const noexcept { return axes_size::value; } /// Total number of bins in the histogram (including underflow/overflow) std::size_t bincount() const noexcept { return storage_.size(); } - /// Sum of all counts in the histogram - element_type sum() const noexcept { - element_type result(0); - for (std::size_t i = 0, n = storage_.size(); i < n; ++i) - result += storage_[i]; - return result; - } - /// Reset bin counters to zero void reset() { storage_ = Storage(bincount_from_axes()); } @@ -246,11 +265,11 @@ public: } const_iterator begin() const noexcept { - return const_iterator(*this, storage_, true); + return const_iterator(*this, storage_, 0); } const_iterator end() const noexcept { - return const_iterator(*this, storage_); + return const_iterator(*this, storage_, storage_.size()); } private: @@ -263,72 +282,139 @@ private: return fc.value; } - template - inline void fill_impl(mpl::false_, mpl::false_, const Args &... args) { + template + void fill_storage_impl(std::size_t idx, detail::weight&& w) { + storage_.add(idx, w); + } + + void fill_storage_impl(std::size_t idx) { + storage_.increase(idx); + } + + template + void fill_impl(detail::dynamic_container_tag, T&&t, Ts&&...ts) { + BOOST_ASSERT_MSG(std::distance(std::begin(t), std::end(t)) == axes_size::value, + "fill container does not match histogram dimension"); std::size_t idx = 0, stride = 1; - xlin<0>(idx, stride, args...); + xlin_iter(mpl::int_(), idx, stride, std::begin(t)); if (stride) { - storage_.increase(idx); + fill_storage_impl(idx, std::forward(ts)...); } } - template - inline void fill_impl(mpl::true_, mpl::false_, const Args &... args) { + template + void fill_impl(detail::static_container_tag, T&&t, Ts&&...ts) { + static_assert(detail::size_of::value == axes_size::value, + "fill container does not match histogram dimension"); std::size_t idx = 0, stride = 1; - typename mpl::deref, detail::is_weight>::type>::type w; - wxlin<0>(idx, stride, w, args...); - if (stride) - storage_.add(idx, w); + xlin_get(mpl::int_(), idx, stride, std::forward(t)); + if (stride) { + fill_storage_impl(idx, std::forward(ts)...); + } } - template - inline void fill_impl(mpl::false_, mpl::true_, const Args &... args) { - // not implemented + template + void fill_impl(detail::no_container_tag, T&&t, Ts&&... ts) { + static_assert(axes_size::value == 1, + "fill argument does not match histogram dimension"); + std::size_t idx = 0, stride = 1; + xlin<0>(idx, stride, std::forward(t)); + if (stride) { + fill_storage_impl(idx, std::forward(ts)...); + } } - template - inline void fill_impl(mpl::true_, mpl::true_, const Args &... args) { - // not implemented + template + const_reference bin_impl(detail::dynamic_container_tag, T && t) const { + BOOST_ASSERT_MSG(std::distance(std::begin(t), std::end(t)) == axes_size::value, + "bin container does not match histogram dimension"); + std::size_t idx = 0, stride = 1; + lin_iter(mpl::int_(), idx, stride, std::begin(t)); + BOOST_ASSERT_MSG(stride > 0, "invalid index"); + return storage_[idx]; + } + + template + const_reference bin_impl(detail::static_container_tag, T&&t) const { + static_assert(detail::size_of::value == axes_size::value, + "bin container does not match histogram dimension"); + std::size_t idx = 0, stride = 1; + lin_get(mpl::int_(), idx, stride, std::forward(t)); + BOOST_ASSERT_MSG(stride > 0, "invalid index"); + return storage_[idx]; + } + + template + const_reference bin_impl(detail::no_container_tag, T&&t) const { + static_assert(axes_size::value == 1, + "bin argument does not match histogram dimension"); + std::size_t idx = 0, stride = 1; + lin<0>(idx, stride, static_cast(t)); + BOOST_ASSERT_MSG(stride > 0, "invalid index"); + return storage_[idx]; + } + + template inline void xlin(std::size_t &, std::size_t &) const {} + + template + inline void xlin(std::size_t &idx, std::size_t &stride, T&&t, + Ts&&... ts) const { + detail::lin(idx, stride, fusion::at_c(axes_), + fusion::at_c(axes_).index(t)); + xlin(idx, stride, std::forward(ts)...); + } + + template + void xlin_iter(mpl::int_<0>, std::size_t &, std::size_t &, Iterator ) const {} + + template + void xlin_iter(mpl::int_, std::size_t &idx, std::size_t &stride, Iterator iter) const { + constexpr unsigned D = axes_size::value - N; + detail::lin(idx, stride, fusion::at_c(axes_), + fusion::at_c(axes_).index(*iter)); + xlin_iter(mpl::int_(), idx, stride, ++iter); } template inline void lin(std::size_t &, std::size_t &) const noexcept {} - template - inline void lin(std::size_t &idx, std::size_t &stride, const First &x, - const Rest &... rest) const noexcept { - detail::lin(idx, stride, fusion::at_c(axes_), static_cast(x)); - lin(idx, stride, rest...); + template + inline void lin(std::size_t &idx, std::size_t &stride, int x, + Ts... ts) const noexcept { + detail::lin(idx, stride, fusion::at_c(axes_), x); + lin<(D+1)>(idx, stride, ts...); } - template inline void xlin(std::size_t &, std::size_t &) const {} + template + void lin_iter(mpl::int_<0>, std::size_t &, std::size_t &, Iterator) const {} - template - inline void xlin(std::size_t &idx, std::size_t &stride, const First &first, - const Rest &... rest) const { - detail::xlin(idx, stride, fusion::at_c(axes_), first); - xlin(idx, stride, rest...); + template + void lin_iter(mpl::int_, std::size_t &idx, std::size_t &stride, + Iterator iter) const { + constexpr unsigned D = axes_size::value - N; + detail::lin(idx, stride, fusion::at_c(axes_), static_cast(*iter)); + lin_iter(mpl::int_<(N-1)>(), idx, stride, ++iter); } - template - inline void wxlin(std::size_t &, std::size_t &, Weight &) const {} + template + void xlin_get(mpl::int_<0>, std::size_t &, std::size_t &, T &&) const {} - // enable_if needed, because gcc thinks the overloads are ambiguous - template - inline typename std::enable_if::value)>::type - wxlin(std::size_t &idx, std::size_t &stride, Weight &w, const First &first, - const Rest &... rest) const { - detail::xlin(idx, stride, fusion::at_c(axes_), first); - wxlin(idx, stride, w, rest...); + template + void xlin_get(mpl::int_, std::size_t &idx, std::size_t &stride, T && t) const { + constexpr unsigned D = detail::size_of::value - N; + detail::lin(idx, stride, fusion::at_c(axes_), + fusion::at_c(axes_).index(std::get(t))); + xlin_get(mpl::int_<(N-1)>(), idx, stride, std::forward(t)); } - template - inline void wxlin(std::size_t &idx, std::size_t &stride, Weight &w, - const detail::weight_t &first, - const Rest &... rest) const { - w = first; - wxlin(idx, stride, w, rest...); + template + void lin_get(mpl::int_<0>, std::size_t &, std::size_t &, T&&) const {} + + template + void lin_get(mpl::int_, std::size_t &idx, std::size_t &stride, T&&t) const { + constexpr unsigned D = detail::size_of::value - N; + detail::lin(idx, stride, fusion::at_c(axes_), static_cast(std::get(t))); + lin_get(mpl::int_<(N-1)>(), idx, stride, std::forward(t)); } struct shape_assign_visitor { diff --git a/include/boost/histogram/storage/adaptive_storage.hpp b/include/boost/histogram/storage/adaptive_storage.hpp index 221935d4..5a579d3c 100644 --- a/include/boost/histogram/storage/adaptive_storage.hpp +++ b/include/boost/histogram/storage/adaptive_storage.hpp @@ -311,10 +311,10 @@ template <> struct radd_visitor : public static_visitor { }; template <> -struct radd_visitor> : public static_visitor { +struct radd_visitor> : public static_visitor { any_array &lhs_any; const std::size_t idx; - const weight_t rhs; + const weight rhs; radd_visitor(any_array &l, const std::size_t i, const double w) : lhs_any(l), idx(i), rhs{w} {} @@ -446,9 +446,9 @@ public: apply_visitor(detail::radd_visitor(buffer_, i, t), buffer_); } - template void add(std::size_t i, const detail::weight_t &w) { + template void add(std::size_t i, const detail::weight &w) { apply_visitor( - detail::radd_visitor>(buffer_, i, w.value), + detail::radd_visitor>(buffer_, i, w.value), buffer_); } diff --git a/include/boost/histogram/storage/operators.hpp b/include/boost/histogram/storage/operators.hpp index a2af99f6..d26c4590 100644 --- a/include/boost/histogram/storage/operators.hpp +++ b/include/boost/histogram/storage/operators.hpp @@ -11,50 +11,16 @@ namespace boost { namespace histogram { -namespace detail { -template -bool equal_impl(std::true_type, std::true_type, const S1 &s1, const S2 &s2) { - for (std::size_t i = 0, n = s1.size(); i < n; ++i) { - auto x1 = s1[i]; - auto x2 = s2[i]; - if (x1.value() != x2.value() || x1.variance() != x2.variance()) - return false; - } - return true; -} - -template -bool equal_impl(std::true_type, std::false_type, const S1 &s1, const S2 &s2) { - for (std::size_t i = 0, n = s1.size(); i < n; ++i) { - auto x1 = s1[i]; - auto x2 = s2[i]; - if (x1.value() != x2 || x1.value() != x1.variance()) - return false; - } - return true; -} - -template -bool equal_impl(std::false_type, std::true_type, const S1 &s1, const S2 &s2) { - return equal_impl(std::true_type(), std::false_type(), s2, s1); -} - -template -bool equal_impl(std::false_type, std::false_type, const S1 &s1, const S2 &s2) { - for (std::size_t i = 0, n = s1.size(); i < n; ++i) - if (s1[i] != s2[i]) - return false; - return true; -} -} // namespace detail template , typename = detail::requires_storage> bool operator==(const S1 &s1, const S2 &s2) noexcept { if (s1.size() != s2.size()) return false; - return detail::equal_impl(detail::has_variance_support_t(), - detail::has_variance_support_t(), s1, s2); + for (std::size_t i = 0, n = s1.size(); i < n; ++i) + if (s1[i] != s2[i]) + return false; + return true; } template , diff --git a/include/boost/histogram/storage/weight_counter.hpp b/include/boost/histogram/storage/weight_counter.hpp index 2b137fc2..b9e60a37 100644 --- a/include/boost/histogram/storage/weight_counter.hpp +++ b/include/boost/histogram/storage/weight_counter.hpp @@ -37,6 +37,7 @@ public: return *this; } + // TODO: explain why this is needed weight_counter &operator+=(const RealType &x) { w += x; w2 += x; @@ -51,7 +52,7 @@ public: } template - weight_counter &operator+=(const detail::weight_t &rhs) { + weight_counter &operator+=(const detail::weight &rhs) { const auto x = static_cast(rhs.value); w += x; w2 += x * x; @@ -127,6 +128,23 @@ bool operator!=(const T &t, const weight_counter &w) { return !(w == t); } +template +weight_counter& operator+(const weight_counter& a, const weight_counter& b) { + weight_counter c = a; + return c += b; +} + +template +weight_counter&& operator+(weight_counter&& a, const weight_counter& b) { + a += b; + return std::move(a); +} + +template +weight_counter&& operator+(const weight_counter& a, weight_counter&& b) { + return operator+(std::move(b), a); +} + } // namespace histogram } // namespace boost diff --git a/test/array_storage_test.cpp b/test/array_storage_test.cpp index a701a7cf..65c669bd 100644 --- a/test/array_storage_test.cpp +++ b/test/array_storage_test.cpp @@ -35,7 +35,7 @@ int main() { c.increase(0); d.increase(0); d.add(1, 5); - d.add(0, weight(2)); + d.add(0, 2); BOOST_TEST_EQ(a[0], 1); BOOST_TEST_EQ(b[0], 1); BOOST_TEST_EQ(c[0], 2); diff --git a/test/detail_test.cpp b/test/detail_test.cpp index aaec9a69..b3cd6aed 100644 --- a/test/detail_test.cpp +++ b/test/detail_test.cpp @@ -5,6 +5,7 @@ // or copy at http://www.boost.org/LICENSE_1_0.txt) #include +#include #include #include #include @@ -17,6 +18,7 @@ #include #include #include +#include using namespace boost::mpl; using namespace boost::histogram::detail; @@ -140,5 +142,29 @@ int main() { true); } + // classify_container + { + using result1 = classify_container_t; + BOOST_TEST_TRAIT_TRUE(( std::is_same )); + + using result1a = classify_container_t; + BOOST_TEST_TRAIT_TRUE(( std::is_same )); + + using result2 = classify_container_t>; + BOOST_TEST_TRAIT_TRUE(( std::is_same )); + + using result2a = classify_container_t&>; + BOOST_TEST_TRAIT_TRUE(( std::is_same )); + + using result3 = classify_container_t>; + BOOST_TEST_TRAIT_TRUE(( std::is_same )); + + using result3a = classify_container_t&>; + BOOST_TEST_TRAIT_TRUE(( std::is_same )); + + using result4 = classify_container_t; + BOOST_TEST_TRAIT_TRUE(( std::is_same )); + } + return boost::report_errors(); } diff --git a/test/histogram_test.cpp b/test/histogram_test.cpp index d52a67c1..f3cf58f4 100644 --- a/test/histogram_test.cpp +++ b/test/histogram_test.cpp @@ -22,6 +22,9 @@ #include #include #include +#include +#include +#include using namespace boost::histogram; using namespace boost::histogram::literals; // to get _c suffix @@ -55,6 +58,16 @@ int expected_moved_from_dim(static_tag, int static_value) { int expected_moved_from_dim(dynamic_tag, int) { return 0; } +template +typename Histogram::element_type sum(const Histogram& h) { + // auto result = typename Histogram::element_type(0); + // for (auto x : h) + // result += x; + // return result; + return std::accumulate(h.begin(), h.end(), + typename Histogram::element_type(0)); +} + template void pass_histogram(boost::histogram::histogram &h) {} @@ -132,7 +145,7 @@ template void run_tests() { { auto h = make_histogram(Type(), axis::integer<>{0, 2}, axis::integer<>{0, 3}); - h.fill(0, 0); + h(0, 0); auto h2 = decltype(h)(h); BOOST_TEST(h2 == h); auto h3 = static_histogram, axis::integer<>>, @@ -144,7 +157,7 @@ template void run_tests() { { auto h = make_histogram(Type(), axis::integer<>(0, 1), axis::integer<>(0, 2)); - h.fill(0, 0); + h(0, 0); auto h2 = decltype(h)(); BOOST_TEST_NE(h, h2); h2 = h; @@ -162,19 +175,19 @@ template void run_tests() { { auto h = make_histogram(Type(), axis::integer<>(0, 1), axis::integer<>(0, 2)); - h.fill(0, 0); + h(0, 0); const auto href = h; decltype(h) h2(std::move(h)); // static axes cannot shrink to zero BOOST_TEST_EQ(h.dim(), expected_moved_from_dim(Type(), 2)); - BOOST_TEST_EQ(h.sum().value(), 0); + BOOST_TEST_EQ(sum(h).value(), 0); BOOST_TEST_EQ(h.bincount(), 0); BOOST_TEST_EQ(h2, href); decltype(h) h3; h3 = std::move(h2); // static axes cannot shrink to zero BOOST_TEST_EQ(h2.dim(), expected_moved_from_dim(Type(), 2)); - BOOST_TEST_EQ(h2.sum().value(), 0); + BOOST_TEST_EQ(sum(h2).value(), 0); BOOST_TEST_EQ(h2.bincount(), 0); BOOST_TEST_EQ(h3, href); } @@ -229,13 +242,13 @@ template void run_tests() { auto d = make_histogram(Type(), axis::regular<>(2, 0, 1)); BOOST_TEST(c != d); BOOST_TEST(d != c); - c.fill(0); + c(0); BOOST_TEST(a != c); BOOST_TEST(c != a); - a.fill(0); + a(0); BOOST_TEST(a == c); BOOST_TEST(c == a); - a.fill(0); + a(0); BOOST_TEST(a != c); BOOST_TEST(c != a); } @@ -243,85 +256,79 @@ template void run_tests() { // d1 { auto h = make_histogram(Type(), axis::integer<>{0, 2}); - h.fill(0); - h.fill(0); - h.fill(-1); - h.fill(10); + h(0); + h(0); + h(-1); + h(10); BOOST_TEST_EQ(h.dim(), 1); BOOST_TEST_EQ(h.axis(0_c).size(), 2); BOOST_TEST_EQ(h.axis(0_c).shape(), 4); - BOOST_TEST_EQ(h.sum(), 4); + BOOST_TEST_EQ(sum(h), 4); - BOOST_TEST_THROWS(h(-2), std::out_of_range); - BOOST_TEST_EQ(h(-1), 1); - BOOST_TEST_EQ(h(0), 2); - BOOST_TEST_EQ(h(1), 0); - BOOST_TEST_EQ(h(2), 1); - BOOST_TEST_THROWS(h(3), std::out_of_range); + BOOST_TEST_EQ(h.bin(-1), 1); + BOOST_TEST_EQ(h.bin(0), 2); + BOOST_TEST_EQ(h.bin(1), 0); + BOOST_TEST_EQ(h.bin(2), 1); } // d1_2 { auto h = make_histogram( Type(), axis::integer<>(0, 2, "", axis::uoflow::off)); - h.fill(0); - h.fill(-0); - h.fill(-1); - h.fill(10); + h(0); + h(-0); + h(-1); + h(10); BOOST_TEST_EQ(h.dim(), 1); BOOST_TEST_EQ(h.axis(0_c).size(), 2); BOOST_TEST_EQ(h.axis(0_c).shape(), 2); - BOOST_TEST_EQ(h.sum(), 2); + BOOST_TEST_EQ(sum(h), 2); - BOOST_TEST_THROWS(h(-1), std::out_of_range); - BOOST_TEST_EQ(h(0), 2); - BOOST_TEST_EQ(h(1), 0); - BOOST_TEST_THROWS(h(2), std::out_of_range); + BOOST_TEST_EQ(h.bin(0), 2); + BOOST_TEST_EQ(h.bin(1), 0); } // d1_3 { auto h = make_histogram( Type(), axis::category({"A", "B"})); - h.fill("A"); - h.fill("B"); - h.fill("D"); - h.fill("E"); + h("A"); + h("B"); + h("D"); + h("E"); BOOST_TEST_EQ(h.dim(), 1); BOOST_TEST_EQ(h.axis(0_c).size(), 2); BOOST_TEST_EQ(h.axis(0_c).shape(), 2); - BOOST_TEST_EQ(h.sum(), 2); + BOOST_TEST_EQ(sum(h), 2); - BOOST_TEST_THROWS(h(-1), std::out_of_range); - BOOST_TEST_EQ(h(0), 1); - BOOST_TEST_EQ(h(1), 1); - BOOST_TEST_THROWS(h(2), std::out_of_range); + BOOST_TEST_EQ(h.bin(0), 1); + BOOST_TEST_EQ(h.bin(1), 1); } // d1w { auto h = make_histogram(Type(), axis::integer<>(0, 2)); - h.fill(-1); - h.fill(0); - h.fill(weight(0.5), 0); - h.fill(1); - h.fill(2, weight(2)); + h(-1); + h(0); + h(weight(0.5), 0); + h(1); + h(weight(2), 2); - BOOST_TEST_EQ(h.sum().value(), 5.5); - BOOST_TEST_EQ(h.sum().variance(), 7.25); + BOOST_TEST_EQ(sum(h).value(), 5.5); + BOOST_TEST_EQ(sum(h).variance(), 7.25); - BOOST_TEST_EQ(h(-1).value(), 1); - BOOST_TEST_EQ(h(0).value(), 1.5); - BOOST_TEST_EQ(h(1).value(), 1); - BOOST_TEST_EQ(h(2).value(), 2); + BOOST_TEST_EQ(h.bin(-1).value(), 1); + BOOST_TEST_EQ(h.bin(0).value(), 1.5); + BOOST_TEST_EQ(h.bin(1).value(), 1); + BOOST_TEST_EQ(h.bin(2).value(), 2); - BOOST_TEST_EQ(h(-1).variance(), 1); - BOOST_TEST_EQ(h(0).variance(), 1.25); - BOOST_TEST_EQ(h(1).variance(), 1); - BOOST_TEST_EQ(h(2).variance(), 4); + BOOST_TEST_EQ(h.bin(-1).variance(), 1); + BOOST_TEST_EQ(h.bin(0).variance(), 1.25); + BOOST_TEST_EQ(h.bin(1).variance(), 1); + BOOST_TEST_EQ(h.bin(2).variance(), 4); } // d2 @@ -329,33 +336,33 @@ template void run_tests() { auto h = make_histogram( Type(), axis::regular<>(2, -1, 1), axis::integer<>(-1, 2, "", axis::uoflow::off)); - h.fill(-1, -1); - h.fill(-1, 0); - h.fill(-1, -10); - h.fill(-10, 0); + h(-1, -1); + h(-1, 0); + h(-1, -10); + h(-10, 0); BOOST_TEST_EQ(h.dim(), 2); BOOST_TEST_EQ(h.axis(0_c).size(), 2); BOOST_TEST_EQ(h.axis(0_c).shape(), 4); BOOST_TEST_EQ(h.axis(1_c).size(), 3); BOOST_TEST_EQ(h.axis(1_c).shape(), 3); - BOOST_TEST_EQ(h.sum(), 3); + BOOST_TEST_EQ(sum(h), 3); - BOOST_TEST_EQ(h(-1, 0), 0); - BOOST_TEST_EQ(h(-1, 1), 1); - BOOST_TEST_EQ(h(-1, 2), 0); + BOOST_TEST_EQ(h.bin(-1, 0), 0); + BOOST_TEST_EQ(h.bin(-1, 1), 1); + BOOST_TEST_EQ(h.bin(-1, 2), 0); - BOOST_TEST_EQ(h(0, 0), 1); - BOOST_TEST_EQ(h(0, 1), 1); - BOOST_TEST_EQ(h(0, 2), 0); + BOOST_TEST_EQ(h.bin(0, 0), 1); + BOOST_TEST_EQ(h.bin(0, 1), 1); + BOOST_TEST_EQ(h.bin(0, 2), 0); - BOOST_TEST_EQ(h(1, 0), 0); - BOOST_TEST_EQ(h(1, 1), 0); - BOOST_TEST_EQ(h(1, 2), 0); + BOOST_TEST_EQ(h.bin(1, 0), 0); + BOOST_TEST_EQ(h.bin(1, 1), 0); + BOOST_TEST_EQ(h.bin(1, 2), 0); - BOOST_TEST_EQ(h(2, 0), 0); - BOOST_TEST_EQ(h(2, 1), 0); - BOOST_TEST_EQ(h(2, 2), 0); + BOOST_TEST_EQ(h.bin(2, 0), 0); + BOOST_TEST_EQ(h.bin(2, 1), 0); + BOOST_TEST_EQ(h.bin(2, 2), 0); } // d2w @@ -363,45 +370,45 @@ template void run_tests() { auto h = make_histogram( Type(), axis::regular<>(2, -1, 1), axis::integer<>(-1, 2, "", axis::uoflow::off)); - h.fill(-1, 0); // -> 0, 1 - h.fill(weight(10), -1, -1); // -> 0, 0 - h.fill(weight(5), -1, -10); // is ignored - h.fill(weight(7), -10, 0); // -> -1, 1 + h(-1, 0); // -> 0, 1 + h(weight(10), -1, -1); // -> 0, 0 + h(weight(5), -1, -10); // is ignored + h(weight(7), -10, 0); // -> -1, 1 - BOOST_TEST_EQ(h.sum().value(), 18); - BOOST_TEST_EQ(h.sum().variance(), 150); + BOOST_TEST_EQ(sum(h).value(), 18); + BOOST_TEST_EQ(sum(h).variance(), 150); - BOOST_TEST_EQ(h(-1, 0).value(), 0); - BOOST_TEST_EQ(h(-1, 1).value(), 7); - BOOST_TEST_EQ(h(-1, 2).value(), 0); + BOOST_TEST_EQ(h.bin(-1, 0).value(), 0); + BOOST_TEST_EQ(h.bin(-1, 1).value(), 7); + BOOST_TEST_EQ(h.bin(-1, 2).value(), 0); - BOOST_TEST_EQ(h(0, 0).value(), 10); - BOOST_TEST_EQ(h(0, 1).value(), 1); - BOOST_TEST_EQ(h(0, 2).value(), 0); + BOOST_TEST_EQ(h.bin(0, 0).value(), 10); + BOOST_TEST_EQ(h.bin(0, 1).value(), 1); + BOOST_TEST_EQ(h.bin(0, 2).value(), 0); - BOOST_TEST_EQ(h(1, 0).value(), 0); - BOOST_TEST_EQ(h(1, 1).value(), 0); - BOOST_TEST_EQ(h(1, 2).value(), 0); + BOOST_TEST_EQ(h.bin(1, 0).value(), 0); + BOOST_TEST_EQ(h.bin(1, 1).value(), 0); + BOOST_TEST_EQ(h.bin(1, 2).value(), 0); - BOOST_TEST_EQ(h(2, 0).value(), 0); - BOOST_TEST_EQ(h(2, 1).value(), 0); - BOOST_TEST_EQ(h(2, 2).value(), 0); + BOOST_TEST_EQ(h.bin(2, 0).value(), 0); + BOOST_TEST_EQ(h.bin(2, 1).value(), 0); + BOOST_TEST_EQ(h.bin(2, 2).value(), 0); - BOOST_TEST_EQ(h(-1, 0).variance(), 0); - BOOST_TEST_EQ(h(-1, 1).variance(), 49); - BOOST_TEST_EQ(h(-1, 2).variance(), 0); + BOOST_TEST_EQ(h.bin(-1, 0).variance(), 0); + BOOST_TEST_EQ(h.bin(-1, 1).variance(), 49); + BOOST_TEST_EQ(h.bin(-1, 2).variance(), 0); - BOOST_TEST_EQ(h(0, 0).variance(), 100); - BOOST_TEST_EQ(h(0, 1).variance(), 1); - BOOST_TEST_EQ(h(0, 2).variance(), 0); + BOOST_TEST_EQ(h.bin(0, 0).variance(), 100); + BOOST_TEST_EQ(h.bin(0, 1).variance(), 1); + BOOST_TEST_EQ(h.bin(0, 2).variance(), 0); - BOOST_TEST_EQ(h(1, 0).variance(), 0); - BOOST_TEST_EQ(h(1, 1).variance(), 0); - BOOST_TEST_EQ(h(1, 2).variance(), 0); + BOOST_TEST_EQ(h.bin(1, 0).variance(), 0); + BOOST_TEST_EQ(h.bin(1, 1).variance(), 0); + BOOST_TEST_EQ(h.bin(1, 2).variance(), 0); - BOOST_TEST_EQ(h(2, 0).variance(), 0); - BOOST_TEST_EQ(h(2, 1).variance(), 0); - BOOST_TEST_EQ(h(2, 2).variance(), 0); + BOOST_TEST_EQ(h.bin(2, 0).variance(), 0); + BOOST_TEST_EQ(h.bin(2, 1).variance(), 0); + BOOST_TEST_EQ(h.bin(2, 2).variance(), 0); } // d3w @@ -412,7 +419,7 @@ template void run_tests() { for (auto i = 0; i < h.axis(0_c).size(); ++i) { for (auto j = 0; j < h.axis(1_c).size(); ++j) { for (auto k = 0; k < h.axis(2_c).size(); ++k) { - h.fill(weight(i + j + k), i, j, k); + h(weight(i + j + k), i, j, k); } } } @@ -420,8 +427,8 @@ template void run_tests() { for (auto i = 0; i < h.axis(0_c).size(); ++i) { for (auto j = 0; j < h.axis(1_c).size(); ++j) { for (auto k = 0; k < h.axis(2_c).size(); ++k) { - BOOST_TEST_EQ(h(i, j, k).value(), i + j + k); - BOOST_TEST_EQ(h(i, j, k).variance(), (i + j + k) * (i + j + k)); + BOOST_TEST_EQ(h.bin(i, j, k).value(), i + j + k); + BOOST_TEST_EQ(h.bin(i, j, k).variance(), (i + j + k) * (i + j + k)); } } } @@ -432,22 +439,22 @@ template void run_tests() { auto a = make_histogram(Type(), axis::integer<>(-1, 2)); auto b = make_histogram>(Type(), axis::integer<>(-1, 2)); - a.fill(-1); - b.fill(1); + a(-1); + b(1); auto c = a; c += b; - BOOST_TEST_EQ(c(-1), 0); - BOOST_TEST_EQ(c(0), 1); - BOOST_TEST_EQ(c(1), 0); - BOOST_TEST_EQ(c(2), 1); - BOOST_TEST_EQ(c(3), 0); + BOOST_TEST_EQ(c.bin(-1), 0); + BOOST_TEST_EQ(c.bin(0), 1); + BOOST_TEST_EQ(c.bin(1), 0); + BOOST_TEST_EQ(c.bin(2), 1); + BOOST_TEST_EQ(c.bin(3), 0); auto d = a; d += b; - BOOST_TEST_EQ(d(-1), 0); - BOOST_TEST_EQ(d(0), 1); - BOOST_TEST_EQ(d(1), 0); - BOOST_TEST_EQ(d(2), 1); - BOOST_TEST_EQ(d(3), 0); + BOOST_TEST_EQ(d.bin(-1), 0); + BOOST_TEST_EQ(d.bin(0), 1); + BOOST_TEST_EQ(d.bin(1), 0); + BOOST_TEST_EQ(d.bin(2), 1); + BOOST_TEST_EQ(d.bin(3), 0); } // add_2 @@ -455,26 +462,26 @@ template void run_tests() { auto a = make_histogram(Type(), axis::integer<>(0, 2)); auto b = make_histogram(Type(), axis::integer<>(0, 2)); - a.fill(0); - BOOST_TEST_EQ(a(0).variance(), 1); - b.fill(1, weight(3)); - BOOST_TEST_EQ(b(1).variance(), 9); + a(0); + BOOST_TEST_EQ(a.bin(0).variance(), 1); + b(weight(3), 1); + BOOST_TEST_EQ(b.bin(1).variance(), 9); auto c = a; c += b; - BOOST_TEST_EQ(c(-1).value(), 0); - BOOST_TEST_EQ(c(0).value(), 1); - BOOST_TEST_EQ(c(0).variance(), 1); - BOOST_TEST_EQ(c(1).value(), 3); - BOOST_TEST_EQ(c(1).variance(), 9); - BOOST_TEST_EQ(c(2).value(), 0); + BOOST_TEST_EQ(c.bin(-1).value(), 0); + BOOST_TEST_EQ(c.bin(0).value(), 1); + BOOST_TEST_EQ(c.bin(0).variance(), 1); + BOOST_TEST_EQ(c.bin(1).value(), 3); + BOOST_TEST_EQ(c.bin(1).variance(), 9); + BOOST_TEST_EQ(c.bin(2).value(), 0); auto d = a; d += b; - BOOST_TEST_EQ(d(-1).value(), 0); - BOOST_TEST_EQ(d(0).value(), 1); - BOOST_TEST_EQ(d(0).variance(), 1); - BOOST_TEST_EQ(d(1).value(), 3); - BOOST_TEST_EQ(d(1).variance(), 9); - BOOST_TEST_EQ(d(2).value(), 0); + BOOST_TEST_EQ(d.bin(-1).value(), 0); + BOOST_TEST_EQ(d.bin(0).value(), 1); + BOOST_TEST_EQ(d.bin(0).variance(), 1); + BOOST_TEST_EQ(d.bin(1).value(), 3); + BOOST_TEST_EQ(d.bin(1).variance(), 9); + BOOST_TEST_EQ(d.bin(2).value(), 0); } // add_3 @@ -483,22 +490,22 @@ template void run_tests() { make_histogram>(Type(), axis::integer<>(-1, 2)); auto b = make_histogram>(Type(), axis::integer<>(-1, 2)); - a.fill(-1); - b.fill(1); + a(-1); + b(1); auto c = a; c += b; - BOOST_TEST_EQ(c(-1), 0); - BOOST_TEST_EQ(c(0), 1); - BOOST_TEST_EQ(c(1), 0); - BOOST_TEST_EQ(c(2), 1); - BOOST_TEST_EQ(c(3), 0); + BOOST_TEST_EQ(c.bin(-1), 0); + BOOST_TEST_EQ(c.bin(0), 1); + BOOST_TEST_EQ(c.bin(1), 0); + BOOST_TEST_EQ(c.bin(2), 1); + BOOST_TEST_EQ(c.bin(3), 0); auto d = a; d += b; - BOOST_TEST_EQ(d(-1), 0); - BOOST_TEST_EQ(d(0), 1); - BOOST_TEST_EQ(d(1), 0); - BOOST_TEST_EQ(d(2), 1); - BOOST_TEST_EQ(d(3), 0); + BOOST_TEST_EQ(d.bin(-1), 0); + BOOST_TEST_EQ(d.bin(0), 1); + BOOST_TEST_EQ(d.bin(1), 0); + BOOST_TEST_EQ(d.bin(2), 1); + BOOST_TEST_EQ(d.bin(3), 0); } // bad_add @@ -508,51 +515,48 @@ template void run_tests() { BOOST_TEST_THROWS(a += b, std::logic_error); } - // bad_index - { - auto a = make_histogram(Type(), axis::integer<>(0, 2)); - BOOST_TEST_THROWS(a(5), std::out_of_range); - } - // functional programming { - auto v = std::vector{0, 1, 2, 3, 4, 5, 6, 7, 8, 9}; - auto h = make_histogram(Type(), axis::integer<>(0, 10)); - std::for_each(v.begin(), v.end(), [&h](int x) { h.fill(weight(2.0), x); }); - BOOST_TEST_EQ(h.sum().value(), 20); + auto v = std::vector{0, 1, 2}; + auto h = std::for_each(v.begin(), v.end(), + make_histogram(Type(), axis::integer<>(0, 3))); + BOOST_TEST_EQ(h.bin(0), 1); + BOOST_TEST_EQ(h.bin(1), 1); + BOOST_TEST_EQ(h.bin(2), 1); + BOOST_TEST_EQ(sum(h), 3); } // operators { auto a = make_histogram(Type(), axis::integer<>(0, 3)); auto b = a; - a.fill(0); - b.fill(1); + a(0); + b(1); auto c = a + b; - BOOST_TEST_EQ(c(0), 1); - BOOST_TEST_EQ(c(1), 1); + BOOST_TEST_EQ(c.bin(0), 1); + BOOST_TEST_EQ(c.bin(1), 1); c += b; - BOOST_TEST_EQ(c(0), 1); - BOOST_TEST_EQ(c(1), 2); + BOOST_TEST_EQ(c.bin(0), 1); + BOOST_TEST_EQ(c.bin(1), 2); auto d = a + b + c; - BOOST_TEST_EQ(d(0), 2); - BOOST_TEST_EQ(d(1), 3); + BOOST_TEST_EQ(d.bin(0), 2); + BOOST_TEST_EQ(d.bin(1), 3); auto e = 3 * a; auto f = b * 2; - BOOST_TEST_EQ(e(0).value(), 3); - BOOST_TEST_EQ(e(1).value(), 0); - BOOST_TEST_EQ(f(0).value(), 0); - BOOST_TEST_EQ(f(1).value(), 2); + BOOST_TEST_EQ(e.bin(0).value(), 3); + BOOST_TEST_EQ(e.bin(1).value(), 0); + BOOST_TEST_EQ(f.bin(0).value(), 0); + BOOST_TEST_EQ(f.bin(1).value(), 2); auto r = a; r += b; r += e; - BOOST_TEST_EQ(r(0).value(), 4); - BOOST_TEST_EQ(r(1).value(), 1); + BOOST_TEST_EQ(r.bin(0).value(), 4); + BOOST_TEST_EQ(r.bin(1).value(), 1); BOOST_TEST_EQ(r, a + b + 3 * a); auto s = r / 4; r /= 4; - BOOST_TEST_EQ(r(0).value(), 1); - BOOST_TEST_EQ(r(1).value(), 0.25); + BOOST_TEST_EQ(r.bin(0).value(), 1); + BOOST_TEST_EQ(r.bin(1).value(), 0.25); BOOST_TEST_EQ(r, s); } @@ -565,7 +569,7 @@ template void run_tests() { axis::regular(3, 1, 100, "lr"), axis::variable<>({0.1, 0.2, 0.3, 0.4, 0.5}, "v"), axis::category<>{A, B, C}, axis::integer<>(0, 2, "i")); - a.fill(0.5, 20, 0.1, 0.25, 1, 0); + a(0.5, 20, 0.1, 0.25, 1, 0); std::string buf; { std::ostringstream os; @@ -597,101 +601,103 @@ template void run_tests() { // histogram_reset { - auto a = make_histogram( + auto h = make_histogram( Type(), axis::integer<>(0, 2, "", axis::uoflow::off)); - a.fill(0); - a.fill(1); - BOOST_TEST_EQ(a(0), 1); - BOOST_TEST_EQ(a(1), 1); - a.reset(); - BOOST_TEST_EQ(a(0), 0); - BOOST_TEST_EQ(a(1), 0); + h(0); + h(1); + BOOST_TEST_EQ(h.bin(0), 1); + BOOST_TEST_EQ(h.bin(1), 1); + BOOST_TEST_EQ(sum(h), 2); + h.reset(); + BOOST_TEST_EQ(h.bin(0), 0); + BOOST_TEST_EQ(h.bin(1), 0); + BOOST_TEST_EQ(sum(h), 0); } // reduce { auto h1 = make_histogram(Type(), axis::integer<>(0, 2), axis::integer<>(0, 3)); - h1.fill(0, 0); - h1.fill(0, 1); - h1.fill(1, 0); - h1.fill(1, 1); - h1.fill(1, 2); + h1(0, 0); + h1(0, 1); + h1(1, 0); + h1(1, 1); + h1(1, 2); auto h1_0 = h1.reduce_to(0_c); BOOST_TEST_EQ(h1_0.dim(), 1); - BOOST_TEST_EQ(h1_0.sum(), 5); - BOOST_TEST_EQ(h1_0(0), 2); - BOOST_TEST_EQ(h1_0(1), 3); + BOOST_TEST_EQ(sum(h1_0), 5); + BOOST_TEST_EQ(h1_0.bin(0), 2); + BOOST_TEST_EQ(h1_0.bin(1), 3); BOOST_TEST_EQ(h1_0.axis()[0].lower(), 0); BOOST_TEST_EQ(h1_0.axis()[1].lower(), 1); BOOST_TEST(axis_equal(Type(), h1_0.axis(), h1.axis(0_c))); auto h1_1 = h1.reduce_to(1_c); BOOST_TEST_EQ(h1_1.dim(), 1); - BOOST_TEST_EQ(h1_1.sum(), 5); - BOOST_TEST_EQ(h1_1(0), 2); - BOOST_TEST_EQ(h1_1(1), 2); - BOOST_TEST_EQ(h1_1(2), 1); + BOOST_TEST_EQ(sum(h1_1), 5); + BOOST_TEST_EQ(h1_1.bin(0), 2); + BOOST_TEST_EQ(h1_1.bin(1), 2); + BOOST_TEST_EQ(h1_1.bin(2), 1); BOOST_TEST(axis_equal(Type(), h1_1.axis(), h1.axis(1_c))); auto h2 = make_histogram(Type(), axis::integer<>(0, 2), axis::integer<>(0, 3), axis::integer<>(0, 4)); - h2.fill(0, 0, 0); - h2.fill(0, 1, 0); - h2.fill(0, 1, 1); - h2.fill(0, 0, 2); - h2.fill(1, 0, 2); + h2(0, 0, 0); + h2(0, 1, 0); + h2(0, 1, 1); + h2(0, 0, 2); + h2(1, 0, 2); auto h2_0 = h2.reduce_to(0_c); BOOST_TEST_EQ(h2_0.dim(), 1); - BOOST_TEST_EQ(h2_0.sum(), 5); - BOOST_TEST_EQ(h2_0(0), 4); - BOOST_TEST_EQ(h2_0(1), 1); + BOOST_TEST_EQ(sum(h2_0), 5); + BOOST_TEST_EQ(h2_0.bin(0), 4); + BOOST_TEST_EQ(h2_0.bin(1), 1); BOOST_TEST(axis_equal(Type(), h2_0.axis(), axis::integer<>(0, 2))); auto h2_1 = h2.reduce_to(1_c); BOOST_TEST_EQ(h2_1.dim(), 1); - BOOST_TEST_EQ(h2_1.sum(), 5); - BOOST_TEST_EQ(h2_1(0), 3); - BOOST_TEST_EQ(h2_1(1), 2); + BOOST_TEST_EQ(sum(h2_1), 5); + BOOST_TEST_EQ(h2_1.bin(0), 3); + BOOST_TEST_EQ(h2_1.bin(1), 2); BOOST_TEST(axis_equal(Type(), h2_1.axis(), axis::integer<>(0, 3))); auto h2_2 = h2.reduce_to(2_c); BOOST_TEST_EQ(h2_2.dim(), 1); - BOOST_TEST_EQ(h2_2.sum(), 5); - BOOST_TEST_EQ(h2_2(0), 2); - BOOST_TEST_EQ(h2_2(1), 1); - BOOST_TEST_EQ(h2_2(2), 2); + BOOST_TEST_EQ(sum(h2_2), 5); + BOOST_TEST_EQ(h2_2.bin(0), 2); + BOOST_TEST_EQ(h2_2.bin(1), 1); + BOOST_TEST_EQ(h2_2.bin(2), 2); BOOST_TEST(axis_equal(Type(), h2_2.axis(), axis::integer<>(0, 4))); auto h2_01 = h2.reduce_to(0_c, 1_c); BOOST_TEST_EQ(h2_01.dim(), 2); - BOOST_TEST_EQ(h2_01.sum(), 5); - BOOST_TEST_EQ(h2_01(0, 0), 2); - BOOST_TEST_EQ(h2_01(0, 1), 2); - BOOST_TEST_EQ(h2_01(1, 0), 1); + BOOST_TEST_EQ(sum(h2_01), 5); + BOOST_TEST_EQ(h2_01.bin(0, 0), 2); + BOOST_TEST_EQ(h2_01.bin(0, 1), 2); + BOOST_TEST_EQ(h2_01.bin(1, 0), 1); BOOST_TEST(axis_equal(Type(), h2_01.axis(0_c), axis::integer<>(0, 2))); BOOST_TEST(axis_equal(Type(), h2_01.axis(1_c), axis::integer<>(0, 3))); auto h2_02 = h2.reduce_to(0_c, 2_c); BOOST_TEST_EQ(h2_02.dim(), 2); - BOOST_TEST_EQ(h2_02.sum(), 5); - BOOST_TEST_EQ(h2_02(0, 0), 2); - BOOST_TEST_EQ(h2_02(0, 1), 1); - BOOST_TEST_EQ(h2_02(0, 2), 1); - BOOST_TEST_EQ(h2_02(1, 2), 1); + BOOST_TEST_EQ(sum(h2_02), 5); + BOOST_TEST_EQ(h2_02.bin(0, 0), 2); + BOOST_TEST_EQ(h2_02.bin(0, 1), 1); + BOOST_TEST_EQ(h2_02.bin(0, 2), 1); + BOOST_TEST_EQ(h2_02.bin(1, 2), 1); BOOST_TEST(axis_equal(Type(), h2_02.axis(0_c), axis::integer<>(0, 2))); BOOST_TEST(axis_equal(Type(), h2_02.axis(1_c), axis::integer<>(0, 4))); auto h2_12 = h2.reduce_to(1_c, 2_c); BOOST_TEST_EQ(h2_12.dim(), 2); - BOOST_TEST_EQ(h2_12.sum(), 5); - BOOST_TEST_EQ(h2_12(0, 0), 1); - BOOST_TEST_EQ(h2_12(1, 0), 1); - BOOST_TEST_EQ(h2_12(1, 1), 1); - BOOST_TEST_EQ(h2_12(0, 2), 2); + BOOST_TEST_EQ(sum(h2_12), 5); + BOOST_TEST_EQ(h2_12.bin(0, 0), 1); + BOOST_TEST_EQ(h2_12.bin(1, 0), 1); + BOOST_TEST_EQ(h2_12.bin(1, 1), 1); + BOOST_TEST_EQ(h2_12.bin(0, 2), 2); BOOST_TEST(axis_equal(Type(), h2_12.axis(0_c), axis::integer<>(0, 3))); BOOST_TEST(axis_equal(Type(), h2_12.axis(1_c), axis::integer<>(0, 4))); } @@ -710,25 +716,25 @@ template void run_tests() { }; auto h = make_histogram(Type(), custom_axis(0, 3)); - h.fill("-10"); - h.fill("0"); - h.fill("1"); - h.fill("9"); + h("-10"); + h("0"); + h("1"); + h("9"); BOOST_TEST_EQ(h.dim(), 1); BOOST_TEST(h.axis() == custom_axis(0, 3)); - BOOST_TEST_EQ(h(0), 1); - BOOST_TEST_EQ(h(1), 1); - BOOST_TEST_EQ(h(2), 0); + BOOST_TEST_EQ(h.bin(0), 1); + BOOST_TEST_EQ(h.bin(1), 1); + BOOST_TEST_EQ(h.bin(2), 0); } // histogram iterator 1D { auto h = make_histogram(Type(), axis::integer<>(0, 3)); const auto &a = h.axis(); - h.fill(0, weight(2)); - h.fill(1); - h.fill(1); + h(weight(2), 0); + h(1); + h(1); auto it = h.begin(); BOOST_TEST_EQ(it.dim(), 1); @@ -745,18 +751,27 @@ template void run_tests() { BOOST_TEST_EQ(it.bin(), a[2]); BOOST_TEST_EQ(*it, 0); ++it; + BOOST_TEST_EQ(it.idx(), 3); + BOOST_TEST_EQ(it.bin(), a[3]); + BOOST_TEST_EQ(*it, 0); + ++it; + BOOST_TEST_EQ(it.idx(), -1); + BOOST_TEST_EQ(it.bin(), a[-1]); + BOOST_TEST_EQ(*it, 0); + ++it; BOOST_TEST(it == h.end()); } // histogram iterator 2D { - auto h = make_histogram(Type(), axis::integer<>(0, 3), - axis::integer<>(2, 4)); + auto h = make_histogram(Type(), axis::integer<>(0, 1), + axis::integer<>(2, 4, "", axis::uoflow::off)); const auto &a0 = h.axis(0_c); const auto &a1 = h.axis(1_c); - h.fill(0, 2, weight(2)); - h.fill(1, 2); - h.fill(1, 3); + h(weight(2), 0, 2); + h(-1, 2); + h(1, 3); + auto it = h.begin(); BOOST_TEST_EQ(it.dim(), 2); @@ -771,40 +786,49 @@ template void run_tests() { BOOST_TEST_EQ(it.idx(1), 0); BOOST_TEST_EQ(it.bin(0_c), a0[1]); BOOST_TEST_EQ(it.bin(1_c), a1[0]); - BOOST_TEST_EQ(*it, 1); + BOOST_TEST_EQ(it->value(), 0); + BOOST_TEST_EQ(it->variance(), 0); ++it; - BOOST_TEST_EQ(it.idx(0), 2); + BOOST_TEST_EQ(it.idx(0), -1); BOOST_TEST_EQ(it.idx(1), 0); - BOOST_TEST_EQ(it.bin(0_c), a0[2]); + BOOST_TEST_EQ(it.bin(0_c), a0[-1]); BOOST_TEST_EQ(it.bin(1_c), a1[0]); - BOOST_TEST_EQ(*it, 0); + BOOST_TEST_EQ(it->value(), 1); + BOOST_TEST_EQ(it->variance(), 1); ++it; BOOST_TEST_EQ(it.idx(0), 0); BOOST_TEST_EQ(it.idx(1), 1); BOOST_TEST_EQ(it.bin(0_c), a0[0]); BOOST_TEST_EQ(it.bin(1_c), a1[1]); - BOOST_TEST_EQ(*it, 0); + BOOST_TEST_EQ(it->value(), 0); + BOOST_TEST_EQ(it->variance(), 0); ++it; BOOST_TEST_EQ(it.idx(0), 1); BOOST_TEST_EQ(it.idx(1), 1); BOOST_TEST_EQ(it.bin(0_c), a0[1]); BOOST_TEST_EQ(it.bin(1_c), a1[1]); - BOOST_TEST_EQ(*it, 1); + BOOST_TEST_EQ(it->value(), 1); + BOOST_TEST_EQ(it->variance(), 1); ++it; - BOOST_TEST_EQ(it.idx(0), 2); + BOOST_TEST_EQ(it.idx(0), -1); BOOST_TEST_EQ(it.idx(1), 1); - BOOST_TEST_EQ(it.bin(0_c), a0[2]); + BOOST_TEST_EQ(it.bin(0_c), a0[-1]); BOOST_TEST_EQ(it.bin(1_c), a1[1]); - BOOST_TEST_EQ(*it, 0); + BOOST_TEST_EQ(it->value(), 0); + BOOST_TEST_EQ(it->variance(), 0); ++it; BOOST_TEST(it == h.end()); + + auto v = sum(h); + BOOST_TEST_EQ(v.value(), 4); + BOOST_TEST_EQ(v.variance(), 6); } // STL compatibility { auto h = make_histogram(Type(), axis::integer<>(0, 3)); for (int i = 0; i < 3; ++i) - h.fill(i); + h(i); auto a = std::vector(); std::partial_sum(h.begin(), h.end(), std::back_inserter(a)); BOOST_TEST_EQ(a[0], 1); @@ -812,6 +836,32 @@ template void run_tests() { BOOST_TEST_EQ(a[2], 3); } + // using STL containers + { + auto h = make_histogram(Type(), + axis::integer<>(0, 2), + axis::regular<>(2, 2, 4)); + // vector in + h(std::vector({0, 2})); + // pair in + h(std::make_pair(1, 3.0)); + + // pair out + BOOST_TEST_EQ(h.bin(std::make_pair(0, 0)), 1); + // tuple out + BOOST_TEST_EQ(h.bin(std::make_tuple(1, 1)), 1); + + // vector in, weights + h(weight(2), std::vector({0, 2})); + // pair in, weights + h(weight(2), std::make_pair(1, 3.0)); + + // vector + BOOST_TEST_EQ(h.bin(std::vector({0, 0})).value(), 3); + // // initializer_list + // BOOST_TEST_EQ(h.bin({1, 1}).variance(), 5); + } + // pass histogram to function { auto h = make_histogram(Type(), axis::integer<>(0, 3)); @@ -842,7 +892,7 @@ template void run_mixed_tests() { axis::integer<>(0, 2)); auto b = make_histogram>(T2{}, axis::regular<>{3, 0, 3}, axis::integer<>(0, 2)); - a.fill(1, 1); + a(1, 1); BOOST_TEST_NE(a, b); b = a; BOOST_TEST_EQ(a, b); @@ -867,34 +917,6 @@ int main() { BOOST_TEST_EQ(h.axis(1), v[1]); } - // using iterator ranges - { - auto h = - make_dynamic_histogram(axis::integer<>(0, 2), axis::integer<>(2, 4)); - auto v = std::vector(2); - auto i = std::array(); - - v = {0, 2}; - h.fill(v.begin(), v.end()); - v = {1, 3}; - h.fill(v.begin(), v.end()); - - i = {{0, 0}}; - BOOST_TEST_EQ(h(i.begin(), i.end()), 1); - i = {{1, 1}}; - BOOST_TEST_EQ(h(i.begin(), i.end()), 1); - - v = {0, 2}; - h.fill(v.begin(), v.end(), weight(2)); - v = {1, 3}; - h.fill(v.begin(), v.end(), weight(2)); - - i = {{0, 0}}; - BOOST_TEST_EQ(h(i.begin(), i.end()).value(), 3); - i = {{1, 1}}; - BOOST_TEST_EQ(h(i.begin(), i.end()).variance(), 5); - } - // axis methods { enum { A, B }; @@ -906,25 +928,25 @@ int main() { { auto h1 = make_dynamic_histogram(axis::integer<>(0, 2), axis::integer<>(0, 3)); - h1.fill(0, 0); - h1.fill(0, 1); - h1.fill(1, 0); - h1.fill(1, 1); - h1.fill(1, 2); + h1(0, 0); + h1(0, 1); + h1(1, 0); + h1(1, 1); + h1(1, 2); auto h1_0 = h1.reduce_to(0); BOOST_TEST_EQ(h1_0.dim(), 1); - BOOST_TEST_EQ(h1_0.sum(), 5); - BOOST_TEST_EQ(h1_0(0), 2); - BOOST_TEST_EQ(h1_0(1), 3); + BOOST_TEST_EQ(sum(h1_0), 5); + BOOST_TEST_EQ(h1_0.bin(0), 2); + BOOST_TEST_EQ(h1_0.bin(1), 3); BOOST_TEST(axis_equal(dynamic_tag(), h1_0.axis(), h1.axis(0_c))); auto h1_1 = h1.reduce_to(1); BOOST_TEST_EQ(h1_1.dim(), 1); - BOOST_TEST_EQ(h1_1.sum(), 5); - BOOST_TEST_EQ(h1_1(0), 2); - BOOST_TEST_EQ(h1_1(1), 2); - BOOST_TEST_EQ(h1_1(2), 1); + BOOST_TEST_EQ(sum(h1_1), 5); + BOOST_TEST_EQ(h1_1.bin(0), 2); + BOOST_TEST_EQ(h1_1.bin(1), 2); + BOOST_TEST_EQ(h1_1.bin(2), 1); BOOST_TEST(axis_equal(dynamic_tag(), h1_1.axis(), h1.axis(1_c))); } @@ -932,9 +954,9 @@ int main() { { auto h = make_dynamic_histogram(axis::integer<>(0, 3)); const auto &a = h.axis(); - h.fill(0, weight(2)); - h.fill(1); - h.fill(1); + h(weight(2), 0); + h(1); + h(1); auto it = h.begin(); BOOST_TEST_EQ(it.dim(), 1); @@ -951,6 +973,14 @@ int main() { BOOST_TEST_EQ(it.bin(0), a[2]); BOOST_TEST_EQ(*it, 0); ++it; + BOOST_TEST_EQ(it.idx(0), 3); + BOOST_TEST_EQ(it.bin(0), a[3]); + BOOST_TEST_EQ(*it, 0); + ++it; + BOOST_TEST_EQ(it.idx(0), -1); + BOOST_TEST_EQ(it.bin(0), a[-1]); + BOOST_TEST_EQ(*it, 0); + ++it; BOOST_TEST(it == h.end()); } diff --git a/test/weight_counter_test.cpp b/test/weight_counter_test.cpp index 6e7d93ef..ced45d00 100644 --- a/test/weight_counter_test.cpp +++ b/test/weight_counter_test.cpp @@ -41,6 +41,10 @@ int main() { BOOST_TEST_EQ(w.variance(), 5); BOOST_TEST_EQ(w, wcount(3, 5)); BOOST_TEST_NE(w, wcount(3)); + + w += wcount(1, 2); + BOOST_TEST_EQ(w.value(), 4); + BOOST_TEST_EQ(w.variance(), 7); } {