diff --git a/CMake/CMakeLists.txt b/CMake/CMakeLists.txt new file mode 100644 index 00000000..9cff97c4 --- /dev/null +++ b/CMake/CMakeLists.txt @@ -0,0 +1,76 @@ +project(histogram) +list(APPEND CMAKE_MODULE_PATH ${CMAKE_SOURCE_DIR}/CMake) + +if(${CMAKE_BUILD_TYPE}) + STRING(TOLOWER ${CMAKE_BUILD_TYPE} CMAKE_BUILD_TYPE) +endif() + +find_package(Boost 1.55 REQUIRED COMPONENTS python iostreams serialization unit_test_framework) +find_package(PythonLibs) +find_package(Numpy) # optional +find_package(ROOT) # only used in one of the tests + +include_directories(include ${Boost_INCLUDE_DIRS} ${PYTHON_INCLUDE_DIRS}) +add_definitions(-DBOOST_TEST_DYN_LINK) # for unit_test_framework +set(LIBRARIES stdc++ ${Boost_LIBRARIES}) + +if(Boost_PYTHON_FOUND) + LIST(APPEND LIBRARIES ${PYTHON_LIBRARIES}) +endif() + +if(NUMPY_FOUND) + include_directories(${NUMPY_INCLUDE_DIR}) + add_definitions(-DUSE_NUMPY) +endif() + +add_library(histogram SHARED + src/axis.cpp + src/histogram_base.cpp + src/nhistogram.cpp + src/nstore.cpp + src/zero_suppression.cpp + # src/whistogram.cpp +) +target_link_libraries(histogram ${LIBRARIES}) + +if(CMAKE_BUILD_TYPE STREQUAL "debug") + message(STATUS "debug mode: optimizations off") +else() + message(STATUS "release mode: optimizations on") + target_compile_options(histogram PUBLIC "-O3 -fomit-frame-pointer -mtune=generic") +endif() + +if(Boost_PYTHON_FOUND) + add_library(pyhistogram MODULE + src/python/module.cpp + src/python/histogram_base.cpp + src/python/nhistogram.cpp + ) + target_link_libraries(pyhistogram histogram ${LIBRARIES}) + set_target_properties(pyhistogram PROPERTIES OUTPUT_NAME "histogram" PREFIX "" SUFFIX ".so") +endif() + +add_executable(sizeof_axis + examples/sizeof_axis.cpp +) +target_link_libraries(sizeof_axis ${LIBRARIES}) + +if(ROOT_FOUND) + add_executable(nhistogram_speed + examples/speed_vs_root.cpp) + target_include_directories(nhistogram_speed PUBLIC ${ROOT_INCLUDE_DIR}) + target_link_libraries(nhistogram_speed histogram ${ROOT_LIBRARIES} ${LIBRARIES}) +endif() + +# tests +enable_testing() +add_executable(zero_suppression_test + test/zero_suppression_test.cpp) +target_link_libraries(zero_suppression_test histogram ${LIBRARIES}) +add_test(zero_suppression_test zero_suppression_test) + +if(Boost_PYTHON_FOUND) + add_custom_target(python_suite_test ALL + cp ${CMAKE_SOURCE_DIR}/test/python_suite_test.py .) + add_test(python_suite_test python_suite_test.py) +endif() diff --git a/CMake/FindNumpy.cmake b/CMake/FindNumpy.cmake new file mode 100644 index 00000000..5c70a24a --- /dev/null +++ b/CMake/FindNumpy.cmake @@ -0,0 +1,24 @@ +# VARIABLES set by this module +# - NUMPY_INCLUDE_DIR +# - NUMPY_FOUND +find_package(PythonInterp) + +set(NUMPY_FOUND FALSE) + +if(PYTHONINTERP_FOUND) + execute_process( + COMMAND ${PYTHON_EXECUTABLE} -c "import numpy; print numpy.get_include()" + OUTPUT_VARIABLE NUMPY_INCLUDE_DIR + OUTPUT_STRIP_TRAILING_WHITESPACE) + + if(EXISTS ${NUMPY_INCLUDE_DIR}) + set(NUMPY_FOUND TRUE) + if (NOT NUMPY_FIND_QUIETLY) + message (STATUS "Found numpy") + endif (NOT NUMPY_FIND_QUIETLY) + endif() +endif() + +if(NOT NUMPY_FOUND AND NUMPY_FIND_REQUIRED) + message (FATAL_ERROR "numpy missing") +endif() diff --git a/CMake/FindROOT.cmake b/CMake/FindROOT.cmake new file mode 100644 index 00000000..c5fef351 --- /dev/null +++ b/CMake/FindROOT.cmake @@ -0,0 +1,159 @@ +# - Finds ROOT instalation +# This module sets up ROOT information +# It defines: +# ROOT_FOUND If the ROOT is found +# ROOT_INCLUDE_DIR PATH to the include directory +# ROOT_LIBRARIES Most common libraries +# ROOT_LIBRARY_DIR PATH to the library directory +# +# Updated by K. Smith (ksmith37@nd.edu) to properly handle +# dependncies in ROOT_GENERATE_DICTIONARY + +find_program(ROOT_CONFIG_EXECUTABLE root-config + PATHS $ENV{ROOTSYS}/bin) + +if(NOT ROOT_CONFIG_EXECUTABLE) + set(ROOT_FOUND FALSE) +else() + set(ROOT_FOUND TRUE) + + execute_process( + COMMAND ${ROOT_CONFIG_EXECUTABLE} --prefix + OUTPUT_VARIABLE ROOTSYS + OUTPUT_STRIP_TRAILING_WHITESPACE) + + execute_process( + COMMAND ${ROOT_CONFIG_EXECUTABLE} --version + OUTPUT_VARIABLE ROOT_VERSION + OUTPUT_STRIP_TRAILING_WHITESPACE) + + execute_process( + COMMAND ${ROOT_CONFIG_EXECUTABLE} --incdir + OUTPUT_VARIABLE ROOT_INCLUDE_DIR + OUTPUT_STRIP_TRAILING_WHITESPACE) + + execute_process( + COMMAND ${ROOT_CONFIG_EXECUTABLE} --libs + OUTPUT_VARIABLE ROOT_LIBRARIES + OUTPUT_STRIP_TRAILING_WHITESPACE) + + #set(ROOT_LIBRARIES ${ROOT_LIBRARIES} -lThread -lMinuit -lHtml -lVMC -lEG -lGeom -lTreePlayer -lXMLIO -lProof) + #set(ROOT_LIBRARIES ${ROOT_LIBRARIES} -lProofPlayer -lMLP -lSpectrum -lEve -lRGL -lGed -lXMLParser -lPhysics) + set(ROOT_LIBRARY_DIR ${ROOTSYS}/lib) +endif() + +include(FindPackageHandleStandardArgs) +find_package_handle_standard_args(ROOT DEFAULT_MSG ROOT_CONFIG_EXECUTABLE + ROOTSYS ROOT_VERSION ROOT_INCLUDE_DIR ROOT_LIBRARIES ROOT_LIBRARY_DIR) + +mark_as_advanced(ROOT_CONFIG_EXECUTABLE) + +include(CMakeParseArguments) +find_program(ROOTCINT_EXECUTABLE rootcint PATHS $ENV{ROOTSYS}/bin) +find_program(GENREFLEX_EXECUTABLE genreflex PATHS $ENV{ROOTSYS}/bin) +find_package(GCCXML) + +#---------------------------------------------------------------------------- +# function ROOT_GENERATE_DICTIONARY( dictionary +# header1 header2 ... +# LINKDEF linkdef1 ... +# OPTIONS opt1...) +function(ROOT_GENERATE_DICTIONARY dictionary) + CMAKE_PARSE_ARGUMENTS(ARG "" "" "LINKDEF;OPTIONS" "" ${ARGN}) + #---Get the list of include directories------------------ + get_directory_property(incdirs INCLUDE_DIRECTORIES) + set(includedirs) + foreach( d ${incdirs}) + set(includedirs ${includedirs} -I${d}) + endforeach() + #---Get the list of header files------------------------- + set(headerfiles) + foreach(fp ${ARG_UNPARSED_ARGUMENTS}) + file(GLOB files ${fp}) + if(files) + foreach(f ${files}) + if(NOT f MATCHES LinkDef) + find_file(headerFile ${f} PATHS ${incdirs}) + set(headerfiles ${headerfiles} ${headerFile}) + unset(headerFile CACHE) + endif() + endforeach() + else() + find_file(headerFile ${fp} PATHS ${incdirs}) + set(headerfiles ${headerfiles} ${headerFile}) + unset(headerFile CACHE) + endif() + endforeach() + #---Get LinkDef.h file------------------------------------ + set(linkdefs) + foreach( f ${ARG_LINKDEF}) + find_file(linkFile ${f} PATHS ${incdirs}) + set(linkdefs ${linkdefs} ${linkFile}) + unset(linkFile CACHE) + endforeach() + #---call rootcint------------------------------------------ + add_custom_command(OUTPUT ${dictionary}.cxx ${dictionary}.h + COMMAND ${ROOTCINT_EXECUTABLE} -cint -f ${dictionary}.cxx + -c ${ARG_OPTIONS} ${includedirs} ${headerfiles} ${linkdefs} + DEPENDS ${headerfiles} ${linkdefs} VERBATIM) +endfunction() + +#---------------------------------------------------------------------------- +# function REFLEX_GENERATE_DICTIONARY(dictionary +# header1 header2 ... +# SELECTION selectionfile ... +# OPTIONS opt1...) +function(REFLEX_GENERATE_DICTIONARY dictionary) + CMAKE_PARSE_ARGUMENTS(ARG "" "" "SELECTION;OPTIONS" "" ${ARGN}) + #---Get the list of header files------------------------- + set(headerfiles) + foreach(fp ${ARG_UNPARSED_ARGUMENTS}) + file(GLOB files ${fp}) + if(files) + foreach(f ${files}) + set(headerfiles ${headerfiles} ${f}) + endforeach() + else() + set(headerfiles ${headerfiles} ${fp}) + endif() + endforeach() + #---Get Selection file------------------------------------ + if(IS_ABSOLUTE ${ARG_SELECTION}) + set(selectionfile ${ARG_SELECTION}) + else() + set(selectionfile ${CMAKE_CURRENT_SOURCE_DIR}/${ARG_SELECTION}) + endif() + #---Get the list of include directories------------------ + get_directory_property(incdirs INCLUDE_DIRECTORIES) + set(includedirs) + foreach( d ${incdirs}) + set(includedirs ${includedirs} -I${d}) + endforeach() + #---Get preprocessor definitions-------------------------- + get_directory_property(defs COMPILE_DEFINITIONS) + foreach( d ${defs}) + set(definitions ${definitions} -D${d}) + endforeach() + #---Nanes and others--------------------------------------- + set(gensrcdict ${dictionary}.cpp) + if(MSVC) + set(gccxmlopts "--gccxmlopt=\"--gccxml-compiler cl\"") + else() + #set(gccxmlopts "--gccxmlopt=\'--gccxml-cxxflags -m64 \'") + set(gccxmlopts) + endif() + #set(rootmapname ${dictionary}Dict.rootmap) + #set(rootmapopts --rootmap=${rootmapname} --rootmap-lib=${libprefix}${dictionary}Dict) + #---Check GCCXML and get path----------------------------- + if(GCCXML) + get_filename_component(gccxmlpath ${GCCXML} PATH) + else() + message(WARNING "GCCXML not found. Install and setup your environment to find 'gccxml' executable") + endif() + #---Actual command---------------------------------------- + add_custom_command(OUTPUT ${gensrcdict} ${rootmapname} + COMMAND ${GENREFLEX_EXECUTABLE} ${headerfiles} -o ${gensrcdict} ${gccxmlopts} ${rootmapopts} --select=${selectionfile} + --gccxmlpath=${gccxmlpath} ${ARG_OPTIONS} ${includedirs} ${definitions} + DEPENDS ${headerfiles} ${selectionfile}) +endfunction() + diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 00000000..8205be45 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,2 @@ +cmake_minimum_required (VERSION 2.8) +include(CMake/CMakeLists.txt) diff --git a/examples/merge.py b/examples/merge.py new file mode 100755 index 00000000..bf907f8e --- /dev/null +++ b/examples/merge.py @@ -0,0 +1,41 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +from icecube.histogram import I3HistogramN +from icecube import dataio, icetray +import numpy as np +import argparse + + +def main(): + parser = argparse.ArgumentParser() + parser.add_argument("input", + nargs="+") + parser.add_argument("-o", + metavar="output", + default="merged.i3.gz", + help="name of output file") + args = parser.parse_args() + + histograms = {} + for fname in args.input: + print "Processing:", fname + frame = dataio.I3File(fname).pop_frame() + for key in frame.keys(): + print " -", key + h = frame[key] + if key not in histograms: + histograms[key] = I3HistogramN(h) + else: + histograms[key] += h + + print "Writing:", args.o + frame = icetray.I3Frame() + for key in histograms: + print " -", key + frame[key] = histograms[key] + f = dataio.I3File(args.o, "w") + f.push(frame) + f.close() + +if __name__ == "__main__": + main() diff --git a/examples/sizeof_axis.cpp b/examples/sizeof_axis.cpp new file mode 100644 index 00000000..8465711a --- /dev/null +++ b/examples/sizeof_axis.cpp @@ -0,0 +1,30 @@ +#include +#include +#include +#include +#include +#include +#include +#include + +int main(int argc, char** argv) { + using namespace boost::histogram; + #define SIZEOF(x) std::printf("%32s: %lu\n", BOOST_STRINGIZE(x), sizeof(x)) + SIZEOF(char); + SIZEOF(int); + SIZEOF(long); + SIZEOF(float); + SIZEOF(double); + SIZEOF(void*); + SIZEOF(std::string); + SIZEOF(std::vector); + SIZEOF(std::valarray); + SIZEOF(boost::container::vector); + SIZEOF(boost::scoped_array); + SIZEOF(regular_axis); + SIZEOF(polar_axis); + SIZEOF(variable_axis); + SIZEOF(category_axis); + SIZEOF(integer_axis); + SIZEOF(axis_type); +} diff --git a/examples/speed_vs_root.cpp b/examples/speed_vs_root.cpp new file mode 100644 index 00000000..eedbf391 --- /dev/null +++ b/examples/speed_vs_root.cpp @@ -0,0 +1,146 @@ +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace std; +using namespace boost::histogram; + +template +struct rng { + boost::random::mt19937 r; + D d; + rng(double a, double b) : d(a, b) {} + double operator()() { return d(r); } +}; + +vector random_array(unsigned n, int type) { + using namespace boost::random; + std::vector result; + switch (type) { + case 0: + std::generate_n(std::back_inserter(result), n, rng >(0.0, 1.0)); + break; + case 1: + std::generate_n(std::back_inserter(result), n, rng >(0.0, 0.3)); + break; + } + return result; +} + +void compare_1d(unsigned n) +{ + vector r = random_array(1000000, 1); + + double best_root = std::numeric_limits::max(); + double best_boost = std::numeric_limits::max(); + for (unsigned k = 0; k < 3; ++k) { + TH1I hroot("", "", 100, 0, 1); + clock_t t = clock(); + for (unsigned i = 0; i < n; ++i) + hroot.Fill(r[i]); + t = clock() - t; + best_root = std::min(best_root, double(t) / CLOCKS_PER_SEC); + + nhistogram h(regular_axis(100, 0, 1, std::string(), true)); + t = clock(); + for (unsigned i = 0; i < n; ++i) + h.fill(r[i]); + t = clock() - t; + best_boost = std::min(best_boost, double(t) / CLOCKS_PER_SEC); + // printf("root %g this %g\n", hroot.GetSum(), h.sum()); + assert(hroot.GetSum() == h.sum()); + } + + printf("1D\n"); + printf("t[root] = %g\n", best_root); + printf("t[boost] = %g\n", best_boost); +} + +void compare_3d(unsigned n) +{ + vector r = random_array(3 * n, 1); + + double best_root = std::numeric_limits::max(); + double best_boost = std::numeric_limits::max(); + for (unsigned k = 0; k < 3; ++k) { + TH3I hroot("", "", 100, 0, 1, 100, 0, 1, 100, 0, 1); + clock_t t = clock(); + for (unsigned i = 0; i < n; ++i) + hroot.Fill(r[3 * i], r[3 * i + 1], r[3 * i + 2]); + t = clock() - t; + best_root = std::min(best_root, double(t) / CLOCKS_PER_SEC); + + nhistogram h(regular_axis(100, 0, 1, "", true), + regular_axis(100, 0, 1, "", true), + regular_axis(100, 0, 1, "", true)); + t = clock(); + for (unsigned i = 0; i < n; ++i) + h.fill(r[3 * i], r[3 * i + 1], r[3 * i + 2]); + t = clock() - t; + best_boost = std::min(best_boost, double(t) / CLOCKS_PER_SEC); + assert(hroot.GetSum() == h.sum()); + } + + printf("3D\n"); + printf("t[root] = %g\n", best_root); + printf("t[boost] = %g\n", best_boost); +} + +void compare_6d(unsigned n) +{ + vector r = random_array(6 * n, 1); + + double best_root = std::numeric_limits::max(); + double best_boost = std::numeric_limits::max(); + for (unsigned k = 0; k < 3; ++k) { + double x[6]; + vector bin(6, 10); + vector min(6, 0); + vector max(6, 1); + THnI hroot("", "", 6, &bin.front(), &min.front(), &max.front()); + + clock_t t = clock(); + for (unsigned i = 0; i < n; ++i) { + for (unsigned k = 0; k < 6; ++k) + x[k] = r[6 * i + k]; + hroot.Fill(x); + } + t = clock() - t; + best_root = std::min(best_root, double(t) / CLOCKS_PER_SEC); + + nhistogram h(regular_axis(10, 0, 1), + regular_axis(10, 0, 1), + regular_axis(10, 0, 1), + regular_axis(10, 0, 1), + regular_axis(10, 0, 1), + regular_axis(10, 0, 1)); + + t = clock(); + for (unsigned i = 0; i < n; ++i) { + for (unsigned k = 0; k < 6; ++k) + x[k] = r[6 * i + k]; + h.fill(x); + } + t = clock() - t; + best_boost = std::min(best_boost, double(t) / CLOCKS_PER_SEC); + } + + printf("6D\n"); + printf("t[root] = %g\n", best_root); + printf("t[boost] = %g\n", best_boost); +} + +int main(int argc, char** argv) { + compare_1d(1000000); + compare_3d(100000); + compare_6d(100000); +} diff --git a/include/boost/histogram.hpp b/include/boost/histogram.hpp new file mode 100644 index 00000000..1adb858a --- /dev/null +++ b/include/boost/histogram.hpp @@ -0,0 +1,7 @@ +#ifndef _HISTOGRAM_I3HISTOGRAM_H_ +#define _HISTOGRAM_I3HISTOGRAM_H_ + +#include +#include + +#endif diff --git a/include/boost/histogram/axis.hpp b/include/boost/histogram/axis.hpp new file mode 100644 index 00000000..005c8bf0 --- /dev/null +++ b/include/boost/histogram/axis.hpp @@ -0,0 +1,252 @@ +#ifndef _BOOST_HISTOGRAM_AXIS_HPP_ +#define _BOOST_HISTOGRAM_AXIS_HPP_ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace boost { +namespace histogram { + +class axis_base { +public: + inline unsigned bins() const { return size_ > 0 ? size_ : -size_; } + inline bool uoflow() const { return size_ > 0; } + const std::string& label() const { return label_; } + void label(const std::string& label) { label_ = label; } + +protected: + axis_base(int, const std::string&, bool); + + axis_base() : size_(0) {} + explicit axis_base(const axis_base&); + axis_base& operator=(const axis_base&); + + bool operator==(const axis_base&) const; + +private: + int size_; + std::string label_; + + friend class serialization::access; + template + void serialize(Archive& ar, unsigned version) + { + using namespace serialization; + ar & size_; + ar & label_; + } +}; + +template +class real_axis { +public: + typedef double value_type; + + double left(int idx) const { + return static_cast(*this)[idx]; + } + + double right(int idx) const { + return static_cast(*this)[idx + 1]; + } + + double center(int idx) const { + return 0.5 * (left(idx) + right(idx)); + } +}; + +// real regular axis (constant bin widths) +class regular_axis: public axis_base, public real_axis { +public: + regular_axis(unsigned n, double min, double max, + const std::string& label = std::string(), + bool uoflow = true); + + regular_axis() {} + explicit regular_axis(const regular_axis&); + regular_axis& operator=(const regular_axis&); + + inline int index(double x) const { + const double z = (x - min_) / range_; + return algorithm::clamp(int(floor(z * bins())), -1, bins()); + } + + double operator[](int idx) const; + bool operator==(const regular_axis&) const; +private: + double min_, range_; + + friend class serialization::access; + template + void serialize(Archive& ar, unsigned version) + { + using namespace serialization; + ar & boost::serialization::base_object(*this); + ar & min_; + ar & range_; + } +}; + +// real polar axis (constant bin widths, wraps around) +class polar_axis: public axis_base, public real_axis { +public: + polar_axis(unsigned n, double start = 0.0, + const std::string& label = std::string()); + + polar_axis() {} + explicit polar_axis(const polar_axis&); + polar_axis& operator=(const polar_axis&); + + inline int index(double x) const { + const double z = (x - start_) / 6.283185307179586; + const int i = int(floor(z * bins())) % bins(); + return i + (i < 0) * bins(); + } + + double operator[](int idx) const; + bool operator==(const polar_axis&) const; +private: + double start_; + + friend class serialization::access; + template + void serialize(Archive& ar, unsigned version) + { + using namespace serialization; + ar & boost::serialization::base_object(*this); + ar & start_; + } +}; + +// real variable axis (varying bin widths) +class variable_axis : public axis_base, public real_axis { +public: + template + variable_axis(const Container& x, + const std::string& label = std::string(), + bool uoflow = true) : + axis_base(x.size() - 1, label, uoflow), + x_(new double[x.size()]) + { + std::copy(x.begin(), x.end(), x_.get()); + std::sort(x_.get(), x_.get() + bins() + 1); + } + + template + variable_axis(const Iterator& begin, const Iterator& end, + const std::string& label = std::string(), + bool uoflow = true) : + axis_base(end - begin - 1, label, uoflow), + x_(new double[end - begin]) + { + std::copy(begin, end, x_.get()); + std::sort(x_.get(), x_.get() + bins() + 1); + } + + variable_axis() {} + explicit variable_axis(const variable_axis&); + variable_axis& operator=(const variable_axis&); + + inline int index(double x) const { + return std::upper_bound(x_.get(), x_.get() + bins() + 1, x) + - x_.get() - 1; + } + + double operator[](int idx) const; + bool operator==(const variable_axis&) const; +private: + boost::scoped_array x_; + + friend class serialization::access; + template + void serialize(Archive& ar, unsigned version) + { + ar & boost::serialization::base_object(*this); + if (Archive::is_loading::value) + x_.reset(new double[bins() + 1]); + ar & serialization::make_array(x_.get(), bins() + 1); + } +}; + +class category_axis { +public: + typedef std::string value_type; + + category_axis(const std::string&); + category_axis(const std::vector&); + + category_axis() {} + explicit category_axis(const category_axis&); + category_axis& operator=(const category_axis&); + + inline unsigned bins() const { return categories_.size(); } + const std::string& operator[](int idx) const + { return categories_[idx]; } + + bool operator==(const category_axis&) const; +private: + std::vector categories_; + + friend class serialization::access; + template + void serialize(Archive& ar, unsigned version) + { + using namespace serialization; + ar & categories_; + } +}; + +class integer_axis: public axis_base { +public: + typedef int value_type; + + integer_axis(int min, int max, + const std::string& label = std::string(), + bool uoflow = true); + + integer_axis() {} + explicit integer_axis(const integer_axis&); + integer_axis& operator=(const integer_axis&); + + inline int index(double x) const + { return algorithm::clamp(rint(x) - min_, -1, bins()); } + int operator[](int idx) const { return min_ + idx; } + + bool operator==(const integer_axis&) const; +private: + int min_; + + friend class serialization::access; + template + void serialize(Archive& ar, unsigned version) + { + using namespace serialization; + ar & boost::serialization::base_object(*this); + ar & min_; + } +}; + +typedef variant< + regular_axis, // most common type + polar_axis, + variable_axis, + category_axis, + integer_axis +> axis_type; + +typedef std::vector axes_type; + +std::ostream& operator<<(std::ostream&, const axis_type&); + +} +} + +#endif diff --git a/include/boost/histogram/detail/zero_suppression.hpp b/include/boost/histogram/detail/zero_suppression.hpp new file mode 100644 index 00000000..83be148b --- /dev/null +++ b/include/boost/histogram/detail/zero_suppression.hpp @@ -0,0 +1,23 @@ +#ifndef _BOOST_HISTOGRAM_DETAIL_ZERO_SUPPRESSION_HPP_ +#define _BOOST_HISTOGRAM_DETAIL_ZERO_SUPPRESSION_HPP_ + +#include +#include + +namespace boost { +namespace histogram { +namespace detail { + +bool +zero_suppression_encode(std::vector& output, std::size_t output_max, + const char* input, std::size_t input_len); + +void +zero_suppression_decode(char* output, std::size_t output_len, + const std::vector& input); + +} +} +} + +#endif diff --git a/include/boost/histogram/histogram_base.hpp b/include/boost/histogram/histogram_base.hpp new file mode 100644 index 00000000..35a87f28 --- /dev/null +++ b/include/boost/histogram/histogram_base.hpp @@ -0,0 +1,121 @@ +#ifndef _BOOST_HISTOGRAM_BASE_HPP_ +#define _BOOST_HISTOGRAM_BASE_HPP_ + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#define BOOST_HISTOGRAM_AXIS_LIMIT 16 + +namespace boost { +namespace histogram { + +// holds an array of axis and computes the internal index +class histogram_base { +public: + typedef uint64_t size_type; + + histogram_base(const histogram_base&); + histogram_base& operator=(const histogram_base&); + ~histogram_base() {} + + unsigned dim() const { return axes_.size(); } + + template + T& axis(unsigned i) + { if (is_same::value) return axes_[i]; + return boost::get(axes_[i]); } + + template + const T& axis(unsigned i) const + { if (is_same::value) return axes_[i]; + return boost::get(axes_[i]); } + + unsigned bins(unsigned i) const { return size_[i]; } + unsigned shape(unsigned i) const { return size_[i] + 2 * uoflow_[i]; } + +protected: + histogram_base() {} + + explicit histogram_base(const axes_type& axes); + + explicit histogram_base(const axis_type& a) : axes_(1, a) { update_buffers(); } + +#define BOOST_HISTOGRAM_BASE_APPEND(z, n, unused) axes_.push_back(a ## n); +#define BOOST_HISTOGRAM_BASE_CTOR(z, n, unused) \ + histogram_base( BOOST_PP_ENUM_PARAMS_Z(z, n, const axis_type& a) ) \ + { \ + axes_.reserve(n); \ + BOOST_PP_REPEAT(n, BOOST_HISTOGRAM_BASE_APPEND, unused) \ + update_buffers(); \ + } + +// generates constructors taking 2 to AXIS_LIMIT arguments +BOOST_PP_REPEAT_FROM_TO(2, BOOST_HISTOGRAM_AXIS_LIMIT, BOOST_HISTOGRAM_BASE_CTOR, nil) + + bool operator==(const histogram_base&) const; + bool operator!=(const histogram_base& o) const + { return !operator==(o); } + + template + inline + size_type pos(const Array& v) const { + int32_t idx[BOOST_HISTOGRAM_AXIS_LIMIT]; + for (unsigned i = 0, n = axes_.size(); i < n; ++i) + idx[i] = apply_visitor(detail::index_visitor(v[i]), axes_[i]); + return linearize(idx); + } + + inline + size_type linearize(const int* idx) const { + size_type stride = 1, k = 0, i = axes_.size(); + while (i) { + --i; + const int size = size_[i]; + const int range = size + 2 * uoflow_[i]; + int j = idx[i]; + // the following three lines work for any overflow setting + j += (j < 0) * (size + 2); // wrap around if j < 0 + if (j >= range) + return size_type(-1); // indicate out of range + k += j * stride; + stride *= range; + } + return k; + } + + // compute the number of fields needed for storage + size_type field_count() const; + +private: + axes_type axes_; + + // internal buffers + int32_t size_[BOOST_HISTOGRAM_AXIS_LIMIT]; + std::bitset uoflow_; + + void update_buffers(); ///< fills size_ and uoflow_ + + friend class serialization::access; + template + void serialize(Archive& ar, unsigned version) + { + using namespace serialization; + ar & axes_; + if (Archive::is_loading::value) + update_buffers(); + } +}; + +} +} + +#endif diff --git a/include/boost/histogram/nhistogram.hpp b/include/boost/histogram/nhistogram.hpp new file mode 100644 index 00000000..c423400e --- /dev/null +++ b/include/boost/histogram/nhistogram.hpp @@ -0,0 +1,111 @@ +#ifndef _BOOST_HISTOGRAM_NHISTOGRAM_HPP_ +#define _BOOST_HISTOGRAM_NHISTOGRAM_HPP_ + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace boost { +namespace histogram { + +class nhistogram : public histogram_base { +public: + nhistogram() {} + + explicit + nhistogram(const axes_type& axes) : + histogram_base(axes), + data_(field_count()) + {} + + explicit + nhistogram(const axis_type& a) : + histogram_base(a), + data_(field_count()) + {} + +#define BOOST_NHISTOGRAM_CTOR(z, n, unused) \ + nhistogram( BOOST_PP_ENUM_PARAMS_Z(z, n, const axis_type& a) ) : \ + histogram_base( BOOST_PP_ENUM_PARAMS_Z(z, n, a) ), \ + data_(field_count()) \ + {} + +// generates constructors taking 2 to AXIS_LIMIT arguments +BOOST_PP_REPEAT_FROM_TO(2, BOOST_HISTOGRAM_AXIS_LIMIT, BOOST_NHISTOGRAM_CTOR, nil) + + double sum() const; + + template + inline + void fill(const Array& v) + { + const size_type k = pos(v); + if (k != uintmax_t(-1)) + data_.increase(k); + } + + void fill(double x, + BOOST_PP_ENUM_PARAMS_WITH_A_DEFAULT(BOOST_PP_DEC(BOOST_HISTOGRAM_AXIS_LIMIT), + double x, 0.0)) + { + const double buffer[BOOST_HISTOGRAM_AXIS_LIMIT] = { + x, BOOST_PP_ENUM_PARAMS(BOOST_PP_DEC(BOOST_HISTOGRAM_AXIS_LIMIT), x) + }; + fill(buffer); + } + + template + inline + size_type operator()(const Array& idx) + const + { return data_.read(linearize(idx)); } + + size_type operator()(int i, + BOOST_PP_ENUM_PARAMS_WITH_A_DEFAULT(BOOST_PP_DEC(BOOST_HISTOGRAM_AXIS_LIMIT), + int i, 0)) + const + { + const int32_t idx[BOOST_HISTOGRAM_AXIS_LIMIT] = { + i, BOOST_PP_ENUM_PARAMS(BOOST_PP_DEC(BOOST_HISTOGRAM_AXIS_LIMIT), i) + }; + return operator()(idx); + } + + bool operator==(const nhistogram& o) const + { return histogram_base::operator==(o) && + data_ == o.data_; } + + nhistogram& operator+=(const nhistogram& o) + { data_ += o.data_; return *this; } + + // needed by boost::python interface to pass internal data to numpy + const nstore& data() const { return data_; } + +private: + nstore data_; + + friend class serialization::access; + template + void serialize(Archive& ar, unsigned version) + { + using namespace serialization; + ar & boost::serialization::base_object(*this); + ar & data_; + } +}; + +nhistogram operator+(const nhistogram& a, const nhistogram& b) { + nhistogram result(a); + return result += b; +} + +} +} + +#endif diff --git a/include/boost/histogram/nstore.hpp b/include/boost/histogram/nstore.hpp new file mode 100644 index 00000000..3b818200 --- /dev/null +++ b/include/boost/histogram/nstore.hpp @@ -0,0 +1,112 @@ +#ifndef _BOOST_HISTOGRAM_NSTORE_HPP_ +#define _BOOST_HISTOGRAM_NSTORE_HPP_ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include // for bad:alloc + +namespace boost { +namespace histogram { + +class nstore { +public: + typedef uint64_t size_type; + typedef unsigned depth_type; + +public: + nstore(); + nstore(const nstore&); + nstore(size_type, depth_type d = sizeof(uint8_t)); + ~nstore() { destroy(); } + + nstore& operator=(const nstore&); + nstore& operator+=(const nstore&); + bool operator==(const nstore&) const; + + uint64_t read(size_type) const; + void write(size_type, uint64_t); + + inline void increase(size_type i) { + switch (depth_) { + #define BOOST_HISTOGRAM_NSTORE_INC(T) \ + case sizeof(T): { \ + T& b = ((T*)buffer_)[i]; \ + if (b == std::numeric_limits::max()) \ + grow(); /* and fall to next case */ \ + else { ++b; break; } \ + } + BOOST_HISTOGRAM_NSTORE_INC(uint8_t); + BOOST_HISTOGRAM_NSTORE_INC(uint16_t); + BOOST_HISTOGRAM_NSTORE_INC(uint32_t); + BOOST_HISTOGRAM_NSTORE_INC(uint64_t); + #undef BOOST_HISTOGRAM_NSTORE_INC + default: assert(!"invalid depth"); + } + } + + const void* buffer() const { return buffer_; } + depth_type depth() const { return depth_; } + +private: + size_type size_; + depth_type depth_; + void* buffer_; + + void create(); + void destroy(); + void grow(); + + uint64_t max_count() const; + + friend class serialization::access; + template + void serialize(Archive& ar, unsigned version) + { + const size_type s = size_; + const depth_type d = depth_; + ar & size_; + ar & depth_; + if (s != size_ || d != depth_) { + // realloc is safe if buffer_ is null + buffer_ = std::realloc(buffer_, size_ * depth_); + } + if (buffer_ == 0 && size_ > 0) + throw std::bad_alloc(); + + if (Archive::is_saving::value) { + std::vector buf; + if (detail::zero_suppression_encode(buf, size_ * depth_, + (char*)buffer_, size_ * depth_)) { + bool is_zero_suppressed = true; + ar & is_zero_suppressed; + ar & buf; + } else { + bool is_zero_suppressed = false; + ar & is_zero_suppressed; + ar & serialization::make_array((char*)buffer_, size_ * depth_); + } + } else { + bool is_zero_suppressed = false; + ar & is_zero_suppressed; + if (is_zero_suppressed) { + std::vector buf; + ar & buf; + detail::zero_suppression_decode((char*)buffer_, size_ * depth_, buf); + } else { + ar & serialization::make_array((char*)buffer_, size_ * depth_); + } + } + } +}; + +} +} + +#endif diff --git a/include/boost/histogram/vector_operators.hpp b/include/boost/histogram/vector_operators.hpp new file mode 100644 index 00000000..84c29666 --- /dev/null +++ b/include/boost/histogram/vector_operators.hpp @@ -0,0 +1,52 @@ +#ifndef _BOOST_HISTOGRAM_VECTOR_OPERATORS_HPP_ +#define _BOOST_HISTOGRAM_VECTOR_OPERATORS_HPP_ + +#include + +namespace boost { +namespace histogram { + + // vectors: generic == + template + bool operator==(const std::vector& a, const std::vector& b) { + if (a.size() != b.size()) + return false; + // pointer arithmetric is faster + const T* da = &a.front(); + const T* db = &b.front(); + for (unsigned i = 0, n = a.size(); i < n; ++i) + if (da[i] != db[i]) + return false; + return true; + } + + // vectors: generic += + template + std::vector& operator+=(std::vector& a, const std::vector& b) { + if (a.size() != b.size()) + throw std::invalid_argument("sizes do not match"); + // pointer arithmetric is faster + T* da = &a.front(); + const T* db = &b.front(); + for (unsigned i = 0, n = a.size(); i < n; ++i) + da[i] += db[i]; + return a; + } + + // vectors: generic -= + template + std::vector& operator-=(std::vector& a, const std::vector& b) { + if (a.size() != b.size()) + throw std::invalid_argument("sizes do not match"); + // pointer arithmetric is faster + T* da = &a.front(); + const T* db = &b.front(); + for (unsigned i = 0, n = a.size(); i < n; ++i) + da[i] -= db[i]; + return a; + } + +} +} + +#endif diff --git a/include/boost/histogram/visitors.hpp b/include/boost/histogram/visitors.hpp new file mode 100644 index 00000000..c249c0cc --- /dev/null +++ b/include/boost/histogram/visitors.hpp @@ -0,0 +1,57 @@ +#ifndef _BOOST_HISTOGARM_VISITORS_HPP_ +#define _BOOST_HISTOGARM_VISITORS_HPP_ + +#include +#include + +namespace boost { +namespace histogram { +namespace detail { + + struct bins_visitor : public static_visitor + { + template + unsigned operator()(const A& a) const { return a.bins(); } + }; + + struct fields_visitor : public static_visitor + { + unsigned operator()(const category_axis& a) const { return a.bins(); } + + template + unsigned operator()(const A& a) const { return a.bins() + 2 * a.uoflow(); } + }; + + struct uoflow_visitor : public static_visitor + { + bool operator()(const category_axis& a) const { return false; } + + template + bool operator()(const A& a) const { return a.uoflow(); } + }; + + struct index_visitor : public static_visitor + { + double x; + index_visitor(double d) : x(d) {} + + int operator()(const category_axis& a) const { return int(x + 0.5); } + + template + int operator()(const A& a) const { return a.index(x); } + }; + + struct cmp_visitor : public static_visitor + { + template + bool operator()(const T& a, const T& b) const + { return a == b; } + + template + bool operator()(const T&, const U&) const { return false; } + }; +} +} +} + +#endif diff --git a/include/boost/histogram/wstore.hpp b/include/boost/histogram/wstore.hpp new file mode 100644 index 00000000..20cf767c --- /dev/null +++ b/include/boost/histogram/wstore.hpp @@ -0,0 +1,110 @@ +#include +#include +#include +#include + +class wstore { +public: + typedef boost::uintmax_t size_type; + struct wcount { double w_, w2_; }; + +private: + size_type size_; + double* buffer_; + + friend class boost::serialization::access; + template + void serialize(Archive & ar, unsigned version) { + using boost::serialization::make_array; + size_type s = size_; + ar & make_nvp("size", size_); + if (s != size_) { + destroy(); + create(); + } + ar & make_nvp("data", \ + make_array(buffer_, size_ * 2)); + } + + void + create() { + buffer_ = (double*)std::calloc(size_, sizeof(wcount)); + } + + void + destroy() { + free(buffer_); + buffer_ = 0; + } + +public: + wstore() + : size_(0), buffer_(0) + {} + + wstore(size_type n) + : size_(n) + { + create(); + } + + wstore(const wstore& other) + : size_(other.size_) + { + create(); + *this = other; + } + + ~wstore() { destroy(); } + + wstore& + operator=(const wstore& other) + { + if (size_ != other.size_) { + destroy(); + size_ = other.size_; + create(); + } + std::memcpy(buffer_, other.buffer_, size_ * sizeof(wcount)); + return *this; + } + + inline + wcount& + operator[](size_type i) + { + return *((wcount*)buffer_ + i); + } + + inline + const wcount& + operator[](size_type i) const + { + return *((wcount*)buffer_ + i); + } + + wstore& + operator+=(const wstore& other) { + if (size_ != other.size_) + throw std::logic_error("sizes do not match"); + + for (size_type i = 0; i < (size_ * 2); ++i) + buffer_[i] += other.buffer_[i]; + return *this; + } + + bool + operator==(const wstore& other) const { + if (size_ != other.size_) + return false; + + for (size_type i = 0; i < (size_ * 2); ++i) + if (buffer_[i] != other.buffer_[i]) + return false; + + return true; + } + + inline size_type size() const { return size_; } + inline const double* buffer() const { return buffer_; } +}; diff --git a/src/axis.cpp b/src/axis.cpp new file mode 100644 index 00000000..92338010 --- /dev/null +++ b/src/axis.cpp @@ -0,0 +1,269 @@ +#include +#include +#include +#include + +namespace boost { +namespace histogram { + +axis_base::axis_base(int size, + const std::string& label, + bool uoflow) : + size_(size), + label_(label) +{ + if (size <= 0) + throw std::logic_error("size must be positive"); + if (!uoflow) size_ = -size_; +} + +axis_base::axis_base(const axis_base& o) : + size_(o.size_), + label_(o.label_) +{} + +axis_base& +axis_base::operator=(const axis_base& o) +{ + size_ = o.size_; + label_ = o.label_; + return *this; +} + +bool axis_base::operator==(const axis_base& o) const +{ return size_ == o.size_ && label_ == o.label_; } + +regular_axis::regular_axis(unsigned n, double min, double max, + const std::string& label, bool uoflow): + axis_base(n, label, uoflow), + min_(min), + range_(max - min) +{ + if (min >= max) + throw std::logic_error("regular_axis: min must be less than max"); +} + +regular_axis::regular_axis(const regular_axis& o) : + axis_base(o), + min_(o.min_), + range_(o.range_) +{} + +regular_axis& +regular_axis::operator=(const regular_axis& o) +{ + axis_base::operator=(o); + min_ = o.min_; + range_ = o.range_; + return *this; +} + +double +regular_axis::operator[](int idx) + const +{ + if (idx < 0) + return -std::numeric_limits::infinity(); + if (idx > bins()) + return std::numeric_limits::infinity(); + const double z = double(idx) / bins(); + return (1.0 - z) * min_ + z * (min_ + range_); +} + +bool +regular_axis::operator==(const regular_axis& o) const +{ + return axis_base::operator==(o) && + min_ == o.min_ && + range_ == o.range_; +} + +polar_axis::polar_axis(unsigned n, double start, + const std::string& label) : + axis_base(n, label, false), + start_(start) +{} + +polar_axis::polar_axis(const polar_axis& o) : + axis_base(o), + start_(o.start_) +{} + +polar_axis& +polar_axis::operator=(const polar_axis& o) +{ + axis_base::operator=(o); + start_ = o.start_; + return *this; +} + +double +polar_axis::operator[](int idx) + const +{ + const double z = double(idx) / bins(); + return z * 6.283185307179586 + start_; +} + +bool +polar_axis::operator==(const polar_axis& o) const +{ + return axis_base::operator==(o) && + start_ == o.start_; +} + +variable_axis::variable_axis(const variable_axis& o) : + axis_base(o), + x_(new double[bins() + 1]) +{ + std::copy(o.x_.get(), o.x_.get() + bins() + 1, x_.get()); +} + +variable_axis& +variable_axis::operator=(const variable_axis& o) +{ + axis_base::operator=(o); + x_.reset(new double[bins() + 1]); + std::copy(o.x_.get(), o.x_.get() + bins() + 1, x_.get()); + return *this; +} + +double +variable_axis::operator[](int idx) + const +{ + if (idx < 0) + return -std::numeric_limits::infinity(); + if (idx > bins()) + return std::numeric_limits::infinity(); + return x_[idx]; +} + +bool +variable_axis::operator==(const variable_axis& o) const +{ + if (!axis_base::operator==(o)) + return false; + for (unsigned i = 0, n = bins() + 1; i < n; ++i) + if (x_[i] != o.x_[i]) + return false; + return true; +} + +category_axis::category_axis(const std::string& s) +{ + std::size_t i = s.find(';'); + categories_.push_back(s.substr(0, i)); + while (i != std::string::npos) { + const std::size_t p = i + 1; + i = s.find(';', p); + categories_.push_back(s.substr(p, i - p)); + } +} + +category_axis::category_axis(const std::vector& c) : + categories_(c) +{} + +category_axis::category_axis(const category_axis& o) : + categories_(o.categories_) +{} + +category_axis& +category_axis::operator=(const category_axis& o) +{ + categories_ = o.categories_; + return *this; +} + +bool +category_axis::operator==(const category_axis& o) const +{ return categories_ == o.categories_; } + +integer_axis::integer_axis(int min, int max, + const std::string& label, + bool uoflow) : + axis_base(max + 1 - min, label, uoflow), + min_(min) +{} + +integer_axis::integer_axis(const integer_axis& a) : + axis_base(a), + min_(a.min_) +{} + +integer_axis& +integer_axis::operator=(const integer_axis& o) +{ + axis_base::operator=(o); + min_ = o.min_; + return *this; +} + +bool +integer_axis::operator==(const integer_axis& o) const +{ + return axis_base::operator==(o) && + min_ == o.min_; +} + +struct stream_visitor : public static_visitor +{ + std::string operator()(const category_axis& a) const { + std::stringstream line; + line << "category_axis("; + for (int i = 0; i < a.bins(); ++i) + line << a[i] << (i == (a.bins() - 1)? ")" : ", "); + return line.str(); + } + + std::string operator()(const integer_axis& a) const { + std::stringstream line; + line << "integer_axis("; + if (a.bins()) + line << a[0] << "," << a[a.bins() - 1]; + line << ")"; + return line.str(); + } + + std::string operator()(const polar_axis& a) const { + std::stringstream line; + line << "polar_axis("; + if (!a.label().empty()) + line << a.label() << ", "; + line << a.bins() << ":"; + for (int i = 0; i <= a.bins(); ++i) + line << " " << a.left(i); + return line.str(); + } + + std::string operator()(const regular_axis& a) const { + std::stringstream line; + line << "regular_axis["; + if (!a.label().empty()) + line << a.label() << ", "; + line << a.bins() << ":"; + for (int i = 0; i <= a.bins(); ++i) + line << " " << a.left(i); + return line.str(); + } + + std::string operator()(const variable_axis& a) const { + std::stringstream line; + line << "variable_axis["; + if (!a.label().empty()) + line << a.label() << ", "; + line << a.bins() << ":"; + for (int i = 0; i <= a.bins(); ++i) + line << " " << a.left(i); + return line.str(); + } +}; + +std::ostream& operator<<(std::ostream& os, const axis_type& a) { + os << apply_visitor(stream_visitor(), a); + return os; +} + +} +} diff --git a/src/histogram_base.cpp b/src/histogram_base.cpp new file mode 100644 index 00000000..a1b3200d --- /dev/null +++ b/src/histogram_base.cpp @@ -0,0 +1,63 @@ +#include +#include +#include + +namespace boost { +namespace histogram { + +histogram_base::histogram_base(const histogram_base& o) : + axes_(o.axes_) +{ + update_buffers(); +} + +histogram_base::histogram_base(const axes_type& axes) : + axes_(axes) +{ + if (axes_.size() > BOOST_HISTOGRAM_AXIS_LIMIT) + throw std::invalid_argument("too many axes"); + update_buffers(); +} + +histogram_base& +histogram_base::operator=(const histogram_base& o) +{ + axes_ = o.axes_; + update_buffers(); + return *this; +} + +bool +histogram_base::operator==(const histogram_base& o) + const +{ + if (axes_.size() != o.axes_.size()) + return false; + for (unsigned i = 0; i < axes_.size(); ++i) { + if (!apply_visitor(detail::cmp_visitor(), axes_[i], o.axes_[i])) + return false; + } + return true; +} + +histogram_base::size_type +histogram_base::field_count() + const +{ + size_type fc = 1; + for (unsigned i = 0, n = axes_.size(); i < n; ++i) + fc *= apply_visitor(detail::fields_visitor(), axes_[i]); + return fc; +} + +void +histogram_base::update_buffers() +{ + for (unsigned i = 0, n = axes_.size(); i < n; ++i) { + size_[i] = apply_visitor(detail::bins_visitor(), axes_[i]); + uoflow_[i] = apply_visitor(detail::uoflow_visitor(), axes_[i]); + } +} + +} +} diff --git a/src/nhistogram.cpp b/src/nhistogram.cpp new file mode 100644 index 00000000..97a5c3f5 --- /dev/null +++ b/src/nhistogram.cpp @@ -0,0 +1,18 @@ +#include +#include + +namespace boost { +namespace histogram { + +double +nhistogram::sum() + const +{ + double result = 0.0; + for (size_type i = 0, n = field_count(); i < n; ++i) + result += data_.read(i); + return result; +} + +} +} diff --git a/src/nstore.cpp b/src/nstore.cpp new file mode 100644 index 00000000..9ecbd81b --- /dev/null +++ b/src/nstore.cpp @@ -0,0 +1,157 @@ +#include +#include +#include +#include // for std::bad_alloc + +namespace boost { +namespace histogram { + +nstore::nstore() : + size_(0), + depth_(0), + buffer_(0) +{} + +nstore::nstore(const nstore& o) : + size_(o.size_), + depth_(o.depth_) +{ + create(); + std::memcpy(buffer_, o.buffer_, size_ * depth_); +} + +nstore::nstore(size_type n, depth_type d) : + size_(n), + depth_(d) +{ + if (d == 0) + throw std::invalid_argument("depth may not be zero"); + if (d > sizeof(uint64_t)) + throw std::invalid_argument("depth > sizeof(uint64_t) is not supported"); + create(); +} + +nstore& +nstore::operator=(const nstore& o) +{ + if (size_ != o.size_ || depth_ != o.depth_) { + destroy(); + size_ = o.size_; + depth_ = o.depth_; + create(); + } + std::memcpy(buffer_, o.buffer_, size_ * depth_); + return *this; +} + +nstore& +nstore::operator+=(const nstore& o) +{ + if (size_ != o.size_) + throw std::logic_error("sizes do not match"); + for (size_type i = 0; i < size_; ++i) + write(i, read(i) + o.read(i)); + return *this; +} + +bool +nstore::operator==(const nstore& o) + const +{ + if (size_ != o.size_ || depth_ != o.depth_) + return false; + return std::memcmp(buffer_, o.buffer_, size_ * depth_) == 0; +} + +uint64_t +nstore::read(size_type i) + const +{ + switch (depth_) { + #define BOOST_HISTOGRAM_NSTORE_READ(T) \ + case sizeof(T): return ((T*)buffer_)[i] + BOOST_HISTOGRAM_NSTORE_READ(uint8_t); + BOOST_HISTOGRAM_NSTORE_READ(uint16_t); + BOOST_HISTOGRAM_NSTORE_READ(uint32_t); + BOOST_HISTOGRAM_NSTORE_READ(uint64_t); + #undef BOOST_HISTOGRAM_NSTORE_READ + default: assert(!"invalid depth"); + } + return 0; +} + +void +nstore::write(size_type i, uint64_t v) +{ + const uint64_t vmax = max_count(); + while (vmax < v) grow(); + + switch (depth_) { + #define BOOST_HISTOGRAM_NSTORE_WRITE(T) \ + case sizeof(T): { \ + ((T*)buffer_)[i] = v; \ + break; \ + } + BOOST_HISTOGRAM_NSTORE_WRITE(uint8_t); + BOOST_HISTOGRAM_NSTORE_WRITE(uint16_t); + BOOST_HISTOGRAM_NSTORE_WRITE(uint32_t); + BOOST_HISTOGRAM_NSTORE_WRITE(uint64_t); + #undef BOOST_HISTOGRAM_NSTORE_WRITE + default: assert(!"invalid depth"); + } +} + +void +nstore::create() +{ + buffer_ = std::calloc(size_, depth_); + if (!buffer_) throw std::bad_alloc(); +} + +void +nstore::destroy() +{ + std::free(buffer_); buffer_ = 0; +} + +void +nstore::grow() +{ + if (depth_ == sizeof(uint64_t)) + throw std::overflow_error("depth > 64 bit is not supported"); + if (depth_ == 0 || buffer_ == 0) + throw std::logic_error("cannot call grow on null buffer"); + buffer_ = std::realloc(buffer_, size_ * 2 * depth_); + if (!buffer_) throw std::bad_alloc(); + switch (depth_) { + #define BOOST_HISTOGRAM_NSTORE_GROW(T0, T1) \ + case sizeof(T0): \ + for (size_type i = size_ - 1; i != size_type(-1); --i) \ + ((T1*)buffer_)[i] = ((T0*)buffer_)[i]; \ + break + BOOST_HISTOGRAM_NSTORE_GROW(uint8_t, uint16_t); + BOOST_HISTOGRAM_NSTORE_GROW(uint16_t, uint32_t); + BOOST_HISTOGRAM_NSTORE_GROW(uint32_t, uint64_t); + #undef BOOST_HISTOGRAM_NSTORE_GROW + } + depth_ *= 2; +} + +uint64_t +nstore::max_count() + const +{ + switch (depth_) { + #define BOOST_HISTOGRAM_NSTORE_CASE(T) \ + case sizeof(T): return std::numeric_limits::max() + BOOST_HISTOGRAM_NSTORE_CASE(uint8_t); + BOOST_HISTOGRAM_NSTORE_CASE(uint16_t); + BOOST_HISTOGRAM_NSTORE_CASE(uint32_t); + BOOST_HISTOGRAM_NSTORE_CASE(uint64_t); + default: assert(!"invalid depth"); + } + return 0; +} + +} +} diff --git a/src/python/histogram_base.cpp b/src/python/histogram_base.cpp new file mode 100644 index 00000000..626c5407 --- /dev/null +++ b/src/python/histogram_base.cpp @@ -0,0 +1,325 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace boost { +namespace histogram { + +namespace { + +python::object +variable_axis_init(python::tuple args, python::dict kwargs) { + using namespace python; + using python::tuple; + + object self = args[0]; + object pyinit = self.attr("__init__"); + + if (len(args) < 2) { + PyErr_SetString(PyExc_TypeError, "require at least two arguments"); + throw_error_already_set(); + } + + std::vector v; + for (int i = 1, n = len(args); i < n; ++i) { + v.push_back(extract(args[i])); + } + + std::string label; + bool uoflow = true; + while (len(kwargs) > 0) { + tuple kv = kwargs.popitem(); + std::string k = extract(kv[0]); + object v = kv[1]; + if (k == "label") + label = extract(v); + else if (k == "uoflow") + uoflow = extract(v); + else { + std::stringstream s; + s << "keyword " << k << " not recognized"; + PyErr_SetString(PyExc_KeyError, s.str().c_str()); + throw_error_already_set(); + } + } + return pyinit(v, label, uoflow); +} + +python::object +category_axis_init(python::tuple args, python::dict kwargs) { + using namespace python; + + object self = args[0]; + object pyinit = self.attr("__init__"); + + if (len(args) == 1) { + PyErr_SetString(PyExc_TypeError, "require at least one argument"); + throw_error_already_set(); + } + + if (len(kwargs) > 0) { + PyErr_SetString(PyExc_TypeError, "unknown keyword argument"); + throw_error_already_set(); + } + + if (len(args) == 2) { + extract es(args[1]); + if (es.check()) + pyinit(es); + else { + PyErr_SetString(PyExc_TypeError, "require one or several string arguments"); + throw_error_already_set(); + } + } + + std::vector c; + for (int i = 1, n = len(args); i < n; ++i) + c.push_back(extract(args[i])); + + return pyinit(c); +} + +// python::object +// regular_axis_data(const regular_axis& self) { +// #if HAVE_NUMPY +// py_intp dims[1] = { self.size() + 1 }; +// PyArrayObject* a = (PyArrayObject*)PyArray_SimpleNew(1, dims, NPY_DOUBLE); +// for (unsigned i = 0; i < self.size() + 1; ++i) { +// double* pv = (double*)PyArray_GETPTR1(a, i); +// *pv = self.left(i); +// } +// return python::object(handle<>((PyObject*)a)); +// #else +// list v; +// for (int i = 0; i < self.size() + 1; ++i) +// v.append(self.left(i)); +// return v; +// #endif +// } + +struct axis_visitor : public static_visitor +{ + template + python::object operator()(const T& t) const { return python::object(T(t)); } +}; + +template +unsigned +axis_len(const T& t) { + return t.bins() + 1; +} + +template <> +unsigned +axis_len(const category_axis& t) { + return t.bins(); +} + +template <> +unsigned +axis_len(const integer_axis& t) { + return t.bins(); +} + +template +typename T::value_type +axis_getitem(const T& t, int i) { + if (i == axis_len(t)) { + PyErr_SetString(PyExc_StopIteration, "no more"); + python::throw_error_already_set(); + } + return t[i]; +} + +template +struct has_index_method +{ + struct yes { char x[1]; }; + struct no { char x[2]; }; + template struct SFINAE {}; + template static yes test( SFINAE* ); + template static no test( ... ); + enum { value = sizeof(test(0)) == sizeof(yes) }; +}; + +std::string escape(const std::string& s) { + std::ostringstream os; + os << '\''; + for (unsigned i = 0; i < s.size(); ++i) { + const char c = s[i]; + if (c == '\'' && (i == 0 || s[i - 1] != '\\')) + os << "\\\'"; + else + os << c; + } + os << '\''; + return os.str(); +} + +std::string axis_repr(const regular_axis& a) { + std::stringstream s; + s << "regular_axis(" << a.bins() << ", " << a[0] << ", " << a[a.bins()]; + if (!a.label().empty()) + s << ", label=" << escape(a.label()); + if (!a.uoflow()) + s << ", uoflow=False"; + s << ")"; + return s.str(); +} + +std::string axis_repr(const polar_axis& a) { + std::stringstream s; + s << "polar_axis(" << a.bins(); + if (a[0] != 0.0) + s << ", " << a[0]; + if (!a.label().empty()) + s << ", label=" << escape(a.label()); + s << ")"; + return s.str(); +} + +std::string axis_repr(const variable_axis& a) { + std::stringstream s; + s << "variable_axis(" << a[0]; + for (int i = 1; i <= a.bins(); ++i) + s << ", " << a.left(i); + if (!a.label().empty()) + s << ", label=" << escape(a.label()); + if (!a.uoflow()) + s << ", uoflow=False"; + s << ")"; + return s.str(); +} + +std::string axis_repr(const category_axis& a) { + std::stringstream s; + s << "category_axis("; + for (int i = 0; i < a.bins(); ++i) + s << escape(a[i]) << (i == (a.bins() - 1)? ")" : ", "); + return s.str(); +} + +std::string axis_repr(const integer_axis& a) { + std::stringstream s; + s << "integer_axis(" << a[0] << ", " << a[a.bins() - 1]; + if (!a.label().empty()) + s << ", label=" << escape(a.label()); + if (!a.uoflow()) + s << ", uoflow=False"; + s << ")"; + return s.str(); +} + +template +struct axis_suite : public python::def_visitor > { + + template + static + typename enable_if_c::value, void>::type + add_axis_index(Class& cl) { + cl.def("index", &U::index); + } + + template + static + typename disable_if_c::value, void>::type + add_axis_index(Class& cl) {} + + template + static + typename enable_if, void>::type + add_axis_label(Class& cl) { + cl.add_property("label", + python::make_function((const std::string&(U::*)() const) &U::label, + python::return_value_policy()), + (void(U::*)(const std::string&)) &U::label); + } + + template + static + typename disable_if, void>::type + add_axis_label(Class& cl) {} + + template + static void + visit(Class& cl) + { + cl.add_property("bins", &T::bins); + add_axis_index(cl); + add_axis_label(cl); + cl.def("__len__", axis_len); + cl.def("__getitem__", axis_getitem); + cl.def("__repr__", (std::string (*)(const T&)) &axis_repr); + cl.def(python::self == python::self); + } +}; + +python::object +histogram_base_axis(const histogram_base& self, unsigned i) +{ + return apply_visitor(axis_visitor(), self.axis(i)); +} + +} // namespace + +void register_histogram_base() { + using namespace python; + using python::arg; + + // used to pass arguments from raw python init to specialized C++ constructors + class_ >("vector_double", no_init); + class_ >("vector_string", no_init); + class_("axes", no_init); + + class_("regular_axis", no_init) + .def(init( + (arg("bin"), arg("min"), arg("max"), + arg("label") = std::string(), + arg("uoflow") = true))) + .def(axis_suite()) + ; + + class_("polar_axis", no_init) + .def(init( + (arg("bin"), arg("start") = 0.0, + arg("label") = std::string()))) + .def(axis_suite()) + ; + + class_("variable_axis", no_init) + .def("__init__", raw_function(variable_axis_init)) + .def(init, std::string, bool>()) + .def(axis_suite()) + ; + + class_("category_axis", no_init) + .def("__init__", raw_function(category_axis_init)) + .def(init()) + .def(init >()) + .def(axis_suite()) + ; + + class_("integer_axis", no_init) + .def(init( + (arg("min"), arg("max"), + arg("label") = std::string(), + arg("uoflow") = true))) + .def(axis_suite()) + ; + + class_("histogram_base", no_init) + .add_property("dim", &histogram_base::dim) + .def("bins", &histogram_base::bins) + .def("shape", &histogram_base::shape) + .def("axis", histogram_base_axis) + ; +} + +} +} diff --git a/src/python/module.cpp b/src/python/module.cpp new file mode 100644 index 00000000..8b579a7a --- /dev/null +++ b/src/python/module.cpp @@ -0,0 +1,22 @@ +#include +#ifdef USE_NUMPY +# define PY_ARRAY_UNIQUE_SYMBOL boost_histogram_ARRAY_API +# define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION +# include +#endif + +namespace boost { +namespace histogram { + void register_histogram_base(); + void register_nhistogram(); +} +} + +BOOST_PYTHON_MODULE(histogram) +{ +#ifdef USE_NUMPY + import_array(); +#endif + boost::histogram::register_histogram_base(); + boost::histogram::register_nhistogram(); +} diff --git a/src/python/nhistogram.cpp b/src/python/nhistogram.cpp new file mode 100644 index 00000000..72642c15 --- /dev/null +++ b/src/python/nhistogram.cpp @@ -0,0 +1,183 @@ +#include +#include +#include +#include +#include +#include +#include "serialization_suite.hpp" +#ifdef USE_NUMPY +# define NO_IMPORT_ARRAY +# define PY_ARRAY_UNIQUE_SYMBOL boost_histogram_ARRAY_API +# define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION +# include +#endif + +namespace boost { +namespace histogram { + +python::object +nhistogram_init(python::tuple args, python::dict kwargs) { + using namespace python; + using python::tuple; + + object self = args[0]; + object pyinit = self.attr("__init__"); + + if (kwargs) { + PyErr_SetString(PyExc_TypeError, "no keyword arguments allowed"); + throw_error_already_set(); + } + + // normal constructor + axes_type axes; + for (unsigned i = 1, n = len(args); i < n; ++i) { + object pa = args[i]; + extract er(pa); + if (er.check()) { axes.push_back(er()); continue; } + extract ep(pa); + if (ep.check()) { axes.push_back(ep()); continue; } + extract ev(pa); + if (ev.check()) { axes.push_back(ev()); continue; } + extract ec(pa); + if (ec.check()) { axes.push_back(ec()); continue; } + extract ei(pa); + if (ei.check()) { axes.push_back(ei()); continue; } + PyErr_SetString(PyExc_TypeError, "require an axis object"); + throw_error_already_set(); + } + return pyinit(axes); +} + +python::dict +nhistogram_array_interface(nhistogram& self) { + using namespace python; + using python::tuple; + using python::make_tuple; + + list shape; + for (unsigned i = 0; i < self.dim(); ++i) + shape.append(self.shape(i)); + dict d; + d["shape"] = tuple(shape); + d["typestr"] = str("(args[0]); + +#ifdef USE_NUMPY + if (nargs == 2) { + object o = args[1]; + if (PySequence_Check(o.ptr())) { + PyArrayObject* a = (PyArrayObject*) + PyArray_FROM_OTF(o.ptr(), NPY_DOUBLE, NPY_ARRAY_IN_ARRAY); + if (!a) { + PyErr_SetString(PyExc_ValueError, "could not convert sequence into array"); + throw_error_already_set(); + } + + npy_intp* dims = PyArray_DIMS(a); + switch (PyArray_NDIM(a)) { + case 1: + if (self.dim() > 1) { + PyErr_SetString(PyExc_ValueError, "array has to be two-dimensional"); + throw_error_already_set(); + } + break; + case 2: + if (self.dim() != dims[1]) + { + PyErr_SetString(PyExc_ValueError, "size of second dimension does not match"); + throw_error_already_set(); + } + break; + default: + PyErr_SetString(PyExc_ValueError, "array has wrong dimension"); + throw_error_already_set(); + } + + for (unsigned i = 0; i < dims[0]; ++i) { + double* v = (double*)PyArray_GETPTR1(a, i); + self.fill(v); + } + + Py_DECREF(a); + return object(); + } + } +#endif + + const unsigned dim = nargs - 1; + if (dim != self.dim()) { + PyErr_SetString(PyExc_TypeError, "wrong number of arguments"); + throw_error_already_set(); + } + + if (kwargs) { + PyErr_SetString(PyExc_TypeError, "no keyword arguments allowed"); + throw_error_already_set(); + } + + double v[BOOST_HISTOGRAM_AXIS_LIMIT]; + for (unsigned i = 0; i < dim; ++i) + v[i] = extract(args[1 + i]); + self.fill(v); + return object(); +} + +unsigned +nhistogram_depth(const nhistogram& self) { + return self.data().depth(); +} + +uint64_t +nhistogram_getitem(const nhistogram& self, python::object oidx) { + using namespace python; + + if (self.dim() == 1) + return self(extract(oidx)()); + + const unsigned dim = len(oidx); + if (dim != self.dim()) { + PyErr_SetString(PyExc_RuntimeError, "wrong number of arguments"); + throw_error_already_set(); + } + + int idx[BOOST_HISTOGRAM_AXIS_LIMIT]; + for (unsigned i = 0; i < dim; ++i) + idx[i] = extract(oidx[i]); + + return self(idx); +} + +void register_nhistogram() +{ + using namespace python; + + class_< + nhistogram, bases, + shared_ptr + >("nhistogram", no_init) + .def("__init__", raw_function(nhistogram_init)) + // shadowed C++ ctors + .def(init()) + .add_property("__array_interface__", nhistogram_array_interface) + .def("fill", raw_function(nhistogram_fill)) + .add_property("depth", nhistogram_depth) + .add_property("sum", &nhistogram::sum) + .def("__getitem__", nhistogram_getitem) + .def(self == self) + .def(self += self) + .def(self + self) + .def_pickle(serialization_suite()) + ; +} + +} +} diff --git a/src/python/serialization_suite.hpp b/src/python/serialization_suite.hpp new file mode 100644 index 00000000..97dcb298 --- /dev/null +++ b/src/python/serialization_suite.hpp @@ -0,0 +1,93 @@ +#ifndef _BOOST_PYTHON_SERIALIZATION_SUITE_HPP_ +#define _BOOST_PYTHON_SERIALIZATION_SUITE_HPP_ + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace boost { + +class str_sink : public iostreams::sink { +public: + str_sink(PyObject** pstr) : + pstr_(pstr), + len_(0), + pos_(0) + { assert(*pstr == 0); } + + std::streamsize write(const char* s, std::streamsize n) + { + if (len_ == 0) { + *pstr_ = PyString_FromStringAndSize(s, n); + if (*pstr_ == 0) { + PyErr_SetString(PyExc_RuntimeError, "cannot allocate memory"); + python::throw_error_already_set(); + } + } else { + if (pos_ + n > len_) { + len_ = pos_ + n; + if (_PyString_Resize(pstr_, len_) == -1) + python::throw_error_already_set(); + } + char* b = PyString_AS_STRING(*pstr_); + std::copy(s, s + n, b + pos_); + } + pos_ += n; + return n; + } + +private: + PyObject** pstr_; + std::streamsize len_; + std::streamsize pos_; +}; + +template +struct serialization_suite : python::pickle_suite +{ + static + python::tuple getstate(python::object obj) + { + PyObject* pobj = 0; + iostreams::stream os(&pobj); + archive::text_oarchive oa(os); + oa << python::extract(obj)(); + os.flush(); + return python::make_tuple(obj.attr("__dict__"), + python::object(python::handle<>(pobj))); + } + + static + void setstate(python::object obj, python::tuple state) + { + if (python::len(state) != 2) { + PyErr_SetObject(PyExc_ValueError, + ("expected 2-item tuple in call to __setstate__; got %s" + % state).ptr()); + python::throw_error_already_set(); + } + + // restore the object's __dict__ + python::dict d = python::extract(obj.attr("__dict__")); + d.update(state[0]); + + // restore the C++ object + python::object o = state[1]; + iostreams::stream + is(python::extract(o)(), python::len(o)); + archive::text_iarchive ia(is); + ia >> python::extract(obj)(); + } + + static + bool getstate_manages_dict() { return true; } +}; + +} + +#endif diff --git a/src/zero_suppression.cpp b/src/zero_suppression.cpp new file mode 100644 index 00000000..0ac233b1 --- /dev/null +++ b/src/zero_suppression.cpp @@ -0,0 +1,63 @@ +#include + +namespace boost { +namespace histogram { +namespace detail { + +bool +zero_suppression_encode(std::vector& output, std::size_t output_max, + const char* input, std::size_t input_len) +{ + #define FILL_OUTPUT { \ + if ((output_max - output.size()) < 2) \ + return false; \ + output.push_back(0); \ + output.push_back(nzero); \ + nzero = 0; \ + } + + const char* input_end = input + input_len; + unsigned nzero = 0; + for (; input != input_end; ++input) { + if (*input != 0) { + if (nzero) + FILL_OUTPUT + if (output.size() == output_max) + return false; + output.push_back(*input); + } + else { + ++nzero; + if (nzero == 256) + FILL_OUTPUT + } + } + if (nzero) + FILL_OUTPUT + return true; +} + +void +zero_suppression_decode(char* output, std::size_t output_len, + const std::vector& input) +{ + const char* inp = &input[0]; + const char* output_end = output + output_len; + while (output != output_end) { + *output = *inp; + if (*inp == 0) { + const char nzero = *(++inp); + for (char j = 1; j != nzero; ++j) { + *(++output) = 0; + if (output == output_end) + return; + } + } + ++inp; + ++output; + } +} + +} +} +} diff --git a/test/python_suite_test.py b/test/python_suite_test.py new file mode 100755 index 00000000..b7e5ead6 --- /dev/null +++ b/test/python_suite_test.py @@ -0,0 +1,504 @@ +#!/usr/bin/env python +import unittest +from math import pi +from histogram import nhistogram, regular_axis, polar_axis, \ + variable_axis, category_axis, integer_axis +import cPickle +import StringIO +try: import numpy # optional, for tests that need numpy +except ImportError: pass + +class test_regular_axis(unittest.TestCase): + + def test_init(self): + regular_axis(1, 1.0, 2.0) + regular_axis(1, 1.0, 2.0, label="ra") + regular_axis(1, 1.0, 2.0, uoflow=False) + regular_axis(1, 1.0, 2.0, label="ra", uoflow=False) + with self.assertRaises(TypeError): + regular_axis() + with self.assertRaises(TypeError): + regular_axis(1) + with self.assertRaises(TypeError): + regular_axis(1, 1.0) + with self.assertRaises(RuntimeError): + regular_axis(0, 1.0, 2.0) + with self.assertRaises(TypeError): + regular_axis("1", 1.0, 2.0) + with self.assertRaises(StandardError): + regular_axis(-1, 1.0, 2.0) + with self.assertRaises(RuntimeError): + regular_axis(1, 2.0, 1.0) + with self.assertRaises(TypeError): + regular_axis(1, 1.0, 2.0, label=0) + with self.assertRaises(TypeError): + regular_axis(1, 1.0, 2.0, label="ra", uoflow="True") + with self.assertRaises(TypeError): + regular_axis(1, 1.0, 2.0, bad_keyword="ra") + a = regular_axis(4, 1.0, 2.0) + self.assertEqual(a, regular_axis(4, 1.0, 2.0)) + self.assertNotEqual(a, regular_axis(3, 1.0, 2.0)) + self.assertNotEqual(a, regular_axis(4, 1.1, 2.0)) + self.assertNotEqual(a, regular_axis(4, 1.0, 2.1)) + + def test_len(self): + a = regular_axis(4, 1.0, 2.0) + self.assertEqual(len(a), 5) + + def test_repr(self): + for s in ("regular_axis(4, 1.1, 2.2)", + "regular_axis(4, 1.1, 2.2, label='ra')", + "regular_axis(4, 1.1, 2.2, uoflow=False)", + "regular_axis(4, 1.1, 2.2, label='ra', uoflow=False)"): + self.assertEqual(str(eval(s)), s) + + def test_getitem(self): + v = [1.0, 1.25, 1.5, 1.75, 2.0] + a = regular_axis(4, 1.0, 2.0) + for i in xrange(5): + self.assertEqual(a[i], v[i]) + + def test_iter(self): + v = [1.0, 1.25, 1.5, 1.75, 2.0] + a = regular_axis(4, 1.0, 2.0) + self.assertEqual([x for x in a], v) + + def test_index(self): + a = regular_axis(4, 1.0, 2.0) + self.assertEqual(a.index(-1), -1) + self.assertEqual(a.index(0.99), -1) + self.assertEqual(a.index(1.0), 0) + self.assertEqual(a.index(1.1), 0) + self.assertEqual(a.index(1.249), 0) + self.assertEqual(a.index(1.250), 1) + self.assertEqual(a.index(1.499), 1) + self.assertEqual(a.index(1.500), 2) + self.assertEqual(a.index(1.749), 2) + self.assertEqual(a.index(1.750), 3) + self.assertEqual(a.index(1.999), 3) + self.assertEqual(a.index(2.000), 4) + self.assertEqual(a.index(2.1), 4) + self.assertEqual(a.index(20), 4) + +class test_polar_axis(unittest.TestCase): + + def test_init(self): + polar_axis(1) + polar_axis(4, 1.0) + polar_axis(4, 1.0, label="pa") + with self.assertRaises(TypeError): + polar_axis() + with self.assertRaises(StandardError): + polar_axis(-1) + with self.assertRaises(TypeError): + polar_axis(4, 1.0, uoflow=True) + with self.assertRaises(TypeError): + polar_axis(1, 1.0, 2.0) + with self.assertRaises(TypeError): + polar_axis(1, 1.0, label=1) + with self.assertRaises(TypeError): + polar_axis("1") + a = polar_axis(4, 1.0) + self.assertEqual(a, polar_axis(4, 1.0)) + self.assertNotEqual(a, polar_axis(2, 1.0)) + self.assertNotEqual(a, polar_axis(4, 0.0)) + + def test_len(self): + self.assertEqual(len(polar_axis(4)), 5) + self.assertEqual(len(polar_axis(4, 1.0)), 5) + + def test_repr(self): + for s in ("polar_axis(4)", + "polar_axis(4, 1)", + "polar_axis(4, 1, label='x')", + "polar_axis(4, label='x')"): + self.assertEqual(str(eval(s)), s) + + def test_getitem(self): + v = [1.0, 1.0 + 0.5 * pi, 1.0 + pi, 1.0 + 1.5 *pi, 1.0 + 2.0 * pi] + a = polar_axis(4, 1.0) + for i in xrange(5): + self.assertEqual(a[i], v[i]) + + def test_iter(self): + a = polar_axis(4, 1.0) + v = [1.0, 1.0 + 0.5 * pi, 1.0 + pi, 1.0 + 1.5 *pi, 1.0 + 2.0 * pi] + self.assertEqual([x for x in a], v) + + def test_index(self): + a = polar_axis(4, 1.0) + d = 0.5 * pi + self.assertEqual(a.index(0.99 - 4*d), 3) + self.assertEqual(a.index(0.99 - 3*d), 0) + self.assertEqual(a.index(0.99 - 2*d), 1) + self.assertEqual(a.index(0.99 - d), 2) + self.assertEqual(a.index(0.99), 3) + self.assertEqual(a.index(1.0), 0) + self.assertEqual(a.index(1.01), 0) + self.assertEqual(a.index(0.99 + d), 0) + self.assertEqual(a.index(1.0 + d), 1) + self.assertEqual(a.index(1.0 + 2*d), 2) + self.assertEqual(a.index(1.0 + 3*d), 3) + self.assertEqual(a.index(1.0 + 4*d), 0) + self.assertEqual(a.index(1.0 + 5*d), 1) + +class test_variable_axis(unittest.TestCase): + + def test_init(self): + variable_axis(0, 1) + variable_axis(1, -1) + variable_axis(0, 1, 2, 3, 4) + variable_axis(0, 1, label="va") + variable_axis(0, 1, uoflow=True) + variable_axis(0, 1, label="va", uoflow=True) + with self.assertRaises(TypeError): + variable_axis() + with self.assertRaises(RuntimeError): + variable_axis(1.0) + with self.assertRaises(TypeError): + variable_axis("1", 2) + with self.assertRaises(TypeError): + regular_axis(1, 1.0, 2.0, bad_keyword="ra") + a = variable_axis(-0.1, 0.2, 0.3) + self.assertEqual(a, variable_axis(-0.1, 0.2, 0.3)) + self.assertNotEqual(a, variable_axis(0, 0.2, 0.3)) + self.assertNotEqual(a, variable_axis(-0.1, 0.1, 0.3)) + self.assertNotEqual(a, variable_axis(-0.1, 0.1)) + + def test_len(self): + self.assertEqual(len(variable_axis(-0.1, 0.2, 0.3)), 3) + + def test_repr(self): + for s in ("variable_axis(-0.1, 0.2)", + "variable_axis(-0.1, 0.2, 0.3)", + "variable_axis(-0.1, 0.2, 0.3, label='va')", + "variable_axis(-0.1, 0.2, 0.3, uoflow=False)", + "variable_axis(-0.1, 0.2, 0.3, label='va', uoflow=False)"): + self.assertEqual(str(eval(s)), s) + + def test_getitem(self): + v = [-0.1, 0.2, 0.3] + a = variable_axis(*v) + for i in xrange(3): + self.assertEqual(a[i], v[i]) + + def test_iter(self): + v = [-0.1, 0.2, 0.3] + a = variable_axis(*v) + self.assertEqual([x for x in a], v) + + def test_index(self): + a = variable_axis(-0.1, 0.2, 0.3) + self.assertEqual(a.index(-10.0), -1) + self.assertEqual(a.index(-0.11), -1) + self.assertEqual(a.index(-0.1), 0) + self.assertEqual(a.index(0.0), 0) + self.assertEqual(a.index(0.19), 0) + self.assertEqual(a.index(0.2), 1) + self.assertEqual(a.index(0.21), 1) + self.assertEqual(a.index(0.29), 1) + self.assertEqual(a.index(0.3), 2) + self.assertEqual(a.index(0.31), 2) + self.assertEqual(a.index(10), 2) + +class test_category_axis(unittest.TestCase): + + def test_init(self): + category_axis("A", "B", "C") + category_axis("A;B;C") + with self.assertRaises(TypeError): + category_axis() + with self.assertRaises(TypeError): + category_axis(1) + with self.assertRaises(TypeError): + category_axis("1", 2) + with self.assertRaises(TypeError): + category_axis("A", "B", "C", label="ca") + with self.assertRaises(TypeError): + category_axis("A", "B", "C", uoflow=True) + self.assertEqual(category_axis("A", "B", "C"), + category_axis("A", "B", "C")) + self.assertEqual(category_axis("A", "B", "C"), + category_axis("A;B;C")) + self.assertEqual(category_axis(";"), category_axis("", "")) + self.assertEqual(category_axis("abc;d"), category_axis("abc", "d")) + self.assertEqual(category_axis("a;bcd"), category_axis("a", "bcd")) + self.assertNotEqual(category_axis("A;B"), category_axis("A;C")) + + def test_len(self): + a = category_axis("A", "B", "C") + self.assertEqual(len(a), 3) + a = category_axis("A;B;C") + self.assertEqual(len(a), 3) + + def test_repr(self): + for s in ("category_axis('A')", + "category_axis('A', 'B')", + "category_axis('A', 'B', 'C')"): + self.assertEqual(str(eval(s)), s) + + def test_getitem(self): + c = "A", "B", "C" + a = category_axis(*c) + for i in xrange(3): + self.assertEqual(a[i], c[i]) + + def test_iter(self): + c = ["A", "B", "C"] + self.assertEqual([x for x in category_axis(*c)], c) + +class test_integer_axis(unittest.TestCase): + + def test_init(self): + integer_axis(-1, 2) + with self.assertRaises(TypeError): + integer_axis() + with self.assertRaises(TypeError): + integer_axis(1) + with self.assertRaises(TypeError): + integer_axis("1", 2) + with self.assertRaises(RuntimeError): + integer_axis(2, -1) + with self.assertRaises(TypeError): + integer_axis(1, 2, 3) + self.assertEqual(integer_axis(-1, 2), integer_axis(-1, 2)) + self.assertNotEqual(integer_axis(-1, 2), integer_axis(-1, 2, label="ia")) + self.assertNotEqual(integer_axis(-1, 2), integer_axis(-1, 2, uoflow=True)) + + def test_len(self): + self.assertEqual(len(integer_axis(-1, 2)), 4) + + def test_repr(self): + for s in ("integer_axis(-1, 1)", + "integer_axis(-1, 1, label='ia')", + "integer_axis(-1, 1, uoflow=False)", + "integer_axis(-1, 1, label='ia', uoflow=False)"): + self.assertEqual(str(eval(s)), s) + + def test_label(self): + self.assertEqual(integer_axis(-1, 2, label="ia").label, "ia") + + def test_getitem(self): + v = [-1, 0, 1, 2] + a = integer_axis(-1, 2) + for i in xrange(4): + self.assertEqual(a[i], v[i]) + + def test_iter(self): + v = [x for x in integer_axis(-1, 2)] + self.assertEqual(v, [-1, 0, 1, 2]) + + def test_index(self): + a = integer_axis(-1, 2) + self.assertEqual(a.index(-3), -1) + self.assertEqual(a.index(-2), -1) + self.assertEqual(a.index(-1), 0) + self.assertEqual(a.index(0), 1) + self.assertEqual(a.index(1), 2) + self.assertEqual(a.index(2), 3) + self.assertEqual(a.index(3), 4) + self.assertEqual(a.index(4), 4) + +class nhistogram_test(unittest.TestCase): + + def test_init(self): + nhistogram(integer_axis(-1, 1)) + with self.assertRaises(TypeError): + nhistogram(1) + with self.assertRaises(TypeError): + nhistogram("bla") + with self.assertRaises(TypeError): + nhistogram([]) + with self.assertRaises(TypeError): + nhistogram(regular_axis) + with self.assertRaises(TypeError): + nhistogram(regular_axis()) + with self.assertRaises(TypeError): + nhistogram([integer_axis(-1, 1)]) + with self.assertRaises(TypeError): + nhistogram(integer_axis(-1, 1), unknown_keyword="nh") + + h = nhistogram(integer_axis(-1, 1)) + self.assertEqual(h.dim, 1) + self.assertEqual(h.axis(0), integer_axis(-1, 1)) + self.assertEqual(h.shape(0), 5) + self.assertEqual(h.depth, 1) + self.assertEqual(nhistogram(integer_axis(-1, 1, uoflow=False)).shape(0), 3) + self.assertNotEqual(h, nhistogram(regular_axis(1, -1, 1))) + self.assertNotEqual(h, nhistogram(integer_axis(-1, 2))) + self.assertNotEqual(h, nhistogram(integer_axis(-1, 1, label="ia"))) + self.assertNotEqual(h, nhistogram(integer_axis(-1, 1, uoflow=True))) + + def test_fill_1d(self): + h0 = nhistogram(integer_axis(-1, 1, uoflow=False)) + h1 = nhistogram(integer_axis(-1, 1, uoflow=True)) + for h in (h0, h1): + with self.assertRaises(StandardError): + h.fill() + with self.assertRaises(StandardError): + h.fill(1, 2) + h.fill(-10) + h.fill(-1) + h.fill(-1) + h.fill(0) + h.fill(1) + h.fill(1) + h.fill(1) + h.fill(10) + self.assertEqual(h0.sum, 6) + self.assertEqual(h0.depth, 1) + self.assertEqual(h0.shape(0), 3) + self.assertEqual(h1.sum, 8) + self.assertEqual(h1.depth, 1) + self.assertEqual(h1.shape(0), 5) + + for h in (h0, h1): + self.assertEqual(h[0], 2) + self.assertEqual(h[1], 1) + self.assertEqual(h[2], 3) + with self.assertRaises(TypeError): + h[0, 1] + + self.assertEqual(h1[-1], 1) + self.assertEqual(h1[3], 1) + + def test_growth(self): + h = nhistogram(integer_axis(-1, 1)) + h.fill(-1) + h.fill(1) + h.fill(1) + for i in xrange(255): + h.fill(0) + self.assertEqual(h.depth, 1) + h.fill(0) + self.assertEqual(h.depth, 2) + for i in xrange(1000-256): + h.fill(0) + self.assertEqual(h[-1], 0) + self.assertEqual(h[0], 1) + self.assertEqual(h[1], 1000) + self.assertEqual(h[2], 2) + self.assertEqual(h[3], 0) + + def test_fill_2d(self): + for uoflow in (False, True): + h = nhistogram(integer_axis(-1, 1, uoflow=uoflow), + regular_axis(4, -2, 2, uoflow=uoflow)) + h.fill(-1, -2) + h.fill(-1, -1) + h.fill(0, 0) + h.fill(0, 1) + h.fill(1, 0) + h.fill(3, -1) + h.fill(0, -3) + with self.assertRaises(StandardError): + h.fill(1) + with self.assertRaises(StandardError): + h.fill(1, 2, 3) + + m = [[1, 1, 0, 0, 0, 0], + [0, 0, 1, 1, 0, 1], + [0, 0, 1, 0, 0, 0], + [0, 1, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0]] + for i in xrange(h.axis(0).bins + 2*uoflow): + for j in xrange(h.axis(1).bins + 2*uoflow): + self.assertEqual(h[i, j], m[i][j]) + + def test_add_2d(self): + for uoflow in (False, True): + h = nhistogram(integer_axis(-1, 1, uoflow=uoflow), + regular_axis(4, -2, 2, uoflow=uoflow)) + h.fill(-1, -2) + h.fill(-1, -1) + h.fill(0, 0) + h.fill(0, 1) + h.fill(1, 0) + h.fill(3, -1) + h.fill(0, -3) + + m = [[1, 1, 0, 0, 0, 0], + [0, 0, 1, 1, 0, 1], + [0, 0, 1, 0, 0, 0], + [0, 1, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0]] + + h2 = h + h + h += h + self.assertEqual(h, h2) + + for i in xrange(h.axis(0).bins + 2*uoflow): + for j in xrange(h.axis(1).bins + 2*uoflow): + self.assertEqual(h[i, j], 2 * m[i][j]) + + def test_pickle(self): + a = nhistogram(category_axis('A', 'B', 'C'), + integer_axis(0, 3, label='ia'), + regular_axis(4, 0.0, 4.0, uoflow=False), + variable_axis(0.0, 1.0, 2.0)) + for i in xrange(a.axis(0).bins): + a.fill(i, 0, 0, 0) + for j in xrange(a.axis(1).bins): + a.fill(i, j, 0, 0) + for k in xrange(a.axis(2).bins): + a.fill(i, j, k, 0) + for l in xrange(a.axis(3).bins): + a.fill(i, j, k, l) + + io = StringIO.StringIO() + cPickle.dump(a, io) + io.seek(0) + b = cPickle.load(io) + self.assertNotEqual(id(a), id(b)) + self.assertEqual(a.dim, b.dim) + self.assertEqual(a.axis(0), b.axis(0)) + self.assertEqual(a.axis(1), b.axis(1)) + self.assertEqual(a.axis(2), b.axis(2)) + self.assertEqual(a.axis(3), b.axis(3)) + self.assertEqual(a.sum, b.sum) + self.assertEqual(a, b) + + @unittest.skipUnless("numpy" in globals(), "requires numpy") + def test_numpy_conversion(self): + a = nhistogram(integer_axis(0, 2, uoflow=False)) + for i in xrange(100): + a.fill(1) + b = numpy.array(a) # a copy + v = numpy.asarray(a) # a view + self.assertEqual(b.dtype, numpy.uint8) + self.assertTrue(numpy.all(b == numpy.array((0, 100, 0)))) + for i in xrange(100): + a.fill(1) + self.assertTrue(numpy.all(b == numpy.array((0, 100, 0)))) + self.assertTrue(numpy.all(v == numpy.array((0, 200, 0)))) + for i in xrange(100): + a.fill(1) + b = numpy.array(a) + self.assertEqual(b.dtype, numpy.uint16) + self.assertTrue(numpy.all(b == numpy.array((0, 300, 0)))) + # view does not follow underlying switch in word size + self.assertFalse(numpy.all(v == b)) + + @unittest.skipUnless("numpy" in globals(), "requires numpy") + def test_fill_with_numpy_array(self): + a = nhistogram(integer_axis(0, 2, uoflow=False)) + a.fill(numpy.array([-1, 0, 1, 2, 1, 4])) + self.assertEqual(a[0], 1) + self.assertEqual(a[1], 2) + self.assertEqual(a[2], 1) + a = nhistogram(integer_axis(0, 1, uoflow=False), + regular_axis(2, 0, 2, uoflow=False)) + a.fill(numpy.array([[-1., -1.], [0., 1.], [1., 0.1]])) + self.assertEqual(a[0, 0], 0) + self.assertEqual(a[0, 1], 1) + self.assertEqual(a[1, 0], 1) + self.assertEqual(a[1, 1], 0) + a = nhistogram(integer_axis(0, 2, uoflow=False)) + a.fill([0, 0, 1, 2]) + a.fill((1, 0, 2, 2)) + self.assertEqual(a[0], 3) + self.assertEqual(a[1], 2) + self.assertEqual(a[2], 3) + +if __name__ == "__main__": + unittest.main() diff --git a/test/zero_suppression_test.cpp b/test/zero_suppression_test.cpp new file mode 100644 index 00000000..acd4425d --- /dev/null +++ b/test/zero_suppression_test.cpp @@ -0,0 +1,154 @@ +#include +#define BOOST_TEST_MODULE zero_suppression_test +#include +#include +#include +#include +using namespace boost::assign; +using boost::histogram::detail::zero_suppression_encode; +using boost::histogram::detail::zero_suppression_decode; +namespace tt = boost::test_tools; + +void +print(const std::vector& v) { + std::cerr << "["; + for (unsigned i = 0; i < v.size(); ++i) { + std::cerr << unsigned(v[i]) << (i < v.size() - 1 ? ", " : "]"); + } + std::cerr << "\n"; +} + +#define EQUAL_VECTORS(a, b) \ + BOOST_CHECK_EQUAL_COLLECTIONS(a.begin(), a.end(), b.begin(), b.end()) + +BOOST_AUTO_TEST_CASE( codec_no_zero ) +{ + std::vector a, b, c(3, 0); + a += 1, 1, 1; + BOOST_CHECK(zero_suppression_encode(b, 3, &a[0], a.size())); + EQUAL_VECTORS(a, b); + zero_suppression_decode(&c[0], c.size(), b); + EQUAL_VECTORS(a, c); +} + +BOOST_AUTO_TEST_CASE( codec_fail ) +{ + std::vector a, b; + a += 1, 1, 1; + BOOST_CHECK(zero_suppression_encode(b, 2, &a[0], a.size()) == false); +} + +BOOST_AUTO_TEST_CASE( codec_empty ) +{ + std::vector a, b, c; + BOOST_CHECK(zero_suppression_encode(b, 3, &a[0], a.size())); + EQUAL_VECTORS(a, b); + zero_suppression_decode(&c[0], c.size(), b); + EQUAL_VECTORS(a, c); +} + +BOOST_AUTO_TEST_CASE( codec_zero_0 ) +{ + std::vector a, b, c, d(2, 0); + a += 0, 1; + c += 0, 1, 1; + BOOST_CHECK(zero_suppression_encode(b, c.size(), &a[0], a.size())); + EQUAL_VECTORS(b, c); + zero_suppression_decode(&d[0], d.size(), b); + EQUAL_VECTORS(a, d); +} + +BOOST_AUTO_TEST_CASE( codec_zero_1 ) +{ + std::vector a, b, c, d(3, 0); + a += 0, 0, 1; + c += 0, 2, 1; + BOOST_CHECK(zero_suppression_encode(b, c.size(), &a[0], a.size())); + EQUAL_VECTORS(b, c); + zero_suppression_decode(&d[0], d.size(), b); + EQUAL_VECTORS(a, d); +} + +BOOST_AUTO_TEST_CASE( codec_zero_2 ) +{ + std::vector a, b, c, d(2, 0); + a += 0, 0; + c += 0, 2; + BOOST_CHECK(zero_suppression_encode(b, c.size(), &a[0], a.size())); + EQUAL_VECTORS(b, c); + zero_suppression_decode(&d[0], d.size(), b); + EQUAL_VECTORS(a, d); +} + +BOOST_AUTO_TEST_CASE( codec_zero_3 ) +{ + std::vector a, b, c, d(1, 0); + a += 0; + c += 0, 1; + BOOST_CHECK(zero_suppression_encode(b, c.size(), &a[0], a.size())); + EQUAL_VECTORS(b, c); + zero_suppression_decode(&d[0], d.size(), b); + EQUAL_VECTORS(a, d); +} + +BOOST_AUTO_TEST_CASE( codec_zero_4 ) +{ + std::vector a, b, c, d(3, 0); + a += 0, 1, 0; + c += 0, 1, 1, 0, 1; + BOOST_CHECK(zero_suppression_encode(b, c.size(), &a[0], a.size())); + EQUAL_VECTORS(b, c); + zero_suppression_decode(&d[0], d.size(), b); + EQUAL_VECTORS(a, d); +} + +BOOST_AUTO_TEST_CASE( codec_zero_5 ) +{ + std::vector a, b, c, d(4, 0); + a += 0, 1, 0, 0; + c += 0, 1, 1, 0, 2; + BOOST_CHECK(zero_suppression_encode(b, c.size(), &a[0], a.size())); + EQUAL_VECTORS(b, c); + zero_suppression_decode(&d[0], d.size(), b); + EQUAL_VECTORS(a, d); +} + +BOOST_AUTO_TEST_CASE( codec_zero_6 ) +{ + std::vector a(255, 0), b, c, d(255, 0); + c += 0, 255; + BOOST_CHECK(zero_suppression_encode(b, c.size(), &a[0], a.size())); + EQUAL_VECTORS(b, c); + zero_suppression_decode(&d[0], d.size(), b); + EQUAL_VECTORS(a, d); +} + +BOOST_AUTO_TEST_CASE( codec_zero_7 ) +{ + std::vector a(256, 0), b, c, d(256, 0); + c += 0, 0; + BOOST_CHECK(zero_suppression_encode(b, c.size(), &a[0], a.size())); + EQUAL_VECTORS(b, c); + zero_suppression_decode(&d[0], d.size(), b); + EQUAL_VECTORS(a, d); +} + +BOOST_AUTO_TEST_CASE( codec_zero_8 ) +{ + std::vector a(257, 0), b, c, d(257, 0); + c += 0, 0, 0, 1; + BOOST_CHECK(zero_suppression_encode(b, c.size(), &a[0], a.size())); + EQUAL_VECTORS(b, c); + zero_suppression_decode(&d[0], d.size(), b); + EQUAL_VECTORS(a, d); +} + +BOOST_AUTO_TEST_CASE( codec_zero_9 ) +{ + std::vector a(1000, 0), b, c, d(1000, 0); + c += 0, 0, 0, 0, 0, 0, 0, 232; + BOOST_CHECK(zero_suppression_encode(b, c.size(), &a[0], a.size())); + EQUAL_VECTORS(b, c); + zero_suppression_decode(&d[0], d.size(), b); + EQUAL_VECTORS(a, d); +}