Direct memory access to accumulators (#194)

* Allow direct memory access, no longer throws error/segfault

* Internalize accumulators

* First working version of conversion

* Fix several bugs, and change the meaning of np.asarray() slightly

* Fix view not returning an editable view

* Normalize names

* Nicer repr and str for views

* Cleanup views with a decorator

* Some small polish

* Fix #195, unlimited and atomic single access

* Adding tests for storages

* README update

* CHANGELOG and renable a test

* Fix unneeded change

* Address review comments

* Move ostream, remove some extra headers

* Fix reprs, add missing tests

* Adding accum[key] access

* Addressing second review

* Rename wsum and wsum2

Author: henryiii

CHANGELOG.md

@@ -3,6 +3,7 @@
 #### User changes
 
 * Histogram and Axis classes now follow PEP 8 naming scheme (`histogram`->`Histogram`, `regular`->`Regular`, etc.) [#192][]
+* You can now view a histogram with accumulators, with property access such as `h.view().value` [#194][]
 * Added axes transforms [#192][]
 * You can now sum over a range with endpoints [#185][]
 * You can now access the functional regular axis directly, `regular_sqrt` becomes `regular.sqrt`, etc.. [#183][]
@@ -12,6 +13,10 @@
 * Signatures are much nicer in Python 3 [#188][]
 
 
+#### Bug fixes:
+* Unlimited and AtomicInt storages now allow single item access [#194][]
+* `.view()` now no longer makes a copy [#194][]
+
 #### Developer changes
 
 * The hist/axis classes are now pure Python, with a C++ object inside [#183][]

README.md

@@ -18,7 +18,6 @@ Python bindings for [Boost::Histogram][] ([source][Boost::Histogram source]), a
 > Join the [discussion on gitter][gitter-link] or [open an issue](https://github.com/scikit-hep/boost-histogram/issues)!
 >
 > #### Known issues (develop):
-> * Non-simple storages do not support `.view()` or the buffer interface; you can access and set one element at a time.
 > * Setting with an array is not yet supported (`h[...] = np.array(...)`).
 
 

boost_histogram/_internal/hist.py

@@ -4,6 +4,7 @@
 
 from .. import _core
 from .axis import _to_axis, Axis
+from .view import _to_view
 from .axistuple import AxesTuple
 from .sig_tools import inject_signature
 from .storage import Double
@@ -53,13 +54,13 @@ def _expand_ellipsis(indexes, rank):
         raise IndexError("an index can only have a single ellipsis ('...')")
 
 
-def _compute_commonindex(hist, index, expand):
+def _compute_commonindex(hist, index, expand_ellipsis):
     # Normalize -> h[i] == h[i,]
     if not isinstance(index, tuple):
         index = (index,)
 
     # Now a list
-    if expand:
+    if expand_ellipsis:
         indexes = _expand_ellipsis(index, hist.rank())
     else:
         indexes = list(index)
@@ -135,9 +136,7 @@ def __repr__(self):
         return self.__class__.__name__ + repr(self._hist)[9:]
 
     def __array__(self):
-        return np.asarray(self._hist)
-        # TODO: .view does not seem to return an editable view
-        #        so we have to use the buffer interface here
+        return self.view()
 
     def __add__(self, other):
         return self.__class__(self._hist + other._hist)
@@ -296,7 +295,7 @@ def view(self, flow=False):
         """
         Return a view into the data, optionally with overflow turned on.
         """
-        return self._hist.view(flow)
+        return _to_view(self._hist.view(flow))
 
     def reset(self):
         """
@@ -334,7 +333,7 @@ def size(self):
 
     def __getitem__(self, index):
 
-        indexes = _compute_commonindex(self._hist, index, expand=True)
+        indexes = _compute_commonindex(self._hist, index, expand_ellipsis=True)
 
         # If this is (now) all integers, return the bin contents
         try:
@@ -406,5 +405,5 @@ def __getitem__(self, index):
             )
 
     def __setitem__(self, index, value):
-        indexes = _compute_commonindex(self._hist, index, expand=True)
+        indexes = _compute_commonindex(self._hist, index, expand_ellipsis=False)
         self._hist._at_set(value, *indexes)

boost_histogram/_internal/view.py

@@ -0,0 +1,120 @@
+from __future__ import absolute_import, division, print_function
+
+from ..accumulators import Mean, WeightedMean, WeightedSum
+
+import numpy as np
+
+
+class View(np.ndarray):
+    __slots__ = ()
+
+    def __getitem__(self, ind):
+        sliced = super(View, self).__getitem__(ind)
+
+        # If the shape is empty, return the parent type
+        if not sliced.shape:
+            return self._PARENT._make(*sliced)
+        # If the dtype has changed, return a normal array (no longer a record)
+        elif sliced.dtype != self.dtype:
+            return np.asarray(sliced)
+        # Otherwise, no change, return the same View type
+        else:
+            return sliced
+
+    def __repr__(self):
+        # Numpy starts the ndarray class name with "array", so we replace it
+        # with our class name
+        return (
+            "{self.__class__.__name__}(\n      ".format(self=self)
+            + repr(self.view(np.ndarray))[6:]
+        )
+
+    def __str__(self):
+        fields = ", ".join(self._FIELDS)
+        return "{self.__class__.__name__}: ({fields})\n{arr}".format(
+            self=self, fields=fields, arr=self.view(np.ndarray)
+        )
+
+
+def fields(*names):
+    """
+    This decorator adds the name to the _FIELDS
+    class property (for printing in reprs), and
+    adds a property that looks like this:
+
+    @property
+    def name(self):
+        return self["name"]
+    @name.setter
+    def name(self, value):
+        self["name"] = value
+    """
+
+    def injector(cls):
+        if hasattr(cls, "_FIELDS"):
+            raise RuntimeError(
+                "{0} already has had a fields decorator applied".format(
+                    self.__class__.__name__
+                )
+            )
+        fields = []
+        for name in names:
+            fields.append(name)
+
+            def fget(self):
+                return self[name]
+
+            def fset(self, value):
+                self[name] = value
+
+            setattr(cls, name, property(fget, fset))
+            cls._FIELDS = tuple(fields)
+        return cls
+
+    return injector
+
+
+@fields("value", "variance")
+class WeightedSumView(View):
+    __slots__ = ()
+    _PARENT = WeightedSum
+
+
+@fields(
+    "sum_of_weights",
+    "sum_of_weights_squared",
+    "value",
+    "sum_of_weighted_deltas_squared",
+)
+class WeightedMeanView(View):
+    __slots__ = ()
+    _PARENT = WeightedMean
+
+    @property
+    def variance(self):
+        return self["sum_of_weighted_deltas_squared"] / (
+            self["sum_of_weights"]
+            - self["sum_of_weights_squared"] / self["sum_of_weights"]
+        )
+
+
+@fields("count", "value", "sum_of_deltas_squared")
+class MeanView(View):
+    __slots__ = ()
+    _PARENT = Mean
+
+    # Variance is a computation
+    @property
+    def variance(self):
+        return self["sum_of_deltas_squared"] / (self["count"] - 1)
+
+
+def _to_view(item, value=False):
+    for cls in View.__subclasses__():
+        if cls._FIELDS == item.dtype.names:
+            ret = item.view(cls)
+            if value and ret.shape:
+                return ret.value
+            else:
+                return ret
+    return item

include/boost/histogram/python/accumulators/mean.hpp

@@ -0,0 +1,96 @@
+// Copyright 2015-2019 Hans Dembinski and Henry Schreiner
+//
+// Distributed under the Boost Software License, version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+//
+// Based on boost/histogram/accumulators/mean.hpp
+// Changes:
+//  * Internal values are public for access from Python
+//  * A special constructor added for construction from Python
+
+#pragma once
+
+#include <boost/core/nvp.hpp>
+#include <boost/histogram/weight.hpp>
+
+namespace accumulators {
+
+/** Calculates mean and variance of sample.
+
+  Uses Welford's incremental algorithm to improve the numerical
+  stability of mean and variance computation.
+*/
+template <class RealType>
+struct mean {
+    mean() = default;
+    mean(const RealType &n, const RealType &mean, const RealType &variance) noexcept
+        : count(n)
+        , value(mean)
+        , sum_of_deltas_squared(variance * (n - 1)) {}
+
+    mean(const RealType &sum,
+         const RealType &mean,
+         const RealType &sum_of_deltas_squared,
+         bool /* Tag to trigger python internal constructor */)
+        : count(sum)
+        , value(mean)
+        , sum_of_deltas_squared(sum_of_deltas_squared) {}
+
+    void operator()(const RealType &x) noexcept {
+        count += static_cast<RealType>(1);
+        const auto delta = x - value;
+        value += delta / count;
+        sum_of_deltas_squared += delta * (x - value);
+    }
+
+    void operator()(const boost::histogram::weight_type<RealType> &w,
+                    const RealType &x) noexcept {
+        count += w.value;
+        const auto delta = x - value;
+        value += w.value * delta / count;
+        sum_of_deltas_squared += w.value * delta * (x - value);
+    }
+
+    template <class T>
+    mean &operator+=(const mean<T> &rhs) noexcept {
+        if(count != 0 || rhs.count != 0) {
+            const auto tmp
+                = value * count + static_cast<RealType>(rhs.value * rhs.count);
+            count += rhs.count;
+            value = tmp / count;
+        }
+        sum_of_deltas_squared += static_cast<RealType>(rhs.sum_of_deltas_squared);
+        return *this;
+    }
+
+    mean &operator*=(const RealType &s) noexcept {
+        value *= s;
+        sum_of_deltas_squared *= s * s;
+        return *this;
+    }
+
+    template <class T>
+    bool operator==(const mean<T> &rhs) const noexcept {
+        return count == rhs.count && value == rhs.value
+               && sum_of_deltas_squared == rhs.sum_of_deltas_squared;
+    }
+
+    template <class T>
+    bool operator!=(const mean<T> &rhs) const noexcept {
+        return !operator==(rhs);
+    }
+
+    RealType variance() const noexcept { return sum_of_deltas_squared / (count - 1); }
+
+    template <class Archive>
+    void serialize(Archive &ar, unsigned) {
+        ar &boost::make_nvp("count", count);
+        ar &boost::make_nvp("value", value);
+        ar &boost::make_nvp("sum_of_deltas_squared", sum_of_deltas_squared);
+    }
+
+    RealType count = 0, value = 0, sum_of_deltas_squared = 0;
+};
+
+} // namespace accumulators