init from private repository

This commit is contained in:
Hans Dembinski 2016-04-06 09:23:37 -04:00
parent e7622c74a3
commit 19a740fa68
27 changed files with 3174 additions and 0 deletions

76
CMake/CMakeLists.txt Normal file
View File

@ -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()

24
CMake/FindNumpy.cmake Normal file
View File

@ -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()

159
CMake/FindROOT.cmake Normal file
View File

@ -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()

2
CMakeLists.txt Normal file
View File

@ -0,0 +1,2 @@
cmake_minimum_required (VERSION 2.8)
include(CMake/CMakeLists.txt)

41
examples/merge.py Executable file
View File

@ -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()

30
examples/sizeof_axis.cpp Normal file
View File

@ -0,0 +1,30 @@
#include <cstdio>
#include <boost/preprocessor.hpp>
#include <boost/histogram/axis.hpp>
#include <boost/container/vector.hpp>
#include <boost/scoped_array.hpp>
#include <vector>
#include <string>
#include <valarray>
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<double>);
SIZEOF(std::valarray<double>);
SIZEOF(boost::container::vector<double>);
SIZEOF(boost::scoped_array<double>);
SIZEOF(regular_axis);
SIZEOF(polar_axis);
SIZEOF(variable_axis);
SIZEOF(category_axis);
SIZEOF(integer_axis);
SIZEOF(axis_type);
}

146
examples/speed_vs_root.cpp Normal file
View File

@ -0,0 +1,146 @@
#include <boost/histogram/nhistogram.hpp>
#include <boost/histogram/axis.hpp>
#include <boost/random.hpp>
#include <TH1I.h>
#include <TH3I.h>
#include <THn.h>
#include <algorithm>
#include <limits>
#include <vector>
#include <ctime>
#include <cstdio>
using namespace std;
using namespace boost::histogram;
template <typename D>
struct rng {
boost::random::mt19937 r;
D d;
rng(double a, double b) : d(a, b) {}
double operator()() { return d(r); }
};
vector<double> random_array(unsigned n, int type) {
using namespace boost::random;
std::vector<double> result;
switch (type) {
case 0:
std::generate_n(std::back_inserter(result), n, rng<uniform_real_distribution<> >(0.0, 1.0));
break;
case 1:
std::generate_n(std::back_inserter(result), n, rng<normal_distribution<> >(0.0, 0.3));
break;
}
return result;
}
void compare_1d(unsigned n)
{
vector<double> r = random_array(1000000, 1);
double best_root = std::numeric_limits<double>::max();
double best_boost = std::numeric_limits<double>::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<double> r = random_array(3 * n, 1);
double best_root = std::numeric_limits<double>::max();
double best_boost = std::numeric_limits<double>::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<double> r = random_array(6 * n, 1);
double best_root = std::numeric_limits<double>::max();
double best_boost = std::numeric_limits<double>::max();
for (unsigned k = 0; k < 3; ++k) {
double x[6];
vector<int> bin(6, 10);
vector<double> min(6, 0);
vector<double> 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);
}

View File

@ -0,0 +1,7 @@
#ifndef _HISTOGRAM_I3HISTOGRAM_H_
#define _HISTOGRAM_I3HISTOGRAM_H_
#include <boost/histogram/nhistogram.hpp>
#include <boost/histogram/whistogram.hpp>
#endif

View File

@ -0,0 +1,252 @@
#ifndef _BOOST_HISTOGRAM_AXIS_HPP_
#define _BOOST_HISTOGRAM_AXIS_HPP_
#include <boost/algorithm/clamp.hpp>
#include <boost/variant.hpp>
#include <boost/scoped_array.hpp>
#include <boost/serialization/access.hpp>
#include <boost/serialization/base_object.hpp>
#include <boost/serialization/array.hpp>
#include <string>
#include <vector>
#include <cmath>
#include <ostream>
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 <class Archive>
void serialize(Archive& ar, unsigned version)
{
using namespace serialization;
ar & size_;
ar & label_;
}
};
template <typename Derived>
class real_axis {
public:
typedef double value_type;
double left(int idx) const {
return static_cast<const Derived&>(*this)[idx];
}
double right(int idx) const {
return static_cast<const Derived&>(*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<regular_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 <class Archive>
void serialize(Archive& ar, unsigned version)
{
using namespace serialization;
ar & boost::serialization::base_object<axis_base>(*this);
ar & min_;
ar & range_;
}
};
// real polar axis (constant bin widths, wraps around)
class polar_axis: public axis_base, public real_axis<polar_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 <class Archive>
void serialize(Archive& ar, unsigned version)
{
using namespace serialization;
ar & boost::serialization::base_object<axis_base>(*this);
ar & start_;
}
};
// real variable axis (varying bin widths)
class variable_axis : public axis_base, public real_axis<variable_axis> {
public:
template <typename Container>
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 <typename Iterator>
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<double> x_;
friend class serialization::access;
template <class Archive>
void serialize(Archive& ar, unsigned version)
{
ar & boost::serialization::base_object<axis_base>(*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<std::string>&);
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<std::string> categories_;
friend class serialization::access;
template <class Archive>
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 <class Archive>
void serialize(Archive& ar, unsigned version)
{
using namespace serialization;
ar & boost::serialization::base_object<axis_base>(*this);
ar & min_;
}
};
typedef variant<
regular_axis, // most common type
polar_axis,
variable_axis,
category_axis,
integer_axis
> axis_type;
typedef std::vector<axis_type> axes_type;
std::ostream& operator<<(std::ostream&, const axis_type&);
}
}
#endif

View File

@ -0,0 +1,23 @@
#ifndef _BOOST_HISTOGRAM_DETAIL_ZERO_SUPPRESSION_HPP_
#define _BOOST_HISTOGRAM_DETAIL_ZERO_SUPPRESSION_HPP_
#include <cstddef>
#include <vector>
namespace boost {
namespace histogram {
namespace detail {
bool
zero_suppression_encode(std::vector<char>& 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<char>& input);
}
}
}
#endif

View File

@ -0,0 +1,121 @@
#ifndef _BOOST_HISTOGRAM_BASE_HPP_
#define _BOOST_HISTOGRAM_BASE_HPP_
#include <boost/histogram/axis.hpp>
#include <boost/histogram/visitors.hpp>
#include <boost/cstdint.hpp>
#include <boost/preprocessor.hpp>
#include <boost/serialization/access.hpp>
#include <boost/serialization/split_member.hpp>
#include <boost/serialization/variant.hpp>
#include <boost/serialization/vector.hpp>
#include <boost/type_traits/is_same.hpp>
#include <bitset>
#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 <typename T>
T& axis(unsigned i)
{ if (is_same<T, axis_type>::value) return axes_[i];
return boost::get<T&>(axes_[i]); }
template <typename T>
const T& axis(unsigned i) const
{ if (is_same<T, axis_type>::value) return axes_[i];
return boost::get<const T&>(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 <typename Array>
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<BOOST_HISTOGRAM_AXIS_LIMIT> uoflow_;
void update_buffers(); ///< fills size_ and uoflow_
friend class serialization::access;
template <class Archive>
void serialize(Archive& ar, unsigned version)
{
using namespace serialization;
ar & axes_;
if (Archive::is_loading::value)
update_buffers();
}
};
}
}
#endif

View File

@ -0,0 +1,111 @@
#ifndef _BOOST_HISTOGRAM_NHISTOGRAM_HPP_
#define _BOOST_HISTOGRAM_NHISTOGRAM_HPP_
#include <boost/histogram/axis.hpp>
#include <boost/histogram/histogram_base.hpp>
#include <boost/histogram/nstore.hpp>
#include <boost/preprocessor.hpp>
#include <boost/serialization/access.hpp>
#include <boost/serialization/base_object.hpp>
#include <boost/shared_ptr.hpp>
#include <ostream>
#include <vector>
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 <typename Array>
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 <typename Array>
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 <class Archive>
void serialize(Archive& ar, unsigned version)
{
using namespace serialization;
ar & boost::serialization::base_object<histogram_base>(*this);
ar & data_;
}
};
nhistogram operator+(const nhistogram& a, const nhistogram& b) {
nhistogram result(a);
return result += b;
}
}
}
#endif

View File

@ -0,0 +1,112 @@
#ifndef _BOOST_HISTOGRAM_NSTORE_HPP_
#define _BOOST_HISTOGRAM_NSTORE_HPP_
#include <boost/serialization/access.hpp>
#include <boost/serialization/array.hpp>
#include <boost/cstdint.hpp>
#include <boost/histogram/detail/zero_suppression.hpp>
#include <cstdlib>
#include <cstring>
#include <limits>
#include <stdexcept>
#include <vector>
#include <new> // 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<T>::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 <class Archive>
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<char> 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<char> buf;
ar & buf;
detail::zero_suppression_decode((char*)buffer_, size_ * depth_, buf);
} else {
ar & serialization::make_array((char*)buffer_, size_ * depth_);
}
}
}
};
}
}
#endif

View File

@ -0,0 +1,52 @@
#ifndef _BOOST_HISTOGRAM_VECTOR_OPERATORS_HPP_
#define _BOOST_HISTOGRAM_VECTOR_OPERATORS_HPP_
#include <stdexcept>
namespace boost {
namespace histogram {
// vectors: generic ==
template<typename T>
bool operator==(const std::vector<T>& a, const std::vector<T>& 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<typename T>
std::vector<T>& operator+=(std::vector<T>& a, const std::vector<T>& 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<typename T>
std::vector<T>& operator-=(std::vector<T>& a, const std::vector<T>& 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

View File

@ -0,0 +1,57 @@
#ifndef _BOOST_HISTOGARM_VISITORS_HPP_
#define _BOOST_HISTOGARM_VISITORS_HPP_
#include <boost/histogram/axis.hpp>
#include <boost/variant/static_visitor.hpp>
namespace boost {
namespace histogram {
namespace detail {
struct bins_visitor : public static_visitor<unsigned>
{
template <typename A>
unsigned operator()(const A& a) const { return a.bins(); }
};
struct fields_visitor : public static_visitor<unsigned>
{
unsigned operator()(const category_axis& a) const { return a.bins(); }
template <typename A>
unsigned operator()(const A& a) const { return a.bins() + 2 * a.uoflow(); }
};
struct uoflow_visitor : public static_visitor<bool>
{
bool operator()(const category_axis& a) const { return false; }
template <typename A>
bool operator()(const A& a) const { return a.uoflow(); }
};
struct index_visitor : public static_visitor<int>
{
double x;
index_visitor(double d) : x(d) {}
int operator()(const category_axis& a) const { return int(x + 0.5); }
template <typename A>
int operator()(const A& a) const { return a.index(x); }
};
struct cmp_visitor : public static_visitor<bool>
{
template <typename T>
bool operator()(const T& a, const T& b) const
{ return a == b; }
template <typename T, typename U>
bool operator()(const T&, const U&) const { return false; }
};
}
}
}
#endif

View File

@ -0,0 +1,110 @@
#include <cstdlib>
#include <cstring>
#include <stdexcept>
#include <boost/serialization/array.hpp>
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 <class Archive>
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_; }
};

269
src/axis.cpp Normal file
View File

@ -0,0 +1,269 @@
#include <boost/histogram/axis.hpp>
#include <limits>
#include <sstream>
#include <stdexcept>
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<double>::infinity();
if (idx > bins())
return std::numeric_limits<double>::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<double>::infinity();
if (idx > bins())
return std::numeric_limits<double>::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<std::string>& 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>
{
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;
}
}
}

63
src/histogram_base.cpp Normal file
View File

@ -0,0 +1,63 @@
#include <boost/histogram/histogram_base.hpp>
#include <boost/histogram/axis.hpp>
#include <stdexcept>
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]);
}
}
}
}

18
src/nhistogram.cpp Normal file
View File

@ -0,0 +1,18 @@
#include <boost/histogram/nhistogram.hpp>
#include <ostream>
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;
}
}
}

157
src/nstore.cpp Normal file
View File

@ -0,0 +1,157 @@
#include <boost/histogram/nstore.hpp>
#include <boost/cstdint.hpp>
#include <stdexcept>
#include <new> // 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<T>::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;
}
}
}

View File

@ -0,0 +1,325 @@
#include <boost/histogram/axis.hpp>
#include <boost/histogram/histogram_base.hpp>
#include <boost/python.hpp>
#include <boost/python/raw_function.hpp>
#include <boost/python/def_visitor.hpp>
#include <boost/utility/enable_if.hpp>
#include <boost/type_traits.hpp>
#include <sstream>
#include <string>
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<double> v;
for (int i = 1, n = len(args); i < n; ++i) {
v.push_back(extract<double>(args[i]));
}
std::string label;
bool uoflow = true;
while (len(kwargs) > 0) {
tuple kv = kwargs.popitem();
std::string k = extract<std::string>(kv[0]);
object v = kv[1];
if (k == "label")
label = extract<std::string>(v);
else if (k == "uoflow")
uoflow = extract<bool>(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<std::string> 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<std::string> c;
for (int i = 1, n = len(args); i < n; ++i)
c.push_back(extract<std::string>(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<python::object>
{
template <typename T>
python::object operator()(const T& t) const { return python::object(T(t)); }
};
template <typename T>
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>
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<typename T>
struct has_index_method
{
struct yes { char x[1]; };
struct no { char x[2]; };
template<typename U, int (U::*)(double) const> struct SFINAE {};
template<typename U> static yes test( SFINAE<U, &U::index>* );
template<typename U> static no test( ... );
enum { value = sizeof(test<T>(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 <class T>
struct axis_suite : public python::def_visitor<axis_suite<T> > {
template <typename Class, typename U>
static
typename enable_if_c<has_index_method<U>::value, void>::type
add_axis_index(Class& cl) {
cl.def("index", &U::index);
}
template <typename Class, typename U>
static
typename disable_if_c<has_index_method<U>::value, void>::type
add_axis_index(Class& cl) {}
template <typename Class, typename U>
static
typename enable_if<is_base_of<axis_base, U>, 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<python::copy_const_reference>()),
(void(U::*)(const std::string&)) &U::label);
}
template <typename Class, typename U>
static
typename disable_if<is_base_of<axis_base, U>, void>::type
add_axis_label(Class& cl) {}
template <class Class>
static void
visit(Class& cl)
{
cl.add_property("bins", &T::bins);
add_axis_index<Class, T>(cl);
add_axis_label<Class, T>(cl);
cl.def("__len__", axis_len<T>);
cl.def("__getitem__", axis_getitem<T>);
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<axis_type>(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_<std::vector<double> >("vector_double", no_init);
class_<std::vector<std::string> >("vector_string", no_init);
class_<axes_type>("axes", no_init);
class_<regular_axis>("regular_axis", no_init)
.def(init<unsigned, double, double, std::string, bool>(
(arg("bin"), arg("min"), arg("max"),
arg("label") = std::string(),
arg("uoflow") = true)))
.def(axis_suite<regular_axis>())
;
class_<polar_axis>("polar_axis", no_init)
.def(init<unsigned, double, std::string>(
(arg("bin"), arg("start") = 0.0,
arg("label") = std::string())))
.def(axis_suite<polar_axis>())
;
class_<variable_axis>("variable_axis", no_init)
.def("__init__", raw_function(variable_axis_init))
.def(init<std::vector<double>, std::string, bool>())
.def(axis_suite<variable_axis>())
;
class_<category_axis>("category_axis", no_init)
.def("__init__", raw_function(category_axis_init))
.def(init<std::string>())
.def(init<std::vector<std::string> >())
.def(axis_suite<category_axis>())
;
class_<integer_axis>("integer_axis", no_init)
.def(init<int, int, std::string, bool>(
(arg("min"), arg("max"),
arg("label") = std::string(),
arg("uoflow") = true)))
.def(axis_suite<integer_axis>())
;
class_<histogram_base>("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)
;
}
}
}

22
src/python/module.cpp Normal file
View File

@ -0,0 +1,22 @@
#include <boost/python/module.hpp>
#ifdef USE_NUMPY
# define PY_ARRAY_UNIQUE_SYMBOL boost_histogram_ARRAY_API
# define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
# include <numpy/arrayobject.h>
#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();
}

183
src/python/nhistogram.cpp Normal file
View File

@ -0,0 +1,183 @@
#include <boost/histogram/axis.hpp>
#include <boost/histogram/nhistogram.hpp>
#include <boost/python.hpp>
#include <boost/python/raw_function.hpp>
#include <boost/foreach.hpp>
#include <boost/shared_ptr.hpp>
#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 <numpy/arrayobject.h>
#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<regular_axis> er(pa);
if (er.check()) { axes.push_back(er()); continue; }
extract<polar_axis> ep(pa);
if (ep.check()) { axes.push_back(ep()); continue; }
extract<variable_axis> ev(pa);
if (ev.check()) { axes.push_back(ev()); continue; }
extract<category_axis> ec(pa);
if (ec.check()) { axes.push_back(ec()); continue; }
extract<integer_axis> 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("<u") + str(self.data().depth());
d["data"] = make_tuple((long long)(self.data().buffer()), false);
return d;
}
python::object
nhistogram_fill(python::tuple args, python::dict kwargs) {
using namespace python;
const unsigned nargs = len(args);
nhistogram& self = extract<nhistogram&>(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<double>(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<int>(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<int>(oidx[i]);
return self(idx);
}
void register_nhistogram()
{
using namespace python;
class_<
nhistogram, bases<histogram_base>,
shared_ptr<nhistogram>
>("nhistogram", no_init)
.def("__init__", raw_function(nhistogram_init))
// shadowed C++ ctors
.def(init<const axes_type&>())
.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<nhistogram>())
;
}
}
}

View File

@ -0,0 +1,93 @@
#ifndef _BOOST_PYTHON_SERIALIZATION_SUITE_HPP_
#define _BOOST_PYTHON_SERIALIZATION_SUITE_HPP_
#include <boost/iostreams/concepts.hpp>
#include <boost/iostreams/stream.hpp>
#include <boost/iostreams/device/array.hpp>
#include <boost/archive/text_iarchive.hpp>
#include <boost/archive/text_oarchive.hpp>
#include <boost/python.hpp>
#include <iosfwd>
#include <algorithm>
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<class T>
struct serialization_suite : python::pickle_suite
{
static
python::tuple getstate(python::object obj)
{
PyObject* pobj = 0;
iostreams::stream<str_sink> os(&pobj);
archive::text_oarchive oa(os);
oa << python::extract<const T&>(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<python::dict>(obj.attr("__dict__"));
d.update(state[0]);
// restore the C++ object
python::object o = state[1];
iostreams::stream<iostreams::array_source>
is(python::extract<const char*>(o)(), python::len(o));
archive::text_iarchive ia(is);
ia >> python::extract<T&>(obj)();
}
static
bool getstate_manages_dict() { return true; }
};
}
#endif

63
src/zero_suppression.cpp Normal file
View File

@ -0,0 +1,63 @@
#include <boost/histogram/detail/zero_suppression.hpp>
namespace boost {
namespace histogram {
namespace detail {
bool
zero_suppression_encode(std::vector<char>& 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<char>& 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;
}
}
}
}
}

504
test/python_suite_test.py Executable file
View File

@ -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()

View File

@ -0,0 +1,154 @@
#include <boost/histogram/detail/zero_suppression.hpp>
#define BOOST_TEST_MODULE zero_suppression_test
#include <boost/test/unit_test.hpp>
#include <boost/test/test_tools.hpp>
#include <boost/assign/std/vector.hpp>
#include <iostream>
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<char>& 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<char> 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<char> 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<char> 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<char> 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<char> 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<char> 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<char> 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<char> 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<char> 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<char> 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<char> 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<char> 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<char> 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);
}