From 016532f011605dec459d65c4fe309e31606dfc1c Mon Sep 17 00:00:00 2001
From: Hans Dembinski <HDembinski@users.noreply.github.com>
Date: Mon, 19 Aug 2019 23:55:43 +0200
Subject: [PATCH] Histogram fill method (#54)

huge refactor, various speed improvements, potential for further improvements and parallelization
---
 benchmark/CMakeLists.txt                      |   4 +-
 benchmark/generator.hpp                       |  11 +-
 benchmark/histogram_filling.cpp               |  86 ++++
 benchmark/histogram_filling_experiments.cpp   | 112 +----
 benchmark/histogram_filling_gsl.cpp           |   2 +
 benchmark/histogram_filling_root.cpp          |   4 +
 doc/Jamfile                                   |   2 +-
 examples/guide_histogram_streaming.cpp        |  37 +-
 include/boost/histogram/detail/args_type.hpp  |   9 +-
 include/boost/histogram/detail/at.hpp         |  65 +++
 include/boost/histogram/detail/attribute.hpp  |  16 -
 include/boost/histogram/detail/axes.hpp       | 122 +++--
 include/boost/histogram/detail/data.hpp       |  48 ++
 include/boost/histogram/detail/detect.hpp     |  96 ++--
 include/boost/histogram/detail/fill.hpp       | 345 ++++++++++++++
 include/boost/histogram/detail/fill_n.hpp     | 291 ++++++++++++
 include/boost/histogram/detail/linearize.hpp  | 421 +++---------------
 .../detail/non_member_container_access.hpp    |  70 +++
 .../boost/histogram/detail/optional_index.hpp |  47 ++
 include/boost/histogram/detail/span.hpp       | 249 +++++++++++
 include/boost/histogram/fwd.hpp               |  14 +-
 include/boost/histogram/histogram.hpp         | 176 ++++++--
 include/boost/histogram/indexed.hpp           |  22 +-
 include/boost/histogram/ostream.hpp           |  11 +-
 test/CMakeLists.txt                           |  19 +-
 test/Jamfile                                  |   5 +-
 test/detail_axes_test.cpp                     |  25 +-
 test/detail_detect_test.cpp                   |   2 +-
 test/detail_linearize_test.cpp                |  64 ---
 test/detail_meta_test.cpp                     | 367 ---------------
 test/detail_misc_test.cpp                     |  57 ++-
 test/detail_static_if_test.cpp                |  33 ++
 test/histogram_custom_axis_test.cpp           | 132 ++++++
 test/histogram_fill_test.cpp                  | 299 +++++++++++++
 test/histogram_growing_test.cpp               |   2 +-
 test/histogram_operators_test.cpp             |  19 +-
 test/histogram_test.cpp                       | 160 ++-----
 test/unlimited_storage_test.cpp               |   3 +
 tools/check_build_system.py                   |   1 +
 39 files changed, 2216 insertions(+), 1232 deletions(-)
 create mode 100644 include/boost/histogram/detail/at.hpp
 delete mode 100644 include/boost/histogram/detail/attribute.hpp
 create mode 100644 include/boost/histogram/detail/data.hpp
 create mode 100644 include/boost/histogram/detail/fill.hpp
 create mode 100644 include/boost/histogram/detail/fill_n.hpp
 create mode 100644 include/boost/histogram/detail/non_member_container_access.hpp
 create mode 100644 include/boost/histogram/detail/optional_index.hpp
 create mode 100644 include/boost/histogram/detail/span.hpp
 delete mode 100644 test/detail_linearize_test.cpp
 delete mode 100644 test/detail_meta_test.cpp
 create mode 100644 test/detail_static_if_test.cpp
 create mode 100644 test/histogram_custom_axis_test.cpp
 create mode 100644 test/histogram_fill_test.cpp

diff --git a/benchmark/CMakeLists.txt b/benchmark/CMakeLists.txt
index 1f227342..77d7c9ee 100644
--- a/benchmark/CMakeLists.txt
+++ b/benchmark/CMakeLists.txt
@@ -47,8 +47,8 @@ boost_fetch(hdembinski/benchmark)
 
 macro(add_benchmark name)
   add_executable(${name} "${name}.cpp")
-  target_compile_options(${name}
-    PRIVATE -DNDEBUG -O3 ${BENCHMARK_FLAGS} -funsafe-math-optimizations)
+  target_compile_options(${name} PRIVATE
+    -DNDEBUG -O3 -march=native ${BENCHMARK_FLAGS} -funsafe-math-optimizations)
   target_link_libraries(${name} PRIVATE Boost::histogram benchmark_main)
 endmacro()
 
diff --git a/benchmark/generator.hpp b/benchmark/generator.hpp
index e45f19d5..f613ced3 100644
--- a/benchmark/generator.hpp
+++ b/benchmark/generator.hpp
@@ -30,20 +30,21 @@ uniform_int init<uniform_int, int>(int n) {
 }
 
 template <class Distribution, std::size_t N = 1 << 15>
-struct generator {
+struct generator : std::array<double, N> {
+  using base_t = std::array<double, N>;
+
   template <class... Ts>
   generator(Ts... ts) {
     std::default_random_engine rng(1);
     auto dis = init<Distribution>(ts...);
-    std::generate(buffer_, buffer_ + N, [&] { return dis(rng); });
+    std::generate(base_t::begin(), base_t::end(), [&] { return dis(rng); });
   }
 
   const double& operator()() {
     ++ptr_;
-    if (ptr_ == buffer_ + N) ptr_ = buffer_;
+    if (ptr_ == base_t::data() + N) ptr_ = base_t::data();
     return *ptr_;
   }
 
-  double buffer_[N];
-  const double* ptr_ = buffer_ - 1;
+  const double* ptr_ = base_t::data() - 1;
 };
diff --git a/benchmark/histogram_filling.cpp b/benchmark/histogram_filling.cpp
index 8259f23f..406b32ba 100644
--- a/benchmark/histogram_filling.cpp
+++ b/benchmark/histogram_filling.cpp
@@ -38,6 +38,15 @@ static void fill_1d(benchmark::State& state) {
   auto h = make_s(Tag(), Storage(), reg(100, 0, 1));
   auto gen = generator<Distribution>();
   for (auto _ : state) benchmark::DoNotOptimize(h(gen()));
+  state.SetItemsProcessed(state.iterations());
+}
+
+template <class Distribution, class Tag, class Storage>
+static void fill_n_1d(benchmark::State& state) {
+  auto h = make_s(Tag(), Storage(), reg(100, 0, 1));
+  auto gen = generator<Distribution>();
+  for (auto _ : state) h.fill(gen);
+  state.SetItemsProcessed(state.iterations() * gen.size());
 }
 
 template <class Distribution, class Tag, class Storage>
@@ -45,6 +54,16 @@ static void fill_2d(benchmark::State& state) {
   auto h = make_s(Tag(), Storage(), reg(100, 0, 1), reg(100, 0, 1));
   auto gen = generator<Distribution>();
   for (auto _ : state) benchmark::DoNotOptimize(h(gen(), gen()));
+  state.SetItemsProcessed(state.iterations() * 2);
+}
+
+template <class Distribution, class Tag, class Storage>
+static void fill_n_2d(benchmark::State& state) {
+  auto h = make_s(Tag(), Storage(), reg(100, 0, 1), reg(100, 0, 1));
+  auto gen = generator<Distribution>();
+  auto v = {gen, gen};
+  for (auto _ : state) h.fill(v);
+  state.SetItemsProcessed(state.iterations() * 2 * gen.size());
 }
 
 template <class Distribution, class Tag, class Storage>
@@ -52,6 +71,16 @@ static void fill_3d(benchmark::State& state) {
   auto h = make_s(Tag(), Storage(), reg(100, 0, 1), reg(100, 0, 1), reg(100, 0, 1));
   auto gen = generator<Distribution>();
   for (auto _ : state) benchmark::DoNotOptimize(h(gen(), gen(), gen()));
+  state.SetItemsProcessed(state.iterations() * 3);
+}
+
+template <class Distribution, class Tag, class Storage>
+static void fill_n_3d(benchmark::State& state) {
+  auto h = make_s(Tag(), Storage(), reg(100, 0, 1), reg(100, 0, 1), reg(100, 0, 1));
+  auto gen = generator<Distribution>();
+  auto v = {gen, gen, gen};
+  for (auto _ : state) h.fill(v);
+  state.SetItemsProcessed(state.iterations() * 3 * gen.size());
 }
 
 template <class Distribution, class Tag, class Storage>
@@ -61,38 +90,95 @@ static void fill_6d(benchmark::State& state) {
   auto gen = generator<Distribution>();
   for (auto _ : state)
     benchmark::DoNotOptimize(h(gen(), gen(), gen(), gen(), gen(), gen()));
+  state.SetItemsProcessed(state.iterations() * 6);
+}
+
+template <class Distribution, class Tag, class Storage>
+static void fill_n_6d(benchmark::State& state) {
+  auto h = make_s(Tag(), Storage(), reg(10, 0, 1), reg(10, 0, 1), reg(10, 0, 1),
+                  reg(10, 0, 1), reg(10, 0, 1), reg(10, 0, 1));
+  auto gen = generator<Distribution>();
+  auto v = {gen, gen, gen, gen, gen, gen};
+  for (auto _ : state) h.fill(v);
+  state.SetItemsProcessed(state.iterations() * 6 * gen.size());
 }
 
 BENCHMARK_TEMPLATE(fill_1d, uniform, static_tag, SStore);
 BENCHMARK_TEMPLATE(fill_1d, uniform, static_tag, DStore);
 BENCHMARK_TEMPLATE(fill_1d, uniform, dynamic_tag, SStore);
 BENCHMARK_TEMPLATE(fill_1d, uniform, dynamic_tag, DStore);
+
+BENCHMARK_TEMPLATE(fill_n_1d, uniform, static_tag, SStore);
+BENCHMARK_TEMPLATE(fill_n_1d, uniform, static_tag, DStore);
+BENCHMARK_TEMPLATE(fill_n_1d, uniform, dynamic_tag, SStore);
+BENCHMARK_TEMPLATE(fill_n_1d, uniform, dynamic_tag, DStore);
+
 BENCHMARK_TEMPLATE(fill_2d, uniform, static_tag, SStore);
 BENCHMARK_TEMPLATE(fill_2d, uniform, static_tag, DStore);
 BENCHMARK_TEMPLATE(fill_2d, uniform, dynamic_tag, SStore);
 BENCHMARK_TEMPLATE(fill_2d, uniform, dynamic_tag, DStore);
+
+BENCHMARK_TEMPLATE(fill_n_2d, uniform, static_tag, SStore);
+BENCHMARK_TEMPLATE(fill_n_2d, uniform, static_tag, DStore);
+BENCHMARK_TEMPLATE(fill_n_2d, uniform, dynamic_tag, SStore);
+BENCHMARK_TEMPLATE(fill_n_2d, uniform, dynamic_tag, DStore);
+
 BENCHMARK_TEMPLATE(fill_3d, uniform, static_tag, SStore);
 BENCHMARK_TEMPLATE(fill_3d, uniform, static_tag, DStore);
 BENCHMARK_TEMPLATE(fill_3d, uniform, dynamic_tag, SStore);
 BENCHMARK_TEMPLATE(fill_3d, uniform, dynamic_tag, DStore);
+
+BENCHMARK_TEMPLATE(fill_n_3d, uniform, static_tag, SStore);
+BENCHMARK_TEMPLATE(fill_n_3d, uniform, static_tag, DStore);
+BENCHMARK_TEMPLATE(fill_n_3d, uniform, dynamic_tag, SStore);
+BENCHMARK_TEMPLATE(fill_n_3d, uniform, dynamic_tag, DStore);
+
 BENCHMARK_TEMPLATE(fill_6d, uniform, static_tag, SStore);
 BENCHMARK_TEMPLATE(fill_6d, uniform, static_tag, DStore);
 BENCHMARK_TEMPLATE(fill_6d, uniform, dynamic_tag, SStore);
 BENCHMARK_TEMPLATE(fill_6d, uniform, dynamic_tag, DStore);
 
+BENCHMARK_TEMPLATE(fill_n_6d, uniform, static_tag, SStore);
+BENCHMARK_TEMPLATE(fill_n_6d, uniform, static_tag, DStore);
+BENCHMARK_TEMPLATE(fill_n_6d, uniform, dynamic_tag, SStore);
+BENCHMARK_TEMPLATE(fill_n_6d, uniform, dynamic_tag, DStore);
+
 BENCHMARK_TEMPLATE(fill_1d, normal, static_tag, SStore);
 BENCHMARK_TEMPLATE(fill_1d, normal, static_tag, DStore);
 BENCHMARK_TEMPLATE(fill_1d, normal, dynamic_tag, SStore);
 BENCHMARK_TEMPLATE(fill_1d, normal, dynamic_tag, DStore);
+
+BENCHMARK_TEMPLATE(fill_n_1d, normal, static_tag, SStore);
+BENCHMARK_TEMPLATE(fill_n_1d, normal, static_tag, DStore);
+BENCHMARK_TEMPLATE(fill_n_1d, normal, dynamic_tag, SStore);
+BENCHMARK_TEMPLATE(fill_n_1d, normal, dynamic_tag, DStore);
+
 BENCHMARK_TEMPLATE(fill_2d, normal, static_tag, SStore);
 BENCHMARK_TEMPLATE(fill_2d, normal, static_tag, DStore);
 BENCHMARK_TEMPLATE(fill_2d, normal, dynamic_tag, SStore);
 BENCHMARK_TEMPLATE(fill_2d, normal, dynamic_tag, DStore);
+
+BENCHMARK_TEMPLATE(fill_n_2d, normal, static_tag, SStore);
+BENCHMARK_TEMPLATE(fill_n_2d, normal, static_tag, DStore);
+BENCHMARK_TEMPLATE(fill_n_2d, normal, dynamic_tag, SStore);
+BENCHMARK_TEMPLATE(fill_n_2d, normal, dynamic_tag, DStore);
+
 BENCHMARK_TEMPLATE(fill_3d, normal, static_tag, SStore);
 BENCHMARK_TEMPLATE(fill_3d, normal, static_tag, DStore);
 BENCHMARK_TEMPLATE(fill_3d, normal, dynamic_tag, SStore);
 BENCHMARK_TEMPLATE(fill_3d, normal, dynamic_tag, DStore);
+
+BENCHMARK_TEMPLATE(fill_n_3d, normal, static_tag, SStore);
+BENCHMARK_TEMPLATE(fill_n_3d, normal, static_tag, DStore);
+BENCHMARK_TEMPLATE(fill_n_3d, normal, dynamic_tag, SStore);
+BENCHMARK_TEMPLATE(fill_n_3d, normal, dynamic_tag, DStore);
+
 BENCHMARK_TEMPLATE(fill_6d, normal, static_tag, SStore);
 BENCHMARK_TEMPLATE(fill_6d, normal, static_tag, DStore);
 BENCHMARK_TEMPLATE(fill_6d, normal, dynamic_tag, SStore);
 BENCHMARK_TEMPLATE(fill_6d, normal, dynamic_tag, DStore);
+
+BENCHMARK_TEMPLATE(fill_n_6d, normal, static_tag, SStore);
+BENCHMARK_TEMPLATE(fill_n_6d, normal, static_tag, DStore);
+BENCHMARK_TEMPLATE(fill_n_6d, normal, dynamic_tag, SStore);
+BENCHMARK_TEMPLATE(fill_n_6d, normal, dynamic_tag, DStore);
diff --git a/benchmark/histogram_filling_experiments.cpp b/benchmark/histogram_filling_experiments.cpp
index 8e76b1c8..6647404f 100644
--- a/benchmark/histogram_filling_experiments.cpp
+++ b/benchmark/histogram_filling_experiments.cpp
@@ -33,16 +33,8 @@ auto make_storage(const U& axes) {
   return std::vector<T>(detail::bincount(axes), 0);
 }
 
-template <class T>
-auto make_strides(const T& axes) {
-  std::vector<std::size_t> strides(detail::axes_rank(axes) + 1, 1);
-  auto sit = strides.begin();
-  detail::for_each_axis(axes, [&](const auto& a) { *++sit *= axis::traits::extent(a); });
-  return strides;
-}
-
 template <class Axes, class Storage, class Tuple>
-void fill_b(const Axes& axes, Storage& storage, const Tuple& t) {
+void fill_a(const Axes& axes, Storage& storage, const Tuple& t) {
   using namespace boost::mp11;
   std::size_t stride = 1, index = 0;
   mp_for_each<mp_iota<mp_size<Tuple>>>([&](auto i) {
@@ -54,60 +46,22 @@ void fill_b(const Axes& axes, Storage& storage, const Tuple& t) {
   ++storage[index];
 }
 
-template <class Axes, class Storage, class Tuple>
-void fill_c(const Axes& axes, const std::size_t* strides, Storage& storage,
-            const Tuple& t) {
-  using namespace boost::mp11;
-  std::size_t index = 0;
-  BOOST_ASSERT(boost::histogram::detail::axes_rank(axes) ==
-               boost::histogram::detail::axes_rank(t));
-  mp_for_each<mp_iota<mp_size<Tuple>>>([&](auto i) {
-    const auto& a = boost::histogram::detail::axis_get<i>(axes);
-    const auto& v = std::get<i>(t);
-    index += (a.index(v) + 1) * *strides++;
-  });
-  ++storage[index];
-}
-
 template <class T, class Distribution>
 static void fill_1d_a(benchmark::State& state) {
   auto axes = std::make_tuple(reg(100, 0, 1));
   generator<Distribution> gen;
   auto storage = make_storage<T>(axes);
-  for (auto _ : state) {
-    const auto i = std::get<0>(axes).index(gen());
-    ++storage[i + 1];
-  }
+  for (auto _ : state) fill_a(axes, storage, std::forward_as_tuple(gen()));
+  state.SetItemsProcessed(state.iterations());
 }
 
 template <class T, class Distribution>
-static void fill_1d_b(benchmark::State& state) {
-  auto axes = std::make_tuple(reg(100, 0, 1));
-  generator<Distribution> gen;
-  auto storage = make_storage<T>(axes);
-  for (auto _ : state) { fill_b(axes, storage, std::forward_as_tuple(gen())); }
-}
-
-template <class T, class Distribution>
-static void fill_1d_c(benchmark::State& state) {
-  auto axes = std::make_tuple(reg(100, 0, 1));
-  generator<Distribution> gen;
-  auto storage = make_storage<T>(axes);
-  auto strides = make_strides(axes);
-  for (auto _ : state) {
-    fill_c(axes, strides.data(), storage, std::forward_as_tuple(gen()));
-  }
-}
-
-template <class T, class Distribution>
-static void fill_1d_c_dyn(benchmark::State& state) {
+static void fill_1d_a_dyn(benchmark::State& state) {
   auto axes = vector_of_variant({reg(100, 0, 1)});
   generator<Distribution> gen;
   auto storage = make_storage<T>(axes);
-  auto strides = make_strides(axes);
-  for (auto _ : state) {
-    fill_c(axes, strides.data(), storage, std::forward_as_tuple(gen()));
-  }
+  for (auto _ : state) fill_a(axes, storage, std::forward_as_tuple(gen()));
+  state.SetItemsProcessed(state.iterations());
 }
 
 template <class T, class Distribution>
@@ -115,63 +69,23 @@ static void fill_2d_a(benchmark::State& state) {
   auto axes = std::make_tuple(reg(100, 0, 1), reg(100, 0, 1));
   generator<Distribution> gen;
   auto storage = make_storage<T>(axes);
-  for (auto _ : state) {
-    const auto i0 = std::get<0>(axes).index(gen());
-    const auto i1 = std::get<1>(axes).index(gen());
-    const auto stride = axis::traits::extent(std::get<0>(axes));
-    ++storage[(i0 + 1) * stride + (i1 + 1)];
-  }
+  for (auto _ : state) fill_a(axes, storage, std::forward_as_tuple(gen(), gen()));
+  state.SetItemsProcessed(state.iterations() * 2);
 }
 
 template <class T, class Distribution>
-static void fill_2d_b(benchmark::State& state) {
-  auto axes = std::make_tuple(reg(100, 0, 1), reg(100, 0, 1));
-  generator<Distribution> gen;
-  auto storage = make_storage<T>(axes);
-  for (auto _ : state) { fill_b(axes, storage, std::forward_as_tuple(gen(), gen())); }
-}
-
-template <class T, class Distribution>
-static void fill_2d_c(benchmark::State& state) {
-  auto axes = std::make_tuple(reg(100, 0, 1), reg(100, 0, 1));
-  generator<Distribution> gen;
-  auto storage = make_storage<T>(axes);
-  auto strides = make_strides(axes);
-  BOOST_ASSERT(strides.size() == 3);
-  BOOST_ASSERT(strides[0] == 1);
-  BOOST_ASSERT(strides[1] == 102);
-  for (auto _ : state) {
-    fill_c(axes, strides.data(), storage, std::forward_as_tuple(gen(), gen()));
-  }
-}
-
-template <class T, class Distribution>
-static void fill_2d_c_dyn(benchmark::State& state) {
+static void fill_2d_a_dyn(benchmark::State& state) {
   auto axes = vector_of_variant({reg(100, 0, 1), reg(100, 0, 1)});
   generator<Distribution> gen;
   auto storage = make_storage<T>(axes);
-  auto strides = make_strides(axes);
-  BOOST_ASSERT(strides.size() == 3);
-  BOOST_ASSERT(strides[0] == 1);
-  BOOST_ASSERT(strides[1] == 102);
-  for (auto _ : state) {
-    fill_c(axes, strides.data(), storage, std::forward_as_tuple(gen(), gen()));
-  }
+  for (auto _ : state) fill_a(axes, storage, std::forward_as_tuple(gen(), gen()));
+  state.SetItemsProcessed(state.iterations() * 2);
 }
 
-BENCHMARK_TEMPLATE(fill_1d_a, int, uniform);
 BENCHMARK_TEMPLATE(fill_1d_a, double, uniform);
-BENCHMARK_TEMPLATE(fill_1d_b, double, uniform);
-BENCHMARK_TEMPLATE(fill_1d_c, double, uniform);
-BENCHMARK_TEMPLATE(fill_1d_c_dyn, double, uniform);
+BENCHMARK_TEMPLATE(fill_1d_a_dyn, double, uniform);
 BENCHMARK_TEMPLATE(fill_2d_a, double, uniform);
-BENCHMARK_TEMPLATE(fill_2d_b, double, uniform);
-BENCHMARK_TEMPLATE(fill_2d_c, double, uniform);
-BENCHMARK_TEMPLATE(fill_2d_c_dyn, double, uniform);
+BENCHMARK_TEMPLATE(fill_2d_a_dyn, double, uniform);
 
 BENCHMARK_TEMPLATE(fill_1d_a, double, normal);
-BENCHMARK_TEMPLATE(fill_1d_b, double, normal);
-BENCHMARK_TEMPLATE(fill_1d_c, double, normal);
 BENCHMARK_TEMPLATE(fill_2d_a, double, normal);
-BENCHMARK_TEMPLATE(fill_2d_b, double, normal);
-BENCHMARK_TEMPLATE(fill_2d_c, double, normal);
diff --git a/benchmark/histogram_filling_gsl.cpp b/benchmark/histogram_filling_gsl.cpp
index 132d6b24..94d73896 100644
--- a/benchmark/histogram_filling_gsl.cpp
+++ b/benchmark/histogram_filling_gsl.cpp
@@ -24,6 +24,7 @@ static void fill_1d(benchmark::State& state) {
   generator<Distribution> gen;
   for (auto _ : state) benchmark::DoNotOptimize(gsl_histogram_increment(h, gen()));
   gsl_histogram_free(h);
+  state.SetItemsProcessed(state.iterations());
 }
 
 template <class Distribution>
@@ -34,6 +35,7 @@ static void fill_2d(benchmark::State& state) {
   for (auto _ : state)
     benchmark::DoNotOptimize(gsl_histogram2d_increment(h, gen(), gen()));
   gsl_histogram2d_free(h);
+  state.SetItemsProcessed(state.iterations() * 2);
 }
 
 BENCHMARK_TEMPLATE(fill_1d, uniform);
diff --git a/benchmark/histogram_filling_root.cpp b/benchmark/histogram_filling_root.cpp
index 05395376..98d63851 100644
--- a/benchmark/histogram_filling_root.cpp
+++ b/benchmark/histogram_filling_root.cpp
@@ -23,6 +23,7 @@ static void fill_1d(benchmark::State& state) {
   TH1I h("", "", 100, 0, 1);
   generator<Distribution> gen;
   for (auto _ : state) benchmark::DoNotOptimize(h.Fill(gen()));
+  state.SetItemsProcessed(state.iterations());
 }
 
 template <class Distribution>
@@ -30,6 +31,7 @@ static void fill_2d(benchmark::State& state) {
   TH2I h("", "", 100, 0, 1, 100, 0, 1);
   generator<Distribution> gen;
   for (auto _ : state) benchmark::DoNotOptimize(h.Fill(gen(), gen()));
+  state.SetItemsProcessed(state.iterations() * 2);
 }
 
 template <class Distribution>
@@ -37,6 +39,7 @@ static void fill_3d(benchmark::State& state) {
   TH3I h("", "", 100, 0, 1, 100, 0, 1, 100, 0, 1);
   generator<Distribution> gen;
   for (auto _ : state) benchmark::DoNotOptimize(h.Fill(gen(), gen(), gen()));
+  state.SetItemsProcessed(state.iterations() * 3);
 }
 
 template <class Distribution>
@@ -50,6 +53,7 @@ static void fill_6d(benchmark::State& state) {
     const double buf[6] = {gen(), gen(), gen(), gen(), gen(), gen()};
     benchmark::DoNotOptimize(h.Fill(buf));
   }
+  state.SetItemsProcessed(state.iterations() * 6);
 }
 
 BENCHMARK_TEMPLATE(fill_1d, uniform);
diff --git a/doc/Jamfile b/doc/Jamfile
index a1ab2604..61278276 100644
--- a/doc/Jamfile
+++ b/doc/Jamfile
@@ -30,7 +30,7 @@ doxygen reference
   <doxygen:param>JAVADOC_AUTOBRIEF=YES
   <doxygen:param>EXCLUDE_SYMBOLS=detail
   <doxygen:param>"PREDEFINED=\"BOOST_HISTOGRAM_DOXYGEN_INVOKED\" \\
-                             \"BOOST_HISTOGRAM_NODISCARD\""
+                             \"BOOST_ATTRIBUTE_NODISCARD\""
 ;
 
 make reference_pp.xml : reference.xml : @doxygen_postprocessing ;
diff --git a/examples/guide_histogram_streaming.cpp b/examples/guide_histogram_streaming.cpp
index a9f56469..eae31865 100644
--- a/examples/guide_histogram_streaming.cpp
+++ b/examples/guide_histogram_streaming.cpp
@@ -17,19 +17,8 @@ int main() {
   using namespace boost::histogram;
   namespace tr = axis::transform;
 
-  auto h =
-      make_histogram(axis::regular<>(2, -1.0, 1.0),
-                     axis::regular<double, tr::log>(2, 1.0, 10.0, "axis 1"),
-                     axis::regular<double, tr::pow, use_default, axis::option::growth_t>(
-                         tr::pow{1.5}, 2, 1.0, 10.0, "axis 2"),
-                     // axis without metadata
-                     axis::circular<double, axis::null_type>(4, 0.0, 360.0),
-                     // axis without under-/overflow bins
-                     axis::variable<double, use_default, axis::option::none_t>(
-                         {-1.0, 0.0, 1.0}, "axis 4"),
-                     axis::category<>({2, 1, 3}, "axis 5"),
-                     axis::category<std::string>({"red", "blue"}, "axis 6"),
-                     axis::integer<>(-1, 1, "axis 7"));
+  auto h = make_histogram(axis::regular<>(2, -1.0, 1.0, "axis 1"),
+                          axis::category<std::string>({"red", "blue"}, "axis 2"));
 
   std::ostringstream os;
   os << h;
@@ -38,14 +27,20 @@ int main() {
 
   assert(os.str() ==
          "histogram(\n"
-         "  regular(2, -1, 1, options=underflow | overflow),\n"
-         "  regular_log(2, 1, 10, metadata=\"axis 1\", options=underflow | overflow),\n"
-         "  regular_pow(2, 1, 10, metadata=\"axis 2\", options=growth, power=1.5),\n"
-         "  regular(4, 0, 360, options=overflow | circular),\n"
-         "  variable(-1, 0, 1, metadata=\"axis 4\", options=none),\n"
-         "  category(2, 1, 3, metadata=\"axis 5\", options=overflow),\n"
-         "  category(\"red\", \"blue\", metadata=\"axis 6\", options=overflow),\n"
-         "  integer(-1, 1, metadata=\"axis 7\", options=underflow | overflow)\n"
+         "  regular(2, -1, 1, metadata=\"axis 1\", options=underflow | overflow),\n"
+         "  category(\"red\", \"blue\", metadata=\"axis 2\", options=overflow),\n"
+         "  0: 0\n"
+         "  1: 0\n"
+         "  2: 0\n"
+         "  3: 0\n"
+         "  4: 0\n"
+         "  5: 0\n"
+         "  6: 0\n"
+         "  7: 0\n"
+         "  8: 0\n"
+         "  9: 0\n"
+         "  10: 0\n"
+         "  11: 0\n"
          ")");
 }
 
diff --git a/include/boost/histogram/detail/args_type.hpp b/include/boost/histogram/detail/args_type.hpp
index a6c442d5..1d5a1765 100644
--- a/include/boost/histogram/detail/args_type.hpp
+++ b/include/boost/histogram/detail/args_type.hpp
@@ -16,7 +16,8 @@
 #if BOOST_WORKAROUND(BOOST_GCC, >= 60000)
 #pragma GCC diagnostic pop
 #endif
-#include <boost/mp11/list.hpp> // mp_pop_front
+#include <boost/mp11/list.hpp>    // mp_pop_front
+#include <boost/mp11/utility.hpp> // mp_if
 #include <tuple>
 #include <type_traits> // is_member_function_pointer
 
@@ -24,9 +25,9 @@ namespace boost {
 namespace histogram {
 namespace detail {
 
-template <class T, class Args = boost::callable_traits::args_t<T>>
-using args_type = std::conditional_t<std::is_member_function_pointer<T>::value,
-                                     mp11::mp_pop_front<Args>, Args>;
+template <class T, class Args = callable_traits::args_t<T>>
+using args_type =
+    mp11::mp_if<std::is_member_function_pointer<T>, mp11::mp_pop_front<Args>, Args>;
 
 template <class T, std::size_t N = 0>
 using arg_type = std::tuple_element_t<N, args_type<T>>;
diff --git a/include/boost/histogram/detail/at.hpp b/include/boost/histogram/detail/at.hpp
new file mode 100644
index 00000000..ea5f7232
--- /dev/null
+++ b/include/boost/histogram/detail/at.hpp
@@ -0,0 +1,65 @@
+// Copyright 2015-2018 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_HISTOGRAM_DETAIL_AT_HPP
+#define BOOST_HISTOGRAM_DETAIL_AT_HPP
+
+#include <boost/assert.hpp>
+#include <boost/histogram/axis/option.hpp>
+#include <boost/histogram/axis/traits.hpp>
+#include <boost/histogram/detail/axes.hpp>
+#include <boost/histogram/detail/linearize.hpp>
+#include <boost/histogram/detail/optional_index.hpp>
+#include <boost/histogram/fwd.hpp>
+#include <boost/mp11.hpp>
+#include <tuple>
+
+namespace boost {
+namespace histogram {
+namespace detail {
+
+template <class A>
+std::size_t linearize_index(optional_index& out, const std::size_t stride, const A& axis,
+                            const axis::index_type i) {
+  // A may be axis or variant, cannot use static option detection here
+  const auto opt = axis::traits::options(axis);
+  const auto shift = opt & axis::option::underflow ? 1 : 0;
+  const auto extent = axis.size() + (opt & axis::option::overflow ? 1 : 0) + shift;
+  // i may be arbitrarily out of range
+  using namespace boost::mp11;
+  linearize(mp_false{}, mp_false{}, out, stride, extent, i + shift);
+  return extent;
+}
+
+template <class A, class... Us>
+optional_index at(const A& axes, const std::tuple<Us...>& args) noexcept {
+  optional_index idx{0};
+  std::size_t stride = 1;
+  using namespace boost::mp11;
+  mp_for_each<mp_iota_c<sizeof...(Us)>>([&](auto i) {
+    stride *= linearize_index(idx, stride, axis_get<i>(axes),
+                              static_cast<axis::index_type>(std::get<i>(args)));
+  });
+  return idx;
+}
+
+template <class A, class U>
+optional_index at(const A& axes, const U& args) noexcept {
+  optional_index idx{0};
+  std::size_t stride = 1;
+  using std::begin;
+  auto it = begin(args);
+  for_each_axis(axes, [&](const auto& a) {
+    stride *= linearize_index(idx, stride, a, static_cast<axis::index_type>(*it++));
+  });
+  return idx;
+}
+
+} // namespace detail
+} // namespace histogram
+} // namespace boost
+
+#endif
diff --git a/include/boost/histogram/detail/attribute.hpp b/include/boost/histogram/detail/attribute.hpp
deleted file mode 100644
index 436b1149..00000000
--- a/include/boost/histogram/detail/attribute.hpp
+++ /dev/null
@@ -1,16 +0,0 @@
-// Copyright 2018 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_HISTOGRAM_DETAIL_ATTRIBUTE_HPP
-#define BOOST_HISTOGRAM_DETAIL_ATTRIBUTE_HPP
-
-#if __cplusplus >= 201703L
-#define BOOST_HISTOGRAM_NODISCARD [[nodiscard]]
-#else
-#define BOOST_HISTOGRAM_NODISCARD
-#endif
-
-#endif
diff --git a/include/boost/histogram/detail/axes.hpp b/include/boost/histogram/detail/axes.hpp
index 1372936d..cb4a1947 100644
--- a/include/boost/histogram/detail/axes.hpp
+++ b/include/boost/histogram/detail/axes.hpp
@@ -71,38 +71,39 @@ decltype(auto) axis_get(const T& axes) {
 }
 
 template <class... Ts>
-decltype(auto) axis_get(std::tuple<Ts...>& axes, unsigned i) {
-  using namespace boost::mp11;
+decltype(auto) axis_get(std::tuple<Ts...>& axes, const unsigned i) {
+  using namespace ::boost::mp11;
   constexpr auto S = sizeof...(Ts);
   using V = mp_unique<axis::variant<Ts*...>>;
   return mp_with_index<S>(i, [&axes](auto i) { return V(&std::get<i>(axes)); });
 }
 
 template <class... Ts>
-decltype(auto) axis_get(const std::tuple<Ts...>& axes, unsigned i) {
-  using namespace boost::mp11;
+decltype(auto) axis_get(const std::tuple<Ts...>& axes, const unsigned i) {
+  using namespace ::boost::mp11;
   constexpr auto S = sizeof...(Ts);
   using V = mp_unique<axis::variant<const Ts*...>>;
   return mp_with_index<S>(i, [&axes](auto i) { return V(&std::get<i>(axes)); });
 }
 
 template <class T>
-decltype(auto) axis_get(T& axes, unsigned i) {
-  return axes.at(i);
+decltype(auto) axis_get(T& axes, const unsigned i) {
+  return axes[i];
 }
 
 template <class T>
-decltype(auto) axis_get(const T& axes, unsigned i) {
-  return axes.at(i);
+decltype(auto) axis_get(const T& axes, const unsigned i) {
+  return axes[i];
 }
 
 template <class... Ts, class... Us>
 bool axes_equal(const std::tuple<Ts...>& ts, const std::tuple<Us...>& us) {
-  return static_if<std::is_same<mp11::mp_list<Ts...>, mp11::mp_list<Us...>>>(
+  using namespace ::boost::mp11;
+  return static_if<std::is_same<mp_list<Ts...>, mp_list<Us...>>>(
       [](const auto& ts, const auto& us) {
-        using N = mp11::mp_size<std::decay_t<decltype(ts)>>;
+        using N = mp_size<std::decay_t<decltype(ts)>>;
         bool equal = true;
-        mp11::mp_for_each<mp11::mp_iota<N>>(
+        mp_for_each<mp_iota<N>>(
             [&](auto I) { equal &= relaxed_equal(std::get<I>(ts), std::get<I>(us)); });
         return equal;
       },
@@ -111,10 +112,10 @@ bool axes_equal(const std::tuple<Ts...>& ts, const std::tuple<Us...>& us) {
 
 template <class T, class... Us>
 bool axes_equal(const T& t, const std::tuple<Us...>& u) {
+  using namespace ::boost::mp11;
   if (t.size() != sizeof...(Us)) return false;
   bool equal = true;
-  mp11::mp_for_each<mp11::mp_iota_c<sizeof...(Us)>>(
-      [&](auto I) { equal &= t[I] == std::get<I>(u); });
+  mp_for_each<mp_iota_c<sizeof...(Us)>>([&](auto I) { equal &= t[I] == std::get<I>(u); });
   return equal;
 }
 
@@ -131,7 +132,8 @@ bool axes_equal(const T& t, const U& u) {
 
 template <class... Ts, class... Us>
 void axes_assign(std::tuple<Ts...>& t, const std::tuple<Us...>& u) {
-  static_if<std::is_same<mp11::mp_list<Ts...>, mp11::mp_list<Us...>>>(
+  using namespace ::boost::mp11;
+  static_if<std::is_same<mp_list<Ts...>, mp_list<Us...>>>(
       [](auto& a, const auto& b) { a = b; },
       [](auto&, const auto&) {
         BOOST_THROW_EXCEPTION(
@@ -142,8 +144,9 @@ void axes_assign(std::tuple<Ts...>& t, const std::tuple<Us...>& u) {
 
 template <class... Ts, class U>
 void axes_assign(std::tuple<Ts...>& t, const U& u) {
-  mp11::mp_for_each<mp11::mp_iota_c<sizeof...(Ts)>>([&](auto I) {
-    using T = mp11::mp_at_c<std::tuple<Ts...>, I>;
+  using namespace ::boost::mp11;
+  mp_for_each<mp_iota_c<sizeof...(Ts)>>([&](auto I) {
+    using T = mp_at_c<std::tuple<Ts...>, I>;
     std::get<I>(t) = axis::get<T>(u[I]);
   });
 }
@@ -152,8 +155,8 @@ template <class T, class... Us>
 void axes_assign(T& t, const std::tuple<Us...>& u) {
   // resize instead of reserve, because t may not be empty and we want exact capacity
   t.resize(sizeof...(Us));
-  mp11::mp_for_each<mp11::mp_iota_c<sizeof...(Us)>>(
-      [&](auto I) { t[I] = std::get<I>(u); });
+  using namespace ::boost::mp11;
+  mp_for_each<mp_iota_c<sizeof...(Us)>>([&](auto I) { t[I] = std::get<I>(u); });
 }
 
 template <typename T, typename U>
@@ -175,30 +178,36 @@ auto make_empty_dynamic_axes(const std::tuple<Ts...>&) {
   return std::vector<mp_if_c<(mp_size<L>::value == 1), mp_first<L>, L>>{};
 }
 
-template <typename T>
+template <class T>
 void axis_index_is_valid(const T& axes, const unsigned N) {
   BOOST_ASSERT_MSG(N < axes_rank(axes), "index out of range");
 }
 
-template <typename F, typename T>
-void for_each_axis_impl(std::true_type, const T& axes, F&& f) {
-  for (const auto& x : axes) { axis::visit(std::forward<F>(f), x); }
+template <class Axes, class V>
+void for_each_axis_impl(std::true_type, Axes&& axes, V&& v) {
+  for (auto&& a : axes) { axis::visit(std::forward<V>(v), a); }
 }
 
-template <typename F, typename T>
-void for_each_axis_impl(std::false_type, const T& axes, F&& f) {
-  for (const auto& x : axes) std::forward<F>(f)(x);
+template <class Axes, class V>
+void for_each_axis_impl(std::false_type, Axes&& axes, V&& v) {
+  for (auto&& a : axes) std::forward<V>(v)(a);
 }
 
-template <typename F, typename T>
-void for_each_axis(const T& axes, F&& f) {
-  using U = mp11::mp_first<T>;
-  for_each_axis_impl(is_axis_variant<U>(), axes, std::forward<F>(f));
+template <class Axes, class V>
+void for_each_axis(Axes&& a, V&& v) {
+  using namespace ::boost::mp11;
+  using T = mp_first<std::decay_t<Axes>>;
+  for_each_axis_impl(is_axis_variant<T>(), std::forward<Axes>(a), std::forward<V>(v));
 }
 
-template <typename F, typename... Ts>
-void for_each_axis(const std::tuple<Ts...>& axes, F&& f) {
-  mp11::tuple_for_each(axes, std::forward<F>(f));
+template <class V, class... Axis>
+void for_each_axis(const std::tuple<Axis...>& a, V&& v) {
+  mp11::tuple_for_each(a, std::forward<V>(v));
+}
+
+template <class V, class... Axis>
+void for_each_axis(std::tuple<Axis...>& a, V&& v) {
+  mp11::tuple_for_each(a, std::forward<V>(v));
 }
 
 template <typename T>
@@ -214,12 +223,12 @@ std::size_t bincount(const T& axes) {
 }
 
 template <class T>
-using tuple_size_t = typename std::tuple_size<T>::type;
+using buffer_size_impl = typename std::tuple_size<T>::type;
 
 template <class T>
 using buffer_size = mp11::mp_eval_or<
-    std::integral_constant<std::size_t, BOOST_HISTOGRAM_DETAIL_AXES_LIMIT>, tuple_size_t,
-    T>;
+    std::integral_constant<std::size_t, BOOST_HISTOGRAM_DETAIL_AXES_LIMIT>,
+    buffer_size_impl, T>;
 
 template <class T, std::size_t N>
 class sub_array : public std::array<T, N> {
@@ -257,6 +266,49 @@ auto make_stack_buffer(const T& t, const Ts&... ts) {
   return stack_buffer<U, T>(axes_rank(t), ts...);
 }
 
+template <class T>
+using has_underflow =
+    decltype(axis::traits::static_options<T>::test(axis::option::underflow));
+
+template <class Axis>
+constexpr auto inclusive_options(Axis&&) {
+  return axis::option::underflow | axis::option::overflow;
+}
+
+template <class... Ts>
+constexpr auto inclusive_options(axis::category<Ts...>&&) {
+  return axis::option::overflow;
+}
+
+template <class T>
+using is_non_inclusive = mp11::mp_not<decltype(
+    axis::traits::static_options<T>::test(inclusive_options(std::declval<T>())))>;
+
+template <class T>
+using is_growing = decltype(axis::traits::static_options<T>::test(axis::option::growth));
+
+template <class T>
+using is_multidim = is_tuple<std::decay_t<arg_type<decltype(&T::index)>>>;
+
+template <template <class> class Trait, class T>
+struct has_special_axis_impl : Trait<T> {};
+
+template <template <class> class Trait, class... Ts>
+struct has_special_axis_impl<Trait, std::tuple<Ts...>> : mp11::mp_or<Trait<Ts>...> {};
+
+template <template <class> class Trait, class... Ts>
+struct has_special_axis_impl<Trait, axis::variant<Ts...>> : mp11::mp_or<Trait<Ts>...> {};
+
+template <template <class> class Trait, class T>
+using has_special_axis = typename has_special_axis_impl<
+    Trait, mp11::mp_if<is_vector_like<T>, mp11::mp_first<T>, T>>::type;
+
+template <class T>
+using has_growing_axis = has_special_axis<is_growing, T>;
+
+template <class T>
+using has_non_inclusive_axis = has_special_axis<is_non_inclusive, T>;
+
 } // namespace detail
 } // namespace histogram
 } // namespace boost
diff --git a/include/boost/histogram/detail/data.hpp b/include/boost/histogram/detail/data.hpp
new file mode 100644
index 00000000..5ef7adcb
--- /dev/null
+++ b/include/boost/histogram/detail/data.hpp
@@ -0,0 +1,48 @@
+// Copyright 2019 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_HISTOGRAM_DETAIL_DATA_HPP
+#define BOOST_HISTOGRAM_DETAIL_DATA_HPP
+
+#include <initializer_list>
+
+namespace boost {
+namespace histogram {
+namespace detail {
+
+#if __cpp_lib_nonmember_container_access >= 201411
+
+using std::data;
+
+#else
+
+template <class C>
+constexpr auto data(C& c) -> decltype(c.data()) {
+  return c.data();
+}
+
+template <class C>
+constexpr auto data(const C& c) -> decltype(c.data()) {
+  return c.data();
+}
+
+template <class T, std::size_t N>
+constexpr T* data(T (&array)[N]) noexcept {
+  return array;
+}
+
+template <class E>
+constexpr const E* data(std::initializer_list<E> il) noexcept {
+  return il.begin();
+}
+
+#endif
+
+} // namespace detail
+} // namespace histogram
+} // namespace boost
+
+#endif // BOOST_HISTOGRAM_DETAIL_DATA_HPP
diff --git a/include/boost/histogram/detail/detect.hpp b/include/boost/histogram/detail/detect.hpp
index c1cf8771..a4b90850 100644
--- a/include/boost/histogram/detail/detect.hpp
+++ b/include/boost/histogram/detail/detect.hpp
@@ -11,6 +11,7 @@
 #include <boost/mp11/algorithm.hpp>
 #include <boost/mp11/function.hpp>
 #include <boost/mp11/utility.hpp>
+#include <boost/variant2/variant.hpp>
 #include <iterator>
 #include <tuple>
 #include <type_traits>
@@ -19,17 +20,17 @@ namespace boost {
 namespace histogram {
 namespace detail {
 
-#define BOOST_HISTOGRAM_DETECT(name, cond)   \
-  template <class T, class = decltype(cond)> \
-  struct name##_impl {};                     \
-  template <class T>                         \
-  struct name : mp11::mp_valid<name##_impl, T>::type {}
+#define BOOST_HISTOGRAM_DETECT(name, cond) \
+  template <class T>                       \
+  using name##_impl = decltype(cond);      \
+  template <class T>                       \
+  using name = typename mp11::mp_valid<name##_impl, T>::type
 
-#define BOOST_HISTOGRAM_DETECT_BINARY(name, cond)     \
-  template <class T, class U, class = decltype(cond)> \
-  struct name##_impl {};                              \
-  template <class T, class U = T>                     \
-  struct name : mp11::mp_valid<name##_impl, T, U>::type {}
+#define BOOST_HISTOGRAM_DETECT_BINARY(name, cond) \
+  template <class T, class U>                     \
+  using name##_impl = decltype(cond);             \
+  template <class T, class U = T>                 \
+  using name = typename mp11::mp_valid<name##_impl, T, U>::type
 
 // metadata has overloads, trying to get pmf in this case always fails
 BOOST_HISTOGRAM_DETECT(has_method_metadata, (std::declval<T&>().metadata()));
@@ -45,12 +46,12 @@ BOOST_HISTOGRAM_DETECT(has_method_lower, &T::lower);
 
 BOOST_HISTOGRAM_DETECT(has_method_value, &T::value);
 
-BOOST_HISTOGRAM_DETECT(has_method_update, (&T::update));
+BOOST_HISTOGRAM_DETECT(has_method_update, &T::update);
 
 // reset has overloads, trying to get pmf in this case always fails
 BOOST_HISTOGRAM_DETECT(has_method_reset, (std::declval<T>().reset(0)));
 
-BOOST_HISTOGRAM_DETECT(has_method_options, (&T::options));
+BOOST_HISTOGRAM_DETECT(has_method_options, &T::options);
 
 BOOST_HISTOGRAM_DETECT(has_allocator, &T::get_allocator);
 
@@ -109,68 +110,77 @@ BOOST_HISTOGRAM_DETECT_BINARY(
 
 BOOST_HISTOGRAM_DETECT(has_threading_support, (T::has_threading_support));
 
-template <typename T>
-struct is_weight_impl : std::false_type {};
+template <class T>
+struct is_weight_impl : mp11::mp_false {};
 
-template <typename T>
-struct is_weight_impl<weight_type<T>> : std::true_type {};
+template <class T>
+struct is_weight_impl<weight_type<T>> : mp11::mp_true {};
 
-template <typename T>
-using is_weight = is_weight_impl<std::decay_t<T>>;
+template <class T>
+using is_weight = is_weight_impl<T>;
 
-template <typename T>
-struct is_sample_impl : std::false_type {};
+template <class T>
+struct is_sample_impl : mp11::mp_false {};
 
-template <typename T>
-struct is_sample_impl<sample_type<T>> : std::true_type {};
+template <class T>
+struct is_sample_impl<sample_type<T>> : mp11::mp_true {};
 
-template <typename T>
-using is_sample = is_sample_impl<std::decay_t<T>>;
+template <class T>
+using is_sample = is_sample_impl<T>;
 
-template <typename T>
+template <class T>
 using is_storage = mp11::mp_and<is_indexable_container<T>, has_method_reset<T>,
                                 has_threading_support<T>>;
 
 template <class T>
 using is_adaptible = mp11::mp_or<is_vector_like<T>, is_array_like<T>, is_map_like<T>>;
 
-template <class T, class _ = std::decay_t<T>,
-          class = std::enable_if_t<(is_storage<_>::value || is_adaptible<_>::value)>>
-struct requires_storage_or_adaptible {};
+template <class T>
+struct is_tuple_impl : mp11::mp_false {};
 
-template <typename T>
-struct is_tuple_impl : std::false_type {};
+template <class... Ts>
+struct is_tuple_impl<std::tuple<Ts...>> : mp11::mp_true {};
 
-template <typename... Ts>
-struct is_tuple_impl<std::tuple<Ts...>> : std::true_type {};
-
-template <typename T>
+template <class T>
 using is_tuple = typename is_tuple_impl<T>::type;
 
-template <typename T>
-struct is_axis_variant_impl : std::false_type {};
+template <class T>
+struct is_variant_impl : mp11::mp_false {};
 
-template <typename... Ts>
-struct is_axis_variant_impl<axis::variant<Ts...>> : std::true_type {};
+template <class... Ts>
+struct is_variant_impl<boost::variant2::variant<Ts...>> : mp11::mp_true {};
 
-template <typename T>
+template <class T>
+using is_variant = typename is_variant_impl<T>::type;
+
+template <class T>
+struct is_axis_variant_impl : mp11::mp_false {};
+
+template <class... Ts>
+struct is_axis_variant_impl<axis::variant<Ts...>> : mp11::mp_true {};
+
+template <class T>
 using is_axis_variant = typename is_axis_variant_impl<T>::type;
 
-template <typename T>
+template <class T>
 using is_any_axis = mp11::mp_or<is_axis<T>, is_axis_variant<T>>;
 
-template <typename T>
+template <class T>
 using is_sequence_of_axis = mp11::mp_and<is_iterable<T>, is_axis<mp11::mp_first<T>>>;
 
-template <typename T>
+template <class T>
 using is_sequence_of_axis_variant =
     mp11::mp_and<is_iterable<T>, is_axis_variant<mp11::mp_first<T>>>;
 
-template <typename T>
+template <class T>
 using is_sequence_of_any_axis =
     mp11::mp_and<is_iterable<T>, is_any_axis<mp11::mp_first<T>>>;
 
 // poor-mans concept checks
+template <class T, class _ = std::decay_t<T>,
+          class = std::enable_if_t<(is_storage<_>::value || is_adaptible<_>::value)>>
+struct requires_storage_or_adaptible {};
+
 template <class T, class = std::enable_if_t<is_iterator<std::decay_t<T>>::value>>
 struct requires_iterator {};
 
diff --git a/include/boost/histogram/detail/fill.hpp b/include/boost/histogram/detail/fill.hpp
new file mode 100644
index 00000000..3724f254
--- /dev/null
+++ b/include/boost/histogram/detail/fill.hpp
@@ -0,0 +1,345 @@
+// Copyright 2015-2018 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_HISTOGRAM_DETAIL_FILL_HPP
+#define BOOST_HISTOGRAM_DETAIL_FILL_HPP
+
+#include <algorithm>
+#include <boost/assert.hpp>
+#include <boost/config/workaround.hpp>
+#include <boost/histogram/axis/traits.hpp>
+#include <boost/histogram/axis/variant.hpp>
+#include <boost/histogram/detail/axes.hpp>
+#include <boost/histogram/detail/linearize.hpp>
+#include <boost/histogram/detail/make_default.hpp>
+#include <boost/histogram/detail/optional_index.hpp>
+#include <boost/histogram/detail/static_if.hpp>
+#include <boost/histogram/detail/tuple_slice.hpp>
+#include <boost/histogram/fwd.hpp>
+#include <boost/mp11.hpp>
+#include <mutex>
+#include <tuple>
+#include <type_traits>
+
+namespace boost {
+namespace histogram {
+namespace detail {
+
+template <class Index, class Axis, class Value>
+std::size_t linearize_growth(Index& o, axis::index_type& shift, const std::size_t stride,
+                             Axis& a, const Value& v) {
+  using O = axis::traits::static_options<Axis>;
+  axis::index_type i;
+  std::tie(i, shift) = axis::traits::update(a, v);
+  linearize(O::test(axis::option::underflow), O::test(axis::option::overflow), o, stride,
+            a.size(), i);
+  return axis::traits::extent(a);
+}
+
+template <class Index, class... Ts, class Value>
+std::size_t linearize_growth(Index& o, axis::index_type& sh, const std::size_t st,
+                             axis::variant<Ts...>& a, const Value& v) {
+  return axis::visit([&](auto& a) { return linearize_growth(o, sh, st, a, v); }, a);
+}
+
+template <class A>
+struct storage_grower {
+  const A& axes_;
+  struct {
+    axis::index_type idx, old_extent;
+    std::size_t new_stride;
+  } data_[buffer_size<A>::value];
+  std::size_t new_size_;
+
+  storage_grower(const A& axes) noexcept : axes_(axes) {}
+
+  void from_shifts(const axis::index_type* shifts) noexcept {
+    auto dit = data_;
+    std::size_t s = 1;
+    for_each_axis(axes_, [&](const auto& a) {
+      const auto n = axis::traits::extent(a);
+      *dit++ = {0, n - std::abs(*shifts++), s};
+      s *= n;
+    });
+    new_size_ = s;
+  }
+
+  // must be extents before any shifts were applied
+  void from_extents(const axis::index_type* old_extents) noexcept {
+    auto dit = data_;
+    std::size_t s = 1;
+    for_each_axis(axes_, [&](const auto& a) {
+      const auto n = axis::traits::extent(a);
+      *dit++ = {0, *old_extents++, s};
+      s *= n;
+    });
+    new_size_ = s;
+  }
+
+  template <class S>
+  void apply(S& storage, const axis::index_type* shifts) {
+    auto new_storage = make_default(storage);
+    new_storage.reset(new_size_);
+    const auto dlast = data_ + axes_rank(axes_) - 1;
+    for (const auto& x : storage) {
+      auto ns = new_storage.begin();
+      auto sit = shifts;
+      auto dit = data_;
+      for_each_axis(axes_, [&](const auto& a) {
+        using opt = axis::traits::static_options<decltype(a)>;
+        if (opt::test(axis::option::underflow)) {
+          if (dit->idx == 0) {
+            // axis has underflow and we are in the underflow bin:
+            // keep storage pointer unchanged
+            ++dit;
+            ++sit;
+            return;
+          }
+        }
+        if (opt::test(axis::option::overflow)) {
+          if (dit->idx == dit->old_extent - 1) {
+            // axis has overflow and we are in the overflow bin:
+            // move storage pointer to corresponding overflow bin position
+            ns += (axis::traits::extent(a) - 1) * dit->new_stride;
+            ++dit;
+            ++sit;
+            return;
+          }
+        }
+        // we are in a normal bin:
+        // move storage pointer to index position, apply positive shifts
+        ns += (dit->idx + std::max(*sit, 0)) * dit->new_stride;
+        ++dit;
+        ++sit;
+      });
+      // assign old value to new location
+      *ns = x;
+      // advance multi-dimensional index
+      dit = data_;
+      ++dit->idx;
+      while (dit != dlast && dit->idx == dit->old_extent) {
+        dit->idx = 0;
+        ++(++dit)->idx;
+      }
+    }
+    storage = std::move(new_storage);
+  }
+};
+
+template <class T, class... Us>
+inline void fill_storage_impl(mp11::mp_false, T&& t, Us&&... args) noexcept {
+  t(std::forward<Us>(args)...);
+}
+
+template <class T>
+inline void fill_storage_impl(mp11::mp_true, T&& t) noexcept {
+  ++t;
+}
+
+template <class T, class U>
+inline void fill_storage_impl(mp11::mp_true, T&& t, U&& w) noexcept {
+  t += w;
+}
+
+template <class T, class... Us>
+inline void fill_storage(T&& t, Us&&... args) noexcept {
+  fill_storage_impl(has_operator_preincrement<std::decay_t<T>>{}, std::forward<T>(t),
+                    std::forward<Us>(args)...);
+}
+
+template <class IW, class IS, class T, class U>
+void fill_storage_parse_args(IW, IS, T&& t, U&& u) noexcept {
+  mp11::tuple_apply(
+      [&](auto&&... args) {
+        fill_storage(std::forward<T>(t), std::get<IW::value>(u).value, args...);
+      },
+      std::get<IS::value>(u).value);
+}
+
+template <class IS, class T, class U>
+void fill_storage_parse_args(mp11::mp_int<-1>, IS, T&& t, U&& u) noexcept {
+  mp11::tuple_apply(
+      [&](const auto&... args) { fill_storage(std::forward<T>(t), args...); },
+      std::get<IS::value>(u).value);
+}
+
+template <class IW, class T, class U>
+void fill_storage_parse_args(IW, mp11::mp_int<-1>, T&& t, U&& u) noexcept {
+  fill_storage(std::forward<T>(t), std::get<IW::value>(u).value);
+}
+
+template <class T, class U>
+void fill_storage_parse_args(mp11::mp_int<-1>, mp11::mp_int<-1>, T&& t, U&&) noexcept {
+  fill_storage(std::forward<T>(t));
+}
+
+template <class L>
+struct args_indices {
+  static constexpr int _size = static_cast<int>(mp11::mp_size<L>::value);
+  static constexpr int _weight = static_cast<int>(mp11::mp_find_if<L, is_weight>::value);
+  static constexpr int _sample = static_cast<int>(mp11::mp_find_if<L, is_sample>::value);
+
+  static constexpr unsigned nargs = _size - (_weight < _size) - (_sample < _size);
+  static constexpr int start =
+      _weight < _size && _sample < _size && (_weight + _sample < 2)
+          ? 2
+          : ((_weight == 0 || _sample == 0) ? 1 : 0);
+  using weight = mp11::mp_int<(_weight < _size ? _weight : -1)>;
+  using sample = mp11::mp_int<(_sample < _size ? _sample : -1)>;
+};
+
+template <int S, int N>
+struct argument_loop {
+  template <class Index, class A, class Args>
+  static void impl(mp11::mp_int<N>, Index&, const std::size_t, A&, const Args&) {}
+
+  template <int I, class Index, class A, class Args>
+  static void impl(mp11::mp_int<I>, Index& o, const std::size_t s, A& ax,
+                   const Args& args) {
+    const auto e = linearize(o, s, axis_get<I>(ax), std::get<(S + I)>(args));
+    impl(mp11::mp_int<(I + 1)>{}, o, s * e, ax, args);
+  }
+
+  template <class Index, class A, class Args>
+  static void apply(Index& o, A& ax, const Args& args) {
+    impl(mp11::mp_int<0>{}, o, 1, ax, args);
+  }
+};
+
+template <int S>
+struct argument_loop<S, 1> {
+  template <class Index, class A, class Args>
+  static void apply(Index& o, A& ax, const Args& args) {
+    linearize(o, 1, axis_get<0>(ax), std::get<S>(args));
+  }
+};
+
+template <class A>
+constexpr unsigned min(const unsigned n) noexcept {
+  constexpr unsigned a = static_cast<unsigned>(buffer_size<A>::value);
+  return a < n ? a : n;
+}
+
+// not growing, only inclusive axes
+template <class S, class A, class Args>
+inline auto fill(mp11::mp_false, mp11::mp_false, S& storage, A& axes, const Args& args) {
+  using pos = args_indices<mp11::mp_transform<std::decay_t, Args>>;
+  std::size_t idx = 0;
+  argument_loop<pos::start, min<A>(pos::nargs)>::apply(idx, axes, args);
+  BOOST_ASSERT(idx < storage.size()); // idx is always valid
+  fill_storage_parse_args(typename pos::weight{}, typename pos::sample{}, storage[idx],
+                          args);
+  return storage.begin() + idx;
+}
+
+// not growing, at least one non-inclusive axis
+template <class S, class A, class Args>
+inline auto fill(mp11::mp_false, mp11::mp_true, S& storage, A& axes, const Args& args) {
+  using pos = args_indices<mp11::mp_transform<std::decay_t, Args>>;
+  optional_index idx{0};
+  argument_loop<pos::start, min<A>(pos::nargs)>::apply(idx, axes, args);
+  if (idx.valid()) {
+    fill_storage_parse_args(typename pos::weight{}, typename pos::sample{}, storage[*idx],
+                            args);
+    return storage.begin() + *idx;
+  }
+  return storage.end();
+}
+
+template <class S, class A, class Args>
+inline auto fill(mp11::mp_false, S& storage, A& axes, const Args& args) {
+  return fill(mp11::mp_false{}, has_non_inclusive_axis<A>{}, storage, axes, args);
+}
+
+template <class S, class A, class Args>
+inline auto fill(mp11::mp_true, S& storage, A& axes, const Args& args) {
+  using pos = args_indices<mp11::mp_transform<std::decay_t, Args>>;
+  axis::index_type shifts[pos::nargs];
+  optional_index idx{0};
+  std::size_t stride = 1;
+  bool update_needed = false;
+  mp11::mp_for_each<mp11::mp_iota_c<min<A>(pos::nargs)>>([&](auto i) {
+    stride *= linearize_growth(idx, shifts[i], stride, axis_get<i>(axes),
+                               std::get<(pos::start + i)>(args));
+    update_needed |= (shifts[i] != 0);
+  });
+  if (update_needed) {
+    storage_grower<A> g(axes);
+    g.from_shifts(shifts);
+    g.apply(storage, shifts);
+  }
+  if (idx.valid()) {
+    fill_storage_parse_args(typename pos::weight{}, typename pos::sample{}, storage[*idx],
+                            args);
+    return storage.begin() + *idx;
+  }
+  return storage.end();
+}
+
+// pack original args tuple into another tuple (which is unpacked later)
+template <int Start, int Size, class IW, class IS, class Args>
+decltype(auto) pack_args(IW, IS, const Args& args) noexcept {
+  return std::make_tuple(std::get<IW::value>(args), std::get<IS::value>(args),
+                         tuple_slice<Start, Size>(args));
+}
+
+template <int Start, int Size, class IW, class Args>
+decltype(auto) pack_args(IW, mp11::mp_int<-1>, const Args& args) noexcept {
+  return std::make_tuple(std::get<IW::value>(args), tuple_slice<Start, Size>(args));
+}
+
+template <int Start, int Size, class IS, class Args>
+decltype(auto) pack_args(mp11::mp_int<-1>, IS, const Args& args) noexcept {
+  return std::make_tuple(std::get<IS::value>(args), tuple_slice<Start, Size>(args));
+}
+
+template <int Start, int Size, class Args>
+decltype(auto) pack_args(mp11::mp_int<-1>, mp11::mp_int<-1>, const Args& args) noexcept {
+  return std::make_tuple(args);
+}
+
+#if BOOST_WORKAROUND(BOOST_MSVC, >= 0)
+#pragma warning(disable : 4702) // fixing warning would reduce code readability a lot
+#endif
+
+template <class S, class A, class Args>
+auto fill(S& storage, A& axes, const Args& args) {
+  using pos = args_indices<mp11::mp_transform<std::decay_t, Args>>;
+  using growing = has_growing_axis<A>;
+
+  // Sometimes we need to pack the tuple into another tuple:
+  // - histogram contains one axis which accepts tuple
+  // - user passes tuple to fill(...)
+  // Tuple is normally unpacked and arguments are processed, this causes pos::nargs > 1.
+  // Now we pack tuple into another tuple so that original tuple is send to axis.
+  // Notes:
+  // - has nice side-effect of making histogram::operator(1, 2) work as well
+  // - cannot detect call signature of axis at compile-time in all configurations
+  //   (axis::variant provides generic call interface and hides concrete
+  //   interface), so we throw at runtime if incompatible argument is passed (e.g.
+  //   3d tuple)
+
+  if (axes_rank(axes) == pos::nargs)
+    return fill(growing{}, storage, axes, args);
+  else if (axes_rank(axes) == 1 && axis::traits::rank(axis_get<0>(axes)) == pos::nargs)
+    return fill(growing{}, storage, axes,
+                pack_args<pos::start, pos::nargs>(typename pos::weight{},
+                                                  typename pos::sample{}, args));
+  else
+    return (BOOST_THROW_EXCEPTION(
+                std::invalid_argument("number of arguments != histogram rank")),
+            storage.end());
+}
+
+#if BOOST_WORKAROUND(BOOST_MSVC, >= 0)
+#pragma warning(default : 4702)
+#endif
+
+} // namespace detail
+} // namespace histogram
+} // namespace boost
+
+#endif
diff --git a/include/boost/histogram/detail/fill_n.hpp b/include/boost/histogram/detail/fill_n.hpp
new file mode 100644
index 00000000..7eb79b1a
--- /dev/null
+++ b/include/boost/histogram/detail/fill_n.hpp
@@ -0,0 +1,291 @@
+// Copyright 2019 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_HISTOGRAM_DETAIL_FILL_N_HPP
+#define BOOST_HISTOGRAM_DETAIL_FILL_N_HPP
+
+#include <algorithm>
+#include <boost/assert.hpp>
+#include <boost/histogram/axis/option.hpp>
+#include <boost/histogram/axis/traits.hpp>
+#include <boost/histogram/detail/axes.hpp>
+#include <boost/histogram/detail/detect.hpp>
+#include <boost/histogram/detail/fill.hpp>
+#include <boost/histogram/detail/linearize.hpp>
+#include <boost/histogram/detail/non_member_container_access.hpp>
+#include <boost/histogram/detail/optional_index.hpp>
+#include <boost/histogram/detail/span.hpp>
+#include <boost/histogram/detail/static_if.hpp>
+#include <boost/histogram/fwd.hpp>
+#include <boost/mp11.hpp>
+#include <boost/throw_exception.hpp>
+#include <boost/variant2/variant.hpp>
+#include <stdexcept>
+#include <type_traits>
+
+namespace boost {
+namespace histogram {
+namespace detail {
+
+template <class Index, class Axis>
+struct indexing_visitor {
+  using index_type = Index;
+  using iterator = index_type*;
+
+  using Opt = axis::traits::static_options<Axis>;
+
+  Axis& axis_;
+  const std::size_t stride_, start_, size_; // start and size of value collection
+  const iterator& begin_;
+  axis::index_type* shift_;
+
+  indexing_visitor(Axis& a, const std::size_t& str, const std::size_t& sta,
+                   const std::size_t& si, const iterator& it, axis::index_type* shift)
+      : axis_(a), stride_(str), start_(sta), size_(si), begin_(it), shift_(shift) {}
+
+  void maybe_shift_previous_indices(iterator it, const axis::index_type shift) {
+    if (shift > 0) {
+      while (it != begin_) *--it += static_cast<std::size_t>(shift) * stride_;
+      *shift_ += shift;
+    }
+  }
+
+  template <class T>
+  void impl(mp11::mp_false, iterator it, const T& t) {
+    constexpr auto u = Opt::test(axis::option::underflow);
+    constexpr auto o = Opt::test(axis::option::overflow);
+    linearize(u, o, *it, stride_, axis_.size(), axis::traits::index(axis_, t));
+  }
+
+  template <class T>
+  void impl(mp11::mp_true, iterator it, const T& t) {
+    constexpr auto u = Opt::test(axis::option::underflow);
+    constexpr auto o = Opt::test(axis::option::overflow);
+    const auto p = axis::traits::update(axis_, t);
+    maybe_shift_previous_indices(it, p.second);
+    linearize(u, o, *it, stride_, axis_.size(), p.first);
+  }
+
+  template <class T>
+  void operator()(const T& t) {
+    static_if<is_iterable<T>>(
+        [this](const auto& t) {
+          // T is iterable, fill N values
+          constexpr auto g = Opt::test(axis::option::growth);
+          const auto tp = data(t) + start_;
+          auto it = begin_;
+          for (const auto& x : make_span(tp, size_)) this->impl(g, it++, x);
+        },
+        [this](const auto& t) {
+          // T is value, fill single value N times
+          constexpr auto g = Opt::test(axis::option::growth);
+          index_type o{0};
+          this->impl(g, &o, t);
+          for (auto&& i : make_span(begin_, size_)) i += o;
+        },
+        t);
+  }
+};
+
+template <class Index, class S, class A, class T>
+void fill_n_indices(Index* indices, const std::size_t start, const std::size_t size,
+                    S& storage, A& axes, const T* viter) {
+  axis::index_type extents[buffer_size<A>::value];
+  axis::index_type shifts[buffer_size<A>::value];
+  for_each_axis(axes, [eit = extents, sit = shifts](const auto& a) mutable {
+    *sit++ = 0;
+    *eit++ = axis::traits::extent(a);
+  });
+
+  std::fill(indices, indices + size, 0); // initialize to zero
+  for_each_axis(axes, [start, size, indices, &viter, stride = static_cast<std::size_t>(1),
+                       shift = shifts](auto& axis) mutable {
+    using Axis = std::decay_t<decltype(axis)>;
+    static_if<is_variant<T>>(
+        [&](const auto& v) {
+          variant2::visit(
+              indexing_visitor<Index, Axis>{axis, stride, start, size, indices, shift},
+              v);
+        },
+        [&](const auto& v) {
+          indexing_visitor<Index, Axis>{axis, stride, start, size, indices, shift}(v);
+        },
+        *viter++);
+    stride *= static_cast<std::size_t>(axis::traits::extent(axis));
+    ++shift;
+  });
+
+  bool update_needed = false;
+  for_each_axis(axes, [&update_needed, eit = extents](const auto& a) mutable {
+    update_needed |= *eit++ != axis::traits::extent(a);
+  });
+  if (update_needed) {
+    storage_grower<A> g(axes);
+    g.from_extents(extents);
+    g.apply(storage, shifts);
+  }
+}
+
+template <class S, class T, class... Us>
+void fill_n_storage(S& s, const std::size_t& idx, const T*&& wptr, std::size_t wsize,
+                    const Us*&&... sptrs) {
+  BOOST_ASSERT(idx < s.size());
+  fill_storage(s[idx], (wsize == 1u ? *wptr : *wptr++), *sptrs++...);
+}
+
+template <class S, class... Ts>
+void fill_n_storage(S& s, const std::size_t& idx, const Ts*&&... sptrs) {
+  BOOST_ASSERT(idx < s.size());
+  fill_storage(s[idx], *sptrs++...);
+}
+
+template <class S, class... Ts>
+void fill_n_storage(S& s, const optional_index& idx, Ts&&... ts) {
+  if (idx.valid()) fill_n_storage(s, *idx, std::forward<Ts>(ts)...);
+}
+
+template <class T>
+std::size_t get_total_size(const T* values, std::size_t vsize) {
+  std::size_t s = 1u;
+  auto vis = [&s](const auto& v) {
+    using U = std::remove_cv_t<std::remove_reference_t<decltype(v)>>;
+    const std::size_t n = static_if<is_iterable<U>>(
+        [](const auto& v) { return size(v); },
+        [](const auto&) { return static_cast<std::size_t>(1); }, v);
+    if (s == 1u)
+      s = n;
+    else if (n > 1u && s != n)
+      throw_exception(std::invalid_argument("spans must have same length"));
+  };
+  for (const auto& v : make_span(values, vsize))
+    static_if<is_iterable<T>>([&vis](const auto& v) { vis(v); },
+                              [&vis](const auto& v) { variant2::visit(vis, v); }, v);
+  return s;
+}
+
+// general Nd treatment
+template <class S, class A, class T, class... Ts>
+void fill_n_nd(S& storage, A& axes, const T* values, std::size_t vsize, Ts&&... rest) {
+  constexpr std::size_t buffer_size = 1ul << 14;
+  using index_type = mp11::mp_if<has_non_inclusive_axis<A>, optional_index, std::size_t>;
+  index_type indices[buffer_size];
+
+  /*
+    Parallelization options for generic case.
+
+    A) Run the whole fill2 method in parallel, each thread fills its own buffer of
+    indices, synchronization (atomics) are needed to synchronize the incrementing of
+    the storage cells. This leads to a lot of congestion for small histograms.
+
+    B) Run only fill_n_indices in parallel, subsections of the indices buffer
+    can be filled by different threads. The final loop that fills the storage runs
+    in the main thread, this requires no synchronization for the storage, cells do
+    not need to support atomic operations.
+
+    C) Like B), then sort the indices in the main thread and fill the
+    storage in parallel, where each thread uses a disjunct set of indices. This
+    should create less congestion and requires no synchronization for the storage.
+
+    Note on C): Let's say we have an axis with 5 bins (with *flow to simplify).
+    Then after filling 10 values, converting to indices and sorting, the index
+    buffer may look like this: 0 0 0 1 2 2 2 4 4 5. Let's use two threads to fill
+    the storage. Still in the main thread, we compute an iterator to the middle of
+    the index buffer and move it to the right until the pointee changes. Now we have
+    two ranges which contain disjunct sets of indices. We pass these ranges to the
+    threads which then fill the storage. Since the threads by construction do not
+    compete to increment the same cell, no further synchronization is required.
+
+    In all cases, growing axes cannot be parallelized.
+  */
+
+  for (std::size_t start = 0; start < vsize; start += buffer_size) {
+    const std::size_t n = std::min(buffer_size, vsize - start);
+    // fill buffer of indices...
+    fill_n_indices(indices, start, n, storage, axes, values);
+    // ...and fill corresponding storage cells
+    for (auto&& idx : make_span(indices, n))
+      fill_n_storage(storage, idx, std::forward<Ts>(rest)...);
+  }
+}
+
+template <class S, class A, class T, class... Us>
+void fill_n_impl(mp11::mp_false, S& storage, A& axes, const T* values, std::size_t vsize,
+                 Us&&... rest) {
+  fill_n_nd(storage, axes, values, vsize, std::forward<Us>(rest)...);
+}
+
+// unpack weight argument, can be iterable or value
+template <class S, class A, class T, class U, class... Us>
+void fill_n_impl(mp11::mp_true, S& storage, A& axes, const T* values, std::size_t vsize,
+                 const weight_type<U>& weights, const Us&... rest) {
+  static_if<is_iterable<std::remove_cv_t<std::remove_reference_t<U>>>>(
+      [&](const auto& w, const auto&... rest) {
+        const auto wsize = size(w);
+        if (vsize != wsize)
+          throw_exception(
+              std::invalid_argument("number of arguments must match histogram rank"));
+        fill_n_impl(mp11::mp_bool<sizeof...(Us)>{}, storage, axes, values, vsize, data(w),
+                    wsize, rest...);
+      },
+      [&](const auto w, const auto&... rest) {
+        fill_n_impl(mp11::mp_bool<sizeof...(Us)>{}, storage, axes, values, vsize, &w,
+                    static_cast<std::size_t>(1), rest...);
+      },
+      weights.value, rest...);
+}
+
+// unpack sample argument after weight was unpacked
+template <class S, class A, class T, class U, class V>
+void fill_n_impl(mp11::mp_true, S& storage, A& axes, const T* values, std::size_t vsize,
+                 const U* wptr, std::size_t wsize, const sample_type<V>& s) {
+  mp11::tuple_apply(
+      [&](const auto&... sargs) {
+        fill_n_impl(mp11::mp_false{}, storage, axes, values, vsize, std::move(wptr),
+                    wsize, data(sargs)...);
+      },
+      s.value);
+}
+
+// unpack sample argument (no weight argument)
+template <class S, class A, class T, class U>
+void fill_n_impl(mp11::mp_true, S& storage, A& axes, const T* values, std::size_t vsize,
+                 const sample_type<U>& samples) {
+  using namespace boost::mp11;
+  tuple_apply(
+      [&](const auto&... sargs) {
+        fill_n_impl(mp_false{}, storage, axes, values, vsize, data(sargs)...);
+      },
+      samples.value);
+}
+
+template <class S, class A, class T, class... Us>
+void fill_n(S& storage, A& axes, const T* values, std::size_t vsize, const Us&... rest) {
+  static_assert(!std::is_pointer<T>::value,
+                "passing iterable of pointers not allowed (cannot determine lengths); "
+                "pass iterable of iterables instead");
+  using namespace boost::mp11;
+  using LUs = mp_list<Us...>;
+  using Unpack = mp_or<mp_count_if<LUs, is_weight>, mp_count_if<LUs, is_sample>>;
+  static_if<mp_or<is_iterable<T>, is_variant<T>>>(
+      [&](const auto& values) {
+        if (axes_rank(axes) != vsize)
+          throw_exception(
+              std::invalid_argument("number of arguments must match histogram rank"));
+        fill_n_impl(Unpack{}, storage, axes, values, get_total_size(values, vsize),
+                    rest...);
+      },
+      [&](const auto& values) {
+        auto s = make_span(values, vsize);
+        fill_n_impl(Unpack{}, storage, axes, &s, vsize, rest...);
+      },
+      values);
+}
+
+} // namespace detail
+} // namespace histogram
+} // namespace boost
+
+#endif // BOOST_HISTOGRAM_DETAIL_FILL_N_HPP
diff --git a/include/boost/histogram/detail/linearize.hpp b/include/boost/histogram/detail/linearize.hpp
index 66e062fe..f84d1c89 100644
--- a/include/boost/histogram/detail/linearize.hpp
+++ b/include/boost/histogram/detail/linearize.hpp
@@ -1,4 +1,4 @@
-// Copyright 2015-2018 Hans Dembinski
+// Copyright 2019 Hans Dembinski
 //
 // Distributed under the Boost Software License, Version 1.0.
 // (See accompanying file LICENSE_1_0.txt
@@ -7,409 +7,102 @@
 #ifndef BOOST_HISTOGRAM_DETAIL_LINEARIZE_HPP
 #define BOOST_HISTOGRAM_DETAIL_LINEARIZE_HPP
 
-#include <algorithm>
 #include <boost/assert.hpp>
+#include <boost/histogram/axis/option.hpp>
 #include <boost/histogram/axis/traits.hpp>
 #include <boost/histogram/axis/variant.hpp>
-#include <boost/histogram/detail/args_type.hpp>
-#include <boost/histogram/detail/axes.hpp>
-#include <boost/histogram/detail/make_default.hpp>
-#include <boost/histogram/detail/static_if.hpp>
-#include <boost/histogram/detail/tuple_slice.hpp>
+#include <boost/histogram/detail/optional_index.hpp>
 #include <boost/histogram/fwd.hpp>
-#include <boost/mp11/algorithm.hpp>
-#include <boost/mp11/function.hpp>
-#include <boost/mp11/integral.hpp>
-#include <boost/mp11/list.hpp>
-#include <boost/mp11/tuple.hpp>
-#include <boost/throw_exception.hpp>
-#include <mutex>
-#include <stdexcept>
-#include <tuple>
-#include <type_traits>
 
 namespace boost {
 namespace histogram {
 namespace detail {
 
-template <class T>
-using has_underflow =
-    decltype(axis::traits::static_options<T>::test(axis::option::underflow));
-
-template <class T>
-using is_growing = decltype(axis::traits::static_options<T>::test(axis::option::growth));
-
-template <class T>
-using is_multidim = is_tuple<std::decay_t<arg_type<decltype(&T::index)>>>;
-
-template <template <class> class Trait, class T>
-struct has_special_axis_impl : Trait<T> {};
-
-template <template <class> class Trait, class... Ts>
-struct has_special_axis_impl<Trait, std::tuple<Ts...>> : mp11::mp_or<Trait<Ts>...> {};
-
-template <template <class> class Trait, class... Ts>
-struct has_special_axis_impl<Trait, axis::variant<Ts...>> : mp11::mp_or<Trait<Ts>...> {};
-
-template <template <class> class Trait, class T>
-using has_special_axis =
-    has_special_axis_impl<Trait, mp11::mp_if<is_vector_like<T>, mp11::mp_first<T>, T>>;
-
-template <class T>
-using has_multidim_axis = has_special_axis<is_multidim, T>;
-
-template <class T>
-using has_growing_axis = has_special_axis<is_growing, T>;
-
-/// Index with an invalid state
-struct optional_index {
-  std::size_t idx = 0;
-  std::size_t stride = 1;
-  operator bool() const { return stride > 0; }
-  std::size_t operator*() const { return idx; }
-};
-
 // no underflow, no overflow
-inline void linearize(std::false_type, std::false_type, optional_index& out,
-                      const axis::index_type size, const axis::index_type i) noexcept {
-  out.idx += i * out.stride;
-  out.stride *= i >= 0 && i < size ? size : 0;
+inline void linearize(mp11::mp_false, mp11::mp_false, optional_index& out,
+                      const std::size_t stride, const axis::index_type size,
+                      const axis::index_type i) noexcept {
+  if (i >= 0 && i < size)
+    out += i * stride;
+  else
+    out = optional_index::invalid;
 }
 
 // no underflow, overflow
-inline void linearize(std::false_type, std::true_type, optional_index& out,
-                      const axis::index_type size, const axis::index_type i) noexcept {
+inline void linearize(mp11::mp_false, mp11::mp_true, optional_index& out,
+                      const std::size_t stride, const axis::index_type size,
+                      const axis::index_type i) noexcept {
   BOOST_ASSERT(i <= size);
-  out.idx += i * out.stride;
-  out.stride *= i >= 0 ? size + 1 : 0;
+  if (i >= 0)
+    out += i * stride;
+  else
+    out = optional_index::invalid;
 }
 
 // underflow, no overflow
-inline void linearize(std::true_type, std::false_type, optional_index& out,
-                      const axis::index_type size, const axis::index_type i) noexcept {
-  // internal index must be shifted by +1 since axis has underflow bin
+inline void linearize(mp11::mp_true, mp11::mp_false, optional_index& out,
+                      const std::size_t stride, const axis::index_type size,
+                      const axis::index_type i) noexcept {
+  // internal index must be shifted by + 1 since axis has underflow bin
   BOOST_ASSERT(i + 1 >= 0);
-  out.idx += (i + 1) * out.stride;
-  out.stride *= i < size ? size + 1 : 0;
+  if (i < size)
+    out += (i + 1) * stride;
+  else
+    out = optional_index::invalid;
 }
 
 // underflow, overflow
-inline void linearize(std::true_type, std::true_type, optional_index& out,
-                      const axis::index_type size, const axis::index_type i) noexcept {
+inline void linearize(mp11::mp_true, mp11::mp_true, optional_index& out,
+                      const std::size_t stride, const axis::index_type size,
+                      const axis::index_type i) noexcept {
   // internal index must be shifted by +1 since axis has underflow bin
   BOOST_ASSERT(i + 1 >= 0);
-  BOOST_ASSERT(i <= size);
-  out.idx += (i + 1) * out.stride;
-  out.stride *= size + 2;
+  BOOST_ASSERT(i < size + 1);
+  out += (i + 1) * stride;
+}
+
+template <class HasUnderflow, class HasOverflow>
+inline void linearize(HasUnderflow, HasOverflow, std::size_t& out,
+                      const std::size_t stride, const axis::index_type size,
+                      const axis::index_type i) noexcept {
+  // internal index must be shifted by +1 if axis has underflow bin
+  BOOST_ASSERT((HasUnderflow::value ? i + 1 : i) >= 0);
+  BOOST_ASSERT(i < (HasOverflow::value ? size + 1 : size));
+  out += (HasUnderflow::value ? i + 1 : i) * stride;
 }
 
 template <class Axis, class Value>
-void linearize_value(optional_index& o, const Axis& a, const Value& v) {
+std::size_t linearize(optional_index& o, const std::size_t s, const Axis& a,
+                      const Value& v) {
   using O = axis::traits::static_options<Axis>;
-  linearize(O::test(axis::option::underflow), O::test(axis::option::overflow), o,
+  linearize(O::test(axis::option::underflow), O::test(axis::option::overflow), o, s,
             a.size(), axis::traits::index(a, v));
-}
-
-template <class... Ts, class Value>
-void linearize_value(optional_index& o, const axis::variant<Ts...>& a, const Value& v) {
-  axis::visit([&o, &v](const auto& a) { linearize_value(o, a, v); }, a);
+  return axis::traits::extent(a);
 }
 
 template <class Axis, class Value>
-axis::index_type linearize_value_growth(optional_index& o, Axis& a, const Value& v) {
+std::size_t linearize(std::size_t& o, const std::size_t s, const Axis& a,
+                      const Value& v) {
   using O = axis::traits::static_options<Axis>;
-  axis::index_type i, s;
-  std::tie(i, s) = axis::traits::update(a, v);
-  linearize(O::test(axis::option::underflow), O::test(axis::option::overflow), o,
-            a.size(), i);
-  return s;
+  linearize(O::test(axis::option::underflow), O::test(axis::option::overflow), o, s,
+            a.size(), axis::traits::index(a, v));
+  return axis::traits::extent(a);
 }
 
 template <class... Ts, class Value>
-axis::index_type linearize_value_growth(optional_index& o, axis::variant<Ts...>& a,
-                                        const Value& v) {
-  return axis::visit([&o, &v](auto& a) { return linearize_value_growth(o, a, v); }, a);
+std::size_t linearize(optional_index& o, const std::size_t s,
+                      const axis::variant<Ts...>& a, const Value& v) {
+  return axis::visit([&o, &s, &v](const auto& a) { return linearize(o, s, a, v); }, a);
 }
 
-template <class A>
-void linearize_index(optional_index& out, const A& axis, const axis::index_type i) {
-  // A may be axis or variant, cannot use static option detection here
-  const auto opt = axis::traits::options(axis);
-  const auto shift = opt & axis::option::underflow ? 1 : 0;
-  const auto extent = axis.size() + (opt & axis::option::overflow ? 1 : 0) + shift;
-  // i may be arbitrarily out of range
-  linearize(std::false_type{}, std::false_type{}, out, extent, i + shift);
-}
-
-template <class S, class A>
-void grow_storage(const A& axes, S& storage, const axis::index_type* shifts) {
-  struct item {
-    axis::index_type idx, old_extent;
-    std::size_t new_stride;
-  } data[buffer_size<A>::value];
-  const auto* sit = shifts;
-  auto dit = data;
-  std::size_t s = 1;
-  for_each_axis(axes, [&](const auto& a) {
-    const auto n = axis::traits::extent(a);
-    *dit++ = {0, n - std::abs(*sit++), s};
-    s *= n;
-  });
-  auto new_storage = make_default(storage);
-  new_storage.reset(bincount(axes));
-  const auto dlast = data + axes_rank(axes) - 1;
-  for (const auto& x : storage) {
-    auto ns = new_storage.begin();
-    sit = shifts;
-    dit = data;
-    for_each_axis(axes, [&](const auto& a) {
-      using opt = axis::traits::static_options<decltype(a)>;
-      if (opt::test(axis::option::underflow)) {
-        if (dit->idx == 0) {
-          // axis has underflow and we are in the underflow bin:
-          // keep storage pointer unchanged
-          ++dit;
-          ++sit;
-          return;
-        }
-      }
-      if (opt::test(axis::option::overflow)) {
-        if (dit->idx == dit->old_extent - 1) {
-          // axis has overflow and we are in the overflow bin:
-          // move storage pointer to corresponding overflow bin position
-          ns += (axis::traits::extent(a) - 1) * dit->new_stride;
-          ++dit;
-          ++sit;
-          return;
-        }
-      }
-      // we are in a normal bin:
-      // move storage pointer to index position, apply positive shifts
-      ns += (dit->idx + (std::max)(*sit, 0)) * dit->new_stride;
-      ++dit;
-      ++sit;
-    });
-    // assign old value to new location
-    *ns = x;
-    // advance multi-dimensional index
-    dit = data;
-    ++dit->idx;
-    while (dit != dlast && dit->idx == dit->old_extent) {
-      dit->idx = 0;
-      ++(++dit)->idx;
-    }
-  }
-  storage = std::move(new_storage);
-}
-
-// histogram has no growing and no multidim axis, axis rank known at compile-time
-template <class S, class... As, class... Us>
-optional_index index(std::false_type, std::false_type, const std::tuple<As...>& axes, S&,
-                     const std::tuple<Us...>& args) {
-  optional_index idx;
-  static_assert(sizeof...(As) == sizeof...(Us), "number of arguments != histogram rank");
-  mp11::mp_for_each<mp11::mp_iota_c<sizeof...(As)>>(
-      [&](auto i) { linearize_value(idx, axis_get<i>(axes), std::get<i>(args)); });
-  return idx;
-}
-
-// histogram has no growing and no multidim axis, axis rank known at run-time
-template <class A, class S, class U>
-optional_index index(std::false_type, std::false_type, const A& axes, S&, const U& args) {
-  constexpr auto nargs = static_cast<unsigned>(std::tuple_size<U>::value);
-  optional_index idx;
-  if (axes_rank(axes) != nargs)
-    BOOST_THROW_EXCEPTION(std::invalid_argument("number of arguments != histogram rank"));
-  constexpr auto nbuf = buffer_size<A>::value;
-  mp11::mp_for_each<mp11::mp_iota_c<(nargs < nbuf ? nargs : nbuf)>>(
-      [&](auto i) { linearize_value(idx, axis_get<i>(axes), std::get<i>(args)); });
-  return idx;
-}
-
-// special case: if histogram::operator()(tuple(1, 2)) is called on 1d histogram
-// with axis that accepts 2d tuple, this should not fail
-// - solution is to forward tuples of size > 1 directly to axis for 1d
-// histograms
-// - has nice side-effect of making histogram::operator(1, 2) work as well
-// - cannot detect call signature of axis at compile-time in all configurations
-//   (axis::variant provides generic call interface and hides concrete
-//   interface), so we throw at runtime if incompatible argument is passed (e.g.
-//   3d tuple)
-
-// histogram has no growing multidim axis, axis rank == 1 known at compile-time
-template <class S, class A, class U>
-optional_index index(std::false_type, std::true_type, const std::tuple<A>& axes, S&,
-                     const U& args) {
-  optional_index idx;
-  linearize_value(idx, axis_get<0>(axes), args);
-  return idx;
-}
-
-// histogram has no growing multidim axis, axis rank > 1 known at compile-time
-template <class S, class U, class A0, class A1, class... As>
-optional_index index(std::false_type, std::true_type,
-                     const std::tuple<A0, A1, As...>& axes, S& s, const U& args) {
-  return index(std::false_type{}, std::false_type{}, axes, s, args);
-}
-
-// histogram has no growing multidim axis, axis rank known at run-time
-template <class A, class S, class U>
-optional_index index(std::false_type, std::true_type, const A& axes, S&, const U& args) {
-  optional_index idx;
-  const auto rank = axes_rank(axes);
-  constexpr auto nargs = static_cast<unsigned>(std::tuple_size<U>::value);
-  if (rank == 1 && nargs > 1)
-    linearize_value(idx, axis_get<0>(axes), args);
-  else {
-    if (rank != nargs)
-      BOOST_THROW_EXCEPTION(
-          std::invalid_argument("number of arguments != histogram rank"));
-    constexpr auto nbuf = buffer_size<A>::value;
-    mp11::mp_for_each<mp11::mp_iota_c<(nargs < nbuf ? nargs : nbuf)>>(
-        [&](auto i) { linearize_value(idx, axis_get<i>(axes), std::get<i>(args)); });
-  }
-  return idx;
-}
-
-// histogram has growing axis
-template <class B, class T, class S, class U>
-optional_index index(std::true_type, B, T& axes, S& storage, const U& args) {
-  optional_index idx;
-  constexpr unsigned nbuf = buffer_size<T>::value;
-  axis::index_type shifts[nbuf];
-  const auto rank = axes_rank(axes);
-  constexpr auto nargs = static_cast<unsigned>(std::tuple_size<U>::value);
-
-  bool update_needed = false;
-  if (rank == 1 && nargs > 1) {
-    shifts[0] = linearize_value_growth(idx, axis_get<0>(axes), args);
-    update_needed = (shifts[0] != 0);
-  } else {
-    if (rank != nargs)
-      BOOST_THROW_EXCEPTION(
-          std::invalid_argument("number of arguments != histogram rank"));
-    mp11::mp_for_each<mp11::mp_iota_c<(nargs < nbuf ? nargs : nbuf)>>([&](auto i) {
-      shifts[i] = linearize_value_growth(idx, axis_get<i>(axes), std::get<i>(args));
-      update_needed |= (shifts[i] != 0);
-    });
-  }
-
-  if (update_needed) grow_storage(axes, storage, shifts);
-  return idx;
-}
-
-template <class U>
-constexpr auto weight_sample_indices() noexcept {
-  if (is_weight<U>::value) return std::make_pair(0, -1);
-  if (is_sample<U>::value) return std::make_pair(-1, 0);
-  return std::make_pair(-1, -1);
-}
-
-template <class U0, class U1, class... Us>
-constexpr auto weight_sample_indices() noexcept {
-  using L = mp11::mp_list<U0, U1, Us...>;
-  const int n = sizeof...(Us) + 1;
-  if (is_weight<mp11::mp_at_c<L, 0>>::value) {
-    if (is_sample<mp11::mp_at_c<L, 1>>::value) return std::make_pair(0, 1);
-    if (is_sample<mp11::mp_at_c<L, n>>::value) return std::make_pair(0, n);
-    return std::make_pair(0, -1);
-  }
-  if (is_sample<mp11::mp_at_c<L, 0>>::value) {
-    if (is_weight<mp11::mp_at_c<L, 1>>::value) return std::make_pair(1, 0);
-    if (is_weight<mp11::mp_at_c<L, n>>::value) return std::make_pair(n, 0);
-    return std::make_pair(-1, 0);
-  }
-  if (is_weight<mp11::mp_at_c<L, n>>::value) {
-    // 0, n already covered
-    if (is_sample<mp11::mp_at_c<L, (n - 1)>>::value) return std::make_pair(n, n - 1);
-    return std::make_pair(n, -1);
-  }
-  if (is_sample<mp11::mp_at_c<L, n>>::value) {
-    // n, 0 already covered
-    if (is_weight<mp11::mp_at_c<L, (n - 1)>>::value) return std::make_pair(n - 1, n);
-    return std::make_pair(-1, n);
-  }
-  return std::make_pair(-1, -1);
-}
-
-template <class T, class U>
-void fill_impl(mp11::mp_int<-1>, mp11::mp_int<-1>, std::false_type, T&& t,
-               const U&) noexcept {
-  t();
-}
-
-template <class T, class U>
-void fill_impl(mp11::mp_int<-1>, mp11::mp_int<-1>, std::true_type, T&& t,
-               const U&) noexcept {
-  ++t;
-}
-
-template <class IW, class T, class U>
-void fill_impl(IW, mp11::mp_int<-1>, std::false_type, T&& t, const U& u) noexcept {
-  t(std::get<IW::value>(u).value);
-}
-
-template <class IW, class T, class U>
-void fill_impl(IW, mp11::mp_int<-1>, std::true_type, T&& t, const U& u) noexcept {
-  t += std::get<IW::value>(u).value;
-}
-
-template <class IS, class T, class U>
-void fill_impl(mp11::mp_int<-1>, IS, std::false_type, T&& t, const U& u) noexcept {
-  mp11::tuple_apply([&t](auto&&... args) { t(args...); }, std::get<IS::value>(u).value);
-}
-
-template <class IW, class IS, class T, class U>
-void fill_impl(IW, IS, std::false_type, T&& t, const U& u) noexcept {
-  mp11::tuple_apply([&](auto&&... args2) { t(std::get<IW::value>(u).value, args2...); },
-                    std::get<IS::value>(u).value);
-}
-
-template <class A, class S, class... Us>
-typename S::iterator fill(A& axes, S& storage, const std::tuple<Us...>& tus) {
-  constexpr auto iws = weight_sample_indices<Us...>();
-  constexpr unsigned n = sizeof...(Us) - (iws.first > -1) - (iws.second > -1);
-  constexpr unsigned i = (iws.first == 0 || iws.second == 0)
-                             ? (iws.first == 1 || iws.second == 1 ? 2 : 1)
-                             : 0;
-  const auto idx = index(has_growing_axis<A>{}, has_multidim_axis<A>{}, axes, storage,
-                         tuple_slice<i, n>(tus));
-  if (idx) {
-    fill_impl(mp11::mp_int<iws.first>{}, mp11::mp_int<iws.second>{},
-              has_operator_preincrement<typename S::value_type>{}, storage[*idx], tus);
-    return storage.begin() + *idx;
-  }
-  return storage.end();
-}
-
-template <class A, class... Us>
-optional_index at(const A& axes, const std::tuple<Us...>& args) {
-  if (axes_rank(axes) != sizeof...(Us))
-    BOOST_THROW_EXCEPTION(std::invalid_argument("number of arguments != histogram rank"));
-  optional_index idx;
-  mp11::mp_for_each<mp11::mp_iota_c<sizeof...(Us)>>([&](auto i) {
-    // axes_get works with static and dynamic axes
-    linearize_index(idx, axis_get<i>(axes),
-                    static_cast<axis::index_type>(std::get<i>(args)));
-  });
-  return idx;
-}
-
-template <class A, class U>
-optional_index at(const A& axes, const U& args) {
-  if (axes_rank(axes) != axes_rank(args))
-    BOOST_THROW_EXCEPTION(std::invalid_argument("number of arguments != histogram rank"));
-  optional_index idx;
-  using std::begin;
-  auto it = begin(args);
-  for_each_axis(axes, [&](const auto& a) {
-    linearize_index(idx, a, static_cast<axis::index_type>(*it++));
-  });
-  return idx;
+template <class... Ts, class Value>
+std::size_t linearize(std::size_t& o, const std::size_t s, const axis::variant<Ts...>& a,
+                      const Value& v) {
+  return axis::visit([&o, &s, &v](const auto& a) { return linearize(o, s, a, v); }, a);
 }
 
 } // namespace detail
 } // namespace histogram
 } // namespace boost
 
-#endif
+#endif // BOOST_HISTOGRAM_DETAIL_LINEARIZE_HPP
diff --git a/include/boost/histogram/detail/non_member_container_access.hpp b/include/boost/histogram/detail/non_member_container_access.hpp
new file mode 100644
index 00000000..729f7b33
--- /dev/null
+++ b/include/boost/histogram/detail/non_member_container_access.hpp
@@ -0,0 +1,70 @@
+// Copyright 2019 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_HISTOGRAM_DETAIL_NON_MEMBER_CONTAINER_ACCESS_HPP
+#define BOOST_HISTOGRAM_DETAIL_NON_MEMBER_CONTAINER_ACCESS_HPP
+
+#if __cpp_lib_nonmember_container_access >= 201411
+
+#include <iterator>
+
+namespace boost {
+namespace histogram {
+namespace detail {
+
+using std::data;
+using std::size;
+
+} // namespace detail
+} // namespace histogram
+} // namespace boost
+
+#else
+
+#include <initializer_list>
+#include <type_traits>
+
+namespace boost {
+namespace histogram {
+namespace detail {
+
+template <class C>
+constexpr auto data(C& c) -> decltype(c.data()) {
+  return c.data();
+}
+
+template <class C>
+constexpr auto data(const C& c) -> decltype(c.data()) {
+  return c.data();
+}
+
+template <class T, std::size_t N>
+constexpr T* data(T (&array)[N]) noexcept {
+  return array;
+}
+
+template <class E>
+constexpr const E* data(std::initializer_list<E> il) noexcept {
+  return il.begin();
+}
+
+template <class C>
+constexpr auto size(const C& c) -> decltype(c.size()) {
+  return c.size();
+}
+
+template <class T, std::size_t N>
+constexpr std::size_t size(const T (&)[N]) noexcept {
+  return N;
+}
+
+} // namespace detail
+} // namespace histogram
+} // namespace boost
+
+#endif
+
+#endif // BOOST_HISTOGRAM_DETAIL_NON_MEMBER_CONTAINER_ACCESS_HPP
diff --git a/include/boost/histogram/detail/optional_index.hpp b/include/boost/histogram/detail/optional_index.hpp
new file mode 100644
index 00000000..cd092409
--- /dev/null
+++ b/include/boost/histogram/detail/optional_index.hpp
@@ -0,0 +1,47 @@
+// Copyright 2019 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_HISTOGRAM_DETAIL_OPTIONAL_INDEX_HPP
+#define BOOST_HISTOGRAM_DETAIL_OPTIONAL_INDEX_HPP
+
+#include <boost/assert.hpp>
+
+namespace boost {
+namespace histogram {
+namespace detail {
+
+// integer with a persistent invalid state, similar to NaN
+struct optional_index {
+  static constexpr auto invalid = ~static_cast<std::size_t>(0);
+
+  optional_index& operator=(const std::size_t x) noexcept {
+    value = x;
+    return *this;
+  }
+
+  optional_index& operator+=(const std::size_t x) noexcept {
+    BOOST_ASSERT(x != invalid);
+    if (value != invalid) { value += x; }
+    return *this;
+  }
+
+  optional_index& operator+=(const optional_index& x) noexcept {
+    if (x.valid()) return operator+=(x.value);
+    value = invalid;
+    return *this;
+  }
+
+  bool valid() const noexcept { return value != invalid; }
+  const std::size_t& operator*() const noexcept { return value; }
+
+  std::size_t value;
+};
+
+} // namespace detail
+} // namespace histogram
+} // namespace boost
+
+#endif
diff --git a/include/boost/histogram/detail/span.hpp b/include/boost/histogram/detail/span.hpp
new file mode 100644
index 00000000..a6a0eb06
--- /dev/null
+++ b/include/boost/histogram/detail/span.hpp
@@ -0,0 +1,249 @@
+// Copyright 2019 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_HISTOGRAM_DETAIL_SPAN_HPP
+#define BOOST_HISTOGRAM_DETAIL_SPAN_HPP
+
+#if __has_include(<span>)
+#include <span>
+
+namespace boost {
+namespace histogram {
+namespace detail {
+using std::span;
+} // namespace detail
+} // namespace histogram
+} // namespace boost
+
+#else
+
+#include <array>
+#include <boost/assert.hpp>
+#include <boost/config.hpp>
+#include <boost/histogram/detail/non_member_container_access.hpp>
+#include <initializer_list>
+#include <iterator>
+#include <type_traits>
+
+namespace boost {
+namespace histogram {
+namespace detail {
+
+static constexpr std::size_t dynamic_extent = ~static_cast<std::size_t>(0);
+
+template <class T, std::size_t N>
+class span_base {
+public:
+  constexpr T* data() noexcept { return begin_; }
+  constexpr const T* data() const noexcept { return begin_; }
+  constexpr std::size_t size() const noexcept { return N; }
+
+protected:
+  constexpr span_base(T* b, std::size_t s) noexcept : begin_(b) {
+    static_assert(N == s, "sizes do not match");
+  }
+  constexpr void set(T* b, std::size_t s) noexcept {
+    begin_ = b;
+    static_assert(N == s, "sizes do not match");
+  }
+
+private:
+  T* begin_;
+};
+
+template <class T>
+class span_base<T, dynamic_extent> {
+public:
+  constexpr T* data() noexcept { return begin_; }
+  constexpr const T* data() const noexcept { return begin_; }
+  constexpr std::size_t size() const noexcept { return size_; }
+
+protected:
+  constexpr span_base(T* b, std::size_t s) noexcept : begin_(b), size_(s) {}
+
+  constexpr void set(T* b, std::size_t s) noexcept {
+    begin_ = b;
+    size_ = s;
+  }
+
+private:
+  T* begin_;
+  std::size_t size_;
+};
+
+template <class T, std::size_t Extent = dynamic_extent>
+class span : public span_base<T, Extent> {
+  using base = span_base<T, Extent>;
+
+public:
+  using element_type = T;
+  using value_type = std::remove_cv_t<T>;
+  using index_type = std::size_t;
+  using difference_type = std::ptrdiff_t;
+  using pointer = T*;
+  using const_pointer = const T*;
+  using reference = T&;
+  using const_reference = const T&;
+  using iterator = pointer;
+  using const_iterator = const_pointer;
+  using reverse_iterator = std::reverse_iterator<iterator>;
+  using const_reverse_iterator = std::reverse_iterator<const_iterator>;
+
+  static constexpr std::size_t extent = Extent;
+
+  template <std::size_t _ = extent,
+            class = std::enable_if_t<(_ == 0 || _ == dynamic_extent)> >
+  constexpr span() noexcept : base(nullptr, 0) {}
+
+  constexpr span(pointer first, pointer last)
+      : span(first, static_cast<std::size_t>(last - first)) {
+    BOOST_ASSERT(extent == dynamic_extent ||
+                 static_cast<difference_type>(extent) == (last - first));
+  }
+
+  constexpr span(pointer ptr, index_type count) : base(ptr, count) {}
+
+  template <std::size_t N>
+  constexpr span(element_type (&arr)[N]) noexcept : span(data(arr), N) {
+    static_assert(extent == dynamic_extent || extent == N, "static sizes do not match");
+  }
+
+  template <std::size_t N,
+            class = std::enable_if_t<(extent == dynamic_extent || extent == N)> >
+  constexpr span(std::array<value_type, N>& arr) noexcept : span(data(arr), N) {}
+
+  template <std::size_t N,
+            class = std::enable_if_t<(extent == dynamic_extent || extent == N)> >
+  constexpr span(const std::array<value_type, N>& arr) noexcept : span(data(arr), N) {}
+
+  template <
+      class Container,
+      class = std::enable_if_t<std::is_convertible<
+          decltype(size(std::declval<Container>()), data(std::declval<Container>())),
+          pointer>::value> >
+  constexpr span(const Container& cont) : span(data(cont), size(cont)) {}
+
+  template <
+      class Container,
+      class = std::enable_if_t<std::is_convertible<
+          decltype(size(std::declval<Container>()), data(std::declval<Container>())),
+          pointer>::value> >
+  constexpr span(Container& cont) : span(data(cont), size(cont)) {}
+
+  template <class U, std::size_t N,
+            class = std::enable_if_t<((extent == dynamic_extent || extent == N) &&
+                                      std::is_convertible<U, element_type>::value)> >
+  constexpr span(const span<U, N>& s) noexcept : span(s.data(), s.size()) {}
+
+  constexpr span(const span& other) noexcept = default;
+
+  constexpr iterator begin() { return base::data(); }
+  constexpr const_iterator begin() const { return base::data(); }
+  constexpr const_iterator cbegin() const { return base::data(); }
+
+  constexpr iterator end() { return base::data() + base::size(); }
+  constexpr const_iterator end() const { return base::data() + base::size(); }
+  constexpr const_iterator cend() const { return base::data() + base::size(); }
+
+  reverse_iterator rbegin() { return reverse_iterator(end()); }
+  const_reverse_iterator rbegin() const { return reverse_iterator(end()); }
+  const_reverse_iterator crbegin() { return reverse_iterator(end()); }
+
+  reverse_iterator rend() { return reverse_iterator(begin()); }
+  const_reverse_iterator rend() const { return reverse_iterator(begin()); }
+  const_reverse_iterator crend() { return reverse_iterator(begin()); }
+
+  constexpr reference front() { *base::data(); }
+  constexpr reference back() { *(base::data() + base::size() - 1); }
+
+  constexpr reference operator[](index_type idx) const { return base::data()[idx]; }
+
+  // constexpr pointer data() const noexcept { return base::data(); }
+  // constexpr std::size_t size() const noexcept { return base::size(); }
+
+  constexpr std::size_t size_bytes() const noexcept {
+    return base::size() * sizeof(element_type);
+  }
+
+#if BOOST_WORKAROUND(BOOST_CLANG, >= 0)
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wc++17-extensions"
+#endif
+
+  BOOST_ATTRIBUTE_NODISCARD constexpr bool empty() const noexcept {
+    return base::size() == 0;
+  }
+
+#if BOOST_WORKAROUND(BOOST_CLANG, >= 0)
+#pragma GCC diagnostic pop
+#endif
+
+  template <std::size_t Count>
+  constexpr span<element_type, Count> first() const {
+    BOOST_ASSERT(Count <= base::size());
+    return span<element_type, Count>(base::data(), Count);
+  }
+
+  constexpr span<element_type, dynamic_extent> first(std::size_t count) const {
+    BOOST_ASSERT(count <= base::size());
+    return span<element_type, dynamic_extent>(base::data(), count);
+  }
+
+  template <std::size_t Count>
+  constexpr span<element_type, Count> last() const {
+    BOOST_ASSERT(Count <= base::size());
+    return span<element_type, Count>(base::data() + base::size() - Count, Count);
+  }
+
+  constexpr span<element_type, dynamic_extent> last(std::size_t count) const {
+    BOOST_ASSERT(count <= base::size());
+    return span<element_type, dynamic_extent>(base::data() + base::size() - count, count);
+  }
+
+  template <std::size_t Offset, std::size_t Count = dynamic_extent>
+  constexpr span<element_type,
+                 (Count != dynamic_extent
+                      ? Count
+                      : (extent != dynamic_extent ? extent - Offset : dynamic_extent))>
+  subspan() const {
+    BOOST_ASSERT(Offset <= base::size());
+    constexpr std::size_t E =
+        (Count != dynamic_extent
+             ? Count
+             : (extent != dynamic_extent ? extent - Offset : dynamic_extent));
+    BOOST_ASSERT(E == dynamic_extent || E <= base::size());
+    return span<element_type, E>(base::data() + Offset,
+                                 Count == dynamic_extent ? base::size() - Offset : Count);
+  }
+
+  constexpr span<element_type, dynamic_extent> subspan(
+      std::size_t offset, std::size_t count = dynamic_extent) const {
+    BOOST_ASSERT(offset <= base::size());
+    const std::size_t s = count == dynamic_extent ? base::size() - offset : count;
+    BOOST_ASSERT(s <= base::size());
+    return span<element_type, dynamic_extent>(base::data() + offset, s);
+  }
+
+private:
+};
+
+#endif
+
+template <class T>
+span<T> make_span(T* begin, T* end) {
+  return span<T>{begin, end};
+}
+
+template <class T>
+span<T> make_span(T* begin, std::size_t size) {
+  return span<T>{begin, size};
+}
+
+} // namespace detail
+} // namespace histogram
+} // namespace boost
+
+#endif
diff --git a/include/boost/histogram/fwd.hpp b/include/boost/histogram/fwd.hpp
index cc042bb7..0dc15fac 100644
--- a/include/boost/histogram/fwd.hpp
+++ b/include/boost/histogram/fwd.hpp
@@ -12,8 +12,9 @@
   Forward declarations, tag types and type aliases.
 */
 
+#include <boost/config.hpp>            // BOOST_ATTRIBUTE_NODISCARD
+#include <boost/config/workaround.hpp> // BOOST_WORKAROUND
 #include <boost/core/use_default.hpp>
-#include <boost/histogram/detail/attribute.hpp> // BOOST_HISTOGRAM_NODISCARD
 #include <vector>
 
 namespace boost {
@@ -119,8 +120,17 @@ using weighted_profile_storage = dense_storage<accumulators::weighted_mean<>>;
 
 #ifndef BOOST_HISTOGRAM_DOXYGEN_INVOKED
 
+#if BOOST_WORKAROUND(BOOST_CLANG, >= 0)
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wc++17-extensions"
+#endif
+
 template <class Axes, class Storage = default_storage>
-class BOOST_HISTOGRAM_NODISCARD histogram;
+class BOOST_ATTRIBUTE_NODISCARD histogram;
+
+#if BOOST_WORKAROUND(BOOST_CLANG, >= 0)
+#pragma GCC diagnostic pop
+#endif
 
 #endif
 } // namespace histogram
diff --git a/include/boost/histogram/histogram.hpp b/include/boost/histogram/histogram.hpp
index 75bc0df4..f7e4d803 100644
--- a/include/boost/histogram/histogram.hpp
+++ b/include/boost/histogram/histogram.hpp
@@ -7,10 +7,13 @@
 #ifndef BOOST_HISTOGRAM_HISTOGRAM_HPP
 #define BOOST_HISTOGRAM_HISTOGRAM_HPP
 
+#include <boost/histogram/detail/at.hpp>
 #include <boost/histogram/detail/axes.hpp>
 #include <boost/histogram/detail/common_type.hpp>
 #include <boost/histogram/detail/compressed_pair.hpp>
-#include <boost/histogram/detail/linearize.hpp>
+#include <boost/histogram/detail/fill.hpp>
+#include <boost/histogram/detail/fill_n.hpp>
+#include <boost/histogram/detail/non_member_container_access.hpp>
 #include <boost/histogram/detail/noop_mutex.hpp>
 #include <boost/histogram/detail/static_if.hpp>
 #include <boost/histogram/fwd.hpp>
@@ -149,41 +152,144 @@ public:
 
   /** Fill histogram with values, an optional weight, and/or a sample.
 
-   Arguments are passed in order to the axis objects. Passing an argument type that is
-   not convertible to the value type accepted by the axis or passing the wrong number
-   of arguments causes a throw of `std::invalid_argument`.
+    Arguments are passed in order to the axis objects. Passing an argument type that is
+    not convertible to the value type accepted by the axis or passing the wrong number
+    of arguments causes a throw of `std::invalid_argument`.
 
-   __Optional weight__
+    __Optional weight__
 
-   An optional weight can be passed as the first or last argument
-   with the [weight](boost/histogram/weight.html) helper function. Compilation fails if
-   the storage elements do not support weights.
+    An optional weight can be passed as the first or last argument
+    with the [weight](boost/histogram/weight.html) helper function. Compilation fails if
+    the storage elements do not support weights.
 
-   __Samples__
+    __Samples__
 
-   If the storage elements accept samples, pass them with the sample helper function
-   in addition to the axis arguments, which can be the first or last argument. The
-   [sample](boost/histogram/sample.html) helper function can pass one or more arguments to
-   the storage element. If samples and weights are used together, they can be passed in
-   any order at the beginning or end of the argument list.
+    If the storage elements accept samples, pass them with the sample helper function
+    in addition to the axis arguments, which can be the first or last argument. The
+    [sample](boost/histogram/sample.html) helper function can pass one or more arguments
+    to the storage element. If samples and weights are used together, they can be passed
+    in any order at the beginning or end of the argument list.
 
-   __Axis with multiple arguments__
+    __Axis with multiple arguments__
 
-   If the histogram contains an axis which accepts a `std::tuple` of arguments, the
-   arguments for that axis need to passed as a `std::tuple`, for example,
-   `std::make_tuple(1.2, 2.3)`. If the histogram contains only this axis and no other,
-   the arguments can be passed directly.
+    If the histogram contains an axis which accepts a `std::tuple` of arguments, the
+    arguments for that axis need to passed as a `std::tuple`, for example,
+    `std::make_tuple(1.2, 2.3)`. If the histogram contains only this axis and no other,
+    the arguments can be passed directly.
   */
-  template <class... Ts>
-  iterator operator()(const Ts&... ts) {
-    return operator()(std::forward_as_tuple(ts...));
+  template <class... Args>
+  iterator operator()(const Args&... args) {
+    return operator()(std::forward_as_tuple(args...));
   }
 
   /// Fill histogram with values, an optional weight, and/or a sample from a `std::tuple`.
   template <class... Ts>
-  iterator operator()(const std::tuple<Ts...>& t) {
+  iterator operator()(const std::tuple<Ts...>& args) {
     std::lock_guard<mutex_type> guard{storage_and_mutex_.second()};
-    return detail::fill(axes_, storage_and_mutex_.first(), t);
+    return detail::fill(storage_and_mutex_.first(), axes_, args);
+  }
+
+  /** Fill histogram with several values at once.
+
+    In general, the argument must be an iterable with the same length as the rank of the
+    histogram. The elements of the iterable may be either identical sub-iterables over
+    values or a variant. The variant may hold sub-iterables over values or single values.
+    All sub-iterables must have the same length.
+
+    Values in sub-iterables are passed to the corresponding histogram axis. A single
+    value is treated as if a sub-iterable was passed that is filled with copies of this
+    value.
+
+    If the histogram has only one axis, a sub-iterable may be passed directly in place of
+    the main iterable.
+
+    @param args iterable as explained in the long description.
+  */
+  template <class Iterable, class = detail::requires_iterable<Iterable>>
+  void fill(const Iterable& args) {
+    std::lock_guard<mutex_type> guard{storage_and_mutex_.second()};
+    detail::fill_n(storage_and_mutex_.first(), axes_, detail::data(args),
+                   detail::size(args));
+  }
+
+  /** Fill histogram with several values and weights at once.
+
+    @param weights single weight or an iterable of weights.
+    @param args iterable of values.
+  */
+  template <class Iterable, class T, class = detail::requires_iterable<Iterable>>
+  void fill(const weight_type<T>& weights, const Iterable& args) {
+    std::lock_guard<mutex_type> guard{storage_and_mutex_.second()};
+    detail::fill_n(storage_and_mutex_.first(), axes_, detail::data(args),
+                   detail::size(args), weights);
+  }
+
+  /** Fill histogram with several values and weights at once.
+
+    @param args iterable of values.
+    @param weights single weight or an iterable of weights.
+  */
+  template <class Iterable, class T, class = detail::requires_iterable<Iterable>>
+  void fill(const Iterable& args, const weight_type<T>& weights) {
+    std::lock_guard<mutex_type> guard{storage_and_mutex_.second()};
+    detail::fill_n(storage_and_mutex_.first(), axes_, detail::data(args),
+                   detail::size(args), weights);
+  }
+
+  /** Fill histogram with several values and samples at once.
+
+    @param samples single sample or an iterable of samples.
+    @param args iterable of values.
+  */
+  template <class Iterable, class T, class = detail::requires_iterable<Iterable>>
+  void fill(const sample_type<T>& samples, const Iterable& args) {
+    std::lock_guard<mutex_type> guard{storage_and_mutex_.second()};
+    detail::fill_n(storage_and_mutex_.first(), axes_, detail::data(args),
+                   detail::size(args), samples);
+  }
+
+  /** Fill histogram with several values and samples at once.
+
+    @param args iterable of values.
+    @param samples single sample or an iterable of samples.
+  */
+  template <class Iterable, class T, class = detail::requires_iterable<Iterable>>
+  void fill(const Iterable& args, const sample_type<T>& samples) {
+    std::lock_guard<mutex_type> guard{storage_and_mutex_.second()};
+    detail::fill_n(storage_and_mutex_.first(), axes_, detail::data(args),
+                   detail::size(args), samples);
+  }
+
+  template <class Iterable, class T, class U, class = detail::requires_iterable<Iterable>>
+  void fill(const sample_type<T>& samples, const weight_type<U>& weights,
+            const Iterable& args) {
+    std::lock_guard<mutex_type> guard{storage_and_mutex_.second()};
+    detail::fill_n(storage_and_mutex_.first(), axes_, detail::data(args),
+                   detail::size(args), weights, samples);
+  }
+
+  template <class Iterable, class T, class U, class = detail::requires_iterable<Iterable>>
+  void fill(const weight_type<T>& weights, const sample_type<U>& samples,
+            const Iterable& args) {
+    std::lock_guard<mutex_type> guard{storage_and_mutex_.second()};
+    detail::fill_n(storage_and_mutex_.first(), axes_, detail::data(args),
+                   detail::size(args), weights, samples);
+  }
+
+  template <class Iterable, class T, class U, class = detail::requires_iterable<Iterable>>
+  void fill(const Iterable& args, const sample_type<T>& samples,
+            const weight_type<U>& weights) {
+    std::lock_guard<mutex_type> guard{storage_and_mutex_.second()};
+    detail::fill_n(storage_and_mutex_.first(), axes_, detail::data(args),
+                   detail::size(args), weights, samples);
+  }
+
+  template <class Iterable, class T, class U, class = detail::requires_iterable<Iterable>>
+  void fill(const Iterable& args, const weight_type<T>& weights,
+            const sample_type<U>& samples) {
+    std::lock_guard<mutex_type> guard{storage_and_mutex_.second()};
+    detail::fill_n(storage_and_mutex_.first(), axes_, detail::data(args),
+                   detail::size(args), weights, samples);
   }
 
   /** Access cell value at integral indices.
@@ -209,10 +315,13 @@ public:
   }
 
   /// Access cell value at integral indices stored in `std::tuple`.
-  template <typename... Indices>
+  template <class... Indices>
   decltype(auto) at(const std::tuple<Indices...>& is) {
+    if (rank() != sizeof...(Indices))
+      BOOST_THROW_EXCEPTION(
+          std::invalid_argument("number of arguments != histogram rank"));
     const auto idx = detail::at(axes_, is);
-    if (!idx)
+    if (!idx.valid())
       BOOST_THROW_EXCEPTION(std::out_of_range("at least one index out of bounds"));
     return storage_and_mutex_.first()[*idx];
   }
@@ -220,8 +329,11 @@ public:
   /// Access cell value at integral indices stored in `std::tuple` (read-only).
   template <typename... Indices>
   decltype(auto) at(const std::tuple<Indices...>& is) const {
+    if (rank() != sizeof...(Indices))
+      BOOST_THROW_EXCEPTION(
+          std::invalid_argument("number of arguments != histogram rank"));
     const auto idx = detail::at(axes_, is);
-    if (!idx)
+    if (!idx.valid())
       BOOST_THROW_EXCEPTION(std::out_of_range("at least one index out of bounds"));
     return storage_and_mutex_.first()[*idx];
   }
@@ -229,8 +341,11 @@ public:
   /// Access cell value at integral indices stored in iterable.
   template <class Iterable, class = detail::requires_iterable<Iterable>>
   decltype(auto) at(const Iterable& is) {
+    if (rank() != detail::axes_rank(is))
+      BOOST_THROW_EXCEPTION(
+          std::invalid_argument("number of arguments != histogram rank"));
     const auto idx = detail::at(axes_, is);
-    if (!idx)
+    if (!idx.valid())
       BOOST_THROW_EXCEPTION(std::out_of_range("at least one index out of bounds"));
     return storage_and_mutex_.first()[*idx];
   }
@@ -238,8 +353,11 @@ public:
   /// Access cell value at integral indices stored in iterable (read-only).
   template <class Iterable, class = detail::requires_iterable<Iterable>>
   decltype(auto) at(const Iterable& is) const {
+    if (rank() != detail::axes_rank(is))
+      BOOST_THROW_EXCEPTION(
+          std::invalid_argument("number of arguments != histogram rank"));
     const auto idx = detail::at(axes_, is);
-    if (!idx)
+    if (!idx.valid())
       BOOST_THROW_EXCEPTION(std::out_of_range("at least one index out of bounds"));
     return storage_and_mutex_.first()[*idx];
   }
diff --git a/include/boost/histogram/indexed.hpp b/include/boost/histogram/indexed.hpp
index 3b4ec823..299d37cc 100644
--- a/include/boost/histogram/indexed.hpp
+++ b/include/boost/histogram/indexed.hpp
@@ -8,8 +8,8 @@
 #define BOOST_HISTOGRAM_INDEXED_HPP
 
 #include <array>
+#include <boost/config.hpp>
 #include <boost/histogram/axis/traits.hpp>
-#include <boost/histogram/detail/attribute.hpp>
 #include <boost/histogram/detail/axes.hpp>
 #include <boost/histogram/detail/iterator_adaptor.hpp>
 #include <boost/histogram/detail/operators.hpp>
@@ -29,6 +29,11 @@ enum class coverage {
   all,   /*!< iterate over all bins, including underflow and overflow */
 };
 
+#if BOOST_WORKAROUND(BOOST_CLANG, >= 0)
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wc++17-extensions"
+#endif
+
 /** Input iterator range over histogram bins with multi-dimensional index.
 
   The iterator returned by begin() can only be incremented. begin() may only be called
@@ -36,7 +41,7 @@ enum class coverage {
   input iterators exist, the other copies become invalid if one of them is incremented.
 */
 template <class Histogram>
-class BOOST_HISTOGRAM_NODISCARD indexed_range {
+class BOOST_ATTRIBUTE_NODISCARD indexed_range {
 private:
   using histogram_type = Histogram;
   static constexpr std::size_t buffer_size =
@@ -54,11 +59,10 @@ public:
 
     Its methods provide access to the current indices and bins and it acts like a pointer
     to the cell value. To interoperate with the algorithms of the standard library, the
-    accessor is implicitly convertible to a cell value. Assignments and comparisons
-    are passed through to the cell. The accessor is coupled to its parent
-    iterator. Moving the parent iterator forward also updates the linked
-    accessor. Accessors are not copyable. They cannot be stored in containers, but
-    range_iterators can be stored.
+    accessor is implicitly convertible to a cell value. Assignments and comparisons are
+    passed through to the cell. The accessor is coupled to its parent iterator. Moving the
+    parent iterator forward also updates the linked accessor. Accessors are not copyable.
+    They cannot be stored in containers, but range_iterators can be stored.
   */
   class accessor : detail::mirrored<accessor, void> {
   public:
@@ -316,6 +320,10 @@ private:
   iterator begin_, end_;
 };
 
+#if BOOST_WORKAROUND(BOOST_CLANG, >= 0)
+#pragma GCC diagnostic pop
+#endif
+
 /** Generates an indexed range of <a
   href="https://en.cppreference.com/w/cpp/named_req/ForwardIterator">forward iterators</a>
   over the histogram cells.
diff --git a/include/boost/histogram/ostream.hpp b/include/boost/histogram/ostream.hpp
index e38a0600..a4ed78d1 100644
--- a/include/boost/histogram/ostream.hpp
+++ b/include/boost/histogram/ostream.hpp
@@ -32,13 +32,10 @@ template <typename CharT, typename Traits, typename A, typename S>
 std::basic_ostream<CharT, Traits>& operator<<(std::basic_ostream<CharT, Traits>& os,
                                               const histogram<A, S>& h) {
   os << "histogram(";
-  unsigned n = 0;
-  h.for_each_axis([&](const auto& a) {
-    if (h.rank() > 1) os << "\n  ";
-    os << a;
-    if (++n < h.rank()) os << ",";
-  });
-  os << (h.rank() > 1 ? "\n)" : ")");
+  h.for_each_axis([&](const auto& a) { os << "\n  " << a << ","; });
+  std::size_t i = 0;
+  for (auto&& x : h) os << "\n  " << i++ << ": " << x;
+  os << (h.rank() ? "\n)" : ")");
   return os;
 }
 
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 315b5d8b..50665ecb 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -2,7 +2,10 @@
 # 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
 
-# keep in sync with Jamfile
+boost_test(SOURCES check_cmake_version.cpp ARGUMENTS ${PROJECT_VERSION}
+  LIBRARIES Boost::core Boost::config)
+
+# keep in sync with Jamfile, this should be automatized...
 boost_test(TYPE compile-fail SOURCES make_histogram_fail0.cpp)
 boost_test(TYPE compile-fail SOURCES make_histogram_fail1.cpp)
 boost_test(TYPE run SOURCES algorithm_project_test.cpp
@@ -43,24 +46,28 @@ boost_test(TYPE run SOURCES detail_make_default_test.cpp
   LIBRARIES Boost::histogram Boost::core)
 boost_test(TYPE run SOURCES detail_misc_test.cpp
   LIBRARIES Boost::histogram Boost::core)
-boost_test(TYPE run SOURCES detail_meta_test.cpp
-  LIBRARIES Boost::histogram Boost::core)
 boost_test(TYPE run SOURCES detail_large_int_test.cpp
   LIBRARIES Boost::histogram Boost::core)
-boost_test(TYPE run SOURCES detail_linearize_test.cpp
-  LIBRARIES Boost::histogram Boost::core)
 boost_test(TYPE run SOURCES detail_iterator_adaptor_test.cpp
   LIBRARIES Boost::histogram Boost::core)
+boost_test(TYPE run SOURCES detail_operators_test.cpp
+  LIBRARIES Boost::histogram Boost::core)
 boost_test(TYPE run SOURCES detail_relaxed_equal_test.cpp
   LIBRARIES Boost::histogram Boost::core)
 boost_test(TYPE run SOURCES detail_replace_default_test.cpp
   LIBRARIES Boost::histogram Boost::core)
 boost_test(TYPE run SOURCES detail_safe_comparison_test.cpp
   LIBRARIES Boost::histogram Boost::core)
+boost_test(TYPE run SOURCES detail_static_if_test.cpp
+  LIBRARIES Boost::histogram Boost::core)
 boost_test(TYPE run SOURCES detail_tuple_slice_test.cpp
   LIBRARIES Boost::histogram Boost::core)
+boost_test(TYPE run SOURCES histogram_custom_axis_test.cpp
+  LIBRARIES Boost::histogram Boost::core)
 boost_test(TYPE run SOURCES histogram_dynamic_test.cpp
   LIBRARIES Boost::histogram Boost::core)
+boost_test(TYPE run SOURCES histogram_fill_test.cpp
+  LIBRARIES Boost::histogram Boost::core)
 boost_test(TYPE run SOURCES histogram_growing_test.cpp
   LIBRARIES Boost::histogram Boost::core)
 boost_test(TYPE run SOURCES histogram_mixed_test.cpp
@@ -105,5 +112,3 @@ endif()
 # boost_test(TYPE run SOURCES histogram_serialization_test.cpp LIBRARIES Boost::histogram Boost::core Boost::serialization)
 # boost_test(TYPE run SOURCES axis_variant_serialization_test.cpp
 #   LIBRARIES Boost::histogram Boost::core Boost::serialization)
-
-boost_test(SOURCES check_cmake_version.cpp ARGUMENTS ${PROJECT_VERSION} LIBRARIES Boost::core Boost::config)
diff --git a/test/Jamfile b/test/Jamfile
index 77f0a735..0d31be10 100644
--- a/test/Jamfile
+++ b/test/Jamfile
@@ -40,17 +40,18 @@ alias cxx14 :
     [ run detail_detect_test.cpp ]
     [ run detail_limits_test.cpp ]
     [ run detail_make_default_test.cpp ]
-    [ run detail_meta_test.cpp ]
     [ run detail_misc_test.cpp ]
     [ run detail_iterator_adaptor_test.cpp ]
     [ run detail_large_int_test.cpp ]
-    [ run detail_linearize_test.cpp ]
     [ run detail_operators_test.cpp ]
     [ run detail_relaxed_equal_test.cpp ]
     [ run detail_replace_default_test.cpp ]
     [ run detail_safe_comparison_test.cpp ]
+    [ run detail_static_if_test.cpp ]
     [ run detail_tuple_slice_test.cpp ]
+    [ run histogram_custom_axis_test.cpp ]
     [ run histogram_dynamic_test.cpp ]
+    [ run histogram_fill_test.cpp ]
     [ run histogram_growing_test.cpp ]
     [ run histogram_mixed_test.cpp ]
     [ run histogram_operators_test.cpp ]
diff --git a/test/detail_axes_test.cpp b/test/detail_axes_test.cpp
index b2a1f781..74a44704 100644
--- a/test/detail_axes_test.cpp
+++ b/test/detail_axes_test.cpp
@@ -14,11 +14,10 @@
 #include <tuple>
 #include <utility>
 #include <vector>
-#include "throw_exception.hpp"
 #include "std_ostream.hpp"
+#include "throw_exception.hpp"
 
 using namespace boost::histogram;
-using namespace boost::histogram::detail;
 
 int main() {
   // dynamic axis_get with tuples
@@ -134,5 +133,27 @@ int main() {
     BOOST_TEST_THROWS(detail::bincount(v), std::overflow_error);
   }
 
+  {
+    struct growing {
+      auto update(int) { return std::make_pair(0, 0); }
+    };
+    using T = growing;
+    using I = axis::integer<>;
+
+    using A = std::tuple<I, T>;
+    using B = std::vector<T>;
+    using C = std::vector<axis::variant<I, T>>;
+    using D = std::tuple<I>;
+    using E = std::vector<I>;
+    using F = std::vector<axis::variant<I>>;
+
+    BOOST_TEST_TRAIT_TRUE((detail::has_growing_axis<A>));
+    BOOST_TEST_TRAIT_TRUE((detail::has_growing_axis<B>));
+    BOOST_TEST_TRAIT_TRUE((detail::has_growing_axis<C>));
+    BOOST_TEST_TRAIT_FALSE((detail::has_growing_axis<D>));
+    BOOST_TEST_TRAIT_FALSE((detail::has_growing_axis<E>));
+    BOOST_TEST_TRAIT_FALSE((detail::has_growing_axis<F>));
+  }
+
   return boost::report_errors();
 }
diff --git a/test/detail_detect_test.cpp b/test/detail_detect_test.cpp
index 138a2ed1..83252f3d 100644
--- a/test/detail_detect_test.cpp
+++ b/test/detail_detect_test.cpp
@@ -21,9 +21,9 @@
 #include <unordered_map>
 #include <utility>
 #include <vector>
+#include "std_ostream.hpp"
 #include "throw_exception.hpp"
 #include "utility_allocator.hpp"
-#include "std_ostream.hpp"
 
 using namespace boost::histogram;
 using namespace boost::histogram::detail;
diff --git a/test/detail_linearize_test.cpp b/test/detail_linearize_test.cpp
deleted file mode 100644
index dd16a7dd..00000000
--- a/test/detail_linearize_test.cpp
+++ /dev/null
@@ -1,64 +0,0 @@
-// Copyright 2015-2017 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)
-
-#include <boost/core/lightweight_test.hpp>
-#include <boost/core/lightweight_test_trait.hpp>
-#include <boost/histogram/axis/integer.hpp>
-#include <boost/histogram/axis/variant.hpp>
-#include <boost/histogram/detail/linearize.hpp>
-#include <tuple>
-#include <utility>
-#include <vector>
-#include "std_ostream.hpp"
-
-using namespace boost::histogram;
-
-int main() {
-  {
-    struct growing {
-      auto update(int) { return std::make_pair(0, 0); }
-    };
-    using T = growing;
-    using I = axis::integer<>;
-
-    using A = std::tuple<I, T>;
-    using B = std::vector<T>;
-    using C = std::vector<axis::variant<I, T>>;
-    using D = std::tuple<I>;
-    using E = std::vector<I>;
-    using F = std::vector<axis::variant<I>>;
-
-    BOOST_TEST_TRAIT_TRUE((detail::has_growing_axis<A>));
-    BOOST_TEST_TRAIT_TRUE((detail::has_growing_axis<B>));
-    BOOST_TEST_TRAIT_TRUE((detail::has_growing_axis<C>));
-    BOOST_TEST_TRAIT_FALSE((detail::has_growing_axis<D>));
-    BOOST_TEST_TRAIT_FALSE((detail::has_growing_axis<E>));
-    BOOST_TEST_TRAIT_FALSE((detail::has_growing_axis<F>));
-  }
-
-  {
-    struct multidim {
-      auto index(const std::tuple<int>&) { return 0; }
-    };
-    using T = multidim;
-    using I = axis::integer<>;
-
-    using A = std::tuple<I, T>;
-    using B = std::vector<T>;
-    using C = std::vector<axis::variant<I, T>>;
-    using D = std::tuple<I>;
-    using E = std::vector<I>;
-    using F = std::vector<axis::variant<I>>;
-
-    BOOST_TEST_TRAIT_TRUE((detail::has_multidim_axis<A>));
-    BOOST_TEST_TRAIT_TRUE((detail::has_multidim_axis<B>));
-    BOOST_TEST_TRAIT_TRUE((detail::has_multidim_axis<C>));
-    BOOST_TEST_TRAIT_FALSE((detail::has_multidim_axis<D>));
-    BOOST_TEST_TRAIT_FALSE((detail::has_multidim_axis<E>));
-    BOOST_TEST_TRAIT_FALSE((detail::has_multidim_axis<F>));
-  }
-  return boost::report_errors();
-}
diff --git a/test/detail_meta_test.cpp b/test/detail_meta_test.cpp
deleted file mode 100644
index 12f61542..00000000
--- a/test/detail_meta_test.cpp
+++ /dev/null
@@ -1,367 +0,0 @@
-// Copyright 2015-2017 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)
-
-#include <array>
-#include <boost/core/lightweight_test.hpp>
-#include <boost/core/lightweight_test_trait.hpp>
-#include <boost/histogram/axis/integer.hpp>
-#include <boost/histogram/axis/regular.hpp>
-#include <boost/histogram/axis/variant.hpp>
-#include <boost/histogram/detail/args_type.hpp>
-#include <boost/histogram/detail/make_default.hpp>
-#include <boost/histogram/histogram.hpp>
-#include <boost/histogram/unlimited_storage.hpp>
-#include <deque>
-#include <iterator>
-#include <map>
-#include <tuple>
-#include <type_traits>
-#include <unordered_map>
-#include <utility>
-#include <vector>
-#include "throw_exception.hpp"
-#include "utility_allocator.hpp"
-#include "std_ostream.hpp"
-
-using namespace boost::histogram;
-using namespace boost::histogram::detail;
-
-int main() {
-  // has_method_value*
-  {
-    struct A {};
-    struct B {
-      A value(int) const { return {}; }
-    };
-    struct C {
-      char value(int) const { return 0; }
-    };
-
-    BOOST_TEST_TRAIT_FALSE((has_method_value<A>));
-    BOOST_TEST_TRAIT_TRUE((has_method_value<B>));
-    BOOST_TEST_TRAIT_TRUE((has_method_value<C>));
-  }
-
-  // has_method_options
-  {
-    struct A {};
-    struct B {
-      void options() {}
-    };
-
-    BOOST_TEST_TRAIT_FALSE((has_method_options<A>));
-    BOOST_TEST_TRAIT_TRUE((has_method_options<B>));
-  }
-
-  // has_method_metadata
-  {
-    struct A {};
-    struct B {
-      void metadata();
-    };
-
-    BOOST_TEST_TRAIT_FALSE((has_method_metadata<A>));
-    BOOST_TEST_TRAIT_TRUE((has_method_metadata<B>));
-  }
-
-  // has_method_update
-  {
-    struct A {};
-    struct B {
-      void update(int) {}
-    };
-    using C = axis::integer<int, axis::null_type, use_default>;
-
-    BOOST_TEST_TRAIT_FALSE((has_method_update<A>));
-    BOOST_TEST_TRAIT_TRUE((has_method_update<B>));
-    BOOST_TEST_TRAIT_TRUE((has_method_update<C>));
-  }
-
-  // has_method_resize
-  {
-    struct A {};
-    using B = std::vector<int>;
-    using C = std::map<int, int>;
-
-    BOOST_TEST_TRAIT_FALSE((has_method_resize<A>));
-    BOOST_TEST_TRAIT_TRUE((has_method_resize<B>));
-    BOOST_TEST_TRAIT_FALSE((has_method_resize<C>));
-  }
-
-  // has_method_size
-  {
-    struct A {};
-    using B = std::vector<int>;
-    using C = std::map<int, int>;
-
-    BOOST_TEST_TRAIT_FALSE((has_method_size<A>));
-    BOOST_TEST_TRAIT_TRUE((has_method_size<B>));
-    BOOST_TEST_TRAIT_TRUE((has_method_size<C>));
-  }
-
-  // has_method_clear
-  {
-    struct A {};
-    using B = std::vector<int>;
-    using C = std::map<int, int>;
-    using D = std::array<int, 10>;
-
-    BOOST_TEST_TRAIT_FALSE((has_method_clear<A>));
-    BOOST_TEST_TRAIT_TRUE((has_method_clear<B>));
-    BOOST_TEST_TRAIT_TRUE((has_method_clear<C>));
-    BOOST_TEST_TRAIT_FALSE((has_method_clear<D>));
-  }
-
-  // has_allocator
-  {
-    struct A {};
-    using B = std::vector<int>;
-    using C = std::map<int, int>;
-    using D = std::array<int, 10>;
-
-    BOOST_TEST_TRAIT_FALSE((has_method_clear<A>));
-    BOOST_TEST_TRAIT_TRUE((has_method_clear<B>));
-    BOOST_TEST_TRAIT_TRUE((has_method_clear<C>));
-    BOOST_TEST_TRAIT_FALSE((has_method_clear<D>));
-  }
-
-  // is_storage
-  {
-    struct A {};
-    using B = std::vector<int>;
-    using C = unlimited_storage<>;
-
-    BOOST_TEST_TRAIT_FALSE((is_storage<A>));
-    BOOST_TEST_TRAIT_FALSE((is_storage<B>));
-    BOOST_TEST_TRAIT_TRUE((is_storage<C>));
-  }
-
-  // is_indexable
-  {
-    struct A {};
-    using B = std::vector<int>;
-    using C = std::map<int, int>;
-    using D = std::map<A, int>;
-
-    BOOST_TEST_TRAIT_FALSE((is_indexable<A>));
-    BOOST_TEST_TRAIT_TRUE((is_indexable<B>));
-    BOOST_TEST_TRAIT_TRUE((is_indexable<C>));
-    BOOST_TEST_TRAIT_FALSE((is_indexable<D>));
-  }
-
-  // is_transform
-  {
-    struct A {};
-    struct B {
-      double forward(double);
-      double inverse(double);
-    };
-
-    BOOST_TEST_TRAIT_FALSE((is_transform<A>));
-    BOOST_TEST_TRAIT_TRUE((is_transform<B>));
-  }
-
-  // is_vector_like
-  {
-    struct A {};
-    using B = std::vector<int>;
-    using C = std::array<int, 10>;
-    using D = std::map<unsigned, int>;
-    using E = std::deque<int>;
-    BOOST_TEST_TRAIT_FALSE((is_vector_like<A>));
-    BOOST_TEST_TRAIT_TRUE((is_vector_like<B>));
-    BOOST_TEST_TRAIT_FALSE((is_vector_like<C>));
-    BOOST_TEST_TRAIT_FALSE((is_vector_like<D>));
-    BOOST_TEST_TRAIT_TRUE((is_vector_like<E>));
-  }
-
-  // is_array_like
-  {
-    struct A {};
-    using B = std::vector<int>;
-    using C = std::array<int, 10>;
-    using D = std::map<unsigned, int>;
-    BOOST_TEST_TRAIT_FALSE((is_array_like<A>));
-    BOOST_TEST_TRAIT_FALSE((is_array_like<B>));
-    BOOST_TEST_TRAIT_TRUE((is_array_like<C>));
-    BOOST_TEST_TRAIT_FALSE((is_array_like<D>));
-  }
-
-  // is_map_like
-  {
-    struct A {};
-    using B = std::vector<int>;
-    using C = std::array<int, 10>;
-    using D = std::map<unsigned, int>;
-    using E = std::unordered_map<unsigned, int>;
-    BOOST_TEST_TRAIT_FALSE((is_map_like<A>));
-    BOOST_TEST_TRAIT_FALSE((is_map_like<B>));
-    BOOST_TEST_TRAIT_FALSE((is_map_like<C>));
-    BOOST_TEST_TRAIT_TRUE((is_map_like<D>));
-    BOOST_TEST_TRAIT_TRUE((is_map_like<E>));
-  }
-
-  // is_axis
-  {
-    struct A {};
-    struct B {
-      int index(double);
-      int size() const;
-    };
-    struct C {
-      int index(double);
-    };
-    struct D {
-      int size();
-    };
-    using E = axis::variant<axis::regular<>>;
-
-    BOOST_TEST_TRAIT_FALSE((is_axis<A>));
-    BOOST_TEST_TRAIT_TRUE((is_axis<B>));
-    BOOST_TEST_TRAIT_FALSE((is_axis<C>));
-    BOOST_TEST_TRAIT_FALSE((is_axis<D>));
-    BOOST_TEST_TRAIT_FALSE((is_axis<E>));
-  }
-
-  // is_iterable
-  {
-    using A = std::vector<int>;
-    using B = int[3];
-    using C = std::initializer_list<int>;
-    BOOST_TEST_TRAIT_FALSE((is_iterable<int>));
-    BOOST_TEST_TRAIT_TRUE((is_iterable<A>));
-    BOOST_TEST_TRAIT_TRUE((is_iterable<B>));
-    BOOST_TEST_TRAIT_TRUE((is_iterable<C>));
-  }
-
-  // is_streamable
-  {
-    struct Foo {};
-    BOOST_TEST_TRAIT_TRUE((is_streamable<int>));
-    BOOST_TEST_TRAIT_TRUE((is_streamable<std::string>));
-    BOOST_TEST_TRAIT_FALSE((is_streamable<Foo>));
-  }
-
-  // is_axis_variant
-  {
-    struct A {};
-    BOOST_TEST_TRAIT_FALSE((is_axis_variant<A>));
-    BOOST_TEST_TRAIT_TRUE((is_axis_variant<axis::variant<>>));
-    BOOST_TEST_TRAIT_TRUE((is_axis_variant<axis::variant<axis::regular<>>>));
-  }
-
-  // args_type, arg_type
-  {
-    struct Foo {
-      static int f1(char);
-      int f2(long) const;
-    };
-
-    BOOST_TEST_TRAIT_SAME(args_type<decltype(&Foo::f1)>, std::tuple<char>);
-    BOOST_TEST_TRAIT_SAME(args_type<decltype(&Foo::f2)>, std::tuple<long>);
-    BOOST_TEST_TRAIT_SAME(arg_type<decltype(&Foo::f1)>, char);
-    BOOST_TEST_TRAIT_SAME(arg_type<decltype(&Foo::f2)>, long);
-  }
-
-  // static_if
-  {
-    struct callable {
-      int operator()() { return 1; };
-    };
-    struct not_callable {};
-    auto fcn = [](auto b, auto x) {
-      return static_if<decltype(b)>([](auto x) { return x(); }, [](auto) { return 2; },
-                                    x);
-    };
-    BOOST_TEST_EQ(fcn(std::true_type(), callable()), 1);
-    BOOST_TEST_EQ(fcn(std::false_type(), not_callable()), 2);
-  }
-
-  // tuple_slice
-  {
-    auto a = std::make_tuple(1, 2, 3, 4);
-    auto b = tuple_slice<1, 2>(a);
-    BOOST_TEST_EQ(b, std::make_tuple(2, 3));
-  }
-
-  // is_sequence_of_axis
-  {
-    using A = std::vector<axis::regular<>>;
-    using B = std::vector<axis::variant<axis::regular<>>>;
-    using C = std::vector<int>;
-    auto v = std::vector<axis::variant<axis::regular<>, axis::integer<>>>();
-    BOOST_TEST_TRAIT_TRUE((is_sequence_of_any_axis<A>));
-    BOOST_TEST_TRAIT_TRUE((is_sequence_of_axis<A>));
-    BOOST_TEST_TRAIT_FALSE((is_sequence_of_axis_variant<A>));
-    BOOST_TEST_TRAIT_TRUE((is_sequence_of_any_axis<B>));
-    BOOST_TEST_TRAIT_TRUE((is_sequence_of_axis_variant<B>));
-    BOOST_TEST_TRAIT_FALSE((is_sequence_of_axis<B>));
-    BOOST_TEST_TRAIT_FALSE((is_sequence_of_any_axis<C>));
-    BOOST_TEST_TRAIT_TRUE((is_sequence_of_any_axis<decltype(v)>));
-  }
-
-  // is_weight
-  {
-    struct A {};
-    using B = int;
-    using C = decltype(weight(1));
-    BOOST_TEST_TRAIT_FALSE((is_weight<A>));
-    BOOST_TEST_TRAIT_FALSE((is_weight<B>));
-    BOOST_TEST_TRAIT_TRUE((is_weight<C>));
-  }
-
-  // is_sample
-  {
-    struct A {};
-    using B = int;
-    using C = decltype(sample(1));
-    using D = decltype(sample(1, 2.0));
-    BOOST_TEST_TRAIT_FALSE((is_sample<A>));
-    BOOST_TEST_TRAIT_FALSE((is_sample<B>));
-    BOOST_TEST_TRAIT_TRUE((is_sample<C>));
-    BOOST_TEST_TRAIT_TRUE((is_sample<D>));
-  }
-
-  // make_default
-  {
-    struct A {};
-    auto a = make_default(A());
-    BOOST_TEST_TRAIT_SAME(decltype(a), A);
-    tracing_allocator_db db;
-    using B = std::vector<int, tracing_allocator<int>>;
-    B b = make_default(B(tracing_allocator<int>(db)));
-    b.resize(100);
-    BOOST_TEST_EQ(db.at<int>().first, 100);
-  }
-
-  // has_operator_equal
-  {
-    struct A {};
-    struct B {
-      bool operator==(const B&) const { return true; }
-    };
-
-    BOOST_TEST_TRAIT_FALSE((has_operator_equal<A, A>));
-    BOOST_TEST_TRAIT_FALSE((has_operator_equal<B, A>));
-    BOOST_TEST_TRAIT_TRUE((has_operator_equal<B, B>));
-    BOOST_TEST_TRAIT_TRUE((has_operator_equal<const B&, const B&>));
-  }
-
-  // has_operator_radd
-  {
-    struct A {};
-    struct B {
-      B& operator+=(const B&) { return *this; }
-    };
-
-    BOOST_TEST_TRAIT_FALSE((has_operator_radd<A, A>));
-    BOOST_TEST_TRAIT_FALSE((has_operator_radd<B, A>));
-    BOOST_TEST_TRAIT_TRUE((has_operator_radd<B, B>));
-    BOOST_TEST_TRAIT_TRUE((has_operator_radd<B&, const B&>));
-  }
-
-  return boost::report_errors();
-}
diff --git a/test/detail_misc_test.cpp b/test/detail_misc_test.cpp
index 7d23fd1e..12f069fc 100644
--- a/test/detail_misc_test.cpp
+++ b/test/detail_misc_test.cpp
@@ -9,6 +9,7 @@
 #include <boost/histogram/accumulators/weighted_sum.hpp>
 #include <boost/histogram/detail/cat.hpp>
 #include <boost/histogram/detail/common_type.hpp>
+#include <boost/histogram/detail/non_member_container_access.hpp>
 #include <boost/histogram/fwd.hpp>
 #include <boost/histogram/literals.hpp>
 #include <boost/histogram/storage_adaptor.hpp>
@@ -16,10 +17,12 @@
 #include "std_ostream.hpp"
 
 using namespace boost::histogram;
+using namespace boost::histogram::detail;
 using namespace boost::histogram::literals;
 
 int main() {
-  BOOST_TEST_EQ(detail::cat("foo", 1, "bar"), "foo1bar");
+  // cat
+  { BOOST_TEST_EQ(cat("foo", 1, "bar"), "foo1bar"); }
 
   // literals
   {
@@ -31,26 +34,44 @@ int main() {
 
   // common_storage
   {
-    BOOST_TEST_TRAIT_SAME(
-        detail::common_storage<unlimited_storage<>, unlimited_storage<>>,
-        unlimited_storage<>);
-    BOOST_TEST_TRAIT_SAME(
-        detail::common_storage<dense_storage<double>, dense_storage<double>>,
-        dense_storage<double>);
-    BOOST_TEST_TRAIT_SAME(
-        detail::common_storage<dense_storage<int>, dense_storage<double>>,
-        dense_storage<double>);
-    BOOST_TEST_TRAIT_SAME(
-        detail::common_storage<dense_storage<double>, dense_storage<int>>,
-        dense_storage<double>);
-    BOOST_TEST_TRAIT_SAME(
-        detail::common_storage<dense_storage<double>, unlimited_storage<>>,
-        dense_storage<double>);
-    BOOST_TEST_TRAIT_SAME(detail::common_storage<dense_storage<int>, unlimited_storage<>>,
+    BOOST_TEST_TRAIT_SAME(common_storage<unlimited_storage<>, unlimited_storage<>>,
                           unlimited_storage<>);
-    BOOST_TEST_TRAIT_SAME(detail::common_storage<dense_storage<double>, weight_storage>,
+    BOOST_TEST_TRAIT_SAME(common_storage<dense_storage<double>, dense_storage<double>>,
+                          dense_storage<double>);
+    BOOST_TEST_TRAIT_SAME(common_storage<dense_storage<int>, dense_storage<double>>,
+                          dense_storage<double>);
+    BOOST_TEST_TRAIT_SAME(common_storage<dense_storage<double>, dense_storage<int>>,
+                          dense_storage<double>);
+    BOOST_TEST_TRAIT_SAME(common_storage<dense_storage<double>, unlimited_storage<>>,
+                          dense_storage<double>);
+    BOOST_TEST_TRAIT_SAME(common_storage<dense_storage<int>, unlimited_storage<>>,
+                          unlimited_storage<>);
+    BOOST_TEST_TRAIT_SAME(common_storage<dense_storage<double>, weight_storage>,
                           weight_storage);
   }
 
+  // size & data
+  {
+    char a[4] = {1, 2, 3, 4};
+    BOOST_TEST_EQ(size(a), 4u);
+    BOOST_TEST_EQ(data(a), a);
+    auto b = {1, 2};
+    BOOST_TEST_EQ(size(b), 2u);
+    BOOST_TEST_EQ(data(b), b.begin());
+    struct C {
+      unsigned size() const { return 3; }
+      int* data() { return buf; }
+      const int* data() const { return buf; }
+      int buf[1];
+    } c;
+    BOOST_TEST_EQ(size(c), 3u);
+    BOOST_TEST_EQ(data(c), c.buf);
+    BOOST_TEST_EQ(data(static_cast<const C&>(c)), c.buf);
+    struct {
+      int size() const { return 5; }
+    } d;
+    BOOST_TEST_EQ(size(d), 5u);
+  }
+
   return boost::report_errors();
 }
diff --git a/test/detail_static_if_test.cpp b/test/detail_static_if_test.cpp
new file mode 100644
index 00000000..d25d7ed7
--- /dev/null
+++ b/test/detail_static_if_test.cpp
@@ -0,0 +1,33 @@
+// Copyright 2015-2017 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)
+
+#include <boost/core/lightweight_test.hpp>
+#include <boost/histogram/detail/static_if.hpp>
+#include "throw_exception.hpp"
+
+using namespace boost::histogram::detail;
+
+int main() {
+  using T = std::true_type;
+  using F = std::false_type;
+
+  // check that branch not taken does not have to compile
+  BOOST_TEST_EQ(static_if<T>([](auto) { return 1; }, [] {}, 0), 1);
+  BOOST_TEST_EQ(static_if<F>([] {}, [](auto) { return 1; }, 0), 1);
+
+  BOOST_TEST_EQ(static_if_c<true>([](auto) { return 1; }, [] {}, 0), 1);
+  BOOST_TEST_EQ(static_if_c<false>([] {}, [](auto) { return 1; }, 0), 1);
+
+  // check that noexcept is correctly propagated
+  auto may_throw = [](auto x) { return x; };
+  auto no_throw = [](auto x) noexcept { return x; };
+
+  // make this work with -fno-exceptions
+  BOOST_TEST_EQ(noexcept(static_if<F>(no_throw, may_throw, 0)), noexcept(may_throw(0)));
+  BOOST_TEST(noexcept(static_if<T>(no_throw, may_throw, 0)));
+
+  return boost::report_errors();
+}
diff --git a/test/histogram_custom_axis_test.cpp b/test/histogram_custom_axis_test.cpp
new file mode 100644
index 00000000..d31b3375
--- /dev/null
+++ b/test/histogram_custom_axis_test.cpp
@@ -0,0 +1,132 @@
+// Copyright 2018-2019 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)
+
+#include <boost/core/lightweight_test.hpp>
+#include <boost/histogram/accumulators/mean.hpp>
+#include <boost/histogram/accumulators/weighted_mean.hpp>
+#include <boost/histogram/axis.hpp>
+#include <boost/histogram/axis/ostream.hpp>
+#include <boost/histogram/histogram.hpp>
+#include <boost/histogram/literals.hpp>
+#include <boost/histogram/make_histogram.hpp>
+#include <boost/histogram/ostream.hpp>
+#include <boost/histogram/storage_adaptor.hpp>
+#include <sstream>
+#include <stdexcept>
+#include <tuple>
+#include <utility>
+#include <vector>
+#include "throw_exception.hpp"
+#include "utility_histogram.hpp"
+
+using namespace boost::histogram;
+
+struct modified_axis : public axis::integer<> {
+  using integer::integer; // inherit ctors of base
+  // customization point: convert argument and call base class
+  auto index(const char* s) const { return integer::index(std::atoi(s)); }
+};
+
+struct minimal {
+  auto index(int x) const { return static_cast<axis::index_type>(x % 2); }
+  auto size() const { return axis::index_type{2}; }
+};
+
+struct axis2d {
+  auto index(const std::tuple<double, double>& x) const {
+    return axis::index_type{std::get<0>(x) == 1 && std::get<1>(x) == 2};
+  }
+  auto size() const { return axis::index_type{2}; }
+};
+
+class axis2d_growing {
+public:
+  auto index(std::tuple<double, double> xy) const {
+    const auto x = std::get<0>(xy);
+    const auto y = std::get<1>(xy);
+    const auto r = std::sqrt(x * x + y * y);
+    return std::min(static_cast<axis::index_type>(r), size());
+  }
+
+  auto update(std::tuple<double, double> xy) {
+    const auto x = std::get<0>(xy);
+    const auto y = std::get<1>(xy);
+    const auto r = std::sqrt(x * x + y * y);
+    const auto n = static_cast<int>(r);
+    const auto old = size_;
+    if (n >= size_) size_ = n + 1;
+    return std::make_pair(n, old - size_);
+  }
+
+  axis::index_type size() const { return size_; }
+
+private:
+  axis::index_type size_ = 0;
+};
+
+template <class Tag>
+void run_tests() {
+  // one 2d axis
+  {
+    auto h = make(Tag(), axis2d());
+    h(1, 2);                  // ok, forwards 2d tuple to axis
+    h(std::make_tuple(1, 2)); // also ok, forwards 2d tuple to axis
+    BOOST_TEST_THROWS(h(1), std::invalid_argument);
+    BOOST_TEST_THROWS(h(1, 2, 3), std::invalid_argument);
+    BOOST_TEST_EQ(h.at(0), 0); // ok, bin access is still 1d
+    BOOST_TEST_EQ(h[std::make_tuple(1)], 2);
+    // also works with weights
+    h(1, 2, weight(2));
+    h(std::make_tuple(weight(3), 1, 2));
+    BOOST_TEST_EQ(h.at(0), 0);
+    BOOST_TEST_EQ(h.at(1), 7);
+
+    auto h2 = make_s(Tag(), profile_storage(), axis2d());
+    h2(1, 2, sample(2));
+    BOOST_TEST_EQ(h2[1].count(), 1);
+    BOOST_TEST_EQ(h2[1].value(), 2);
+
+    auto h3 = make_s(Tag(), weighted_profile_storage(), axis2d());
+    h3(1, 2, weight(3), sample(2));
+    BOOST_TEST_EQ(h3[1].sum_of_weights(), 3);
+    BOOST_TEST_EQ(h3[1].value(), 2);
+  }
+
+  // several axes, one 2d
+  {
+    auto h = make(Tag(), modified_axis(0, 3), minimal(), axis2d());
+    h("0", 1, std::make_tuple(1.0, 2.0));
+    h("1", 2, std::make_tuple(2.0, 1.0));
+
+    BOOST_TEST_EQ(h.rank(), 3);
+    BOOST_TEST_EQ(h.at(0, 0, 0), 0);
+    BOOST_TEST_EQ(h.at(0, 1, 1), 1);
+    BOOST_TEST_EQ(h.at(1, 0, 0), 1);
+  }
+
+  // growing axis
+  {
+    auto h = make_s(Tag{}, std::vector<int>{}, axis2d_growing{});
+    BOOST_TEST_EQ(h.size(), 0);
+    h(0, 0);
+    BOOST_TEST_EQ(h.size(), 1);
+    h(1, 0);
+    h(0, 1);
+    BOOST_TEST_EQ(h.size(), 2);
+    h(10, 0);
+    BOOST_TEST_EQ(h.size(), 11);
+    BOOST_TEST_EQ(h[0], 1);
+    BOOST_TEST_EQ(h[1], 2);
+    BOOST_TEST_EQ(h[10], 1);
+    BOOST_TEST_THROWS(h(0), std::invalid_argument);
+  }
+}
+
+int main() {
+  run_tests<static_tag>();
+  run_tests<dynamic_tag>();
+  return boost::report_errors();
+}
diff --git a/test/histogram_fill_test.cpp b/test/histogram_fill_test.cpp
new file mode 100644
index 00000000..7ad2025c
--- /dev/null
+++ b/test/histogram_fill_test.cpp
@@ -0,0 +1,299 @@
+// Copyright 2019 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)
+
+#include <boost/config.hpp>
+#include <boost/core/lightweight_test.hpp>
+#include <boost/histogram/accumulators.hpp>
+#include <boost/histogram/accumulators/ostream.hpp>
+#include <boost/histogram/algorithm/sum.hpp>
+#include <boost/histogram/axis.hpp>
+#include <boost/histogram/axis/ostream.hpp>
+#include <boost/histogram/histogram.hpp>
+#include <boost/histogram/literals.hpp>
+#include <boost/histogram/make_histogram.hpp>
+#include <boost/histogram/ostream.hpp>
+#include <boost/histogram/storage_adaptor.hpp>
+#include <boost/variant2/variant.hpp>
+#include <sstream>
+#include <stdexcept>
+#include <tuple>
+#include <utility>
+#include <vector>
+#include "throw_exception.hpp"
+#include "utility_histogram.hpp"
+
+using namespace boost::histogram;
+using namespace boost::histogram::algorithm;
+using namespace boost::histogram::literals; // to get _c suffix
+using boost::variant2::variant;
+
+using in = axis::integer<int, axis::null_type>;
+using in0 = axis::integer<int, axis::null_type, axis::option::none_t>;
+using ing = axis::integer<int, axis::null_type,
+                          decltype(axis::option::growth | axis::option::underflow |
+                                   axis::option::overflow)>;
+
+struct axis2d {
+  auto size() const { return axis::index_type{2}; }
+
+  auto index(const std::tuple<double, double>& xy) const {
+    const auto x = std::get<0>(xy);
+    const auto y = std::get<1>(xy);
+    const auto r = std::sqrt(x * x + y * y);
+    return std::min(static_cast<axis::index_type>(r), size());
+  }
+
+  friend std::ostream& operator<<(std::ostream& os, const axis2d&) {
+    os << "axis2d()";
+    return os;
+  }
+};
+
+template <class Tag>
+void run_tests() {
+  constexpr auto ndata = 1 << 20;
+
+  // 1D simple
+  {
+    auto h = make(Tag(), in{1, 3});
+    auto h3 = h;
+    auto h2 = h;
+
+    const int data[1][4] = {{0, 1, 2, 3}};
+    for (auto&& x : data[0]) h(x);
+    h2.fill(data[0]); // uses 1D specialization
+    h3.fill(data);    // uses generic form
+
+    BOOST_TEST_EQ(sum(h), 4);
+    BOOST_TEST_EQ(h, h2);
+    BOOST_TEST_EQ(h, h3);
+
+    BOOST_TEST_EQ(h2[-1], 1);
+    BOOST_TEST_EQ(h2[0], 1);
+    BOOST_TEST_EQ(h2[1], 1);
+    BOOST_TEST_EQ(h2[2], 1);
+
+    BOOST_TEST_EQ(h3[-1], 1);
+    BOOST_TEST_EQ(h3[0], 1);
+    BOOST_TEST_EQ(h3[1], 1);
+    BOOST_TEST_EQ(h3[2], 1);
+
+#ifndef BOOST_NO_EXCEPTIONS
+    int bad1[2][4];
+    std::vector<std::array<int, 4>> bad2;
+    BOOST_TEST_THROWS(h.fill(bad1), std::invalid_argument);
+    BOOST_TEST_THROWS(h.fill(bad2), std::invalid_argument);
+#endif
+  }
+
+  // 2D simple
+  {
+    auto h = make(Tag(), in{1, 3}, in0{1, 5});
+    auto h2 = h;
+
+    const std::array<int, 4> x = {0, 1, 2, 3};
+    const std::array<int, 4> y = {1, 2, 3, 4};
+    for (int i = 0; i < 4; ++i) h(x[i], y[i]);
+
+    const auto xy = {x, y};
+    h2.fill(xy);
+
+    BOOST_TEST_EQ(h, h2);
+
+    BOOST_TEST_EQ(h2.at(-1, 0), 1);
+    BOOST_TEST_EQ(h2.at(0, 1), 1);
+    BOOST_TEST_EQ(h2.at(1, 2), 1);
+    BOOST_TEST_EQ(h2.at(2, 3), 1);
+
+    int xy2[2][4] = {{0, 1, 2, 3}, {1, 2, 3, 4}};
+    h2.fill(xy2);
+
+    BOOST_TEST_EQ(h2.at(-1, 0), 2);
+    BOOST_TEST_EQ(h2.at(0, 1), 2);
+    BOOST_TEST_EQ(h2.at(1, 2), 2);
+    BOOST_TEST_EQ(h2.at(2, 3), 2);
+
+    BOOST_TEST_THROWS(h2.fill(std::array<std::vector<int>, 2>(
+                          {std::vector<int>(2), std::vector<int>(3)})),
+                      std::invalid_argument);
+  }
+
+  // 1D variant and large input collection and weight
+  {
+    auto h = make(Tag(), in{1, 3});
+    auto h2 = h;
+
+    constexpr auto n = 1 << 20;
+    const auto x1 = 1.0;
+    std::vector<double> x2(n);
+    for (unsigned i = 0; i < n; ++i) x2[i] = i % 4;
+
+    h(x1);
+    for (auto&& xi : x2) h(xi);
+
+    using V = variant<double, std::vector<double>>;
+    std::vector<V> v(1);
+    v[0] = x1;
+    h2.fill(v);
+    v[0] = x2;
+    h2.fill(v);
+
+    BOOST_TEST_EQ(h, h2);
+
+    std::vector<double> w(n);
+    for (unsigned i = 0; i < n; ++i) {
+      w[i] = i + 1;
+      h(weight(w[i]), x2[i]);
+    }
+    h2.fill(weight(w), v);
+
+    for (unsigned i = 0; i < n; ++i) { h(weight(2), x2[i]); }
+    h2.fill(weight(2), v);
+
+    BOOST_TEST_EQ(h, h2);
+
+    w.resize(n - 1);
+    BOOST_TEST_THROWS(h2.fill(v, weight(w)), std::invalid_argument);
+  }
+
+  // 2D variant and large input collection and weight
+  {
+    auto h = make(Tag(), in{1, 3}, in0{1, 5});
+
+    std::vector<double> v(ndata);
+    for (unsigned i = 0; i < ndata; ++i) v[i] = 1 + i % 4;
+
+    using V = variant<double, std::vector<double>>;
+    V xy[2];
+
+    {
+      xy[0] = 3;
+      xy[1] = v;
+      auto h1 = h;
+      auto h2 = h;
+      for (auto&& vi : v) h1(3, vi);
+      h2.fill(xy);
+      BOOST_TEST_EQ(h1, h2);
+    }
+
+    {
+      xy[0] = v;
+      xy[1] = 3;
+      auto h1 = h;
+      auto h2 = h;
+      for (auto&& vi : v) h1(vi, 3);
+      h2.fill(xy);
+      BOOST_TEST_EQ(h1, h2);
+    }
+
+    {
+      xy[0] = 3;
+      xy[1] = v;
+      auto h1 = h;
+      auto h2 = h;
+      for (unsigned i = 0; i < ndata; ++i) { h1(3, v[i], weight(v[i])); }
+      h2.fill(xy, weight(v));
+      BOOST_TEST_EQ(h1, h2);
+    }
+  }
+
+  // 1D growing
+  {
+    auto h = make(Tag(), ing());
+    auto h2 = h;
+    std::vector<int> x;
+    x.reserve(ndata + 2);
+    for (unsigned i = 0; i < ndata; ++i) { x.push_back(i % 4); }
+    x.push_back(-10);
+    x.push_back(10);
+
+    for (auto&& xi : x) h(xi);
+    h2.fill(x);
+
+    BOOST_TEST_EQ(h.size(), 23);
+    BOOST_TEST_EQ(h2.size(), 23);
+    BOOST_TEST_EQ(sum(h), sum(h2));
+    BOOST_TEST_EQ(h, h2);
+  }
+
+  // 2D growing with weights
+  {
+    auto h = make(Tag(), in(1, 3), ing());
+    auto h2 = h;
+
+    std::vector<int> xy[2];
+    for (auto&& v : xy) {
+      v.reserve(ndata + 2);
+      v.push_back(-10);
+      for (unsigned i = 0; i < ndata; ++i) v.push_back(i % 4);
+      v.push_back(10);
+    }
+
+    for (unsigned i = 0; i < ndata + 2; ++i) h(xy[0][i], xy[1][i]);
+    h2.fill(xy);
+
+    BOOST_TEST_EQ(h.size(), 4 * 23);
+    BOOST_TEST_EQ(h2.size(), 4 * 23);
+    BOOST_TEST_EQ(sum(h), sum(h2));
+    BOOST_TEST_EQ(h, h2);
+  }
+
+  // 1D with samples
+  {
+    auto h = make_s(Tag(), profile_storage(), in(1, 3));
+    auto h2 = h;
+    std::vector<double> x;
+    x.reserve(ndata);
+    for (unsigned i = 0; i < ndata; ++i) { x.push_back(static_cast<double>(i % 4)); }
+
+    for (auto&& xi : x) h(xi, sample(xi));
+    h2.fill(x, sample(x));
+
+    BOOST_TEST_EQ(h, h2);
+  }
+
+  // 2D with samples and weights
+  {
+    auto h = make_s(Tag(), weighted_profile_storage(), in(1, 3), in0(1, 3));
+    auto h2 = h;
+    std::vector<double> x;
+    x.reserve(ndata);
+    for (unsigned i = 0; i < ndata; ++i) x.push_back(static_cast<double>(i % 4));
+
+    using V = variant<double, std::vector<double>>;
+    std::array<V, 2> v;
+    v[0] = x;
+    v[1] = 3;
+
+    for (auto&& xi : x) h(xi, 3, sample(xi), weight(xi));
+    h2.fill(v, sample(x), weight(x));
+
+    BOOST_TEST_EQ(h, h2);
+  }
+
+  // axis2d
+  {
+    auto h = make(Tag(), axis2d{});
+    auto h2 = h;
+
+    std::vector<std::tuple<double, double>> xy;
+    xy.reserve(ndata);
+    for (unsigned i = 0; i < ndata; ++i)
+      xy.emplace_back(static_cast<double>(i % 4), static_cast<double>(i % 4));
+
+    for (auto&& xyi : xy) h(xyi);
+    h2.fill(xy);
+
+    BOOST_TEST_EQ(h, h2);
+  }
+}
+
+int main() {
+  run_tests<static_tag>();
+  run_tests<dynamic_tag>();
+
+  return boost::report_errors();
+}
diff --git a/test/histogram_growing_test.cpp b/test/histogram_growing_test.cpp
index 52857991..2ab9a9c9 100644
--- a/test/histogram_growing_test.cpp
+++ b/test/histogram_growing_test.cpp
@@ -12,9 +12,9 @@
 #include <boost/histogram/ostream.hpp>
 #include <string>
 #include <utility>
+#include "std_ostream.hpp"
 #include "throw_exception.hpp"
 #include "utility_histogram.hpp"
-#include "std_ostream.hpp"
 
 using namespace boost::histogram;
 
diff --git a/test/histogram_operators_test.cpp b/test/histogram_operators_test.cpp
index 544c3634..3499ce3f 100644
--- a/test/histogram_operators_test.cpp
+++ b/test/histogram_operators_test.cpp
@@ -168,19 +168,17 @@ void run_tests() {
   // histogram_ostream
   {
     auto a = make(Tag(), axis::regular<>(3, -1, 1, "r"));
-    std::ostringstream os1;
-    os1 << a;
-    BOOST_TEST_EQ(os1.str(), std::string("histogram(regular(3, -1, 1, metadata=\"r\", "
-                                         "options=underflow | overflow))"));
-
-    auto b = make(Tag(), axis::regular<>(3, -1, 1, "r"), axis::integer<>(0, 2, "i"));
-    std::ostringstream os2;
-    os2 << b;
+    std::ostringstream os;
+    os << a;
     BOOST_TEST_EQ(
-        os2.str(),
+        os.str(),
         std::string("histogram(\n"
                     "  regular(3, -1, 1, metadata=\"r\", options=underflow | overflow),\n"
-                    "  integer(0, 2, metadata=\"i\", options=underflow | overflow)\n"
+                    "  0: 0\n"
+                    "  1: 0\n"
+                    "  2: 0\n"
+                    "  3: 0\n"
+                    "  4: 0\n"
                     ")"));
   }
 }
@@ -190,6 +188,7 @@ int main() {
   run_tests<dynamic_tag>();
 
   {
+    // cannot make empty static histogram
     auto h = histogram<std::vector<axis::regular<>>>();
     std::ostringstream os;
     os << h;
diff --git a/test/histogram_test.cpp b/test/histogram_test.cpp
index bd39948e..54585002 100644
--- a/test/histogram_test.cpp
+++ b/test/histogram_test.cpp
@@ -20,16 +20,16 @@
 #include <utility>
 #include <vector>
 #include "is_close.hpp"
+#include "std_ostream.hpp"
 #include "throw_exception.hpp"
 #include "utility_allocator.hpp"
 #include "utility_axis.hpp"
 #include "utility_histogram.hpp"
-#include "std_ostream.hpp"
 
 using namespace boost::histogram;
 using namespace boost::histogram::literals; // to get _c suffix
 
-template <typename A, typename S>
+template <class A, class S>
 void pass_histogram(boost::histogram::histogram<A, S>& h) {
   BOOST_TEST_EQ(h.at(0), 0);
   BOOST_TEST_EQ(h.at(1), 1);
@@ -37,7 +37,7 @@ void pass_histogram(boost::histogram::histogram<A, S>& h) {
   BOOST_TEST_EQ(h.axis(0_c), axis::integer<>(0, 3));
 }
 
-template <typename Tag>
+template <class Tag>
 void run_tests() {
   // init_1
   {
@@ -51,41 +51,6 @@ void run_tests() {
   }
 
   // init_2
-  {
-    auto h = make(Tag(), axis::regular<>{3, -1, 1}, axis::integer<>{-1, 3});
-    BOOST_TEST_EQ(h.rank(), 2);
-    BOOST_TEST_EQ(h.size(), 30);
-    BOOST_TEST_EQ(h.axis(0_c).size(), 3);
-    BOOST_TEST_EQ(h.axis(1_c).size(), 4);
-    auto h2 = make_s(Tag(), std::vector<unsigned>(), axis::regular<>{3, -1, 1},
-                     axis::integer<>{-1, 3});
-    BOOST_TEST_EQ(h2, h);
-  }
-
-  // init_3
-  {
-    auto h = make(Tag(), axis::regular<>{3, -1, 1}, axis::integer<>{-1, 2},
-                  axis::circular<>{2, 0, 360});
-    BOOST_TEST_EQ(h.rank(), 3);
-    BOOST_TEST_EQ(h.size(), 5 * 5 * 3);
-    auto h2 = make_s(Tag(), std::vector<unsigned>(), axis::regular<>{3, -1, 1},
-                     axis::integer<>{-1, 2}, axis::circular<>{2, 0, 360});
-    BOOST_TEST_EQ(h2, h);
-  }
-
-  // init_4
-  {
-    auto h = make(Tag(), axis::regular<>{3, -1, 1}, axis::integer<>{-1, 2},
-                  axis::circular<>{2, 0, 360}, axis::variable<>{-1, 0, 1});
-    BOOST_TEST_EQ(h.rank(), 4);
-    BOOST_TEST_EQ(h.size(), 5 * 5 * 3 * 4);
-    auto h2 = make_s(Tag(), std::vector<unsigned>(), axis::regular<>{3, -1, 1},
-                     axis::integer<>{-1, 2}, axis::circular<>{2, 0, 360},
-                     axis::variable<>{-1, 0, 1});
-    BOOST_TEST_EQ(h2, h);
-  }
-
-  // init_5
   {
     auto h = make(Tag(), axis::regular<>{3, -1, 1}, axis::integer<>{-1, 2},
                   axis::circular<>{2, 0, 360}, axis::variable<>{-1, 0, 1},
@@ -141,7 +106,7 @@ void run_tests() {
 
   // axis methods
   {
-    auto a = make(Tag(), axis::regular<>(1, 1, 2, "foo"));
+    auto a = make(Tag(), axis::integer<double>(1, 2, "foo"));
     BOOST_TEST_EQ(a.axis().size(), 1);
     BOOST_TEST_EQ(a.axis().bin(0).lower(), 1);
     BOOST_TEST_EQ(a.axis().bin(0).upper(), 2);
@@ -149,7 +114,7 @@ void run_tests() {
     unsafe_access::axis(a, 0).metadata() = "bar";
     BOOST_TEST_EQ(a.axis().metadata(), "bar");
 
-    auto b = make(Tag(), axis::regular<>(1, 1, 2, "foo"), axis::integer<>(1, 3));
+    auto b = make(Tag(), axis::integer<double>(1, 2, "foo"), axis::integer<>(1, 3));
 
     // check static access
     BOOST_TEST_EQ(b.axis(0_c).size(), 1);
@@ -174,14 +139,13 @@ void run_tests() {
     unsafe_access::axis(b, 0).metadata() = "baz";
     BOOST_TEST_EQ(b.axis(0).metadata(), "baz");
 
-    enum class C { A = 3, B = 5 };
-    auto c = make(Tag(), axis::category<C>({C::A, C::B}));
+    auto c = make(Tag(), axis::category<>({1, 2}));
     BOOST_TEST_EQ(c.axis().size(), 2);
     unsafe_access::axis(c, 0).metadata() = "foo";
     BOOST_TEST_EQ(c.axis().metadata(), "foo");
     // need to cast here for this to work with Tag == dynamic_tag, too
-    auto ca = axis::get<axis::category<C>>(c.axis());
-    BOOST_TEST(ca.bin(0) == C::A);
+    auto ca = axis::get<axis::category<>>(c.axis());
+    BOOST_TEST(ca.bin(0) == 1);
   }
 
   // equal_compare
@@ -209,10 +173,9 @@ void run_tests() {
     BOOST_TEST(c != a);
   }
 
-  // d1
+  // 1D
   {
-    auto h = make_s(Tag(), std::vector<unsigned>(),
-                    axis::integer<double, axis::null_type>{0, 2});
+    auto h = make(Tag(), axis::integer<int, axis::null_type>{0, 2});
     h(0);
     auto i = h(0);
     BOOST_TEST(i == h.begin() + 1); // +1 because of underflow
@@ -231,7 +194,7 @@ void run_tests() {
     BOOST_TEST_EQ(h.at(2), 1);
   }
 
-  // d1_2
+  // 1D no *flow
   {
     auto h = make(Tag(), axis::integer<int, axis::null_type, axis::option::none_t>(0, 2));
     h(0);
@@ -250,13 +213,13 @@ void run_tests() {
     BOOST_TEST_EQ(h.at(1), 0);
   }
 
-  // d1_3
+  // 1D category axis
   {
-    auto h = make(Tag(), axis::category<std::string>({"A", "B"}));
-    h("A");
-    h("B");
-    h("D");
-    h("E");
+    auto h = make(Tag(), axis::category<>({1, 2}));
+    h(1);
+    h(2);
+    h(4);
+    h(5);
 
     BOOST_TEST_EQ(h.rank(), 1);
     BOOST_TEST_EQ(h.axis().size(), 2);
@@ -267,10 +230,9 @@ void run_tests() {
     BOOST_TEST_EQ(h.at(2), 2); // overflow bin
   }
 
-  // d1 weight
+  // 1D weight
   {
-    auto h =
-        make_s(Tag(), std::vector<accumulators::weighted_sum<>>(), axis::integer<>(0, 2));
+    auto h = make_s(Tag(), weight_storage(), axis::integer<>(0, 2));
     h(-1);
     h(0);
     h(weight(0.5), 0);
@@ -290,10 +252,9 @@ void run_tests() {
     BOOST_TEST_EQ(h[2].variance(), 4);
   }
 
-  // d1 mean
+  // 1D profile
   {
-    auto h =
-        make_s(Tag(), std::vector<accumulators::mean<double>>(), axis::integer<>(0, 2));
+    auto h = make_s(Tag(), profile_storage(), axis::integer<>(0, 2));
 
     h(0, sample(1));
     h(0, sample(2));
@@ -310,10 +271,9 @@ void run_tests() {
     BOOST_TEST_EQ(h[1].variance(), 1);
   }
 
-  // d1 weighted mean
+  // 1D weighted profile
   {
-    auto h = make_s(Tag(), std::vector<accumulators::weighted_mean<double>>(),
-                    axis::integer<>(0, 2));
+    auto h = make_s(Tag(), weighted_profile_storage(), axis::integer<>(0, 2));
 
     h(0, sample(1));
     h(sample(1), 0);
@@ -333,9 +293,9 @@ void run_tests() {
     BOOST_TEST_EQ(h[1].value(), 2.5);
   }
 
-  // d2
+  // 2D
   {
-    auto h = make(Tag(), axis::regular<>(2, -1, 1),
+    auto h = make(Tag(), axis::integer<>(-1, 1),
                   axis::integer<int, axis::null_type, axis::option::none_t>(-1, 2));
     h(-1, -1);
     h(-1, 0);
@@ -364,10 +324,9 @@ void run_tests() {
     BOOST_TEST_EQ(h.at(2, 2), 0);
   }
 
-  // d2w
+  // 2D weight
   {
-    auto h = make_s(Tag(), std::vector<accumulators::weighted_sum<>>(),
-                    axis::regular<>(2, -1, 1),
+    auto h = make_s(Tag(), weight_storage(), axis::integer<>(-1, 1),
                     axis::integer<int, axis::null_type, axis::option::none_t>(-1, 2));
     h(-1, 0);              // -> 0, 1
     h(weight(10), -1, -1); // -> 0, 0
@@ -410,15 +369,13 @@ void run_tests() {
     BOOST_TEST_EQ(h.at(2, 2).variance(), 0);
   }
 
-  // d3w
+  // 3D weight
   {
-    auto h = make_s(Tag(), std::vector<accumulators::weighted_sum<>>(),
-                    axis::integer<>(0, 3), axis::integer<>(0, 4), axis::integer<>(0, 5));
-    for (auto i = 0; i < h.axis(0_c).size(); ++i) {
-      for (auto j = 0; j < h.axis(1_c).size(); ++j) {
-        for (auto k = 0; k < h.axis(2_c).size(); ++k) { h(i, j, k, weight(i + j + k)); }
-      }
-    }
+    auto h = make_s(Tag(), weight_storage(), axis::integer<>(0, 3), axis::integer<>(0, 4),
+                    axis::integer<>(0, 5));
+    for (auto i = 0; i < h.axis(0_c).size(); ++i)
+      for (auto j = 0; j < h.axis(1_c).size(); ++j)
+        for (auto k = 0; k < h.axis(2_c).size(); ++k) h(i, j, k, weight(i + j + k));
 
     for (auto i = 0; i < h.axis(0_c).size(); ++i) {
       for (auto j = 0; j < h.axis(1_c).size(); ++j) {
@@ -464,40 +421,10 @@ void run_tests() {
     BOOST_TEST_EQ(algorithm::sum(h), 0);
   }
 
-  // custom axes
-  {
-    struct modified_axis : public axis::integer<> {
-      using integer::integer; // inherit ctors of base
-      // customization point: convert argument and call base class
-      auto index(const char* s) const { return integer::index(std::atoi(s)); }
-    };
-
-    struct minimal {
-      auto index(int x) const { return static_cast<axis::index_type>(x % 2); }
-      auto size() const { return axis::index_type{2}; }
-    };
-
-    struct axis2d {
-      auto index(const std::tuple<double, double>& x) const {
-        return axis::index_type{std::get<0>(x) == 1 && std::get<1>(x) == 2};
-      }
-      auto size() const { return axis::index_type{2}; }
-    };
-
-    auto h = make(Tag(), modified_axis(0, 3), minimal(), axis2d());
-    h("0", 1, std::make_tuple(1.0, 2.0));
-    h("1", 2, std::make_tuple(2.0, 1.0));
-
-    BOOST_TEST_EQ(h.rank(), 3);
-    BOOST_TEST_EQ(h.at(0, 0, 0), 0);
-    BOOST_TEST_EQ(h.at(0, 1, 1), 1);
-    BOOST_TEST_EQ(h.at(1, 0, 0), 1);
-  }
-
   // using containers or input and output
   {
-    auto h = make_s(Tag(), std::vector<accumulators::weighted_sum<>>(),
-                    axis::integer<>(0, 2), axis::regular<>(2, 2, 4));
+    auto h = make_s(Tag(), weight_storage(), axis::integer<>(0, 2),
+                    axis::integer<double>(2, 4));
     // tuple in
     h(std::make_tuple(0, 2.0));
     h(std::make_tuple(1, 3.0));
@@ -535,23 +462,6 @@ void run_tests() {
     auto h1 = make(Tag(), axis::integer<>(0, 2));
     h1(std::make_tuple(0));                      // as if one had passed 0 directly
     BOOST_TEST_EQ(h1.at(std::make_tuple(0)), 1); // as if one had passed 0 directly
-
-    struct axis2d {
-      auto index(std::tuple<int, int> x) const {
-        return axis::index_type{std::get<0>(x) == 1 && std::get<1>(x) == 2};
-      }
-      auto size() const { return axis::index_type{2}; }
-    };
-    auto h2 = make(Tag(), axis2d());
-    h2(std::make_tuple(1, 2));  // ok, forwards 2d tuple to axis
-    BOOST_TEST_EQ(h2.at(0), 0); // ok, bin access is still 1d
-    BOOST_TEST_EQ(h2[std::make_tuple(1)], 1);
-    // passing two arguments directly also works
-    h2(1, 2);
-    // also works with weights
-    h2(1, 2, weight(2));
-    h2(std::make_tuple(weight(3), 1, 2));
-    BOOST_TEST_EQ(h2.at(1), 7);
   }
 
   // bad bin access
diff --git a/test/unlimited_storage_test.cpp b/test/unlimited_storage_test.cpp
index 48704fea..1da2601e 100644
--- a/test/unlimited_storage_test.cpp
+++ b/test/unlimited_storage_test.cpp
@@ -7,6 +7,7 @@
 #include <algorithm>
 #include <boost/core/lightweight_test.hpp>
 #include <boost/core/lightweight_test_trait.hpp>
+#include <boost/histogram/detail/detect.hpp>
 #include <boost/histogram/storage_adaptor.hpp>
 #include <boost/histogram/unlimited_storage.hpp>
 #include <boost/histogram/unsafe_access.hpp>
@@ -406,6 +407,8 @@ int main() {
     BOOST_TEST_EQ(unsafe_access::unlimited_storage_buffer(d).type, 5);
     BOOST_TEST_EQ(d[0], -2);
     BOOST_TEST_EQ(d[1], 2);
+
+    BOOST_TEST_TRAIT_TRUE((detail::has_operator_preincrement<decltype(d[0])>));
   }
 
   // iterators
diff --git a/tools/check_build_system.py b/tools/check_build_system.py
index f18d8354..dfbf18bb 100755
--- a/tools/check_build_system.py
+++ b/tools/check_build_system.py
@@ -21,6 +21,7 @@ for dir in ("test", "examples"):
         run = set(re.findall("([a-zA-Z0-9_]+\.cpp)", open(filename).read()))
 
         diff = cpp - run
+        diff.discard("check_cmake_version.cpp") # ignore
 
         if diff:
             print("NOT TESTED in %s\n  " % filename +