diff --git a/CMakeLists.txt b/CMakeLists.txt
index c803d2ec..c70f736f 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -81,6 +81,7 @@ endfunction()
 
 compiled_test(test/adaptive_storage_test.cpp)
 compiled_test(test/algorithm_project_test.cpp)
+compiled_test(test/algorithm_reduce_test.cpp)
 compiled_test(test/algorithm_sum_test.cpp)
 compiled_test(test/axis_regular_test.cpp)
 compiled_test(test/axis_circular_test.cpp)
@@ -92,7 +93,6 @@ compiled_test(test/detail_test.cpp)
 compiled_test(test/histogram_dynamic_test.cpp)
 compiled_test(test/histogram_mixed_test.cpp)
 compiled_test(test/histogram_test.cpp)
-compiled_test(test/index_mapper_test.cpp)
 compiled_test(test/internal_accumulators_test.cpp)
 compiled_test(test/meta_test.cpp)
 compiled_test(test/storage_adaptor_test.cpp)
diff --git a/include/boost/histogram/algorithm/project.hpp b/include/boost/histogram/algorithm/project.hpp
index 47291082..b2d6da8f 100644
--- a/include/boost/histogram/algorithm/project.hpp
+++ b/include/boost/histogram/algorithm/project.hpp
@@ -8,22 +8,17 @@
 #define BOOST_HISTOGRAM_ALGORITHM_PROJECT_HPP
 
 #include <algorithm>
-#include <boost/assert.hpp>
 #include <boost/histogram/detail/axes.hpp>
 #include <boost/histogram/detail/index_mapper.hpp>
-#include <boost/histogram/detail/is_set.hpp>
 #include <boost/histogram/histogram_fwd.hpp>
 #include <boost/histogram/unsafe_access.hpp>
 #include <boost/mp11.hpp>
 #include <stdexcept>
-#include <tuple>
 
 namespace boost {
 namespace histogram {
 namespace algorithm {
 
-// TODO: make generic reduce, which can sum over axes, shrink, rebin
-
 /**
   Returns a lower-dimensional histogram, summing over removed axes.
 
@@ -49,59 +44,70 @@ auto project(const histogram<A, S>& h, mp11::mp_size_t<I> n, Ns... ns) {
   std::size_t s = 1;
   h.for_each_axis([&](const auto& a) {
     const auto n = axis::traits::extend(a);
-    im.ntotal *= n;
-    iter->first = s;
+    im.total *= n;
+    iter->stride[0] = s;
     s *= n;
-    iter->second = 0;
+    iter->stride[1] = 0;
     ++iter;
   });
 
   s = 1;
   mp11::mp_for_each<LN>([&](auto J) {
-    im[J].second = s;
+    im[J].stride[1] = s;
     s *= axis::traits::extend(detail::axis_get<J>(axes));
   });
 
-  do {
-    const auto x = unsafe_access::storage(h)[im.first];
-    unsafe_access::storage(r_h).add(im.second, x);
-  } while (im.next());
+  im(unsafe_access::storage(r_h), unsafe_access::storage(h));
   return r_h;
 }
 
 /**
   Returns a lower-dimensional histogram, summing over removed axes.
 
-  This version accepts an iterator range that represents the indices which are kept.
+  This version accepts an iterable range that represents the indices which are kept.
 */
-template <typename A, typename S, typename Iterator,
-          typename = detail::requires_axis_vector<A>,
-          typename = detail::requires_iterator<Iterator>>
-auto project(const histogram<A, S>& h, Iterator begin, Iterator end) {
-  BOOST_ASSERT_MSG(detail::is_set(begin, end), "indices must be unique");
+template <typename A, typename S, typename C, typename = detail::requires_axis_vector<A>,
+          typename = detail::requires_iterable<C>>
+auto project(const histogram<A, S>& h, C c) {
   using H = histogram<A, S>;
 
+  auto begin = std::begin(c);
+  auto end = std::end(c);
+
   const auto& axes = unsafe_access::axes(h);
-  auto r_axes = typename H::axes_type(axes.get_allocator());
+  auto r_axes = detail::static_if<detail::has_allocator<A>>(
+      [](const auto& axes) {
+        using T = detail::unqual<decltype(axes)>;
+        return T(axes.get_allocator());
+      },
+      [](const auto& axes) {
+        using T = detail::unqual<decltype(axes)>;
+        return T();
+      },
+      axes);
   r_axes.reserve(std::distance(begin, end));
 
   detail::index_mapper im(h.rank());
   auto iter = im.begin();
-  std::size_t s = 1;
+  std::size_t stride = 1;
   h.for_each_axis([&](const auto& a) {
     const auto n = axis::traits::extend(a);
-    im.ntotal *= n;
-    iter->first = s;
-    s *= n;
-    iter->second = 0;
+    im.total *= n;
+    iter->stride[0] = stride;
+    stride *= n;
+    iter->stride[1] = 0;
     ++iter;
   });
 
-  s = 1;
+  stride = 1;
   for (auto it = begin; it != end; ++it) {
     r_axes.emplace_back(axes[*it]);
-    im[*it].second = s;
-    s *= axis::traits::extend(axes[*it]);
+    auto& stride_ref = im[*it].stride[1];
+    if (stride_ref)
+      throw std::invalid_argument("indices must be unique");
+    else
+      stride_ref = stride;
+    stride *= axis::traits::extend(axes[*it]);
   }
 
   auto r_h = H(std::move(r_axes),
@@ -109,10 +115,7 @@ auto project(const histogram<A, S>& h, Iterator begin, Iterator end) {
                    [&h](auto) { return S(unsafe_access::storage(h).get_allocator()); },
                    [](auto) { return S(); }, 0));
 
-  do {
-    const auto x = unsafe_access::storage(h)[im.first];
-    unsafe_access::storage(r_h).add(im.second, x);
-  } while (im.next());
+  im(unsafe_access::storage(r_h), unsafe_access::storage(h));
   return r_h;
 }
 
diff --git a/include/boost/histogram/algorithm/reduce.hpp b/include/boost/histogram/algorithm/reduce.hpp
new file mode 100644
index 00000000..eff3a5b0
--- /dev/null
+++ b/include/boost/histogram/algorithm/reduce.hpp
@@ -0,0 +1,147 @@
+// 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_ALGORITHM_REDUCE_HPP
+#define BOOST_HISTOGRAM_ALGORITHM_REDUCE_HPP
+
+#include <boost/assert.hpp>
+#include <boost/container/static_vector.hpp>
+#include <boost/histogram/detail/axes.hpp>
+#include <boost/histogram/detail/index_mapper.hpp>
+#include <boost/histogram/detail/meta.hpp>
+#include <boost/histogram/histogram_fwd.hpp>
+#include <boost/histogram/unsafe_access.hpp>
+#include <cmath>
+#include <limits>
+#include <stdexcept>
+#include <type_traits>
+
+namespace boost {
+namespace histogram {
+namespace algorithm {
+
+struct reduce_option_type {
+  unsigned iaxis = 0;
+  bool has_range = false;
+  double lower, upper;
+  unsigned merge = 0;
+
+  reduce_option_type() noexcept = default;
+
+  reduce_option_type(unsigned i, bool b, double l, double u, unsigned m)
+      : iaxis(i), has_range(b), lower(l), upper(u), merge(m) {
+    if (has_range && l == u) throw std::invalid_argument("lower != upper required");
+    if (merge == 0) throw std::invalid_argument("merge > 0 required");
+  }
+
+  operator bool() const noexcept { return merge; }
+};
+
+reduce_option_type shrink(unsigned iaxis, double lower, double upper) {
+  return {iaxis, true, lower, upper, 1};
+}
+
+reduce_option_type shrink_and_rebin(unsigned iaxis, double lower, double upper,
+                                    unsigned merge) {
+  return {iaxis, true, lower, upper, merge};
+}
+
+reduce_option_type rebin(unsigned iaxis, unsigned merge) {
+  return {iaxis, false, 0.0, 0.0, merge};
+}
+
+template <typename A, typename S, typename C, typename = detail::requires_iterable<C>>
+histogram<A, S> reduce(const histogram<A, S>& h, const C& c) {
+  auto options =
+      boost::container::static_vector<reduce_option_type, axis::limit>(h.rank());
+  for (const auto& o : c) {
+    auto& opt_ref = options[o.iaxis];
+    if (opt_ref) throw std::invalid_argument("indices must be unique");
+    opt_ref.lower = o.lower;
+    opt_ref.upper = o.upper;
+    opt_ref.merge = o.merge;
+  }
+
+  auto r_axes = detail::make_empty_axes(unsafe_access::axes(h));
+
+  detail::index_mapper_reduce im(h.rank());
+  auto im_iter = im.begin();
+  std::size_t stride[2] = {1, 1};
+  unsigned iaxis = 0;
+  h.for_each_axis([&](const auto& a) {
+    using T = detail::unqual<decltype(a)>;
+
+    const auto n = axis::traits::extend(a);
+    im.total *= n;
+    im_iter->stride[0] = stride[0];
+    stride[0] *= n;
+    auto set_flow = [im_iter](int i, const auto& a) {
+      switch (axis::traits::options(a)) {
+        case axis::option_type::overflow:
+          im_iter->overflow[i] = a.size();
+          im_iter->underflow[i] = -1;
+          break;
+        case axis::option_type::underflow_and_overflow:
+          im_iter->overflow[i] = a.size();
+          im_iter->underflow[i] = a.size() + 1;
+          break;
+        default: im_iter->underflow[i] = -1; im_iter->overflow[i] = -1;
+      };
+    };
+    set_flow(0, a);
+
+    const auto& opt = options[iaxis];
+    unsigned begin = 0, end = a.size(), merge = 1;
+    if (opt) {
+      merge = opt.merge;
+      if (opt.lower < opt.upper) {
+        while (begin != end && a.value(begin) < opt.lower) ++begin;
+        while (end != begin && a.value(end - 1) >= opt.upper) --end;
+      } else if (opt.upper < opt.lower) {
+        while (begin != end && a.value(begin) > opt.lower) ++begin;
+        while (end != begin && a.value(end - 1) <= opt.upper) --end;
+      }
+      end -= (end - begin) % merge;
+      auto a2 = T(a, begin, end, merge);
+      axis::get<T>(detail::axis_get(r_axes, iaxis)) = a2;
+      im_iter->stride[1] = stride[1];
+      stride[1] *= axis::traits::extend(a2);
+      set_flow(1, a2);
+    } else {
+      axis::get<T>(detail::axis_get(r_axes, iaxis)) = a;
+      im_iter->stride[1] = stride[1];
+      stride[1] *= axis::traits::extend(a);
+      set_flow(1, a);
+    }
+    im_iter->begin = begin;
+    im_iter->end = end;
+    im_iter->merge = merge;
+    ++im_iter;
+    ++iaxis;
+  });
+
+  auto h_r = histogram<A, S>(
+      std::move(r_axes),
+      detail::static_if<detail::has_allocator<S>>(
+          [&h](auto) { return S(unsafe_access::storage(h).get_allocator()); },
+          [](auto) { return S(); }, 0));
+
+  im(unsafe_access::storage(h_r), unsafe_access::storage(h));
+  return h_r;
+}
+
+template <typename A, typename S, typename... Ts>
+histogram<A, S> reduce(const histogram<A, S>& h, const reduce_option_type& t,
+                       Ts&&... ts) {
+  // this must be in one line, because any of the ts could be a temporary
+  return reduce(h, std::initializer_list<reduce_option_type>{t, ts...});
+}
+
+} // namespace algorithm
+} // namespace histogram
+} // namespace boost
+
+#endif
diff --git a/include/boost/histogram/axis/base.hpp b/include/boost/histogram/axis/base.hpp
index 1e79ddae..e1854cd3 100644
--- a/include/boost/histogram/axis/base.hpp
+++ b/include/boost/histogram/axis/base.hpp
@@ -44,7 +44,7 @@ public:
   void serialize(Archive&, unsigned);
 
 protected:
-  base(unsigned n, metadata_type&& m, option_type opt)
+  base(unsigned n, metadata_type m, option_type opt)
       : size_meta_(n, std::move(m)), opt_(opt) {
     if (size() == 0) { throw std::invalid_argument("bins > 0 required"); }
     const auto max_index =
diff --git a/include/boost/histogram/axis/integer.hpp b/include/boost/histogram/axis/integer.hpp
index cbd100d1..37af2185 100644
--- a/include/boost/histogram/axis/integer.hpp
+++ b/include/boost/histogram/axis/integer.hpp
@@ -35,19 +35,26 @@ class integer : public base<MetaData>, public iterator_mixin<integer<IntType, Me
   using bin_view =
       std::conditional_t<std::is_integral<value_type>::value, value_bin_view<integer>,
                          interval_bin_view<integer>>;
+  using index_type = std::conditional_t<std::is_integral<value_type>::value, int, double>;
 
 public:
   /** Construct over semi-open integer interval [start, stop).
    *
    * \param start    first integer of covered range.
-   * \param stop      one past last integer of covered range.
+   * \param stop     one past last integer of covered range.
    * \param metadata description of the axis.
    * \param options  extra bin options.
    */
   integer(value_type start, value_type stop, metadata_type m = metadata_type(),
           option_type o = option_type::underflow_and_overflow)
-      : base_type(static_cast<unsigned>(stop - start), std::move(m), o), min_(start) {
-    if (start >= stop) { throw std::invalid_argument("start < stop required"); }
+      : base_type(static_cast<unsigned>(stop - start > 0 ? stop - start : 0),
+                  std::move(m), o)
+      , min_(start) {}
+
+  /// Constructor used by algorithm::reduce to shrink and rebin.
+  integer(const integer& src, unsigned begin, unsigned end, unsigned merge)
+      : base_type(end - begin, src.metadata(), src.options()), min_(src.min_ + begin) {
+    if (merge > 1) { throw std::invalid_argument("cannot merge bins for integer axis"); }
   }
 
   integer() = default;
@@ -60,7 +67,7 @@ public:
   }
 
   /// Returns axis value for index.
-  value_type value(value_type i) const noexcept {
+  value_type value(index_type i) const noexcept {
     if (i < 0) {
       return detail::static_if<std::is_integral<value_type>>(
           [](auto) { return std::numeric_limits<value_type>::min(); },
diff --git a/include/boost/histogram/axis/polymorphic_bin.hpp b/include/boost/histogram/axis/polymorphic_bin.hpp
index 0c409715..6022c189 100644
--- a/include/boost/histogram/axis/polymorphic_bin.hpp
+++ b/include/boost/histogram/axis/polymorphic_bin.hpp
@@ -7,6 +7,8 @@
 #ifndef BOOST_HISTOGRAM_AXIS_POLYMORPHIC_BIN_HPP
 #define BOOST_HISTOGRAM_AXIS_POLYMORPHIC_BIN_HPP
 
+#include <cmath>
+#include <stdexcept>
 #include <tuple>
 
 namespace boost {
@@ -64,7 +66,7 @@ public:
 
   explicit operator int() const noexcept { return idx_; }
 
-  bool is_discrete() const noexcept { return std::get<1>(data_) == std::get<0>(data_); }
+  bool is_discrete() const noexcept { return std::isnan(std::get<2>(data_)); }
 
 private:
   const int idx_;
diff --git a/include/boost/histogram/axis/regular.hpp b/include/boost/histogram/axis/regular.hpp
index 830b5bfa..a9d56fa6 100644
--- a/include/boost/histogram/axis/regular.hpp
+++ b/include/boost/histogram/axis/regular.hpp
@@ -83,8 +83,8 @@ class regular : public base<MetaData>,
   using transform_type = Transform;
   using external_type = detail::return_type<decltype(&Transform::inverse)>;
   using internal_type = detail::return_type<decltype(&Transform::forward)>;
-  static_assert(!std::is_integral<internal_type>::value,
-                "type returned by forward transform cannot be integral");
+  static_assert(std::is_floating_point<internal_type>::value,
+                "forward transform must return floating point number");
   using metadata_type = MetaData;
 
 public:
@@ -102,7 +102,7 @@ public:
       : base_type(n, std::move(m), o)
       , transform_type(std::move(trans))
       , min_(this->forward(start))
-      , delta_((this->forward(stop) - this->forward(start)) / n) {
+      , delta_((this->forward(stop) - min_) / base_type::size()) {
     if (!std::isfinite(min_) || !std::isfinite(delta_))
       throw std::invalid_argument("forward transform of start or stop invalid");
     if (delta_ == 0)
@@ -121,6 +121,15 @@ public:
           option_type o = option_type::underflow_and_overflow)
       : regular({}, n, start, stop, std::move(m), o) {}
 
+  /// Constructor used by algorithm::reduce to shrink and rebin (not for users).
+  regular(const regular& src, unsigned begin, unsigned end, unsigned merge)
+      : base_type((end - begin) / merge, src.metadata(), src.options())
+      , transform_type(src.transform())
+      , min_(this->forward(src.value(begin)))
+      , delta_(src.delta_ * merge) {
+    BOOST_ASSERT((end - begin) % merge == 0);
+  }
+
   regular() = default;
 
   /// Returns instance of the transform type
diff --git a/include/boost/histogram/axis/variant.hpp b/include/boost/histogram/axis/variant.hpp
index bb079960..ee907311 100644
--- a/include/boost/histogram/axis/variant.hpp
+++ b/include/boost/histogram/axis/variant.hpp
@@ -41,7 +41,7 @@ struct get_polymorphic_bin_data
           using Arg = detail::unqual<detail::arg_type<detail::unqual<decltype(a)>>>;
           const auto x = a.value(idx);
           return detail::static_if<std::is_integral<Arg>>(
-              [&](auto) { return T(x, x, 0); },
+              [&](auto) { return T(x, 0, std::numeric_limits<double>::quiet_NaN()); },
               [&](auto) { return T(x, a.value(idx + 1), a.value(idx + 0.5)); }, 0);
         },
         [](const auto&) -> T {
@@ -304,20 +304,8 @@ const T* get(const variant<Us...>* v) {
 
 // pass-through version for generic programming, if U is axis instead of variant
 template <typename T, typename U, typename = detail::requires_axis<detail::unqual<U>>>
-T& get(U& u) {
-  return static_cast<T&>(u);
-}
-
-// pass-through version for generic programming, if U is axis instead of variant
-template <typename T, typename U, typename = detail::requires_axis<detail::unqual<U>>>
-T&& get(U&& u) {
-  return static_cast<T&&>(u);
-}
-
-// pass-through version for generic programming, if U is axis instead of variant
-template <typename T, typename U, typename = detail::requires_axis<detail::unqual<U>>>
-const T& get(const U& u) {
-  return static_cast<const T&>(u);
+auto get(U&& u) -> detail::copy_qualifiers<U, T> {
+  return std::forward<U>(u);
 }
 
 // pass-through version for generic programming, if U is axis instead of variant
diff --git a/include/boost/histogram/detail/axes.hpp b/include/boost/histogram/detail/axes.hpp
index 843c9084..74c19fe8 100644
--- a/include/boost/histogram/detail/axes.hpp
+++ b/include/boost/histogram/detail/axes.hpp
@@ -126,15 +126,14 @@ void axes_assign(T& t, const U& u) {
   t.assign(u.begin(), u.end());
 }
 
-template <typename T, std::size_t N = std::tuple_size<T>::value>
-constexpr std::size_t axes_size(const T&) noexcept {
-  return N;
-}
-
-// static to fix gcc warning about mangled names changing in C++17
-template <typename T, typename = decltype(&T::size)>
-static std::size_t axes_size(const T& axes) noexcept {
-  return axes.size();
+template <typename T>
+constexpr std::size_t axes_size(const T& axes) noexcept {
+  return static_if<has_fixed_size<unqual<T>>>(
+      [](const auto& a) {
+        using U = unqual<decltype(a)>;
+        return std::tuple_size<U>::value;
+      },
+      [&](const auto& a) { return a.size(); }, axes);
 }
 
 template <typename T>
@@ -152,6 +151,22 @@ void for_each_axis(const T& axes, F&& f) {
   for (const auto& x : axes) { axis::visit(std::forward<F>(f), x); }
 }
 
+template <typename T>
+auto make_empty_axes(const T& t) {
+  auto r = T();
+  static_if<is_vector_like<T>>([&](auto) { r.reserve(t.size()); }, [](auto) {}, 0);
+  for_each_axis(t, [&r](const auto& a) {
+    using U = unqual<decltype(a)>;
+    r.emplace_back(U());
+  });
+  return r;
+}
+
+template <typename... Ts>
+auto make_empty_axes(const std::tuple<Ts...>&) {
+  return std::tuple<Ts...>();
+}
+
 template <typename T>
 std::size_t bincount(const T& axes) {
   std::size_t n = 1;
diff --git a/include/boost/histogram/detail/index_mapper.hpp b/include/boost/histogram/detail/index_mapper.hpp
index 05752cc5..055217e9 100644
--- a/include/boost/histogram/detail/index_mapper.hpp
+++ b/include/boost/histogram/detail/index_mapper.hpp
@@ -10,29 +10,75 @@
 #include <boost/container/static_vector.hpp>
 #include <boost/histogram/histogram_fwd.hpp>
 #include <cstddef>
+#include <functional>
 
 namespace boost {
 namespace histogram {
 namespace detail {
+
+struct index_mapper_item {
+  std::size_t stride[2];
+};
+
 class index_mapper
-    : public boost::container::static_vector<std::pair<std::size_t, std::size_t>,
-                                             axis::limit> {
+    : public boost::container::static_vector<index_mapper_item, axis::limit> {
 public:
-  std::size_t first = 0, second = 0, ntotal = 1;
+  std::size_t total = 1;
 
-  using static_vector<std::pair<std::size_t, std::size_t>, axis::limit>::static_vector;
+  index_mapper(unsigned dim) : static_vector(dim) {}
 
-  bool next() {
-    ++first;
-    second = 0;
-    auto f = first;
-    for (auto it = end(); it != begin(); --it) {
-      const auto& d = *(it - 1);
-      // compiler usually optimizes div & mod into one div
-      second += f / d.first * d.second;
-      f %= d.first;
+  template <typename T, typename U>
+  void operator()(T& dst, const U& src) {
+    for (std::size_t i = 0; i < total; ++i) {
+      std::size_t j = 0;
+      auto imod = i;
+      for (auto it = end(); it != begin(); --it) {
+        const auto& d = *(it - 1);
+        // compiler usually optimizes div & mod into one div
+        const auto k = imod / d.stride[0];
+        imod %= d.stride[0];
+        j += k * d.stride[1];
+      }
+      dst.add(j, src[i]);
+    }
+  }
+};
+
+struct index_mapper_reduce_item {
+  std::size_t stride[2];
+  int underflow[2], overflow[2], begin, end, merge;
+};
+
+class index_mapper_reduce
+    : public boost::container::static_vector<index_mapper_reduce_item, axis::limit> {
+public:
+  std::size_t total = 1;
+
+  index_mapper_reduce(unsigned dim) : static_vector(dim) {}
+
+  template <typename T, typename U>
+  void operator()(T& dst, const U& src) {
+    for (std::size_t i = 0; i < total; ++i) {
+      std::size_t j = 0;
+      auto imod = i;
+      bool drop = false;
+      for (auto it = end(); it != begin(); --it) {
+        const auto& d = *(it - 1);
+        // compiler usually optimizes div & mod into one div
+        int k = imod / d.stride[0];
+        imod %= d.stride[0];
+        if (k < d.begin || k == d.underflow[0]) {
+          k = d.underflow[1];
+        } else if (k >= d.end || k == d.overflow[0]) {
+          k = d.overflow[1];
+        } else {
+          k = (k - d.begin) / d.merge;
+        }
+        drop |= k < 0;
+        j += k * d.stride[1];
+      }
+      if (!drop) dst.add(j, src[i]);
     }
-    return first < ntotal;
   }
 };
 } // namespace detail
diff --git a/include/boost/histogram/detail/is_set.hpp b/include/boost/histogram/detail/is_set.hpp
index 182b8b15..57fc7dcd 100644
--- a/include/boost/histogram/detail/is_set.hpp
+++ b/include/boost/histogram/detail/is_set.hpp
@@ -7,10 +7,10 @@
 #ifndef BOOST_HISTOGRAM_DETAIL_IS_SET_HPP
 #define BOOST_HISTOGRAM_DETAIL_IS_SET_HPP
 
-#include <algorithm>
+#include <boost/container/flat_set.hpp>
 #include <boost/container/static_vector.hpp>
 #include <boost/histogram/detail/meta.hpp>
-#include <boost/histogram/histogram_fwd.hpp>
+#include <boost/histogram/histogram_fwd.hpp> // for axis::limit
 
 namespace boost {
 namespace histogram {
@@ -18,10 +18,9 @@ namespace detail {
 template <class Iterator>
 bool is_set(Iterator begin, Iterator end) {
   using T = iterator_value_type<Iterator>;
-  boost::container::static_vector<T, axis::limit> v(begin, end);
-  std::sort(v.begin(), v.end());
-  auto end2 = std::unique(v.begin(), v.end());
-  return static_cast<std::size_t>(std::distance(v.begin(), end2)) == v.size();
+  using C = boost::container::static_vector<T, axis::limit>;
+  boost::container::flat_set<T, std::less<T>, C> s(begin, end);
+  return static_cast<std::size_t>(std::distance(begin, end)) == s.size();
 }
 } // namespace detail
 } // namespace histogram
diff --git a/include/boost/histogram/detail/meta.hpp b/include/boost/histogram/detail/meta.hpp
index 0245259f..e34a5dbb 100644
--- a/include/boost/histogram/detail/meta.hpp
+++ b/include/boost/histogram/detail/meta.hpp
@@ -176,11 +176,10 @@ BOOST_HISTOGRAM_MAKE_SFINAE(is_indexable_container, (std::declval<T&>()[0], &T::
                                                      std::begin(std::declval<T&>()),
                                                      std::end(std::declval<T&>())));
 
-BOOST_HISTOGRAM_MAKE_SFINAE(is_tuple, (std::get<0>(std::declval<T&>())));
-
 BOOST_HISTOGRAM_MAKE_SFINAE(is_equal_comparable,
                             (std::declval<T&>() == std::declval<T&>()));
 
+// is_axis is false for axis::variant, because operator() is templated
 BOOST_HISTOGRAM_MAKE_SFINAE(is_axis, (&T::size, &T::operator()));
 
 BOOST_HISTOGRAM_MAKE_SFINAE(is_iterable, (std::begin(std::declval<T&>()),
@@ -191,6 +190,17 @@ BOOST_HISTOGRAM_MAKE_SFINAE(is_streamable,
 
 BOOST_HISTOGRAM_MAKE_SFINAE(is_incrementable, (++std::declval<T&>()));
 
+BOOST_HISTOGRAM_MAKE_SFINAE(has_fixed_size, (std::tuple_size<T>::value));
+
+template <typename T>
+struct is_tuple_impl : std::false_type {};
+
+template <typename... Ts>
+struct is_tuple_impl<std::tuple<Ts...>> : std::true_type {};
+
+template <typename T>
+using is_tuple = typename is_tuple_impl<T>::type;
+
 template <typename T>
 struct is_axis_variant_impl : std::false_type {};
 
diff --git a/include/boost/histogram/histogram.hpp b/include/boost/histogram/histogram.hpp
index ff0ee819..70222c8b 100644
--- a/include/boost/histogram/histogram.hpp
+++ b/include/boost/histogram/histogram.hpp
@@ -11,7 +11,6 @@
 #include <boost/assert.hpp>
 #include <boost/histogram/arithmetic_operators.hpp>
 #include <boost/histogram/detail/axes.hpp>
-#include <boost/histogram/detail/index_mapper.hpp>
 #include <boost/histogram/detail/meta.hpp>
 #include <boost/histogram/histogram_fwd.hpp>
 #include <boost/histogram/iterator.hpp>
diff --git a/test/algorithm_project_test.cpp b/test/algorithm_project_test.cpp
index 07422702..6636bed5 100644
--- a/test/algorithm_project_test.cpp
+++ b/test/algorithm_project_test.cpp
@@ -150,7 +150,7 @@ int main() {
     std::vector<int> x;
 
     x = {0};
-    auto hx = project(h, x.begin(), x.end());
+    auto hx = project(h, x);
     BOOST_TEST_EQ(hx.rank(), 1);
     BOOST_TEST_EQ(sum(hx), 6);
     BOOST_TEST_EQ(hx.at(0), 2);
@@ -158,7 +158,7 @@ int main() {
     BOOST_TEST(hx.axis() == h.axis(0_c));
 
     x = {1};
-    auto hy = project(h, x.begin(), x.end());
+    auto hy = project(h, x);
     BOOST_TEST_EQ(hy.rank(), 1);
     BOOST_TEST_EQ(sum(hy), 6);
     BOOST_TEST_EQ(hy.at(0), 2);
@@ -167,7 +167,7 @@ int main() {
     BOOST_TEST(hy.axis() == h.axis(1_c));
 
     x = {1, 0};
-    auto hyx = project(h, x.begin(), x.end());
+    auto hyx = project(h, x);
     BOOST_TEST_EQ(hyx.rank(), 2);
     BOOST_TEST_EQ(hyx.axis(0_c), h.axis(1_c));
     BOOST_TEST_EQ(hyx.axis(1_c), h.axis(0_c));
@@ -177,6 +177,9 @@ int main() {
     BOOST_TEST_EQ(hyx.at(0, 1), 1);
     BOOST_TEST_EQ(hyx.at(1, 1), 1);
     BOOST_TEST_EQ(hyx.at(2, 1), 2);
+
+    x = {0, 0};
+    BOOST_TEST_THROWS(project(h, x), std::invalid_argument);
   }
 
   return boost::report_errors();
diff --git a/test/algorithm_reduce_test.cpp b/test/algorithm_reduce_test.cpp
new file mode 100644
index 00000000..fa5575bc
--- /dev/null
+++ b/test/algorithm_reduce_test.cpp
@@ -0,0 +1,136 @@
+// 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)
+
+#include <boost/core/lightweight_test.hpp>
+#include <boost/histogram/algorithm/reduce.hpp>
+#include <boost/histogram/algorithm/sum.hpp>
+#include <boost/histogram/axis/ostream_operators.hpp>
+#include <boost/histogram/axis/regular.hpp>
+#include <boost/histogram/ostream_operators.hpp>
+#include <vector>
+#include "utility_histogram.hpp"
+
+using namespace boost::histogram;
+using namespace boost::histogram::algorithm;
+
+template <typename Tag>
+void run_tests() {
+  using regular = axis::regular<axis::transform::identity<>, axis::null_type>;
+  {
+    auto h = make(Tag(), regular(4, 1, 5), regular(3, -1, 2));
+    h(1, -1);
+    h(1, 0);
+    h(2, 0);
+    h(2, 1);
+    h(2, 1);
+    h(3, -1);
+    h(3, 1);
+    h(4, 1);
+    h(4, 1);
+    h(4, 1);
+
+    /*
+      matrix layout:
+      x ->
+    y 1 0 1 0
+    | 1 1 0 0
+    v 0 2 1 3
+    */
+
+    // should do nothing, index order does not matter
+    auto hr = reduce(h, shrink(1, -1, 2), rebin(0, 1));
+    BOOST_TEST_EQ(hr.rank(), 2);
+    BOOST_TEST_EQ(sum(hr), 10);
+    BOOST_TEST_EQ(hr.axis(0).size(), h.axis(0).size());
+    BOOST_TEST_EQ(hr.axis(1).size(), h.axis(1).size());
+    BOOST_TEST_EQ(hr.axis(0)[0].lower(), 1);
+    BOOST_TEST_EQ(hr.axis(0)[3].upper(), 5);
+    BOOST_TEST_EQ(hr.axis(1)[0].lower(), -1);
+    BOOST_TEST_EQ(hr.axis(1)[2].upper(), 2);
+    BOOST_TEST_EQ(hr, h);
+
+    // not allowed: repeated indices
+    BOOST_TEST_THROWS(reduce(h, rebin(0, 2), rebin(0, 2)), std::invalid_argument);
+    BOOST_TEST_THROWS(reduce(h, shrink(1, 0, 2), shrink(1, 0, 2)), std::invalid_argument);
+    // not allowed: shrink with lower == upper
+    BOOST_TEST_THROWS(reduce(h, shrink(0, 0, 0)), std::invalid_argument);
+    // not allowed: shrink axis to zero size
+    BOOST_TEST_THROWS(reduce(h, shrink(0, 10, 11)), std::invalid_argument);
+    // not allowed: rebin with zero merge
+    BOOST_TEST_THROWS(reduce(h, rebin(0, 0)), std::invalid_argument);
+
+    hr = reduce(h, shrink(0, 2, 4));
+    BOOST_TEST_EQ(hr.rank(), 2);
+    BOOST_TEST_EQ(sum(hr), 10);
+    BOOST_TEST_EQ(hr.axis(0).size(), 2);
+    BOOST_TEST_EQ(hr.axis(1).size(), h.axis(1).size());
+    BOOST_TEST_EQ(hr.axis(0)[0].lower(), 2);
+    BOOST_TEST_EQ(hr.axis(0)[1].upper(), 4);
+    BOOST_TEST_EQ(hr.axis(1)[0].lower(), -1);
+    BOOST_TEST_EQ(hr.axis(1)[2].upper(), 2);
+    BOOST_TEST_EQ(hr.at(-1, 0), 1); // underflow
+    BOOST_TEST_EQ(hr.at(0, 0), 0);
+    BOOST_TEST_EQ(hr.at(1, 0), 1);
+    BOOST_TEST_EQ(hr.at(2, 0), 0); // overflow
+    BOOST_TEST_EQ(hr.at(-1, 1), 1);
+    BOOST_TEST_EQ(hr.at(0, 1), 1);
+    BOOST_TEST_EQ(hr.at(1, 1), 0);
+    BOOST_TEST_EQ(hr.at(2, 1), 0);
+    BOOST_TEST_EQ(hr.at(-1, 2), 0);
+    BOOST_TEST_EQ(hr.at(0, 2), 2);
+    BOOST_TEST_EQ(hr.at(1, 2), 1);
+    BOOST_TEST_EQ(hr.at(2, 2), 3);
+
+    /*
+      matrix layout:
+      x ->
+    y 1 0 1 0
+    | 1 1 0 0
+    v 0 2 1 3
+    */
+
+    hr = reduce(h, shrink_and_rebin(0, 2, 5, 2), rebin(1, 3));
+    BOOST_TEST_EQ(hr.rank(), 2);
+    BOOST_TEST_EQ(sum(hr), 10);
+    BOOST_TEST_EQ(hr.axis(0).size(), 1);
+    BOOST_TEST_EQ(hr.axis(1).size(), 1);
+    BOOST_TEST_EQ(hr.axis(0)[0].lower(), 2);
+    BOOST_TEST_EQ(hr.axis(0)[0].upper(), 4);
+    BOOST_TEST_EQ(hr.axis(1)[0].lower(), -1);
+    BOOST_TEST_EQ(hr.axis(1)[0].upper(), 2);
+    BOOST_TEST_EQ(hr.at(-1, 0), 2); // underflow
+    BOOST_TEST_EQ(hr.at(0, 0), 5);
+    BOOST_TEST_EQ(hr.at(1, 0), 3); // overflow
+
+    std::vector<reduce_option_type> opts{{shrink_and_rebin(0, 2, 5, 2), rebin(1, 3)}};
+    auto hr2 = reduce(h, opts);
+    BOOST_TEST_EQ(hr2, hr);
+  }
+
+  // rebin on integer axis must fail
+  {
+    auto h = make(Tag(), axis::integer<>(1, 4));
+    BOOST_TEST_THROWS(reduce(h, rebin(0, 2)), std::invalid_argument);
+  }
+
+  // reduce on axis with inverted range
+  {
+    auto h = make(Tag(), regular(4, 2, -2));
+    auto hr = reduce(h, shrink(0, 1, -1));
+    BOOST_TEST_EQ(hr.axis().size(), 2);
+    BOOST_TEST_EQ(hr.axis()[0].lower(), 1);
+    BOOST_TEST_EQ(hr.axis()[1].upper(), -1);
+  }
+
+  // reduce does not work with arguments not convertible to double
+}
+
+int main() {
+  run_tests<static_tag>();
+  run_tests<dynamic_tag>();
+
+  return boost::report_errors();
+}
diff --git a/test/axis_integer_test.cpp b/test/axis_integer_test.cpp
index f5df1e04..4b7a0c30 100644
--- a/test/axis_integer_test.cpp
+++ b/test/axis_integer_test.cpp
@@ -71,5 +71,23 @@ int main() {
   // iterators
   { test_axis_iterator(axis::integer<>(0, 4, ""), 0, 4); }
 
+  // shrink and rebin
+  {
+    using A = axis::integer<>;
+    auto a = A(0, 5);
+    auto b = A(a, 1, 4, 1);
+    BOOST_TEST_EQ(b.size(), 3);
+    BOOST_TEST_EQ(b.value(0), 1);
+    BOOST_TEST_EQ(b.value(3), 4);
+    auto c = A(a, 0, 4, 1);
+    BOOST_TEST_EQ(c.size(), 4);
+    BOOST_TEST_EQ(c.value(0), 0);
+    BOOST_TEST_EQ(c.value(4), 4);
+    auto e = A(a, 1, 4, 1);
+    BOOST_TEST_EQ(e.size(), 3);
+    BOOST_TEST_EQ(e.value(0), 1);
+    BOOST_TEST_EQ(e.value(3), 4);
+  }
+
   return boost::report_errors();
 }
diff --git a/test/axis_regular_test.cpp b/test/axis_regular_test.cpp
index 6dd3e764..929ac05b 100644
--- a/test/axis_regular_test.cpp
+++ b/test/axis_regular_test.cpp
@@ -153,5 +153,23 @@ int main() {
     test(a[0], "[0, 0.5)");
   }
 
+  // shrink and rebin
+  {
+    using A = axis::regular<>;
+    auto a = A(5, 0, 5);
+    auto b = A(a, 1, 4, 1);
+    BOOST_TEST_EQ(b.size(), 3);
+    BOOST_TEST_EQ(b.value(0), 1);
+    BOOST_TEST_EQ(b.value(3), 4);
+    auto c = A(a, 0, 4, 2);
+    BOOST_TEST_EQ(c.size(), 2);
+    BOOST_TEST_EQ(c.value(0), 0);
+    BOOST_TEST_EQ(c.value(2), 4);
+    auto e = A(a, 1, 5, 2);
+    BOOST_TEST_EQ(e.size(), 2);
+    BOOST_TEST_EQ(e.value(0), 1);
+    BOOST_TEST_EQ(e.value(2), 5);
+  }
+
   return boost::report_errors();
 }
diff --git a/test/axis_variant_test.cpp b/test/axis_variant_test.cpp
index 2de93512..aff2b249 100644
--- a/test/axis_variant_test.cpp
+++ b/test/axis_variant_test.cpp
@@ -13,6 +13,7 @@
 #include <boost/histogram/axis/regular.hpp>
 #include <boost/histogram/axis/variable.hpp>
 #include <boost/histogram/axis/variant.hpp>
+#include <functional> // for std::ref
 #include <sstream>
 #include <string>
 #include <type_traits>
@@ -35,6 +36,8 @@ int main() {
     BOOST_TEST_EQ(a(10), 2);
     BOOST_TEST_EQ(a[-1].lower(), -std::numeric_limits<double>::infinity());
     BOOST_TEST_EQ(a[a.size()].upper(), std::numeric_limits<double>::infinity());
+    BOOST_TEST_EQ(a[-10].lower(), -std::numeric_limits<double>::infinity());
+    BOOST_TEST_EQ(a[a.size() + 10].upper(), std::numeric_limits<double>::infinity());
     BOOST_TEST_EQ(a.metadata(), "int");
     BOOST_TEST_EQ(a.options(), axis::option_type::underflow_and_overflow);
 
@@ -200,13 +203,12 @@ int main() {
 
   // vector of axes with custom allocators
   {
-    struct null {};
     using M = std::vector<char, tracing_allocator<char>>;
     using T1 = axis::regular<axis::transform::identity<>, M>;
-    using T2 = axis::circular<double, null>;
-    using T3 = axis::variable<double, tracing_allocator<double>, null>;
-    using T4 = axis::integer<int, null>;
-    using T5 = axis::category<long, tracing_allocator<long>, null>;
+    using T2 = axis::circular<double, axis::null_type>;
+    using T3 = axis::variable<double, tracing_allocator<double>, axis::null_type>;
+    using T4 = axis::integer<int, axis::null_type>;
+    using T5 = axis::category<long, tracing_allocator<long>, axis::null_type>;
     using axis_type = axis::variant<T1, T2, T3, T4, T5>; // no heap allocation
     using axes_type = std::vector<axis_type, tracing_allocator<axis_type>>;
 
@@ -258,5 +260,19 @@ int main() {
   // iterators
   test_axis_iterator(axis::variant<axis::regular<>>(axis::regular<>(5, 0, 1)), 0, 5);
 
+  // variant of references
+  {
+    using A = axis::integer<int, axis::null_type>;
+    using VARef = axis::variant<A&>;
+    auto a = A(1, 5);
+    VARef ref(a);
+    BOOST_TEST_EQ(ref.size(), 4);
+    BOOST_TEST_EQ(ref.value(0), 1);
+    // change a through ref
+    axis::get<A>(ref) = A(7, 14);
+    BOOST_TEST_EQ(a.size(), 7);
+    BOOST_TEST_EQ(a.value(0), 7);
+  }
+
   return boost::report_errors();
 }
diff --git a/test/detail_test.cpp b/test/detail_test.cpp
index 5c6b915f..85c86ecb 100644
--- a/test/detail_test.cpp
+++ b/test/detail_test.cpp
@@ -13,6 +13,7 @@
 #include <boost/histogram/axis/variable.hpp>
 #include <boost/histogram/detail/axes.hpp>
 #include <boost/histogram/detail/cat.hpp>
+#include <boost/histogram/detail/is_set.hpp>
 #include <tuple>
 #include <vector>
 #include "utility_meta.hpp"
@@ -59,6 +60,16 @@ int main() {
     BOOST_TEST_NOT(detail::axes_equal(std_vector3, tuple3));
   }
 
+  // axes_size
+  {
+    std::tuple<int, int> a;
+    std::vector<int> b(3);
+    std::array<int, 4> c;
+    BOOST_TEST_EQ(detail::axes_size(a), 2);
+    BOOST_TEST_EQ(detail::axes_size(b), 3);
+    BOOST_TEST_EQ(detail::axes_size(c), 4);
+  }
+
   // sequence assign
   {
     enum { A, B, C, D };
@@ -107,5 +118,13 @@ int main() {
                   std::make_tuple(a0, a1, a2));
   }
 
+  // is_set
+  {
+    std::vector<int> v({3, 1, 2});
+    BOOST_TEST(detail::is_set(v.begin(), v.end()));
+    v = {3, 1, 2, 1};
+    BOOST_TEST_NOT(detail::is_set(v.begin(), v.end()));
+  }
+
   return boost::report_errors();
 }
diff --git a/test/index_mapper_test.cpp b/test/index_mapper_test.cpp
deleted file mode 100644
index bc591602..00000000
--- a/test/index_mapper_test.cpp
+++ /dev/null
@@ -1,70 +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/detail/index_mapper.hpp>
-#include <vector>
-
-using namespace boost::histogram::detail;
-
-int main() {
-  // index_mapper 1
-  {
-    // shape: 2, 3; 2, 0
-    index_mapper m(2);
-    m[0] = std::make_pair(1, 1);
-    m[1] = std::make_pair(2, 0);
-    m.ntotal = 6;
-    BOOST_TEST_EQ(m.first, 0);
-    BOOST_TEST_EQ(m.second, 0);
-    BOOST_TEST_EQ(m.next(), true);
-    BOOST_TEST_EQ(m.first, 1);
-    BOOST_TEST_EQ(m.second, 1);
-    BOOST_TEST_EQ(m.next(), true);
-    BOOST_TEST_EQ(m.first, 2);
-    BOOST_TEST_EQ(m.second, 0);
-    BOOST_TEST_EQ(m.next(), true);
-    BOOST_TEST_EQ(m.first, 3);
-    BOOST_TEST_EQ(m.second, 1);
-    BOOST_TEST_EQ(m.next(), true);
-    BOOST_TEST_EQ(m.first, 4);
-    BOOST_TEST_EQ(m.second, 0);
-    BOOST_TEST_EQ(m.next(), true);
-    BOOST_TEST_EQ(m.first, 5);
-    BOOST_TEST_EQ(m.second, 1);
-    BOOST_TEST_EQ(m.next(), false);
-  }
-
-  // index_mapper 2
-  {
-    // shape: 2, 3; 0, 3
-    index_mapper m(2);
-    m[0] = std::make_pair(1, 0);
-    m[1] = std::make_pair(2, 1);
-    m.ntotal = 6;
-    BOOST_TEST_EQ(m.first, 0);
-    BOOST_TEST_EQ(m.second, 0);
-    BOOST_TEST_EQ(m.next(), true);
-    BOOST_TEST_EQ(m.first, 1);
-    BOOST_TEST_EQ(m.second, 0);
-    BOOST_TEST_EQ(m.next(), true);
-    BOOST_TEST_EQ(m.first, 2);
-    BOOST_TEST_EQ(m.second, 1);
-    BOOST_TEST_EQ(m.next(), true);
-    BOOST_TEST_EQ(m.first, 3);
-    BOOST_TEST_EQ(m.second, 1);
-    BOOST_TEST_EQ(m.next(), true);
-    BOOST_TEST_EQ(m.first, 4);
-    BOOST_TEST_EQ(m.second, 2);
-    BOOST_TEST_EQ(m.next(), true);
-    BOOST_TEST_EQ(m.first, 5);
-    BOOST_TEST_EQ(m.second, 2);
-    BOOST_TEST_EQ(m.next(), false);
-  }
-
-  return boost::report_errors();
-}