diff --git a/build/CMakeLists.txt b/build/CMakeLists.txt
index b59685f2..4cb1ea58 100644
--- a/build/CMakeLists.txt
+++ b/build/CMakeLists.txt
@@ -82,10 +82,10 @@ if(BUILD_CHECKS)
   target_link_libraries(speed_cpp_hs_ss ${LIBRARIES})
   target_compile_definitions(speed_cpp_hs_ss PUBLIC HISTOGRAM_TYPE=1)
 
-  add_executable(speed_cpp_hs_ds
+  add_executable(speed_cpp_hs_sd
                  ../test/check/speed_cpp.cpp)
-  target_link_libraries(speed_cpp_hs_ds ${LIBRARIES})
-  target_compile_definitions(speed_cpp_hs_ds PUBLIC HISTOGRAM_TYPE=2)
+  target_link_libraries(speed_cpp_hs_sd ${LIBRARIES})
+  target_compile_definitions(speed_cpp_hs_sd PUBLIC HISTOGRAM_TYPE=2)
 
   add_executable(speed_cpp_hd_ss
                  ../test/check/speed_cpp.cpp)
diff --git a/include/boost/histogram.hpp b/include/boost/histogram.hpp
index 9e7e9835..3c9d22bb 100644
--- a/include/boost/histogram.hpp
+++ b/include/boost/histogram.hpp
@@ -7,11 +7,9 @@
 #ifndef BOOST_HISTOGRAM_HPP_
 #define BOOST_HISTOGRAM_HPP_
 
-#include <boost/histogram/histogram.hpp>
-#include <boost/histogram/axis.hpp>
-#include <boost/histogram/visitors.hpp>
-#include <boost/histogram/dynamic_storage.hpp>
-#include <boost/histogram/static_storage.hpp>
+#include <boost/histogram/dynamic_histogram.hpp>
+#include <boost/histogram/static_histogram.hpp>
+#include <boost/histogram/utility.hpp>
 
 /**
  * \file boost/histogram.hpp
diff --git a/include/boost/histogram/detail/axis_visitor.hpp b/include/boost/histogram/detail/axis_visitor.hpp
new file mode 100644
index 00000000..729353b4
--- /dev/null
+++ b/include/boost/histogram/detail/axis_visitor.hpp
@@ -0,0 +1,95 @@
+// Copyright 2015-2016 Hans Dembinski
+//
+// 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)
+
+#ifndef _BOOST_HISTOGARM_AXIS_VISITOR_HPP_
+#define _BOOST_HISTOGARM_AXIS_VISITOR_HPP_
+
+#include <boost/variant/static_visitor.hpp>
+
+namespace boost {
+namespace histogram {
+namespace detail {
+
+  struct bins : public static_visitor<int>
+  {
+    template <typename A>
+    int operator()(const A& a) const { return a.bins(); }
+  };
+
+  struct shape : public static_visitor<int>
+  {
+    template <typename A>
+    int operator()(const A& a) const { return a.shape(); }
+  };
+
+  struct uoflow : public static_visitor<bool>
+  {
+    template <typename A>
+    bool operator()(const A& a) const { return a.uoflow(); }
+  };
+
+  template <typename V>
+  struct index : public static_visitor<int>
+  {
+    const V v;
+    index(const V x) : v(x) {}
+    template <typename A>
+    int operator()(const A& a) const { return a.index(v); }
+  };
+
+  struct cmp_axis : 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; }
+  };
+
+  struct linearize : public static_visitor<void> {
+    int j;
+    std::size_t out = 0, stride = 1;
+
+    void set(int in) { j = in; }
+
+    template <typename A>
+    void operator()(const A& a) {
+      // the following is highly optimized code that runs in a hot loop;
+      // please measure the performance impact of changes
+      const int uoflow = a.uoflow();
+      // set stride to zero if 'j' is not in range,
+      // this communicates the out-of-range condition to the caller
+      stride *= (j >= -uoflow) * (j < (a.bins() + uoflow));
+      j += (j < 0) * (a.bins() + 2); // wrap around if in < 0
+      out += j * stride;
+      #pragma GCC diagnostic ignored "-Wstrict-overflow"
+      stride *= a.shape();
+    }
+  };
+
+  struct linearize_x : public static_visitor<void> {
+    double x;
+    std::size_t out = 0, stride = 1;
+
+    void set(double in) { x = in; }
+
+    template <typename A>
+    void
+    operator()(const A& a) {
+      // j is guaranteed to be in range [-1, bins]
+      int j = a.index(x);
+      j += (j < 0) * (a.bins() + 2); // wrap around if j < 0
+      out += j * stride;
+      #pragma GCC diagnostic ignored "-Wstrict-overflow"
+      stride *= (j < a.shape()) * a.shape(); // stride == 0 indicates out-of-range
+    }
+  };
+
+}
+}
+}
+
+#endif
diff --git a/include/boost/histogram/detail/utility.hpp b/include/boost/histogram/detail/utility.hpp
index 2056fd31..e060875b 100644
--- a/include/boost/histogram/detail/utility.hpp
+++ b/include/boost/histogram/detail/utility.hpp
@@ -9,56 +9,26 @@
 
 #include <string>
 #include <type_traits>
-#include <boost/variant/static_visitor.hpp>
 
 namespace boost {
 namespace histogram {
 namespace detail {
 
+    // assumes that Storage1 and Storage2 have equal size
+    template <typename Storage1, typename Storage2>
+    bool storage_content_equal(const Storage1& s1, const Storage2& s2)
+    {
+      for (unsigned i = 0, n = s1.size(); i < n; ++i)
+        if (s1.value(i) != s2.value(i) || s1.variance(i) != s2.variance(i))
+          return false;
+      return true;
+    }
+
     template <typename T,
               typename = decltype(*std::declval<T&>()),
               typename = decltype(++std::declval<T&>())>
     struct is_iterator {};
 
-    struct linearize : public static_visitor<void> {
-        int j;
-        std::size_t out = 0, stride = 1;
-
-        void set(int in) { j = in; }
-
-        template <typename A>
-        void operator()(const A& a) {
-          // the following is highly optimized code that runs in a hot loop;
-          // please measure the performance impact of changes
-          const int uoflow = a.uoflow();
-          // set stride to zero if 'j' is not in range,
-          // this communicates the out-of-range condition to the caller
-          stride *= (j >= -uoflow) * (j < (a.bins() + uoflow));
-          j += (j < 0) * (a.bins() + 2); // wrap around if in < 0
-          out += j * stride;
-          #pragma GCC diagnostic ignored "-Wstrict-overflow"
-          stride *= a.shape();
-        }
-    };
-
-    struct linearize_x : public static_visitor<void> {
-        double x;
-        std::size_t out = 0, stride = 1;
-
-        void set(double in) { x = in; }
-
-        template <typename A>
-        void
-        operator()(const A& a) {
-          // j is guaranteed to be in range [-1, bins]
-          int j = a.index(x);
-          j += (j < 0) * (a.bins() + 2); // wrap around if j < 0
-          out += j * stride;
-          #pragma GCC diagnostic ignored "-Wstrict-overflow"
-          stride *= (j < a.shape()) * a.shape(); // stride == 0 indicates out-of-range
-        }
-    };
-
     inline
     std::string escape(const std::string& s) {
         std::string os;
diff --git a/include/boost/histogram/dynamic_histogram.hpp b/include/boost/histogram/dynamic_histogram.hpp
index f0de4619..ae8a0762 100644
--- a/include/boost/histogram/dynamic_histogram.hpp
+++ b/include/boost/histogram/dynamic_histogram.hpp
@@ -12,11 +12,12 @@
 #include <boost/mpl/vector.hpp>
 #include <boost/mpl/empty.hpp>
 #include <boost/variant.hpp>
-#include <boost/variant/static_visitor.hpp>
 #include <boost/histogram/axis.hpp>
 #include <boost/histogram/static_storage.hpp>
 #include <boost/histogram/dynamic_storage.hpp>
 #include <boost/histogram/utility.hpp>
+#include <boost/histogram/detail/utility.hpp>
+#include <boost/histogram/detail/axis_visitor.hpp>
 #include <cstddef>
 #include <array>
 #include <vector>
@@ -99,16 +100,15 @@ public:
     return *this;
   } 
 
-  template <typename OtherStorage>
-  bool operator==(const dynamic_histogram<OtherStorage, Axes>& other) const
+  template <typename OtherStorage, typename OtherAxes>
+  bool operator==(const dynamic_histogram<OtherStorage, OtherAxes>& other) const
   {
     if (dim() != other.dim())
       return false;
-    for (decltype(dim()) i = 0, n = dim(); i < n; ++i)
-      if (!(axes_[i] == other.axes_[i]))
-        return false;
-    if (!(storage_ == other.storage_))
-        return false;
+    if (!axes_equal_to(other.axes_))
+      return false;
+    if (!detail::storage_content_equal(storage_, other.storage_))
+      return false;
     return true;
   }
 
@@ -119,7 +119,7 @@ public:
       throw std::logic_error("dimensions of histograms differ");
     if (size() != other.size())
       throw std::logic_error("sizes of histograms differ");
-    if (axes_ != other.axes_)
+    if (!axes_equal_to(other.axes_))
       throw std::logic_error("axes of histograms differ");
     storage_ += other.storage_;
     return *this;
@@ -136,8 +136,9 @@ public:
       storage_.increase(lin.out);
   }
 
-  template <typename Iterator>
-  void fill_iter(Iterator values_begin, Iterator values_end)
+  template <typename Iterator,
+            typename = detail::is_iterator<Iterator>>
+  void fill(Iterator values_begin, Iterator values_end)
   {
     BOOST_ASSERT_MSG(std::distance(values_begin, values_end) == dim(),
                      "number of arguments does not match histogram dimension");
@@ -150,6 +151,8 @@ public:
   template <typename... Args>
   void wfill(Args... args)
   {
+    static_assert(std::is_same<Storage, dynamic_storage>::value,
+                  "wfill only supported for dynamic_storage");
     BOOST_ASSERT_MSG((sizeof...(args) - 1) == dim(),
                      "number of arguments does not match histogram dimension");
     detail::linearize_x lin;
@@ -158,9 +161,12 @@ public:
       storage_.increase(lin.out, w);
   }
 
-  template <typename Iterator>
-  void wfill_iter(Iterator values_begin, Iterator values_end, double w)
+  template <typename Iterator,
+            typename = detail::is_iterator<Iterator>>
+  void wfill(Iterator values_begin, Iterator values_end, double w)
   {
+    static_assert(std::is_same<Storage, dynamic_storage>::value,
+                  "wfill only supported for dynamic_storage");
     BOOST_ASSERT_MSG(std::distance(values_begin, values_end) == dim(),
                      "number of arguments does not match histogram dimension");
     detail::linearize_x lin;
@@ -181,8 +187,9 @@ public:
     return storage_.value(lin.out);
   }
 
-  template <typename Iterator>
-  value_t value_iter(Iterator indices_begin, Iterator indices_end) const
+  template <typename Iterator,
+            typename = detail::is_iterator<Iterator>>
+  value_t value(Iterator indices_begin, Iterator indices_end) const
   {
     BOOST_ASSERT_MSG(std::distance(indices_begin, indices_end) == dim(),
                      "number of arguments does not match histogram dimension");
@@ -206,8 +213,9 @@ public:
     return storage_.variance(lin.out);
   }
 
-  template <typename Iterator>
-  value_t variance_iter(Iterator indices_begin, Iterator indices_end) const
+  template <typename Iterator,
+            typename = detail::is_iterator<Iterator>>
+  value_t variance(Iterator indices_begin, Iterator indices_end) const
   {
     BOOST_ASSERT_MSG(std::distance(indices_begin, indices_end) == dim(),
                      "number of arguments does not match histogram dimension");
@@ -253,6 +261,16 @@ private:
     return fc;
   }
 
+  template <typename OtherAxes>
+  bool axes_equal_to(const OtherAxes& other_axes) const
+  {
+    detail::cmp_axis ca;
+    for (unsigned i = 0; i < dim(); ++i)
+      if (!apply_visitor(ca, axes_[i], other_axes[i]))
+        return false;
+    return true;
+  }
+
   template <typename Linearize, typename First, typename... Rest>
   void index_impl(Linearize& lin, First first, Rest... rest) const
   {
@@ -326,6 +344,25 @@ operator+(const dynamic_histogram<StoragePolicyA, Axes>& a,
   return tmp;
 }
 
+
+template <typename... Axes>
+inline
+dynamic_histogram<>
+make_dynamic_histogram(const Axes&... axes)
+{
+  return dynamic_histogram<>(axes...);
+}
+
+
+template <typename Storage, typename... Axes>
+inline
+dynamic_histogram<Storage>
+make_dynamic_histogram_with(const Axes&... axes)
+{
+  return dynamic_histogram<Storage>(axes...);
+}
+
+
 }
 }
 
diff --git a/include/boost/histogram/dynamic_storage.hpp b/include/boost/histogram/dynamic_storage.hpp
index 77f35ed3..794635f3 100644
--- a/include/boost/histogram/dynamic_storage.hpp
+++ b/include/boost/histogram/dynamic_storage.hpp
@@ -45,7 +45,6 @@ class dynamic_storage {
   using buffer_t = detail::buffer_t;
 
 public:
-  struct storage_tag {};
   using value_t = double;
   using variance_t = double;
 
diff --git a/include/boost/histogram/static_histogram.hpp b/include/boost/histogram/static_histogram.hpp
index 457e1bd2..e19bc152 100644
--- a/include/boost/histogram/static_histogram.hpp
+++ b/include/boost/histogram/static_histogram.hpp
@@ -14,6 +14,8 @@
 #include <boost/fusion/include/mpl.hpp>
 #include <boost/fusion/sequence.hpp>
 #include <boost/fusion/include/sequence.hpp>
+#include <boost/fusion/sequence/comparison.hpp>
+#include <boost/fusion/include/comparison.hpp>
 #include <boost/fusion/algorithm.hpp>
 #include <boost/fusion/include/algorithm.hpp>
 #include <boost/fusion/container/vector/convert.hpp>
@@ -22,6 +24,7 @@
 #include <boost/histogram/static_storage.hpp>
 #include <boost/histogram/dynamic_storage.hpp>
 #include <boost/histogram/detail/utility.hpp>
+#include <boost/histogram/detail/axis_visitor.hpp>
 #include <type_traits>
 
 namespace boost {
@@ -33,14 +36,11 @@ class static_histogram
   static_assert(!mpl::empty<Axes>::value, "at least one axis required");
 
 public:
-  using axes_t = typename mpl::transform<
-      typename fusion::result_of::as_vector<Axes>::type,
-      std::add_const<mpl::_1>
-    >::type;
+  using axes_t = typename fusion::result_of::as_vector<Axes>::type;
   using value_t = typename Storage::value_t;
   using variance_t = typename Storage::variance_t;
 
-  static_histogram() = delete; // makes no sense
+  static_histogram() = default;
 
   template <typename... Axes1>
   explicit
@@ -61,12 +61,12 @@ public:
 
   static_histogram& operator=(const static_histogram& other)
   {
-    if (this != &other) {
+    if (static_cast<const void*>(this) != static_cast<const void*>(&other)) {
       axes_ = other.axes_;
       storage_ = other.storage_;
     }
     return *this;
-  }
+  } 
 
   static_histogram& operator=(static_histogram&& other)
   {
@@ -75,19 +75,45 @@ public:
     return *this;
   } 
 
+  template <typename OtherStorage>
+  static_histogram(const static_histogram<OtherStorage, Axes>& other) :
+    axes_(other.axes_), storage_(other.storage_)
+  {}
+
+  template <typename OtherStorage>
+  static_histogram(static_histogram<OtherStorage, Axes>&& other) :
+    axes_(std::move(other.axes_)),
+    storage_(std::move(other.storage_))
+  {}
+
+  template <typename OtherStorage>
+  static_histogram& operator=(const static_histogram<OtherStorage, Axes>& other)
+  {
+    if (static_cast<const void*>(this) != static_cast<const void*>(&other)) {
+      axes_ = other.axes_;
+      storage_ = other.storage_;
+    }
+    return *this;
+  }
+
+  template <typename OtherStorage>
+  static_histogram& operator=(static_histogram<OtherStorage, Axes>&& other)
+  {
+    axes_ = std::move(other.axes_);
+    storage_ = std::move(other.storage_);
+    return *this;
+  }
+
   template <typename OtherStorage>
   bool operator==(const static_histogram<OtherStorage, Axes>& other) const
   {
-    if (axes_ != other.axes_)
+    if (!axes_equal_to(other.axes_))
       return false;
-    for (decltype(size()) i = 0, n = size(); i < n; ++i)
-      if (storage_.value(i) != other.storage_.value(i))
-        return false;
-    return true;
+    return detail::storage_content_equal(storage_, other.storage_);
   }
 
   template <typename OtherStorage, typename OtherAxes>
-  bool operator==(const static_histogram<OtherStorage, OtherAxes>& other) const
+  constexpr bool operator==(const static_histogram<OtherStorage, OtherAxes>& other)
   {
     return false;
   }
@@ -95,7 +121,7 @@ public:
   template <typename OtherStorage>
   static_histogram& operator+=(const static_histogram<OtherStorage, Axes>& other)
   {
-    if (axes_ != other.axes_)
+    if (!axes_equal_to(other.axes_))
       throw std::logic_error("axes of histograms differ");
     storage_ += other.storage_;
     return *this;
@@ -105,52 +131,55 @@ public:
   void fill(Vs... values)
   {
     static_assert(sizeof...(values) == dim(),
-                  "number of arguments does not match static_histogram dimension");
+                  "number of arguments does not match histogram dimension");
     detail::linearize_x lin;
     index_impl(lin, values...);
     if (lin.stride)
       storage_.increase(lin.out);
   }
 
-  // template <typename... Args>
-  // void wfill(Args... args)
-  // {
-  //   BOOST_ASSERT_MSG(sizeof...(args) == (dim() + 1),
-  //                    "number of arguments does not match static_histogram dimension");
-  //   linearize_x lin;
-  //   const double w = windex_impl(lin, args...);
-  //   if (lin.stride)
-  //     storage_.increase(lin.out, w);
-  // }
+  template <typename... Args>
+  void wfill(Args... args)
+  {
+    static_assert(std::is_same<Storage, dynamic_storage>::value,
+                  "wfill only supported for dynamic_storage");
+    static_assert(sizeof...(args) == (dim() + 1),
+                  "number of arguments does not match histogram dimension");
+    detail::linearize_x lin;
+    const double w = windex_impl(lin, args...);
+    if (lin.stride)
+      storage_.increase(lin.out, w);
+  }
 
-  // template <typename... Args>
-  // value_t value(Args... args) const
-  // {
-  //   BOOST_ASSERT_MSG(sizeof...(args) == dim(),
-  //                    "number of arguments does not match static_histogram dimension");
-  //   linearize lin;
-  //   index_impl(lin, args...);
-  //   if (lin.stride == 0)
-  //     throw std::out_of_range("invalid index");
-  //   return storage_.value(lin.out);
-  // }
+  template <typename... Args>
+  value_t value(Args... args) const
+  {
+    static_assert(sizeof...(args) == dim(),
+                  "number of arguments does not match histogram dimension");
+    detail::linearize lin;
+    index_impl(lin, args...);
+    if (lin.stride == 0)
+      throw std::out_of_range("invalid index");
+    return storage_.value(lin.out);
+  }
 
-  // template <typename... Args>
-  // variance_t variance(Args... args) const
-  // {
-  //   BOOST_ASSERT_MSG(sizeof...(args) == dim(),
-  //                    "number of arguments does not match static_histogram dimension");
-  //   linearize lin;
-  //   index_impl(lin, args...);
-  //   if (lin.stride == 0)
-  //     throw std::out_of_range("invalid index");
-  //   return storage_.variance(lin.out);
-  // }
+  template <typename... Args>
+  variance_t variance(Args... args) const
+  {
+    static_assert(sizeof...(args) == dim(),
+                  "number of arguments does not match histogram dimension");
+    detail::linearize lin;
+    index_impl(lin, args...);
+    if (lin.stride == 0)
+      throw std::out_of_range("invalid index");
+    return storage_.variance(lin.out);
+  }
 
-  /// Number of axes (dimensions) of static_histogram
-  constexpr unsigned dim() const { return fusion::size(axes_); }
+  /// Number of axes (dimensions) of histogram
+  static constexpr unsigned dim()
+  { return fusion::result_of::size<Axes>::type::value; }
 
-  /// Total number of bins in the static_histogram (including underflow/overflow)
+  /// Total number of bins in the histogram (including underflow/overflow)
   std::size_t size() const { return storage_.size(); }
 
   /// Access to internal data (used by Python bindings)
@@ -159,7 +188,7 @@ public:
   /// Memory used by a bin in bytes
   unsigned depth() const { return storage_.depth(); }
 
-  /// Sum of all counts in the static_histogram
+  /// Sum of all counts in the histogram
   double sum() const
   {
     double result = 0.0;
@@ -169,7 +198,9 @@ public:
   }
 
   template <unsigned N>
-  typename fusion::result_of::at_c<axes_t, N>::type&
+  typename std::add_const<
+    typename fusion::result_of::value_at_c<axes_t, N>::type
+  >::type&
   axis() const
   { return fusion::at_c<N>(axes_); }
 
@@ -190,16 +221,52 @@ private:
     return fc.value;
   }
 
-  // template <typename Linearize, typename First, typename... Rest>
-  // void index_impl(Linearize& lin, First first, Rest... rest) const
-  // {
-  //   lin.set(first);
-  //   lin(std::get<(dim() - sizeof...(Rest) - 1)>(axes_));
-  //   index_impl(lin, rest...);
-  // }
+  template <typename OtherAxes>
+  bool axes_equal_to(const OtherAxes& other_axes) const {
+    return false;
+  }
 
-  // template <typename Linearize>
-  // void index_impl(Linearize&) const {} // stop recursion
+  bool axes_equal_to(const axes_t& other_axes) const {
+    return axes_ == other_axes;
+  }
+
+  template <typename Linearize, typename First, typename... Rest>
+  void index_impl(Linearize& lin, First first, Rest... rest) const
+  {
+    lin.set(first);
+    lin(fusion::at_c<(dim() - sizeof...(Rest) - 1)>(axes_));
+    index_impl(lin, rest...);
+  }
+
+  template <typename Linearize>
+  void index_impl(Linearize&) const {} // stop recursion
+
+  template <typename Linearize, typename First, typename... Rest>
+  double windex_impl(Linearize& lin, First first, Rest... rest) const
+  {
+    lin.set(first);
+    lin(fusion::at_c<(dim() - sizeof...(Rest))>(axes_));
+    return windex_impl(lin, rest...);
+  }
+
+  template <typename Linearize, typename Last>
+  double windex_impl(Linearize& lin, Last last) const
+  {
+    return last;
+  }
+
+  template <typename Linearize, typename Iterator>
+  void iter_args_impl(Linearize& lin,
+                      Iterator& args_begin,
+                      const Iterator& args_end) const {
+    auto axes_iter = axes_.begin();
+    while (args_begin != args_end) {
+      lin.set(*args_begin);
+      apply_visitor(lin, *axes_iter);
+      ++args_begin;
+      ++axes_iter;
+    }
+  }
 
   template <typename OtherStorage, typename OtherAxes>
   friend class static_histogram;
@@ -209,6 +276,34 @@ private:
 };
 
 
+// when adding histograms with different storage, use storage with more capacity as return type
+template <typename StoragePolicyA,
+          typename StoragePolicyB,
+          typename Axes>
+inline
+static_histogram<
+  typename std::conditional<
+    (std::numeric_limits<typename StoragePolicyA::value_t>::max() >
+     std::numeric_limits<typename StoragePolicyB::value_t>::max()),
+    StoragePolicyA, StoragePolicyB
+  >::type,
+  Axes
+>
+operator+(const static_histogram<StoragePolicyA, Axes>& a,
+          const static_histogram<StoragePolicyB, Axes>& b)
+{
+  static_histogram<
+    typename std::conditional<
+      (std::numeric_limits<typename StoragePolicyA::value_t>::max() >
+       std::numeric_limits<typename StoragePolicyB::value_t>::max()),
+      StoragePolicyA, StoragePolicyB>::type,
+    Axes
+  > tmp = a;
+  tmp += b;
+  return tmp;
+}
+
+
 /// default static type factory
 template <typename... Axes>
 inline
@@ -218,6 +313,7 @@ make_static_histogram(const Axes&... axes)
   return static_histogram<dynamic_storage, mpl::vector<Axes...>>(axes...);
 }
 
+
 /// static type factory with variable storage type
 template <typename Storage, typename... Axes>
 inline
diff --git a/include/boost/histogram/static_storage.hpp b/include/boost/histogram/static_storage.hpp
index d09030d6..d54a8e86 100644
--- a/include/boost/histogram/static_storage.hpp
+++ b/include/boost/histogram/static_storage.hpp
@@ -17,11 +17,12 @@ namespace histogram {
     using buffer_t = detail::buffer_t;
 
   public:
-    struct storage_tag {};
     using value_t = T;
     using variance_t = T;
 
-    explicit static_storage(std::size_t n=0) : data_(n * sizeof(T)) {}
+    static_storage() = default;
+
+    explicit static_storage(std::size_t n) : data_(n * sizeof(T)) {}
 
     static_storage(const static_storage<T>& other) :
       data_(other.data_)
@@ -40,7 +41,7 @@ namespace histogram {
     }
 
     template <typename OtherStorage,
-              typename = typename OtherStorage::storage_tag>
+              typename = typename OtherStorage::value_t>
     static_storage(const OtherStorage& other) :
       data_(other.size() * sizeof(T))
     {
@@ -49,7 +50,7 @@ namespace histogram {
     }
 
     template <typename OtherStorage,
-              typename = typename OtherStorage::storage_tag>
+              typename = typename OtherStorage::value_t>
     static_storage(OtherStorage&& other) :
       data_(other.size() * sizeof(T))
     {
@@ -58,8 +59,7 @@ namespace histogram {
       other = OtherStorage();
     }
 
-    template <typename OtherStorage,
-              typename = typename OtherStorage::storage_tag>
+    template <typename OtherStorage>
     static_storage& operator=(const OtherStorage& other)
     {
       if (static_cast<const void*>(this) != static_cast<const void*>(&other)) {
@@ -70,8 +70,7 @@ namespace histogram {
       return *this;
     }
 
-    template <typename OtherStorage,
-              typename = typename OtherStorage::storage_tag>
+    template <typename OtherStorage>
     static_storage& operator=(OtherStorage&& other)
     {
       if (static_cast<void*>(this) != static_cast<void*>(&other)) {
@@ -90,8 +89,7 @@ namespace histogram {
     value_t value(std::size_t i) const { return data_.at<T>(i); }
     variance_t variance(std::size_t i) const { return data_.at<T>(i); }
 
-    template <typename OtherStorage,
-              typename = typename OtherStorage::storage_tag>
+    template <typename OtherStorage>
     void operator+=(const OtherStorage& other)
     {
       for (std::size_t i = 0, n = size(); i < n; ++i)
diff --git a/include/boost/histogram/utility.hpp b/include/boost/histogram/utility.hpp
index daaa366c..a9c64dc9 100644
--- a/include/boost/histogram/utility.hpp
+++ b/include/boost/histogram/utility.hpp
@@ -2,36 +2,31 @@
 #define BOOST_HISTOGRAM_UTILITY_HPP_
 
 #include <boost/variant.hpp>
-#include <boost/histogram/visitors.hpp>
+#include <boost/histogram/detail/axis_visitor.hpp>
 
 namespace boost {
 namespace histogram {
 
-template <typename Storage1,
-          typename Storage2,
-          typename = typename Storage1::storage_tag,
-          typename = typename Storage2::storage_tag>
-bool operator==(const Storage1& s1, const Storage2& s2)
-{
-    if (s1.size() != s2.size())
-        return false;
-    for (std::size_t i = 0, n = s1.size(); i < n; ++i)
-        if (s1.value(i) != s2.value(i) || s1.variance(i) != s2.variance(i))
-            return false;
-    return true;
-}
-
 template <typename A>
 int bins(const A& a) { return a.bins(); }
 
 template <typename... Axes>
-int bins(const boost::variant<Axes...>& a) { return boost::apply_visitor(visitor::bins(), a); }
+int bins(const boost::variant<Axes...>& a)
+{ return apply_visitor(detail::bins(), a); }
 
 template <typename A>
 int shape(const A& a) { return a.shape(); }
 
 template <typename... Axes>
-int shape(const boost::variant<Axes...>& a) { return boost::apply_visitor(visitor::shape(), a); }
+int shape(const boost::variant<Axes...>& a)
+{ return apply_visitor(detail::shape(), a); }
+
+template <typename A, typename V>
+int index(const A& a, const V v) { return a.index(v); }
+
+template <typename... Axes, typename V>
+int index(const boost::variant<Axes...>& a, const V v)
+{ return apply_visitor(detail::index<V>(v), a); }
 
 }
 }
diff --git a/include/boost/histogram/visitors.hpp b/include/boost/histogram/visitors.hpp
deleted file mode 100644
index 488a4322..00000000
--- a/include/boost/histogram/visitors.hpp
+++ /dev/null
@@ -1,52 +0,0 @@
-// Copyright 2015-2016 Hans Dembinski
-//
-// 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)
-
-#ifndef _BOOST_HISTOGARM_VISITORS_HPP_
-#define _BOOST_HISTOGARM_VISITORS_HPP_
-
-#include <boost/variant/static_visitor.hpp>
-
-namespace boost {
-namespace histogram {
-namespace visitor {
-
-  struct bins : public static_visitor<int>
-  {
-    template <typename A>
-    int operator()(const A& a) const { return a.bins(); }
-  };
-
-  struct shape : public static_visitor<int>
-  {
-    template <typename A>
-    int operator()(const A& a) const { return a.shape(); }
-  };
-
-  struct uoflow : public static_visitor<bool>
-  {
-    template <typename A>
-    bool operator()(const A& a) const { return a.uoflow(); }
-  };
-
-  struct index : public static_visitor<int>
-  {
-    template <typename A, typename V>
-    int operator()(const A& a, const V v) const { return a.index(v); }
-  };
-
-  struct cmp : 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
diff --git a/src/histogram.cpp b/src/histogram.cpp
index 3bb9ac7c..fbaf458a 100644
--- a/src/histogram.cpp
+++ b/src/histogram.cpp
@@ -143,7 +143,7 @@ histogram_fill(python::tuple args, python::dict kwargs) {
           for (unsigned i = 0; i < dims[0]; ++i) {
             double* v = reinterpret_cast<double*>(PyArray_GETPTR1(a, i) );
             double* w = reinterpret_cast<double*>(PyArray_GETPTR1(aw, i));
-            self.wfill_iter(v, v+self.dim(), *w);
+            self.wfill(v, v+self.dim(), *w);
           }
 
           Py_DECREF(aw);
@@ -154,7 +154,7 @@ histogram_fill(python::tuple args, python::dict kwargs) {
       } else {
         for (unsigned i = 0; i < dims[0]; ++i) {
           double* v = reinterpret_cast<double*>(PyArray_GETPTR1(a, i));
-          self.fill_iter(v, v+self.dim());
+          self.fill(v, v+self.dim());
         }
       }
 
@@ -175,11 +175,11 @@ histogram_fill(python::tuple args, python::dict kwargs) {
     v[i] = extract<double>(args[1 + i]);
 
   if (ow.is_none()) {
-    self.fill_iter(v, v+self.dim());
+    self.fill(v, v+self.dim());
 
   } else {
     const double w = extract<double>(ow);
-    self.wfill_iter(v, v+self.dim(), w);
+    self.wfill(v, v+self.dim(), w);
   }
 
   return object();
@@ -204,7 +204,7 @@ histogram_value(python::tuple args, python::dict kwargs) {
   for (unsigned i = 0; i < self.dim(); ++i)
     idx[i] = extract<int>(args[1 + i]);
 
-  return object(self.value_iter(idx + 0, idx + self.dim()));
+  return object(self.value(idx + 0, idx + self.dim()));
 }
 
 python::object
@@ -226,7 +226,7 @@ histogram_variance(python::tuple args, python::dict kwargs) {
   for (unsigned i = 0; i < self.dim(); ++i)
     idx[i] = extract<int>(args[1 + i]);
 
-  return object(self.variance_iter(idx + 0, idx + self.dim()));
+  return object(self.variance(idx + 0, idx + self.dim()));
 }
 
 class histogram_access {
diff --git a/test/check/speed_cpp.cpp b/test/check/speed_cpp.cpp
index a5a1a2a5..af802e8e 100644
--- a/test/check/speed_cpp.cpp
+++ b/test/check/speed_cpp.cpp
@@ -4,6 +4,7 @@
 // (See accompanying file LICENSE_1_0.txt
 // or copy at http://www.boost.org/LICENSE_1_0.txt)
 
+#include <boost/histogram/static_histogram.hpp>
 #include <boost/histogram/dynamic_histogram.hpp>
 #include <boost/histogram/axis.hpp>
 
@@ -15,6 +16,7 @@
 #include <cstdio>
 
 using namespace boost::histogram;
+namespace mpl = boost::mpl; 
 
 std::vector<double> random_array(unsigned n, int type) {
   std::vector<double> result(n);
@@ -109,39 +111,95 @@ int main() {
     printf("1D\n");
     printf("t[boost]   = %.3f\n",
 #if HISTOGRAM_TYPE == 1
-      compare_1d<dynamic::histogram<static_storage<int>>>(12000000, itype)
+      compare_1d<
+        static_histogram<
+          static_storage<int>,
+          mpl::vector<regular_axis>
+        >
+      >(12000000, itype)
 #elif HISTOGRAM_TYPE == 2
-      compare_1d<dynamic::histogram<dynamic_storage>>(12000000, itype)
+      compare_1d<
+        static_histogram<
+          dynamic_storage,
+          mpl::vector<regular_axis>
+        >
+      >(12000000, itype)
 #elif HISTOGRAM_TYPE == 3
-      compare_1d<dynamic::histogram<static_storage<int>>>(12000000, itype)
+      compare_1d<
+        dynamic_histogram<
+          static_storage<int>
+        >
+      >(12000000, itype)
 #elif HISTOGRAM_TYPE == 4
-      compare_1d<dynamic::histogram<dynamic_storage>>(12000000, itype)
+      compare_1d<
+        dynamic_histogram<
+          dynamic_storage
+        >
+      >(12000000, itype)
 #endif
     );
 
     printf("3D\n");
     printf("t[boost]   = %.3f\n",
 #if HISTOGRAM_TYPE == 1
-      compare_3d<dynamic::histogram<static_storage<int>>>(4000000, itype)
+      compare_3d<
+        static_histogram<
+          static_storage<int>,
+          mpl::vector<regular_axis, regular_axis, regular_axis>
+        >
+      >(4000000, itype)
 #elif HISTOGRAM_TYPE == 2
-      compare_3d<dynamic::histogram<dynamic_storage>>(4000000, itype)
+      compare_3d<
+        static_histogram<
+          dynamic_storage,
+          mpl::vector<regular_axis, regular_axis, regular_axis>
+        >
+      >(4000000, itype)
 #elif HISTOGRAM_TYPE == 3
-      compare_3d<dynamic::histogram<static_storage<int>>>(4000000, itype)
+      compare_3d<
+        dynamic_histogram<
+          static_storage<int>
+        >
+      >(4000000, itype)
 #elif HISTOGRAM_TYPE == 4
-      compare_3d<dynamic::histogram<dynamic_storage>>(4000000, itype)
+      compare_3d<
+        dynamic_histogram<
+          dynamic_storage
+        >
+      >(4000000, itype)
 #endif
     );
 
     printf("6D\n");
     printf("t[boost]   = %.3f\n",
 #if HISTOGRAM_TYPE == 1
-      compare_6d<dynamic::histogram<static_storage<int>>>(2000000, itype)
+      compare_6d<
+        static_histogram<
+          static_storage<int>,
+          mpl::vector<regular_axis, regular_axis, regular_axis,
+                      regular_axis, regular_axis, regular_axis>
+        >
+      >(2000000, itype)
 #elif HISTOGRAM_TYPE == 2
-      compare_6d<dynamic::histogram<dynamic_storage>>(2000000, itype)
+      compare_6d<
+        static_histogram<
+          dynamic_storage,
+          mpl::vector<regular_axis, regular_axis, regular_axis,
+                      regular_axis, regular_axis, regular_axis>
+        >
+      >(2000000, itype)
 #elif HISTOGRAM_TYPE == 3
-      compare_6d<dynamic::histogram<static_storage<int>>>(2000000, itype)
+      compare_6d<
+        dynamic_histogram<
+          static_storage<int>
+        >
+      >(2000000, itype)
 #elif HISTOGRAM_TYPE == 4
-      compare_6d<dynamic::histogram<dynamic_storage>>(2000000, itype)
+      compare_6d<
+        dynamic_histogram<
+          dynamic_storage
+        >
+      >(2000000, itype)
 #endif
     );
   }
diff --git a/test/dynamic_histogram_test.cpp b/test/dynamic_histogram_test.cpp
index 158f59e4..b59e9f93 100644
--- a/test/dynamic_histogram_test.cpp
+++ b/test/dynamic_histogram_test.cpp
@@ -80,25 +80,25 @@ BOOST_AUTO_TEST_CASE(init_4)
 
 BOOST_AUTO_TEST_CASE(init_5)
 {
-    auto h = dynamic_histogram<dynamic_storage>(regular_axis{3, -1, 1},
-                                        integer_axis{-1, 1},
-                                        polar_axis{3},
-                                        variable_axis{-1, 0, 1},
-                                        category_axis{"A", "B", "C"});
+    auto h = make_dynamic_histogram(regular_axis{3, -1, 1},
+                                    integer_axis{-1, 1},
+                                    polar_axis{3},
+                                    variable_axis{-1, 0, 1},
+                                    category_axis{"A", "B", "C"});
     BOOST_CHECK_EQUAL(h.dim(), 5);
     BOOST_CHECK_EQUAL(h.size(), 900);
-    auto h2 = dynamic_histogram<dynamic_storage>(regular_axis{3, -1, 1},
-                                         integer_axis{-1, 1},
-                                         polar_axis{3},
-                                         variable_axis{-1, 0, 1},
-                                         category_axis{"A", "B", "C"});
+    auto h2 = make_dynamic_histogram(regular_axis{3, -1, 1},
+                                     integer_axis{-1, 1},
+                                     polar_axis{3},
+                                     variable_axis{-1, 0, 1},
+                                     category_axis{"A", "B", "C"});
     BOOST_CHECK(h2 == h);
 }
 
 BOOST_AUTO_TEST_CASE(copy_ctor)
 {
-    auto h = dynamic_histogram<dynamic_storage>(integer_axis(0, 1),
-                                        integer_axis(0, 2));
+    auto h = make_dynamic_histogram_with<dynamic_storage>(integer_axis(0, 1),
+                                                          integer_axis(0, 2));
     h.fill(0, 0);
     auto h2 = decltype(h)(h);
     BOOST_CHECK(h2 == h);
@@ -108,8 +108,8 @@ BOOST_AUTO_TEST_CASE(copy_ctor)
 
 BOOST_AUTO_TEST_CASE(copy_assign)
 {
-    auto h = dynamic_histogram<dynamic_storage>(integer_axis(0, 1),
-                                        integer_axis(0, 2));
+    auto h = make_dynamic_histogram_with<dynamic_storage>(integer_axis(0, 1),
+                                                          integer_axis(0, 2));
     h.fill(0, 0);
     auto h2 = decltype(h)();
     BOOST_CHECK(!(h == h2));
@@ -125,8 +125,8 @@ BOOST_AUTO_TEST_CASE(copy_assign)
 
 BOOST_AUTO_TEST_CASE(move_ctor)
 {
-    auto h = dynamic_histogram<dynamic_storage>(integer_axis(0, 1),
-                                        integer_axis(0, 2));
+    auto h = make_dynamic_histogram(integer_axis(0, 1),
+                                    integer_axis(0, 2));
     h.fill(0, 0);
     const auto href = h;
     auto h2 = decltype(h)(std::move(h));
@@ -143,8 +143,8 @@ BOOST_AUTO_TEST_CASE(move_ctor)
 
 BOOST_AUTO_TEST_CASE(move_assign)
 {
-    auto h = dynamic_histogram<dynamic_storage>(integer_axis(0, 1),
-                                        integer_axis(0, 2));
+    auto h = make_dynamic_histogram(integer_axis(0, 1),
+                                    integer_axis(0, 2));
     h.fill(0.0, 0.0);
     auto href = h;
     auto h2 = decltype(h)();
@@ -189,7 +189,7 @@ BOOST_AUTO_TEST_CASE(equal_compare)
 
 BOOST_AUTO_TEST_CASE(d1)
 {
-    auto h = dynamic_histogram<dynamic_storage>(integer_axis(0, 1));
+    auto h = make_dynamic_histogram(integer_axis(0, 1));
     h.fill(0);
     h.fill(0);
     h.fill(-1);
@@ -217,7 +217,7 @@ BOOST_AUTO_TEST_CASE(d1)
 
 BOOST_AUTO_TEST_CASE(d1_2)
 {
-    auto h = dynamic_histogram<dynamic_storage>(integer_axis(0, 1, "", false));
+    auto h = make_dynamic_histogram(integer_axis(0, 1, "", false));
     h.fill(0);
     h.fill(-0);
     h.fill(-1);
@@ -241,7 +241,7 @@ BOOST_AUTO_TEST_CASE(d1_2)
 
 BOOST_AUTO_TEST_CASE(d1w)
 {
-    auto h = dynamic_histogram<dynamic_storage>(regular_axis(2, -1, 1));
+    auto h = make_dynamic_histogram(regular_axis(2, -1, 1));
     h.fill(0);
     h.wfill(-1.0, 2.0);
     h.fill(-1.0);
@@ -263,8 +263,8 @@ BOOST_AUTO_TEST_CASE(d1w)
 
 BOOST_AUTO_TEST_CASE(d2)
 {
-    auto h = dynamic_histogram<dynamic_storage>(regular_axis(2, -1, 1),
-                                        integer_axis(-1, 1, std::string(), false));
+    auto h = make_dynamic_histogram(regular_axis(2, -1, 1),
+                                    integer_axis(-1, 1, std::string(), false));
     h.fill(-1, -1);
     h.fill(-1, 0);
     h.fill(-1, -10);
@@ -312,8 +312,8 @@ BOOST_AUTO_TEST_CASE(d2)
 
 BOOST_AUTO_TEST_CASE(d2w)
 {
-    auto h = dynamic_histogram<dynamic_storage>(regular_axis(2, -1, 1),
-                                        integer_axis(-1, 1, std::string(), false));
+    auto h = make_dynamic_histogram(regular_axis(2, -1, 1),
+                                    integer_axis(-1, 1, std::string(), false));
     h.fill(-1, 0);       // -> 0, 1
     h.wfill(-1, -1, 10); // -> 0, 0
     h.wfill(-1, -10, 5); // is ignored
@@ -356,9 +356,9 @@ BOOST_AUTO_TEST_CASE(d2w)
 
 BOOST_AUTO_TEST_CASE(d3w)
 {
-    auto h = dynamic_histogram<dynamic_storage>(integer_axis(0, 3),
-                                                integer_axis(0, 4),
-                                                integer_axis(0, 5));
+    auto h = make_dynamic_histogram(integer_axis(0, 3),
+                                    integer_axis(0, 4),
+                                    integer_axis(0, 5));
     for (auto i = 0; i < bins(h.axis(0)); ++i)
         for (auto j = 0; j < bins(h.axis(1)); ++j)
             for (auto k = 0; k < bins(h.axis(2)); ++k)
@@ -374,17 +374,17 @@ BOOST_AUTO_TEST_CASE(d3w)
 
 BOOST_AUTO_TEST_CASE(add_0)
 {
-    auto a = dynamic_histogram<dynamic_storage>(integer_axis(-1, 1));
-    auto b = dynamic_histogram<dynamic_storage>(regular_axis(3, -1, 1));
-    auto c = dynamic_histogram<dynamic_storage>(regular_axis(3, -1.1, 1));
+    auto a = make_dynamic_histogram(integer_axis(-1, 1));
+    auto b = make_dynamic_histogram(regular_axis(3, -1, 1));
+    auto c = make_dynamic_histogram(regular_axis(3, -1.1, 1));
     BOOST_CHECK_THROW(a + b, std::logic_error);
     BOOST_CHECK_THROW(b + c, std::logic_error);
 }
 
 BOOST_AUTO_TEST_CASE(add_1)
 {
-    auto a = dynamic_histogram<dynamic_storage>(integer_axis(-1, 1));
-    auto b = dynamic_histogram<static_storage<int>>(integer_axis(-1, 1));
+    auto a = make_dynamic_histogram(integer_axis(-1, 1));
+    auto b = make_dynamic_histogram_with<static_storage<int>>(integer_axis(-1, 1));
     a.fill(-1);
     b.fill(1);
     auto c = a;
@@ -404,8 +404,8 @@ BOOST_AUTO_TEST_CASE(add_1)
 
 BOOST_AUTO_TEST_CASE(add_2)
 {
-    auto a = dynamic_histogram<dynamic_storage>(integer_axis(-1, 1));
-    auto b = dynamic_histogram<dynamic_storage>(integer_axis(-1, 1));
+    auto a = make_dynamic_histogram(integer_axis(-1, 1));
+    auto b = make_dynamic_histogram(integer_axis(-1, 1));
 
     a.fill(0);
     b.wfill(-1, 3);
@@ -426,8 +426,8 @@ BOOST_AUTO_TEST_CASE(add_2)
 
 BOOST_AUTO_TEST_CASE(add_3)
 {
-    auto a = dynamic_histogram<static_storage<char>>(integer_axis(-1, 1));
-    auto b = dynamic_histogram<static_storage<int>>(integer_axis(-1, 1));
+    auto a = make_dynamic_histogram_with<static_storage<char>>(integer_axis(-1, 1));
+    auto b = make_dynamic_histogram_with<static_storage<int>>(integer_axis(-1, 1));
     a.fill(-1);
     b.fill(1);
     auto c = a;
@@ -452,7 +452,7 @@ BOOST_AUTO_TEST_CASE(doc_example_0)
     namespace bh = boost::histogram;
     // create 1d-histogram with 10 equidistant bins from -1.0 to 2.0,
     // with axis of histogram labeled as "x"
-    auto h = bh::dynamic_histogram<>(bh::regular_axis(10, -1.0, 2.0, "x"));
+    auto h = bh::make_dynamic_histogram(bh::regular_axis(10, -1.0, 2.0, "x"));
 
     // fill histogram with data
     h.fill(-1.5); // put in underflow bin
@@ -467,8 +467,8 @@ BOOST_AUTO_TEST_CASE(doc_example_0)
 
     std::ostringstream os1;
     // access histogram counts
-    for (int i = -1; i <= bins(h.axis(0)); ++i) {
-        const auto& a = boost::get<bh::regular_axis>(h.axis(0));
+    const auto& a = boost::get<bh::regular_axis>(h.axis(0));
+    for (int i = -1, n = bins(h.axis(0)) + 1; i < n; ++i) {
         os1 << "bin " << i
             << " x in [" << a[i] << ", " << a[i+1] << "): "
             << h.value(i) << " +/- " << std::sqrt(h.variance(i))
diff --git a/test/dynamic_storage_test.cpp b/test/dynamic_storage_test.cpp
index de312915..82014bc2 100644
--- a/test/dynamic_storage_test.cpp
+++ b/test/dynamic_storage_test.cpp
@@ -10,11 +10,20 @@
 #include <boost/histogram/dynamic_storage.hpp>
 #include <boost/histogram/static_storage.hpp>
 #include <boost/histogram/utility.hpp>
+#include <boost/histogram/detail/utility.hpp>
 #include <sstream>
 #include <string>
 #include <limits>
 using namespace boost::histogram;
 
+template <typename Storage1, typename Storage2>
+bool operator==(const Storage1& a, const Storage2& b)
+{
+    if (a.size() != b.size())
+        return false;
+    return boost::histogram::detail::storage_content_equal(a, b);
+}
+
 BOOST_AUTO_TEST_CASE(equal_operator)
 {
     dynamic_storage a(1), b(1), c(1), d(2);
diff --git a/test/static_histogram_test.cpp b/test/static_histogram_test.cpp
index 0adec298..f7bdc68a 100644
--- a/test/static_histogram_test.cpp
+++ b/test/static_histogram_test.cpp
@@ -16,6 +16,27 @@
 #include <sstream>
 
 using namespace boost::histogram;
+namespace mpl = boost::mpl;
+
+BOOST_AUTO_TEST_CASE(init_0)
+{
+    auto h = static_histogram<
+      dynamic_storage,
+      mpl::vector<integer_axis>
+    >();
+    BOOST_CHECK_EQUAL(h.dim(), 1);
+    BOOST_CHECK_EQUAL(h.size(), 0);
+    auto h2 = static_histogram<
+      static_storage<int>,
+      mpl::vector<integer_axis>
+    >();
+    BOOST_CHECK(h2 == h);
+    auto h3 = static_histogram<
+      dynamic_storage,
+      mpl::vector<regular_axis>
+    >();
+    BOOST_CHECK(!(h3 == h));
+}
 
 BOOST_AUTO_TEST_CASE(init_1)
 {
@@ -77,407 +98,421 @@ BOOST_AUTO_TEST_CASE(init_5)
                                                          category_axis{"A", "B", "C"});
     BOOST_CHECK_EQUAL(h.dim(), 5);
     BOOST_CHECK_EQUAL(h.size(), 900);
-    auto h2 = make_static_histogram_with<dynamic_storage>(regular_axis{3, -1, 1},
-                                                          integer_axis{-1, 1},
-                                                          polar_axis{3},
-                                                          variable_axis{-1, 0, 1},
-                                                          category_axis{"A", "B", "C"});
+    auto h2 = make_static_histogram_with<static_storage<int>>(regular_axis{3, -1, 1},
+                                                              integer_axis{-1, 1},
+                                                              polar_axis{3},
+                                                              variable_axis{-1, 0, 1},
+                                                              category_axis{"A", "B", "C"});
     BOOST_CHECK(h2 == h);
 }
 
-// BOOST_AUTO_TEST_CASE(copy_ctor)
-// {
-//     auto h = histogram<dynamic_storage>(integer_axis(0, 1),
-//                                         integer_axis(0, 2));
-//     h.fill(0, 0);
-//     auto h2 = decltype(h)(h);
-//     BOOST_CHECK(h2 == h);
-//     auto h3 = histogram<static_storage<int>>(h);
-//     BOOST_CHECK(h3 == h);
-// }
+BOOST_AUTO_TEST_CASE(copy_ctor)
+{
+    auto h = make_static_histogram(integer_axis(0, 1),
+                                   integer_axis(0, 2));
+    h.fill(0, 0);
+    auto h2 = decltype(h)(h);
+    BOOST_CHECK(h2 == h);
+    auto h3 = static_histogram<
+      static_storage<int>,
+      mpl::vector<integer_axis, integer_axis>
+    >(h);
+    BOOST_CHECK(h3 == h);
+}
 
-// BOOST_AUTO_TEST_CASE(copy_assign)
-// {
-//     auto h = histogram<dynamic_storage>(integer_axis(0, 1),
-//                                         integer_axis(0, 2));
-//     h.fill(0, 0);
-//     auto h2 = decltype(h)();
-//     BOOST_CHECK(!(h == h2));
-//     h2 = h;
-//     BOOST_CHECK(h == h2);
-//     // test self-assign
-//     h2 = h2;
-//     BOOST_CHECK(h == h2);
-//     auto h3 = histogram<static_storage<int>>();
-//     h3 = h;
-//     BOOST_CHECK(h == h3);
-// }
+BOOST_AUTO_TEST_CASE(copy_assign)
+{
+    auto h = make_static_histogram(integer_axis(0, 1),
+                                   integer_axis(0, 2));
+    h.fill(0, 0);
+    auto h2 = decltype(h)();
+    BOOST_CHECK(!(h == h2));
+    h2 = h;
+    BOOST_CHECK(h == h2);
+    // test self-assign
+    h2 = h2;
+    BOOST_CHECK(h == h2);
+    auto h3 = static_histogram<
+      static_storage<int>,
+      mpl::vector<integer_axis, integer_axis>
+    >();
+    h3 = h;
+    BOOST_CHECK(h == h3);
+}
 
-// BOOST_AUTO_TEST_CASE(move_ctor)
-// {
-//     auto h = histogram<dynamic_storage>(integer_axis(0, 1),
-//                                         integer_axis(0, 2));
-//     h.fill(0, 0);
-//     const auto href = h;
-//     auto h2 = decltype(h)(std::move(h));
-//     BOOST_CHECK_EQUAL(h.dim(), 0);
-//     BOOST_CHECK_EQUAL(h.sum(), 0);
-//     BOOST_CHECK_EQUAL(h.size(), 0);
-//     BOOST_CHECK(h2 == href);
-//     auto h3 = histogram<static_storage<int>>(std::move(h2));
-//     BOOST_CHECK_EQUAL(h2.dim(), 0);
-//     BOOST_CHECK_EQUAL(h2.sum(), 0);
-//     BOOST_CHECK_EQUAL(h2.size(), 0);
-//     BOOST_CHECK(h3 == href);
-// }
+BOOST_AUTO_TEST_CASE(move_ctor)
+{
+    auto h = make_static_histogram(integer_axis(0, 1),
+                                   integer_axis(0, 2));
+    h.fill(0, 0);
+    const auto href = h;
+    auto h2 = decltype(h)(std::move(h));
+    BOOST_CHECK_EQUAL(h.dim(), 2);
+    BOOST_CHECK_EQUAL(h.sum(), 0);
+    BOOST_CHECK_EQUAL(h.size(), 0);
+    BOOST_CHECK(h2 == href);
+    auto h3 = static_histogram<
+      static_storage<int>,
+      mpl::vector<integer_axis, integer_axis>
+    >(std::move(h2));
+    BOOST_CHECK_EQUAL(h2.dim(), 2);
+    BOOST_CHECK_EQUAL(h2.sum(), 0);
+    BOOST_CHECK_EQUAL(h2.size(), 0);
+    BOOST_CHECK(h3 == href);
+}
 
-// BOOST_AUTO_TEST_CASE(move_assign)
-// {
-//     auto h = histogram<dynamic_storage>(integer_axis(0, 1),
-//                                         integer_axis(0, 2));
-//     h.fill(0.0, 0.0);
-//     auto href = h;
-//     auto h2 = decltype(h)();
-//     h2 = std::move(h);
-//     BOOST_CHECK(h2 == href);
-//     BOOST_CHECK_EQUAL(h.dim(), 0);
-//     BOOST_CHECK_EQUAL(h.sum(), 0);
-//     BOOST_CHECK_EQUAL(h.size(), 0);
-//     auto h3 = histogram<static_storage<int>>();
-//     h3 = std::move(h2);
-//     BOOST_CHECK_EQUAL(h2.dim(), 0);
-//     BOOST_CHECK_EQUAL(h2.sum(), 0);
-//     BOOST_CHECK_EQUAL(h2.size(), 0);
-//     BOOST_CHECK(h3 == href);
-// }
+BOOST_AUTO_TEST_CASE(move_assign)
+{
+    auto h = make_static_histogram(integer_axis(0, 1),
+                                   integer_axis(0, 2));
+    h.fill(0.0, 0.0);
+    auto href = h;
+    auto h2 = decltype(h)();
+    h2 = std::move(h);
+    BOOST_CHECK(h2 == href);
+    BOOST_CHECK_EQUAL(h.dim(), 2);
+    BOOST_CHECK_EQUAL(h.sum(), 0);
+    BOOST_CHECK_EQUAL(h.size(), 0);
+    auto h3 = static_histogram<
+      static_storage<int>,
+      mpl::vector<integer_axis, integer_axis>
+    >();
+    h3 = std::move(h2);
+    BOOST_CHECK_EQUAL(h2.dim(), 2);
+    BOOST_CHECK_EQUAL(h2.sum(), 0);
+    BOOST_CHECK_EQUAL(h2.size(), 0);
+    BOOST_CHECK(h3 == href);
+}
 
-// BOOST_AUTO_TEST_CASE(equal_compare)
-// {
-//     using H = histogram<dynamic_storage>;
-//     auto a = H(integer_axis(0, 1));
-//     auto b = H(integer_axis(0, 1), integer_axis(0, 2));
-//     BOOST_CHECK(!(a == b));
-//     BOOST_CHECK(!(b == a));
-//     auto c = H(integer_axis(0, 1));
-//     BOOST_CHECK(!(b == c));
-//     BOOST_CHECK(!(c == b));
-//     BOOST_CHECK(a == c);
-//     BOOST_CHECK(c == a);
-//     auto d = H(regular_axis(2, 0, 1));
-//     BOOST_CHECK(!(c == d));
-//     BOOST_CHECK(!(d == c));
-//     c.fill(0);
-//     BOOST_CHECK(!(a == c));
-//     BOOST_CHECK(!(c == a));
-//     a.fill(0);
-//     BOOST_CHECK(a == c);
-//     BOOST_CHECK(c == a);
-//     a.fill(0);
-//     BOOST_CHECK(!(a == c));
-//     BOOST_CHECK(!(c == a));
-// }
+BOOST_AUTO_TEST_CASE(equal_compare)
+{
+    auto a = make_static_histogram(integer_axis(0, 1));
+    auto b = make_static_histogram(integer_axis(0, 1), integer_axis(0, 2));
+    BOOST_CHECK(!(a == b));
+    BOOST_CHECK(!(b == a));
+    auto c = make_static_histogram(integer_axis(0, 1));
+    BOOST_CHECK(!(b == c));
+    BOOST_CHECK(!(c == b));
+    BOOST_CHECK(a == c);
+    BOOST_CHECK(c == a);
+    auto d = make_static_histogram(regular_axis(2, 0, 1));
+    BOOST_CHECK(!(c == d));
+    BOOST_CHECK(!(d == c));
+    c.fill(0);
+    BOOST_CHECK(!(a == c));
+    BOOST_CHECK(!(c == a));
+    a.fill(0);
+    BOOST_CHECK(a == c);
+    BOOST_CHECK(c == a);
+    a.fill(0);
+    BOOST_CHECK(!(a == c));
+    BOOST_CHECK(!(c == a));
+}
 
-// BOOST_AUTO_TEST_CASE(d1)
-// {
-//     auto h = histogram<dynamic_storage>(integer_axis(0, 1));
-//     h.fill(0);
-//     h.fill(0);
-//     h.fill(-1);
-//     h.fill(10);
+BOOST_AUTO_TEST_CASE(d1)
+{
+    auto h = make_static_histogram(integer_axis(0, 1));
+    h.fill(0);
+    h.fill(0);
+    h.fill(-1);
+    h.fill(10);
 
-//     BOOST_CHECK_EQUAL(h.dim(), 1);
-//     BOOST_CHECK_EQUAL(bins(h.axis(0)), 2);
-//     BOOST_CHECK_EQUAL(shape(h.axis(0)), 4);
-//     BOOST_CHECK_EQUAL(h.sum(), 4);
+    BOOST_CHECK_EQUAL(h.dim(), 1);
+    BOOST_CHECK_EQUAL(bins(h.axis<0>()), 2);
+    BOOST_CHECK_EQUAL(shape(h.axis<0>()), 4);
+    BOOST_CHECK_EQUAL(h.sum(), 4);
 
-//     BOOST_CHECK_THROW(h.value(-2), std::out_of_range);
-//     BOOST_CHECK_EQUAL(h.value(-1), 1.0);
-//     BOOST_CHECK_EQUAL(h.value(0), 2.0);
-//     BOOST_CHECK_EQUAL(h.value(1), 0.0);
-//     BOOST_CHECK_EQUAL(h.value(2), 1.0);
-//     BOOST_CHECK_THROW(h.value(3), std::out_of_range);
+    BOOST_CHECK_THROW(h.value(-2), std::out_of_range);
+    BOOST_CHECK_EQUAL(h.value(-1), 1.0);
+    BOOST_CHECK_EQUAL(h.value(0), 2.0);
+    BOOST_CHECK_EQUAL(h.value(1), 0.0);
+    BOOST_CHECK_EQUAL(h.value(2), 1.0);
+    BOOST_CHECK_THROW(h.value(3), std::out_of_range);
 
-//     BOOST_CHECK_THROW(h.variance(-2), std::out_of_range);
-//     BOOST_CHECK_EQUAL(h.variance(-1), 1.0);
-//     BOOST_CHECK_EQUAL(h.variance(0), 2.0);
-//     BOOST_CHECK_EQUAL(h.variance(1), 0.0);
-//     BOOST_CHECK_EQUAL(h.variance(2), 1.0);
-//     BOOST_CHECK_THROW(h.variance(3), std::out_of_range);
-// }
+    BOOST_CHECK_THROW(h.variance(-2), std::out_of_range);
+    BOOST_CHECK_EQUAL(h.variance(-1), 1.0);
+    BOOST_CHECK_EQUAL(h.variance(0), 2.0);
+    BOOST_CHECK_EQUAL(h.variance(1), 0.0);
+    BOOST_CHECK_EQUAL(h.variance(2), 1.0);
+    BOOST_CHECK_THROW(h.variance(3), std::out_of_range);
+}
 
-// BOOST_AUTO_TEST_CASE(d1_2)
-// {
-//     auto h = histogram<dynamic_storage>(integer_axis(0, 1, "", false));
-//     h.fill(0);
-//     h.fill(-0);
-//     h.fill(-1);
-//     h.fill(10);
+BOOST_AUTO_TEST_CASE(d1_2)
+{
+    auto h = make_static_histogram(integer_axis(0, 1, "", false));
+    h.fill(0);
+    h.fill(-0);
+    h.fill(-1);
+    h.fill(10);
 
-//     BOOST_CHECK_EQUAL(h.dim(), 1);
-//     BOOST_CHECK_EQUAL(bins(h.axis(0)), 2);
-//     BOOST_CHECK_EQUAL(shape(h.axis(0)), 2);
-//     BOOST_CHECK_EQUAL(h.sum(), 2);
+    BOOST_CHECK_EQUAL(h.dim(), 1);
+    BOOST_CHECK_EQUAL(bins(h.axis<0>()), 2);
+    BOOST_CHECK_EQUAL(shape(h.axis<0>()), 2);
+    BOOST_CHECK_EQUAL(h.sum(), 2);
 
-//     BOOST_CHECK_THROW(h.value(-1), std::out_of_range);
-//     BOOST_CHECK_EQUAL(h.value(0), 2.0);
-//     BOOST_CHECK_EQUAL(h.value(1), 0.0);
-//     BOOST_CHECK_THROW(h.value(2), std::out_of_range);
+    BOOST_CHECK_THROW(h.value(-1), std::out_of_range);
+    BOOST_CHECK_EQUAL(h.value(0), 2.0);
+    BOOST_CHECK_EQUAL(h.value(1), 0.0);
+    BOOST_CHECK_THROW(h.value(2), std::out_of_range);
 
-//     BOOST_CHECK_THROW(h.variance(-1), std::out_of_range);
-//     BOOST_CHECK_EQUAL(h.variance(0), 2.0);
-//     BOOST_CHECK_EQUAL(h.variance(1), 0.0);
-//     BOOST_CHECK_THROW(h.variance(2), std::out_of_range);
-// }
+    BOOST_CHECK_THROW(h.variance(-1), std::out_of_range);
+    BOOST_CHECK_EQUAL(h.variance(0), 2.0);
+    BOOST_CHECK_EQUAL(h.variance(1), 0.0);
+    BOOST_CHECK_THROW(h.variance(2), std::out_of_range);
+}
 
-// BOOST_AUTO_TEST_CASE(d1w)
-// {
-//     auto h = histogram<dynamic_storage>(regular_axis(2, -1, 1));
-//     h.fill(0);
-//     h.wfill(-1.0, 2.0);
-//     h.fill(-1.0);
-//     h.fill(-2.0);
-//     h.wfill(10.0, 5.0);
+BOOST_AUTO_TEST_CASE(d1w)
+{
+    auto h = make_static_histogram(regular_axis(2, -1, 1));
+    h.fill(0);
+    h.wfill(-1.0, 2.0);
+    h.fill(-1.0);
+    h.fill(-2.0);
+    h.wfill(10.0, 5.0);
 
-//     BOOST_CHECK_EQUAL(h.sum(), 10);
+    BOOST_CHECK_EQUAL(h.sum(), 10);
 
-//     BOOST_CHECK_EQUAL(h.value(-1), 1.0);
-//     BOOST_CHECK_EQUAL(h.value(0), 3.0);
-//     BOOST_CHECK_EQUAL(h.value(1), 1.0);
-//     BOOST_CHECK_EQUAL(h.value(2), 5.0);
+    BOOST_CHECK_EQUAL(h.value(-1), 1.0);
+    BOOST_CHECK_EQUAL(h.value(0), 3.0);
+    BOOST_CHECK_EQUAL(h.value(1), 1.0);
+    BOOST_CHECK_EQUAL(h.value(2), 5.0);
 
-//     BOOST_CHECK_EQUAL(h.variance(-1), 1.0);
-//     BOOST_CHECK_EQUAL(h.variance(0), 5.0);
-//     BOOST_CHECK_EQUAL(h.variance(1), 1.0);
-//     BOOST_CHECK_EQUAL(h.variance(2), 25.0);
-// }
+    BOOST_CHECK_EQUAL(h.variance(-1), 1.0);
+    BOOST_CHECK_EQUAL(h.variance(0), 5.0);
+    BOOST_CHECK_EQUAL(h.variance(1), 1.0);
+    BOOST_CHECK_EQUAL(h.variance(2), 25.0);
+}
 
-// BOOST_AUTO_TEST_CASE(d2)
-// {
-//     auto h = histogram<dynamic_storage>(regular_axis(2, -1, 1),
-//                                         integer_axis(-1, 1, std::string(), false));
-//     h.fill(-1, -1);
-//     h.fill(-1, 0);
-//     h.fill(-1, -10);
-//     h.fill(-10, 0);
+BOOST_AUTO_TEST_CASE(d2)
+{
+    auto h = make_static_histogram(regular_axis(2, -1, 1),
+                                   integer_axis(-1, 1, std::string(), false));
+    h.fill(-1, -1);
+    h.fill(-1, 0);
+    h.fill(-1, -10);
+    h.fill(-10, 0);
 
-//     BOOST_CHECK_EQUAL(h.dim(), 2);
-//     BOOST_CHECK_EQUAL(bins(h.axis(0)), 2);
-//     BOOST_CHECK_EQUAL(shape(h.axis(0)), 4);
-//     BOOST_CHECK_EQUAL(bins(h.axis(1)), 3);
-//     BOOST_CHECK_EQUAL(shape(h.axis(1)), 3);
-//     BOOST_CHECK_EQUAL(h.sum(), 3);
+    BOOST_CHECK_EQUAL(h.dim(), 2);
+    BOOST_CHECK_EQUAL(bins(h.axis<0>()), 2);
+    BOOST_CHECK_EQUAL(shape(h.axis<0>()), 4);
+    BOOST_CHECK_EQUAL(bins(h.axis<1>()), 3);
+    BOOST_CHECK_EQUAL(shape(h.axis<1>()), 3);
+    BOOST_CHECK_EQUAL(h.sum(), 3);
 
-//     BOOST_CHECK_EQUAL(h.value(-1, 0), 0.0);
-//     BOOST_CHECK_EQUAL(h.value(-1, 1), 1.0);
-//     BOOST_CHECK_EQUAL(h.value(-1, 2), 0.0);
+    BOOST_CHECK_EQUAL(h.value(-1, 0), 0.0);
+    BOOST_CHECK_EQUAL(h.value(-1, 1), 1.0);
+    BOOST_CHECK_EQUAL(h.value(-1, 2), 0.0);
 
-//     BOOST_CHECK_EQUAL(h.value(0, 0), 1.0);
-//     BOOST_CHECK_EQUAL(h.value(0, 1), 1.0);
-//     BOOST_CHECK_EQUAL(h.value(0, 2), 0.0);
+    BOOST_CHECK_EQUAL(h.value(0, 0), 1.0);
+    BOOST_CHECK_EQUAL(h.value(0, 1), 1.0);
+    BOOST_CHECK_EQUAL(h.value(0, 2), 0.0);
 
-//     BOOST_CHECK_EQUAL(h.value(1, 0), 0.0);
-//     BOOST_CHECK_EQUAL(h.value(1, 1), 0.0);
-//     BOOST_CHECK_EQUAL(h.value(1, 2), 0.0);
+    BOOST_CHECK_EQUAL(h.value(1, 0), 0.0);
+    BOOST_CHECK_EQUAL(h.value(1, 1), 0.0);
+    BOOST_CHECK_EQUAL(h.value(1, 2), 0.0);
 
-//     BOOST_CHECK_EQUAL(h.value(2, 0), 0.0);
-//     BOOST_CHECK_EQUAL(h.value(2, 1), 0.0);
-//     BOOST_CHECK_EQUAL(h.value(2, 2), 0.0);
+    BOOST_CHECK_EQUAL(h.value(2, 0), 0.0);
+    BOOST_CHECK_EQUAL(h.value(2, 1), 0.0);
+    BOOST_CHECK_EQUAL(h.value(2, 2), 0.0);
 
-//     BOOST_CHECK_EQUAL(h.variance(-1, 0), 0.0);
-//     BOOST_CHECK_EQUAL(h.variance(-1, 1), 1.0);
-//     BOOST_CHECK_EQUAL(h.variance(-1, 2), 0.0);
+    BOOST_CHECK_EQUAL(h.variance(-1, 0), 0.0);
+    BOOST_CHECK_EQUAL(h.variance(-1, 1), 1.0);
+    BOOST_CHECK_EQUAL(h.variance(-1, 2), 0.0);
 
-//     BOOST_CHECK_EQUAL(h.variance(0, 0), 1.0);
-//     BOOST_CHECK_EQUAL(h.variance(0, 1), 1.0);
-//     BOOST_CHECK_EQUAL(h.variance(0, 2), 0.0);
+    BOOST_CHECK_EQUAL(h.variance(0, 0), 1.0);
+    BOOST_CHECK_EQUAL(h.variance(0, 1), 1.0);
+    BOOST_CHECK_EQUAL(h.variance(0, 2), 0.0);
 
-//     BOOST_CHECK_EQUAL(h.variance(1, 0), 0.0);
-//     BOOST_CHECK_EQUAL(h.variance(1, 1), 0.0);
-//     BOOST_CHECK_EQUAL(h.variance(1, 2), 0.0);
+    BOOST_CHECK_EQUAL(h.variance(1, 0), 0.0);
+    BOOST_CHECK_EQUAL(h.variance(1, 1), 0.0);
+    BOOST_CHECK_EQUAL(h.variance(1, 2), 0.0);
 
-//     BOOST_CHECK_EQUAL(h.variance(2, 0), 0.0);
-//     BOOST_CHECK_EQUAL(h.variance(2, 1), 0.0);
-//     BOOST_CHECK_EQUAL(h.variance(2, 2), 0.0);
-// }
+    BOOST_CHECK_EQUAL(h.variance(2, 0), 0.0);
+    BOOST_CHECK_EQUAL(h.variance(2, 1), 0.0);
+    BOOST_CHECK_EQUAL(h.variance(2, 2), 0.0);
+}
 
-// BOOST_AUTO_TEST_CASE(d2w)
-// {
-//     auto h = histogram<dynamic_storage>(regular_axis(2, -1, 1),
-//                                         integer_axis(-1, 1, std::string(), false));
-//     h.fill(-1, 0);       // -> 0, 1
-//     h.wfill(-1, -1, 10); // -> 0, 0
-//     h.wfill(-1, -10, 5); // is ignored
-//     h.wfill(-10, 0, 7);  // -> -1, 1
+BOOST_AUTO_TEST_CASE(d2w)
+{
+    auto h = make_static_histogram(regular_axis(2, -1, 1),
+                                   integer_axis(-1, 1, std::string(), false));
+    h.fill(-1, 0);       // -> 0, 1
+    h.wfill(-1, -1, 10); // -> 0, 0
+    h.wfill(-1, -10, 5); // is ignored
+    h.wfill(-10, 0, 7);  // -> -1, 1
 
-//     BOOST_CHECK_EQUAL(h.sum(), 18);
+    BOOST_CHECK_EQUAL(h.sum(), 18);
 
-//     BOOST_CHECK_EQUAL(h.value(-1, 0), 0.0);
-//     BOOST_CHECK_EQUAL(h.value(-1, 1), 7.0);
-//     BOOST_CHECK_EQUAL(h.value(-1, 2), 0.0);
+    BOOST_CHECK_EQUAL(h.value(-1, 0), 0.0);
+    BOOST_CHECK_EQUAL(h.value(-1, 1), 7.0);
+    BOOST_CHECK_EQUAL(h.value(-1, 2), 0.0);
 
-//     BOOST_CHECK_EQUAL(h.value(0, 0), 10.0);
-//     BOOST_CHECK_EQUAL(h.value(0, 1), 1.0);
-//     BOOST_CHECK_EQUAL(h.value(0, 2), 0.0);
+    BOOST_CHECK_EQUAL(h.value(0, 0), 10.0);
+    BOOST_CHECK_EQUAL(h.value(0, 1), 1.0);
+    BOOST_CHECK_EQUAL(h.value(0, 2), 0.0);
 
-//     BOOST_CHECK_EQUAL(h.value(1, 0), 0.0);
-//     BOOST_CHECK_EQUAL(h.value(1, 1), 0.0);
-//     BOOST_CHECK_EQUAL(h.value(1, 2), 0.0);
+    BOOST_CHECK_EQUAL(h.value(1, 0), 0.0);
+    BOOST_CHECK_EQUAL(h.value(1, 1), 0.0);
+    BOOST_CHECK_EQUAL(h.value(1, 2), 0.0);
 
-//     BOOST_CHECK_EQUAL(h.value(2, 0), 0.0);
-//     BOOST_CHECK_EQUAL(h.value(2, 1), 0.0);
-//     BOOST_CHECK_EQUAL(h.value(2, 2), 0.0);
+    BOOST_CHECK_EQUAL(h.value(2, 0), 0.0);
+    BOOST_CHECK_EQUAL(h.value(2, 1), 0.0);
+    BOOST_CHECK_EQUAL(h.value(2, 2), 0.0);
 
-//     BOOST_CHECK_EQUAL(h.variance(-1, 0), 0.0);
-//     BOOST_CHECK_EQUAL(h.variance(-1, 1), 49.0);
-//     BOOST_CHECK_EQUAL(h.variance(-1, 2), 0.0);
+    BOOST_CHECK_EQUAL(h.variance(-1, 0), 0.0);
+    BOOST_CHECK_EQUAL(h.variance(-1, 1), 49.0);
+    BOOST_CHECK_EQUAL(h.variance(-1, 2), 0.0);
 
-//     BOOST_CHECK_EQUAL(h.variance(0, 0), 100.0);
-//     BOOST_CHECK_EQUAL(h.variance(0, 1), 1.0);
-//     BOOST_CHECK_EQUAL(h.variance(0, 2), 0.0);
+    BOOST_CHECK_EQUAL(h.variance(0, 0), 100.0);
+    BOOST_CHECK_EQUAL(h.variance(0, 1), 1.0);
+    BOOST_CHECK_EQUAL(h.variance(0, 2), 0.0);
 
-//     BOOST_CHECK_EQUAL(h.variance(1, 0), 0.0);
-//     BOOST_CHECK_EQUAL(h.variance(1, 1), 0.0);
-//     BOOST_CHECK_EQUAL(h.variance(1, 2), 0.0);
+    BOOST_CHECK_EQUAL(h.variance(1, 0), 0.0);
+    BOOST_CHECK_EQUAL(h.variance(1, 1), 0.0);
+    BOOST_CHECK_EQUAL(h.variance(1, 2), 0.0);
 
-//     BOOST_CHECK_EQUAL(h.variance(2, 0), 0.0);
-//     BOOST_CHECK_EQUAL(h.variance(2, 1), 0.0);
-//     BOOST_CHECK_EQUAL(h.variance(2, 2), 0.0);
-// }
+    BOOST_CHECK_EQUAL(h.variance(2, 0), 0.0);
+    BOOST_CHECK_EQUAL(h.variance(2, 1), 0.0);
+    BOOST_CHECK_EQUAL(h.variance(2, 2), 0.0);
+}
 
-// BOOST_AUTO_TEST_CASE(d3w)
-// {
-//     auto h = histogram<dynamic_storage>(integer_axis(0, 3),
-//                                         integer_axis(0, 4),
-//                                         integer_axis(0, 5));
-//     for (auto i = 0; i < bins(h.axis(0)); ++i)
-//         for (auto j = 0; j < bins(h.axis(1)); ++j)
-//             for (auto k = 0; k < bins(h.axis(2)); ++k)
-//     {
-//         h.wfill(i, j, k, i+j+k);
-//     }
+BOOST_AUTO_TEST_CASE(d3w)
+{
+    auto h = make_static_histogram(integer_axis(0, 3),
+                                   integer_axis(0, 4),
+                                   integer_axis(0, 5));
+    for (auto i = 0; i < bins(h.axis<0>()); ++i)
+        for (auto j = 0; j < bins(h.axis<1>()); ++j)
+            for (auto k = 0; k < bins(h.axis<2>()); ++k)
+    {
+        h.wfill(i, j, k, i+j+k);
+    }
 
-//     for (auto i = 0; i < bins(h.axis(0)); ++i)
-//         for (auto j = 0; j < bins(h.axis(1)); ++j)
-//             for (auto k = 0; k < bins(h.axis(2)); ++k)
-//         BOOST_CHECK_EQUAL(h.value(i, j, k), i+j+k);
-// }
+    for (auto i = 0; i < bins(h.axis<0>()); ++i)
+        for (auto j = 0; j < bins(h.axis<1>()); ++j)
+            for (auto k = 0; k < bins(h.axis<2>()); ++k)
+        BOOST_CHECK_EQUAL(h.value(i, j, k), i+j+k);
+}
 
-// BOOST_AUTO_TEST_CASE(add_0)
-// {
-//     auto a = histogram<dynamic_storage>(integer_axis(-1, 1));
-//     auto b = histogram<dynamic_storage>(regular_axis(3, -1, 1));
-//     auto c = histogram<dynamic_storage>(regular_axis(3, -1.1, 1));
-//     BOOST_CHECK_THROW(a + b, std::logic_error);
-//     BOOST_CHECK_THROW(b + c, std::logic_error);
-// }
+BOOST_AUTO_TEST_CASE(add_1)
+{
+    auto a = make_static_histogram_with<
+      dynamic_storage
+    >(integer_axis(-1, 1));
+    auto b = make_static_histogram_with<
+      static_storage<int>
+    >(integer_axis(-1, 1));
+    a.fill(-1);
+    b.fill(1);
+    auto c = a;
+    c += b;
+    BOOST_CHECK_EQUAL(c.value(-1), 0);
+    BOOST_CHECK_EQUAL(c.value(0), 1);
+    BOOST_CHECK_EQUAL(c.value(1), 0);
+    BOOST_CHECK_EQUAL(c.value(2), 1);
+    BOOST_CHECK_EQUAL(c.value(3), 0);
+    auto d = a + b;
+    BOOST_CHECK_EQUAL(d.value(-1), 0);
+    BOOST_CHECK_EQUAL(d.value(0), 1);
+    BOOST_CHECK_EQUAL(d.value(1), 0);
+    BOOST_CHECK_EQUAL(d.value(2), 1);
+    BOOST_CHECK_EQUAL(d.value(3), 0);
+}
 
-// BOOST_AUTO_TEST_CASE(add_1)
-// {
-//     auto a = histogram<dynamic_storage>(integer_axis(-1, 1));
-//     auto b = histogram<static_storage<int>>(integer_axis(-1, 1));
-//     a.fill(-1);
-//     b.fill(1);
-//     auto c = a;
-//     c += b;
-//     BOOST_CHECK_EQUAL(c.value(-1), 0);
-//     BOOST_CHECK_EQUAL(c.value(0), 1);
-//     BOOST_CHECK_EQUAL(c.value(1), 0);
-//     BOOST_CHECK_EQUAL(c.value(2), 1);
-//     BOOST_CHECK_EQUAL(c.value(3), 0);
-//     auto d = a + b;
-//     BOOST_CHECK_EQUAL(d.value(-1), 0);
-//     BOOST_CHECK_EQUAL(d.value(0), 1);
-//     BOOST_CHECK_EQUAL(d.value(1), 0);
-//     BOOST_CHECK_EQUAL(d.value(2), 1);
-//     BOOST_CHECK_EQUAL(d.value(3), 0);
-// }
+BOOST_AUTO_TEST_CASE(add_2)
+{
+    auto a = make_static_histogram_with<
+      dynamic_storage
+    >(integer_axis(-1, 1));
+    auto b = make_static_histogram_with<
+      dynamic_storage
+    >(integer_axis(-1, 1));
 
-// BOOST_AUTO_TEST_CASE(add_2)
-// {
-//     auto a = histogram<dynamic_storage>(integer_axis(-1, 1));
-//     auto b = histogram<dynamic_storage>(integer_axis(-1, 1));
+    a.fill(0);
+    b.wfill(-1, 3);
+    auto c = a;
+    c += b;
+    BOOST_CHECK_EQUAL(c.value(-1), 0);
+    BOOST_CHECK_EQUAL(c.value(0), 3);
+    BOOST_CHECK_EQUAL(c.value(1), 1);
+    BOOST_CHECK_EQUAL(c.value(2), 0);
+    BOOST_CHECK_EQUAL(c.value(3), 0);    
+    auto d = a + b;
+    BOOST_CHECK_EQUAL(d.value(-1), 0);
+    BOOST_CHECK_EQUAL(d.value(0), 3);
+    BOOST_CHECK_EQUAL(d.value(1), 1);
+    BOOST_CHECK_EQUAL(d.value(2), 0);
+    BOOST_CHECK_EQUAL(d.value(3), 0);
+}
 
-//     a.fill(0);
-//     b.wfill(-1, 3);
-//     auto c = a;
-//     c += b;
-//     BOOST_CHECK_EQUAL(c.value(-1), 0);
-//     BOOST_CHECK_EQUAL(c.value(0), 3);
-//     BOOST_CHECK_EQUAL(c.value(1), 1);
-//     BOOST_CHECK_EQUAL(c.value(2), 0);
-//     BOOST_CHECK_EQUAL(c.value(3), 0);    
-//     auto d = a + b;
-//     BOOST_CHECK_EQUAL(d.value(-1), 0);
-//     BOOST_CHECK_EQUAL(d.value(0), 3);
-//     BOOST_CHECK_EQUAL(d.value(1), 1);
-//     BOOST_CHECK_EQUAL(d.value(2), 0);
-//     BOOST_CHECK_EQUAL(d.value(3), 0);
-// }
+BOOST_AUTO_TEST_CASE(add_3)
+{
+    auto a = make_static_histogram_with<
+      static_storage<char>
+    >(integer_axis(-1, 1));
+    auto b = make_static_histogram_with<
+      static_storage<int>
+    >(integer_axis(-1, 1));
+    a.fill(-1);
+    b.fill(1);
+    auto c = a;
+    c += b;
+    BOOST_CHECK_EQUAL(c.depth(), sizeof(char));
+    BOOST_CHECK_EQUAL(c.value(-1), 0);
+    BOOST_CHECK_EQUAL(c.value(0), 1);
+    BOOST_CHECK_EQUAL(c.value(1), 0);
+    BOOST_CHECK_EQUAL(c.value(2), 1);
+    BOOST_CHECK_EQUAL(c.value(3), 0);
+    auto d = a + b;
+    BOOST_CHECK_EQUAL(d.depth(), sizeof(int));
+    BOOST_CHECK_EQUAL(d.value(-1), 0);
+    BOOST_CHECK_EQUAL(d.value(0), 1);
+    BOOST_CHECK_EQUAL(d.value(1), 0);
+    BOOST_CHECK_EQUAL(d.value(2), 1);
+    BOOST_CHECK_EQUAL(d.value(3), 0);
+}
 
-// BOOST_AUTO_TEST_CASE(add_3)
-// {
-//     auto a = histogram<static_storage<char>>(integer_axis(-1, 1));
-//     auto b = histogram<static_storage<int>>(integer_axis(-1, 1));
-//     a.fill(-1);
-//     b.fill(1);
-//     auto c = a;
-//     c += b;
-//     BOOST_CHECK_EQUAL(c.depth(), sizeof(char));
-//     BOOST_CHECK_EQUAL(c.value(-1), 0);
-//     BOOST_CHECK_EQUAL(c.value(0), 1);
-//     BOOST_CHECK_EQUAL(c.value(1), 0);
-//     BOOST_CHECK_EQUAL(c.value(2), 1);
-//     BOOST_CHECK_EQUAL(c.value(3), 0);
-//     auto d = a + b;
-//     BOOST_CHECK_EQUAL(d.depth(), sizeof(int));
-//     BOOST_CHECK_EQUAL(d.value(-1), 0);
-//     BOOST_CHECK_EQUAL(d.value(0), 1);
-//     BOOST_CHECK_EQUAL(d.value(1), 0);
-//     BOOST_CHECK_EQUAL(d.value(2), 1);
-//     BOOST_CHECK_EQUAL(d.value(3), 0);
-// }
+BOOST_AUTO_TEST_CASE(doc_example_0)
+{
+    namespace bh = boost::histogram;
+    // create 1d-histogram with 10 equidistant bins from -1.0 to 2.0,
+    // with axis of histogram labeled as "x"
+    auto h = bh::make_static_histogram(bh::regular_axis(10, -1.0, 2.0, "x"));
 
-// BOOST_AUTO_TEST_CASE(doc_example_0)
-// {
-//     namespace bh = boost::histogram;
-//     // create 1d-histogram with 10 equidistant bins from -1.0 to 2.0,
-//     // with axis of histogram labeled as "x"
-//     auto h = bh::dynamic::histogram<>(bh::regular_axis(10, -1.0, 2.0, "x"));
+    // fill histogram with data
+    h.fill(-1.5); // put in underflow bin
+    h.fill(-1.0); // included in first bin, bin interval is semi-open
+    h.fill(-0.5);
+    h.fill(1.1);
+    h.fill(0.3);
+    h.fill(1.7);
+    h.fill(2.0);  // put in overflow bin, bin interval is semi-open
+    h.fill(20.0); // put in overflow bin
+    h.wfill(0.1, 5.0); // fill with a weighted entry, weight is 5.0
 
-//     // fill histogram with data
-//     h.fill(-1.5); // put in underflow bin
-//     h.fill(-1.0); // included in first bin, bin interval is semi-open
-//     h.fill(-0.5);
-//     h.fill(1.1);
-//     h.fill(0.3);
-//     h.fill(1.7);
-//     h.fill(2.0);  // put in overflow bin, bin interval is semi-open
-//     h.fill(20.0); // put in overflow bin
-//     h.wfill(0.1, 5.0); // fill with a weighted entry, weight is 5.0
+    std::ostringstream os1;
+    // access histogram counts
+    const auto& a = h.axis<0>();
+    for (int i = -1, n = bins(a) + 1; i < n; ++i) {
+        os1 << "bin " << i
+            << " x in [" << a[i] << ", " << a[i+1] << "): "
+            << h.value(i) << " +/- " << std::sqrt(h.variance(i))
+            << "\n";
+    }
 
-//     std::ostringstream os1;
-//     // access histogram counts
-//     for (int i = -1; i <= bins(h.axis(0)); ++i) {
-//         const auto& a = boost::get<bh::regular_axis>(h.axis(0));
-//         os1 << "bin " << i
-//             << " x in [" << a[i] << ", " << a[i+1] << "): "
-//             << h.value(i) << " +/- " << std::sqrt(h.variance(i))
-//             << "\n";
-//     }
+    std::ostringstream os2;
+    os2 << "bin -1 x in [-inf, -1): 1 +/- 1\n"
+           "bin 0 x in [-1, -0.7): 1 +/- 1\n"
+           "bin 1 x in [-0.7, -0.4): 1 +/- 1\n"
+           "bin 2 x in [-0.4, -0.1): 0 +/- 0\n"
+           "bin 3 x in [-0.1, 0.2): 5 +/- 5\n"
+           "bin 4 x in [0.2, 0.5): 1 +/- 1\n"
+           "bin 5 x in [0.5, 0.8): 0 +/- 0\n"
+           "bin 6 x in [0.8, 1.1): 0 +/- 0\n"
+           "bin 7 x in [1.1, 1.4): 1 +/- 1\n"
+           "bin 8 x in [1.4, 1.7): 0 +/- 0\n"
+           "bin 9 x in [1.7, 2): 1 +/- 1\n"
+           "bin 10 x in [2, inf): 2 +/- 1.41421\n";
 
-//     std::ostringstream os2;
-//     os2 << "bin -1 x in [-inf, -1): 1 +/- 1\n"
-//            "bin 0 x in [-1, -0.7): 1 +/- 1\n"
-//            "bin 1 x in [-0.7, -0.4): 1 +/- 1\n"
-//            "bin 2 x in [-0.4, -0.1): 0 +/- 0\n"
-//            "bin 3 x in [-0.1, 0.2): 5 +/- 5\n"
-//            "bin 4 x in [0.2, 0.5): 1 +/- 1\n"
-//            "bin 5 x in [0.5, 0.8): 0 +/- 0\n"
-//            "bin 6 x in [0.8, 1.1): 0 +/- 0\n"
-//            "bin 7 x in [1.1, 1.4): 1 +/- 1\n"
-//            "bin 8 x in [1.4, 1.7): 0 +/- 0\n"
-//            "bin 9 x in [1.7, 2): 1 +/- 1\n"
-//            "bin 10 x in [2, inf): 2 +/- 1.41421\n";
-
-//     BOOST_CHECK_EQUAL(os1.str(), os2.str());
-// }
+    BOOST_CHECK_EQUAL(os1.str(), os2.str());
+}