sparse_shape.h

Go to the documentation of this file.

/*
 *  This file is a part of TiledArray.
 *  Copyright (C) 2013  Virginia Tech
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 *  Justus Calvin
 *  Department of Chemistry, Virginia Tech
 *
 *  sparse_shape.h
 *  Jul 9, 2013
 *
 */
 
#ifndef TILEDARRAY_SPARSE_SHAPE_H__INCLUDED
#define TILEDARRAY_SPARSE_SHAPE_H__INCLUDED
 
#include <TiledArray/tensor.h>
#include <TiledArray/tensor/shift_wrapper.h>
#include <TiledArray/tensor/tensor_interface.h>
#include <TiledArray/tiled_range.h>
#include <TiledArray/val_array.h>
#include <typeinfo>
 
namespace TiledArray {
 
 
template <typename T>
class SparseShape {
 public:
  typedef SparseShape<T> SparseShape_;  
  typedef T value_type;                 
  using index1_type = TA_1INDEX_TYPE;
  static_assert(TiledArray::detail::is_scalar_v<T>,
                "SparseShape<T> only supports scalar numeric types for T");
  typedef typename Tensor<value_type>::size_type size_type;  
 
 private:
  // T must be a numeric type
  static_assert(std::is_floating_point<T>::value,
                "SparseShape template type T must be a floating point type");
 
  // Internal typedefs
  typedef detail::ValArray<value_type> vector_type;
 
  Tensor<value_type> tile_norms_;  
  mutable std::unique_ptr<Tensor<value_type>> tile_norms_unscaled_ =
      nullptr;  
  std::shared_ptr<vector_type>
      size_vectors_;  
  size_type zero_tile_count_;    
  static value_type threshold_;  
 
  template <typename Op>
  static vector_type recursive_outer_product(
      const vector_type* const size_vectors, const unsigned int dim,
      const Op& op) {
    vector_type result;
 
    if (dim == 1u) {
      // Construct a modified copy of size_vector[0]
      result = op(*size_vectors);
    } else {
      // Compute split the range and compute the outer products
      const unsigned int middle = (dim >> 1u) + (dim & 1u);
      const vector_type left =
          recursive_outer_product(size_vectors, middle, op);
      const vector_type right =
          recursive_outer_product(size_vectors + middle, dim - middle, op);
 
      // Compute the outer product of left and right
 
      result = vector_type(left.size() * right.size());
      result.outer_fill(left, right,
                        [](const value_type left, const value_type right) {
                          return left * right;
                        });
    }
 
    return result;
  }
 
  enum class ScaleBy { Volume, InverseVolume };
 
 
  template <ScaleBy ScaleBy_, bool Screen = true>
  static size_type scale_tile_norms(
      Tensor<T>& tile_norms,
      const vector_type* MADNESS_RESTRICT const size_vectors) {
    const unsigned int dim = tile_norms.range().rank();
    const value_type threshold = threshold_;
    madness::AtomicInt zero_tile_count;
    zero_tile_count = 0;
 
    if (dim == 1u) {
      // This is the easy case where the data is a vector and can be
      // normalized directly.
      math::inplace_vector_op(
          [threshold, &zero_tile_count](value_type& norm,
                                        const value_type size) {
            if (ScaleBy_ == ScaleBy::Volume)
              norm *= size;
            else
              norm /= size;
            if (Screen && norm < threshold) {
              norm = value_type(0);
              ++zero_tile_count;
            }
          },
          size_vectors[0].size(), tile_norms.data(), size_vectors[0].data());
    } else {
      // Here the normalization constants are computed and multiplied by the
      // norm data using a recursive, outer algorithm. This is done to
      // minimize temporary memory requirements, memory bandwidth, and work.
 
      auto noop = [](const vector_type& size_vector) -> const vector_type& {
        return size_vector;
      };
      auto inv_vec_op = [](const vector_type& size_vector) {
        return vector_type(size_vector, [](const value_type size) {
          return value_type(1) / size;
        });
      };
 
      // Compute the left and right outer products
      const unsigned int middle = (dim >> 1u) + (dim & 1u);
      const vector_type left =
          ScaleBy_ == ScaleBy::Volume
              ? recursive_outer_product(size_vectors, middle, noop)
              : recursive_outer_product(size_vectors, middle, inv_vec_op);
      const vector_type right =
          ScaleBy_ == ScaleBy::Volume
              ? recursive_outer_product(size_vectors + middle, dim - middle,
                                        noop)
              : recursive_outer_product(size_vectors + middle, dim - middle,
                                        inv_vec_op);
 
      math::outer(
          left.size(), right.size(), left.data(), right.data(),
          tile_norms.data(),
          [threshold, &zero_tile_count](value_type& norm, const value_type x,
                                        const value_type y) {
            norm *= x * y;
            if (Screen && norm < threshold) {
              norm = value_type(0);
              ++zero_tile_count;
            }
          });
    }
 
    return Screen ? zero_tile_count : 0;
  }
 
  static std::shared_ptr<vector_type> initialize_size_vectors(
      const TiledRange& trange) {
    // Allocate memory for size vectors
    const unsigned int dim = trange.tiles_range().rank();
    std::shared_ptr<vector_type> size_vectors(
        new vector_type[dim], std::default_delete<vector_type[]>());
 
    // Initialize the size vectors
    for (unsigned int i = 0ul; i != dim; ++i) {
      const size_type n = trange.data()[i].tiles_range().second -
                          trange.data()[i].tiles_range().first;
 
      size_vectors.get()[i] =
          vector_type(n, &(*trange.data()[i].begin()),
                      [](const TiledRange1::range_type& tile) {
                        return value_type(tile.second - tile.first);
                      });
    }
 
    return size_vectors;
  }
 
  std::shared_ptr<vector_type> perm_size_vectors(
      const Permutation& perm) const {
    const unsigned int n = tile_norms_.range().rank();
 
    // Allocate memory for the contracted size vectors
    std::shared_ptr<vector_type> result_size_vectors(
        new vector_type[n], std::default_delete<vector_type[]>());
 
    // Initialize the size vectors
    for (unsigned int i = 0u; i < n; ++i) {
      const unsigned int perm_i = perm[i];
      result_size_vectors.get()[perm_i] = size_vectors_.get()[i];
    }
 
    return result_size_vectors;
  }
 
  decltype(zero_tile_count_) compute_zero_tile_count() {
    decltype(zero_tile_count_) zero_tile_count = 0;
    for (auto&& n : tile_norms_) {
      if (n < threshold()) {
        ++zero_tile_count;
      }
    }
    return zero_tile_count;
  }
 
  SparseShape(const Tensor<T>& tile_norms,
              const std::shared_ptr<vector_type>& size_vectors,
              const size_type zero_tile_count)
      : tile_norms_(tile_norms),
        size_vectors_(size_vectors),
        zero_tile_count_(zero_tile_count) {}
 
 public:
 
  SparseShape() : tile_norms_(), size_vectors_(), zero_tile_count_(0ul) {}
 
 
  SparseShape(const value_type& tile_norm, const TiledRange& trange)
      : tile_norms_(trange.tiles_range(),
                    (tile_norm < threshold_ ? 0 : tile_norm)),
        size_vectors_(initialize_size_vectors(trange)),
        zero_tile_count_(tile_norm < threshold_ ? trange.tiles_range().area()
                                                : 0ul) {}
 
 
  SparseShape(const Tensor<value_type>& tile_norms, const TiledRange& trange,
              bool do_not_scale = false)
      : tile_norms_(tile_norms.clone()),
        size_vectors_(initialize_size_vectors(trange)),
        zero_tile_count_(0ul) {
    TA_ASSERT(!tile_norms_.empty());
    TA_ASSERT(tile_norms_.range() == trange.tiles_range());
 
    if (!do_not_scale) {
      zero_tile_count_ = scale_tile_norms<ScaleBy::InverseVolume>(
          tile_norms_, size_vectors_.get());
    } else {
      zero_tile_count_ = compute_zero_tile_count();
    }
  }
 
 
  template <typename SparseNormSequence,
            typename = std::enable_if_t<
                TiledArray::detail::has_member_function_begin_anyreturn<
                    std::decay_t<SparseNormSequence>>::value &&
                TiledArray::detail::has_member_function_end_anyreturn<
                    std::decay_t<SparseNormSequence>>::value>>
  SparseShape(const SparseNormSequence& tile_norms, const TiledRange& trange,
              bool do_not_scale = false)
      : tile_norms_(trange.tiles_range(), value_type(0)),
        size_vectors_(initialize_size_vectors(trange)),
        zero_tile_count_(trange.tiles_range().volume()) {
    const auto dim = tile_norms_.range().rank();
    for (const auto& pair_idx_norm : tile_norms) {
      auto compute_tile_volume = [dim, this, pair_idx_norm]() -> uint64_t {
        uint64_t tile_volume = 1;
        for (size_t d = 0; d != dim; ++d)
          tile_volume *= size_vectors_.get()[d].at(pair_idx_norm.first[d]);
        return tile_volume;
      };
      auto norm_per_element =
          do_not_scale ? pair_idx_norm.second
                       : (pair_idx_norm.second / compute_tile_volume());
      if (norm_per_element >= threshold()) {
        tile_norms_[pair_idx_norm.first] = norm_per_element;
        --zero_tile_count_;
      }
    }
  }
 
 
  SparseShape(World& world, const Tensor<value_type>& tile_norms,
              const TiledRange& trange, bool do_not_scale = false)
      : tile_norms_(tile_norms.clone()),
        size_vectors_(initialize_size_vectors(trange)),
        zero_tile_count_(0ul) {
    TA_ASSERT(!tile_norms_.empty());
    TA_ASSERT(tile_norms_.range() == trange.tiles_range());
 
    // reduce norm data from all processors
    world.gop.max(tile_norms_.data(), tile_norms_.size());
 
    if (!do_not_scale) {
      zero_tile_count_ = scale_tile_norms<ScaleBy::InverseVolume>(
          tile_norms_, size_vectors_.get());
      ;
    } else {
      zero_tile_count_ = compute_zero_tile_count();
    }
  }
 
 
  template <typename SparseNormSequence>
  SparseShape(World& world, const SparseNormSequence& tile_norms,
              const TiledRange& trange)
      : SparseShape(tile_norms, trange) {
    world.gop.max(tile_norms_.data(), tile_norms_.size());
    zero_tile_count_ = compute_zero_tile_count();
  }
 
 
  SparseShape(const SparseShape<T>& other)
      : tile_norms_(other.tile_norms_),
        tile_norms_unscaled_(
            other.tile_norms_unscaled_
                ? std::make_unique<decltype(tile_norms_)>(
                      other.tile_norms_unscaled_.get()->clone())
                : nullptr),
        size_vectors_(other.size_vectors_),
        zero_tile_count_(other.zero_tile_count_) {}
 
 
  SparseShape<T>& operator=(const SparseShape<T>& other) {
    tile_norms_ = other.tile_norms_;
    tile_norms_unscaled_ = other.tile_norms_unscaled_
                               ? std::make_unique<decltype(tile_norms_)>(
                                     other.tile_norms_unscaled_.get()->clone())
                               : nullptr;
    size_vectors_ = other.size_vectors_;
    zero_tile_count_ = other.zero_tile_count_;
    return *this;
  }
 
 
  bool validate(const Range& range) const {
    if (tile_norms_.empty()) return false;
    return (range == tile_norms_.range());
  }
 
 
  template <typename Index>
  bool is_zero(const Index& i) const {
    TA_ASSERT(!tile_norms_.empty());
    return tile_norms_[i] < threshold_;
  }
 
 
  static constexpr bool is_dense() { return false; }
 
 
  float sparsity() const {
    TA_ASSERT(!tile_norms_.empty());
    return float(zero_tile_count_) / float(tile_norms_.size());
  }
 
 
  static value_type threshold() { return threshold_; }
 
 
  static void threshold(const value_type thresh) { threshold_ = thresh; }
 
 
  template <typename Index>
  value_type operator[](const Index& index) const {
    TA_ASSERT(!tile_norms_.empty());
    return tile_norms_[index];
  }
 
 
  template <typename Op>
  SparseShape_ transform(Op&& op) const {
    Tensor<T> new_norms = op(tile_norms_);
    madness::AtomicInt zero_tile_count;
    zero_tile_count = 0;
 
    const value_type threshold = threshold_;
    auto apply_threshold = [threshold, &zero_tile_count](value_type& norm) {
      TA_ASSERT(norm >= value_type(0));
      if (norm < threshold) {
        norm = value_type(0);
        ++zero_tile_count;
      }
    };
 
    math::inplace_vector_op(apply_threshold, new_norms.range().volume(),
                            new_norms.data());
 
    return SparseShape_(std::move(new_norms), size_vectors_, zero_tile_count);
  }
 
 
  const Tensor<value_type>& data() const { return tile_norms_; }
 
 
  const Tensor<value_type>& tile_norms() const {
    if (tile_norms_unscaled_ == nullptr) {
      tile_norms_unscaled_ =
          std::make_unique<decltype(tile_norms_)>(tile_norms_.clone());
      [[maybe_unused]] auto should_be_zero =
          scale_tile_norms<ScaleBy::Volume, false>(*tile_norms_unscaled_,
                                                   size_vectors_.get());
      TA_ASSERT(should_be_zero == 0);
    }
    return *(tile_norms_unscaled_.get());
  }
 
 
  bool empty() const { return tile_norms_.empty(); }
 
 
  SparseShape_ mask(const SparseShape_& mask_shape) const {
    TA_ASSERT(!tile_norms_.empty());
    TA_ASSERT(!mask_shape.empty());
    TA_ASSERT(tile_norms_.range() == mask_shape.tile_norms_.range());
 
    const value_type threshold = threshold_;
    madness::AtomicInt zero_tile_count;
    zero_tile_count = zero_tile_count_;
    auto op = [threshold, &zero_tile_count](value_type left,
                                            const value_type right) {
      if (left >= threshold && right < threshold) {
        left = value_type(0);
        ++zero_tile_count;
      }
 
      return left;
    };
 
    Tensor<value_type> result_tile_norms =
        tile_norms_.binary(mask_shape.tile_norms_, op);
 
    return SparseShape_(result_tile_norms, size_vectors_, zero_tile_count);
  }
 
  // clang-format off
 
  // clang-format on
  template <typename Index1, typename Index2,
            typename = std::enable_if_t<detail::is_integral_range_v<Index1> &&
                                        detail::is_integral_range_v<Index2>>>
  SparseShape update_block(const Index1& lower_bound, const Index2& upper_bound,
                           const SparseShape& other) const {
    Tensor<value_type> result_tile_norms = tile_norms_.clone();
 
    auto result_tile_norms_blk =
        result_tile_norms.block(lower_bound, upper_bound);
    const value_type threshold = threshold_;
    madness::AtomicInt zero_tile_count;
    zero_tile_count = zero_tile_count_;
    result_tile_norms_blk.inplace_binary(
        other.tile_norms_,
        [threshold, &zero_tile_count](value_type& l, const value_type r) {
          // Update the zero tile count for the result
          if ((l < threshold) && (r >= threshold))
            ++zero_tile_count;
          else if ((l >= threshold) && (r < threshold))
            --zero_tile_count;
 
          // Update the tile norm value
          l = r;
        });
 
    return SparseShape_(result_tile_norms, size_vectors_, zero_tile_count);
  }
 
  // clang-format off
 
  // clang-format on
  template <typename Index1, typename Index2,
            typename = std::enable_if_t<std::is_integral_v<Index1> &&
                                        std::is_integral_v<Index2>>>
  SparseShape update_block(const std::initializer_list<Index1>& lower_bound,
                           const std::initializer_list<Index2>& upper_bound,
                           const SparseShape& other) const {
    return update_block<std::initializer_list<Index1>,
                        std::initializer_list<Index2>>(lower_bound, upper_bound,
                                                       other);
  }
 
  // clang-format off
 
  // clang-format on
  template <typename PairRange,
            typename = std::enable_if_t<detail::is_gpair_range_v<PairRange>>>
  SparseShape update_block(const PairRange& bounds,
                           const SparseShape& other) const {
    Tensor<value_type> result_tile_norms = tile_norms_.clone();
 
    auto result_tile_norms_blk = result_tile_norms.block(bounds);
    const value_type threshold = threshold_;
    madness::AtomicInt zero_tile_count;
    zero_tile_count = zero_tile_count_;
    result_tile_norms_blk.inplace_binary(
        other.tile_norms_,
        [threshold, &zero_tile_count](value_type& l, const value_type r) {
          // Update the zero tile count for the result
          if ((l < threshold) && (r >= threshold))
            ++zero_tile_count;
          else if ((l >= threshold) && (r < threshold))
            --zero_tile_count;
 
          // Update the tile norm value
          l = r;
        });
 
    return SparseShape_(result_tile_norms, size_vectors_, zero_tile_count);
  }
 
  // clang-format off
 
  // clang-format on
  template <typename Index,
            typename = std::enable_if_t<std::is_integral_v<Index>>>
  SparseShape update_block(
      const std::initializer_list<std::initializer_list<Index>>& bounds,
      const SparseShape& other) const {
    return update_block<std::initializer_list<std::initializer_list<Index>>>(
        bounds, other);
  }
 
 
  inline bool operator==(const SparseShape<T>& other) const {
    bool equal = this->zero_tile_count_ == other.zero_tile_count_;
    if (equal) {
      const unsigned int dim = tile_norms_.range().rank();
      for (unsigned d = 0; d != dim && equal; ++d) {
        equal =
            equal && (size_vectors_.get()[d] == other.size_vectors_.get()[d]);
      }
      if (equal) {
        equal = equal && (tile_norms_ == other.tile_norms_);
      }
    }
    return equal;
  }
 
 private:
 
  template <typename Index1, typename Index2,
            typename = std::enable_if_t<detail::is_integral_range_v<Index1> &&
                                        detail::is_integral_range_v<Index2>>>
  std::shared_ptr<vector_type> block_range(const Index1& lower_bound,
                                           const Index2& upper_bound) const {
    // Get the number dimensions of the shape
    const auto rank = tile_norms_.range().rank();
    std::shared_ptr<vector_type> size_vectors(
        new vector_type[rank], std::default_delete<vector_type[]>());
 
    int d = 0;
    using std::begin;
    using std::end;
    auto lower_it = begin(lower_bound);
    auto upper_it = begin(upper_bound);
    const auto lower_end = end(lower_bound);
    const auto upper_end = end(upper_bound);
    for (; lower_it != lower_end && upper_it != upper_end;
         ++d, ++lower_it, ++upper_it) {
      // Get the new range size
      const auto lower_d = *lower_it;
      const auto upper_d = *upper_it;
      const auto extent_d = upper_d - lower_d;
 
      // Check that the input indices are in range
      TA_ASSERT(lower_d >= tile_norms_.range().lobound(d));
      TA_ASSERT(lower_d < upper_d);
      TA_ASSERT(upper_d <= tile_norms_.range().upbound(d));
 
      // Construct the size vector for rank i
      size_vectors.get()[d] =
          vector_type(extent_d, size_vectors_.get()[d].data() + lower_d);
    }
    TA_ASSERT(d == rank);
 
    return size_vectors;
  }
 
 
  template <typename PairRange,
            typename = std::enable_if_t<detail::is_gpair_range_v<PairRange>>>
  std::shared_ptr<vector_type> block_range(const PairRange& bounds) const {
    // Get the number dimensions of the shape
    const auto rank = tile_norms_.range().rank();
    std::shared_ptr<vector_type> size_vectors(
        new vector_type[rank], std::default_delete<vector_type[]>());
 
    int d = 0;
    for (auto&& bound_d : bounds) {
      // Get the new range size
      const auto lower_d = detail::at(bound_d, 0);
      const auto upper_d = detail::at(bound_d, 1);
      const auto extent_d = upper_d - lower_d;
 
      // Check that the input indices are in range
      TA_ASSERT(lower_d >= tile_norms_.range().lobound(d));
      TA_ASSERT(lower_d < upper_d);
      TA_ASSERT(upper_d <= tile_norms_.range().upbound(d));
 
      // Construct the size vector for rank i
      size_vectors.get()[d] =
          vector_type(extent_d, size_vectors_.get()[d].data() + lower_d);
 
      ++d;
    }
    TA_ASSERT(d == rank);
 
    return size_vectors;
  }
 
  template <typename Op>
  static SparseShape_ make_block(
      const std::shared_ptr<vector_type>& size_vectors,
      const TensorConstView<value_type>& block_view, const Op& op) {
    // Copy the data from arg to result
    const value_type threshold = threshold_;
    madness::AtomicInt zero_tile_count;
    zero_tile_count = 0;
    auto copy_op = [threshold, &zero_tile_count, &op](
                       value_type& MADNESS_RESTRICT result,
                       const value_type arg) {
      result = op(arg);
      if (result < threshold) {
        ++zero_tile_count;
        result = value_type(0);
      }
    };
 
    // Construct the result norms tensor
    Tensor<value_type> result_norms(Range(block_view.range().extent()));
    result_norms.inplace_binary(shift(block_view), copy_op);
 
    return SparseShape(result_norms, size_vectors, zero_tile_count);
  }
 
 public:
 
  template <typename Index1, typename Index2,
            typename = std::enable_if_t<detail::is_integral_range_v<Index1> &&
                                        detail::is_integral_range_v<Index2>>>
  SparseShape block(const Index1& lower_bound,
                    const Index2& upper_bound) const {
    return make_block(block_range(lower_bound, upper_bound),
                      tile_norms_.block(lower_bound, upper_bound),
                      [](auto&& arg) { return arg; });
  }
 
 
  template <typename Index1, typename Index2,
            typename = std::enable_if_t<std::is_integral_v<Index1> &&
                                        std::is_integral_v<Index2>>>
  SparseShape block(const std::initializer_list<Index1>& lower_bound,
                    const std::initializer_list<Index2>& upper_bound) const {
    return this
        ->block<std::initializer_list<Index1>, std::initializer_list<Index2>>(
            lower_bound, upper_bound);
  }
 
 
  template <typename PairRange,
            typename = std::enable_if_t<detail::is_gpair_range_v<PairRange>>>
  SparseShape block(const PairRange& bounds) const {
    return make_block(block_range(bounds), tile_norms_.block(bounds),
                      [](auto&& arg) { return arg; });
  }
 
 
  template <typename Index,
            typename = std::enable_if_t<std::is_integral_v<Index>>>
  SparseShape block(
      const std::initializer_list<std::initializer_list<Index>>& bounds) const {
    return make_block(block_range(bounds), tile_norms_.block(bounds),
                      [](auto&& arg) { return arg; });
  }
 
 
  template <typename Index1, typename Index2, typename Scalar,
            typename = std::enable_if_t<detail::is_integral_range_v<Index1> &&
                                        detail::is_integral_range_v<Index2> &&
                                        detail::is_numeric_v<Scalar>>>
  SparseShape block(const Index1& lower_bound, const Index2& upper_bound,
                    const Scalar factor) const {
    const value_type abs_factor = to_abs_factor(factor);
    return make_block(block_range(lower_bound, upper_bound),
                      tile_norms_.block(lower_bound, upper_bound),
                      [&abs_factor](auto&& arg) { return abs_factor * arg; });
  }
 
 
  template <typename Index1, typename Index2, typename Scalar,
            typename = std::enable_if_t<std::is_integral_v<Index1> &&
                                        std::is_integral_v<Index2> &&
                                        detail::is_numeric_v<Scalar>>>
  SparseShape block(const std::initializer_list<Index1>& lower_bound,
                    const std::initializer_list<Index2>& upper_bound,
                    const Scalar factor) const {
    return this->block<std::initializer_list<Index1>,
                       std::initializer_list<Index2>, Scalar>(
        lower_bound, upper_bound, factor);
  }
 
 
  template <typename PairRange, typename Scalar,
            typename = std::enable_if_t<detail::is_numeric_v<Scalar> &&
                                        detail::is_gpair_range_v<PairRange>>>
  SparseShape block(const PairRange& bounds, const Scalar factor) const {
    const value_type abs_factor = to_abs_factor(factor);
    return make_block(block_range(bounds), tile_norms_.block(bounds),
                      [&abs_factor](auto&& arg) { return abs_factor * arg; });
  }
 
 
  template <typename Index, typename Scalar,
            typename = std::enable_if_t<detail::is_numeric_v<Scalar> &&
                                        std::is_integral_v<Index>>>
  SparseShape block(
      const std::initializer_list<std::initializer_list<Index>>& bounds,
      const Scalar factor) const {
    const value_type abs_factor = to_abs_factor(factor);
    return make_block(block_range(bounds), tile_norms_.block(bounds),
                      [&abs_factor](auto&& arg) { return abs_factor * arg; });
  }
 
 
  template <typename Index1, typename Index2,
            typename = std::enable_if_t<detail::is_integral_range_v<Index1> &&
                                        detail::is_integral_range_v<Index2>>>
  SparseShape block(const Index1& lower_bound, const Index2& upper_bound,
                    const Permutation& perm) const {
    return block(lower_bound, upper_bound).perm(perm);
  }
 
 
  template <typename Index1, typename Index2,
            typename = std::enable_if_t<std::is_integral_v<Index1> &&
                                        std::is_integral_v<Index2>>>
  SparseShape block(const std::initializer_list<Index1>& lower_bound,
                    const std::initializer_list<Index2>& upper_bound,
                    const Permutation& perm) const {
    return block(lower_bound, upper_bound).perm(perm);
  }
 
 
  template <typename PairRange,
            typename = std::enable_if_t<detail::is_gpair_range_v<PairRange>>>
  SparseShape block(const PairRange& bounds, const Permutation& perm) const {
    return block(bounds).perm(perm);
  }
 
 
  template <typename Index,
            typename = std::enable_if_t<std::is_integral_v<Index>>>
  SparseShape block(
      const std::initializer_list<std::initializer_list<Index>>& bounds,
      const Permutation& perm) const {
    return block(bounds).perm(perm);
  }
 
 
  template <typename Index1, typename Index2, typename Scalar,
            typename = std::enable_if_t<detail::is_integral_range_v<Index1> &&
                                        detail::is_integral_range_v<Index2> &&
                                        detail::is_numeric_v<Scalar>>>
  SparseShape block(const Index1& lower_bound, const Index2& upper_bound,
                    const Scalar factor, const Permutation& perm) const {
    return block(lower_bound, upper_bound, factor).perm(perm);
  }
 
 
  template <typename Index1, typename Index2, typename Scalar,
            typename = std::enable_if_t<std::is_integral_v<Index1> &&
                                        std::is_integral_v<Index2> &&
                                        detail::is_numeric_v<Scalar>>>
  SparseShape block(const std::initializer_list<Index1>& lower_bound,
                    const std::initializer_list<Index2>& upper_bound,
                    const Scalar factor, const Permutation& perm) const {
    return block(lower_bound, upper_bound, factor).perm(perm);
  }
 
 
  template <typename PairRange, typename Scalar,
            typename = std::enable_if_t<detail::is_numeric_v<Scalar> &&
                                        detail::is_gpair_range_v<PairRange>>>
  SparseShape block(const PairRange& bounds, const Scalar factor,
                    const Permutation& perm) const {
    const value_type abs_factor = to_abs_factor(factor);
    return make_block(block_range(bounds), tile_norms_.block(bounds),
                      [&abs_factor](auto&& arg) { return abs_factor * arg; })
        .perm(perm);
  }
 
 
  template <typename Index, typename Scalar,
            typename = std::enable_if_t<detail::is_numeric_v<Scalar> &&
                                        std::is_integral_v<Index>>>
  SparseShape block(
      const std::initializer_list<std::initializer_list<Index>>& bounds,
      const Scalar factor, const Permutation& perm) const {
    const value_type abs_factor = to_abs_factor(factor);
    return make_block(block_range(bounds), tile_norms_.block(bounds),
                      [&abs_factor](auto&& arg) { return abs_factor * arg; })
        .perm(perm);
  }
 
 
  SparseShape_ perm(const Permutation& perm) const {
    return SparseShape_(tile_norms_.permute(perm), perm_size_vectors(perm),
                        zero_tile_count_);
  }
 
 
  template <typename Scalar,
            typename = std::enable_if_t<detail::is_numeric_v<Scalar>>>
  SparseShape_ scale(const Scalar factor) const {
    TA_ASSERT(!tile_norms_.empty());
    const value_type threshold = threshold_;
    const value_type abs_factor = to_abs_factor(factor);
    madness::AtomicInt zero_tile_count;
    zero_tile_count = 0;
    auto op = [threshold, &zero_tile_count, abs_factor](value_type value) {
      value *= abs_factor;
      if (value < threshold) {
        value = value_type(0);
        ++zero_tile_count;
      }
      return value;
    };
 
    Tensor<value_type> result_tile_norms = tile_norms_.unary(op);
 
    return SparseShape_(result_tile_norms, size_vectors_, zero_tile_count);
  }
 
 
  template <typename Factor>
  SparseShape_ scale(const Factor factor, const Permutation& perm) const {
    TA_ASSERT(!tile_norms_.empty());
    const value_type threshold = threshold_;
    const value_type abs_factor = to_abs_factor(factor);
    madness::AtomicInt zero_tile_count;
    zero_tile_count = 0;
    auto op = [threshold, &zero_tile_count, abs_factor](value_type value) {
      value *= abs_factor;
      if (value < threshold) {
        value = value_type(0);
        ++zero_tile_count;
      }
      return value;
    };
 
    Tensor<value_type> result_tile_norms = tile_norms_.unary(op, perm);
 
    return SparseShape_(result_tile_norms, perm_size_vectors(perm),
                        zero_tile_count);
  }
 
 
  SparseShape_ add(const SparseShape_& other) const {
    TA_ASSERT(!tile_norms_.empty());
    const value_type threshold = threshold_;
    madness::AtomicInt zero_tile_count;
    zero_tile_count = 0;
    auto op = [threshold, &zero_tile_count](value_type left,
                                            const value_type right) {
      left += right;
      if (left < threshold) {
        left = value_type(0);
        ++zero_tile_count;
      }
      return left;
    };
 
    Tensor<value_type> result_tile_norms =
        tile_norms_.binary(other.tile_norms_, op);
 
    return SparseShape_(result_tile_norms, size_vectors_, zero_tile_count);
  }
 
 
  SparseShape_ add(const SparseShape_& other, const Permutation& perm) const {
    TA_ASSERT(!tile_norms_.empty());
    const value_type threshold = threshold_;
    madness::AtomicInt zero_tile_count;
    zero_tile_count = 0;
    auto op = [threshold, &zero_tile_count](value_type left,
                                            const value_type right) {
      left += right;
      if (left < threshold) {
        left = value_type(0);
        ++zero_tile_count;
      }
      return left;
    };
 
    Tensor<value_type> result_tile_norms =
        tile_norms_.binary(other.tile_norms_, op, perm);
 
    return SparseShape_(result_tile_norms, perm_size_vectors(perm),
                        zero_tile_count);
  }
 
 
  template <typename Factor>
  SparseShape_ add(const SparseShape_& other, const Factor factor) const {
    TA_ASSERT(!tile_norms_.empty());
    const value_type threshold = threshold_;
    const value_type abs_factor = to_abs_factor(factor);
    madness::AtomicInt zero_tile_count;
    zero_tile_count = 0;
    auto op = [threshold, &zero_tile_count, abs_factor](
                  value_type left, const value_type right) {
      left += right;
      left *= abs_factor;
      if (left < threshold) {
        left = value_type(0);
        ++zero_tile_count;
      }
      return left;
    };
 
    Tensor<value_type> result_tile_norms =
        tile_norms_.binary(other.tile_norms_, op);
 
    return SparseShape_(result_tile_norms, size_vectors_, zero_tile_count);
  }
 
 
  template <typename Factor>
  SparseShape_ add(const SparseShape_& other, const Factor factor,
                   const Permutation& perm) const {
    TA_ASSERT(!tile_norms_.empty());
    const value_type threshold = threshold_;
    const value_type abs_factor = to_abs_factor(factor);
    madness::AtomicInt zero_tile_count;
    zero_tile_count = 0;
    auto op = [threshold, &zero_tile_count, abs_factor](
                  value_type left, const value_type right) {
      left += right;
      left *= abs_factor;
      if (left < threshold) {
        left = value_type(0);
        ++zero_tile_count;
      }
      return left;
    };
 
    Tensor<value_type> result_tile_norms =
        tile_norms_.binary(other.tile_norms_, op, perm);
 
    return SparseShape_(result_tile_norms, perm_size_vectors(perm),
                        zero_tile_count);
  }
 
  SparseShape_ add(value_type value) const {
    TA_ASSERT(!tile_norms_.empty());
    const value_type threshold = threshold_;
    madness::AtomicInt zero_tile_count;
    zero_tile_count = 0;
 
    Tensor<T> result_tile_norms(tile_norms_.range());
 
    value = std::abs(value);
    const unsigned int dim = tile_norms_.range().rank();
    const vector_type* MADNESS_RESTRICT const size_vectors =
        size_vectors_.get();
 
    if (dim == 1u) {
      auto add_const_op = [threshold, &zero_tile_count, value](
                              value_type norm, const value_type size) {
        norm += value / std::sqrt(size);
        if (norm < threshold) {
          norm = 0;
          ++zero_tile_count;
        }
        return norm;
      };
 
      // This is the easy case where the data is a vector and can be
      // normalized directly.
      math::vector_op(add_const_op, size_vectors[0].size(),
                      result_tile_norms.data(), tile_norms_.data(),
                      size_vectors[0].data());
 
    } else {
      // Here the normalization constants are computed and multiplied by the
      // norm data using a recursive, outer algorithm. This is done to
      // minimize temporary memory requirements, memory bandwidth, and work.
 
      auto inv_sqrt_vec_op = [](const vector_type size_vector) {
        return vector_type(size_vector, [](const value_type size) {
          return value_type(1) / std::sqrt(size);
        });
      };
 
      // Compute the left and right outer products
      const unsigned int middle = (dim >> 1u) + (dim & 1u);
      const vector_type left =
          recursive_outer_product(size_vectors, middle, inv_sqrt_vec_op);
      const vector_type right = recursive_outer_product(
          size_vectors + middle, dim - middle, inv_sqrt_vec_op);
 
      math::outer_fill(
          left.size(), right.size(), left.data(), right.data(),
          tile_norms_.data(), result_tile_norms.data(),
          [threshold, &zero_tile_count, value](
              value_type& norm, const value_type x, const value_type y) {
            norm += value * x * y;
            if (norm < threshold) {
              norm = value_type(0);
              ++zero_tile_count;
            }
          });
    }
 
    return SparseShape_(result_tile_norms, size_vectors_, zero_tile_count);
  }
 
  SparseShape_ add(const value_type value, const Permutation& perm) const {
    // TODO: Optimize this function so it does the permutation at the same
    // time as the addition.
    return add(value).perm(perm);
  }
 
  SparseShape_ subt(const SparseShape_& other) const { return add(other); }
 
  SparseShape_ subt(const SparseShape_& other, const Permutation& perm) const {
    return add(other, perm);
  }
 
  template <typename Factor>
  SparseShape_ subt(const SparseShape_& other, const Factor factor) const {
    return add(other, factor);
  }
 
  template <typename Factor>
  SparseShape_ subt(const SparseShape_& other, const Factor factor,
                    const Permutation& perm) const {
    return add(other, factor, perm);
  }
 
  SparseShape_ subt(const value_type value) const { return add(value); }
 
  SparseShape_ subt(const value_type value, const Permutation& perm) const {
    return add(value, perm);
  }
 
  SparseShape_ mult(const SparseShape_& other) const {
    // TODO: Optimize this function so that the tensor arithmetic and
    // scale_tile_norms operations are performed in one step instead of two.
 
    TA_ASSERT(!tile_norms_.empty());
    Tensor<T> result_tile_norms = tile_norms_.mult(other.tile_norms_);
    const size_type zero_tile_count = scale_tile_norms<ScaleBy::Volume>(
        result_tile_norms, size_vectors_.get());
 
    return SparseShape_(result_tile_norms, size_vectors_, zero_tile_count);
  }
 
  SparseShape_ mult(const SparseShape_& other, const Permutation& perm) const {
    // TODO: Optimize this function so that the tensor arithmetic and
    // scale_tile_norms operations are performed in one step instead of two.
 
    TA_ASSERT(!tile_norms_.empty());
    Tensor<T> result_tile_norms = tile_norms_.mult(other.tile_norms_, perm);
    std::shared_ptr<vector_type> result_size_vector = perm_size_vectors(perm);
    const size_type zero_tile_count = scale_tile_norms<ScaleBy::Volume>(
        result_tile_norms, result_size_vector.get());
 
    return SparseShape_(result_tile_norms, result_size_vector, zero_tile_count);
  }
 
  template <typename Factor>
  SparseShape_ mult(const SparseShape_& other, const Factor factor) const {
    // TODO: Optimize this function so that the tensor arithmetic and
    // scale_tile_norms operations are performed in one step instead of two.
 
    TA_ASSERT(!tile_norms_.empty());
    const value_type abs_factor = to_abs_factor(factor);
    Tensor<T> result_tile_norms =
        tile_norms_.mult(other.tile_norms_, abs_factor);
    const size_type zero_tile_count = scale_tile_norms<ScaleBy::Volume>(
        result_tile_norms, size_vectors_.get());
 
    return SparseShape_(result_tile_norms, size_vectors_, zero_tile_count);
  }
 
  template <typename Factor>
  SparseShape_ mult(const SparseShape_& other, const Factor factor,
                    const Permutation& perm) const {
    // TODO: Optimize this function so that the tensor arithmetic and
    // scale_tile_norms operations are performed in one step instead of two.
 
    TA_ASSERT(!tile_norms_.empty());
    const value_type abs_factor = to_abs_factor(factor);
    Tensor<T> result_tile_norms =
        tile_norms_.mult(other.tile_norms_, abs_factor, perm);
    std::shared_ptr<vector_type> result_size_vector = perm_size_vectors(perm);
    const size_type zero_tile_count = scale_tile_norms<ScaleBy::Volume>(
        result_tile_norms, result_size_vector.get());
 
    return SparseShape_(result_tile_norms, result_size_vector, zero_tile_count);
  }
 
  template <typename Factor>
  SparseShape_ gemm(const SparseShape_& other, const Factor factor,
                    const math::GemmHelper& gemm_helper) const {
    TA_ASSERT(!tile_norms_.empty());
 
    const value_type abs_factor = to_abs_factor(factor);
    const value_type threshold = threshold_;
    madness::AtomicInt zero_tile_count;
    zero_tile_count = 0;
    using integer = TiledArray::math::blas::integer;
    integer M = 0, N = 0, K = 0;
    gemm_helper.compute_matrix_sizes(M, N, K, tile_norms_.range(),
                                     other.tile_norms_.range());
 
    // Allocate memory for the contracted size vectors
    std::shared_ptr<vector_type> result_size_vectors(
        new vector_type[gemm_helper.result_rank()],
        std::default_delete<vector_type[]>());
 
    // Initialize the result size vectors
    unsigned int x = 0ul;
    for (unsigned int i = gemm_helper.left_outer_begin();
         i < gemm_helper.left_outer_end(); ++i, ++x)
      result_size_vectors.get()[x] = size_vectors_.get()[i];
    for (unsigned int i = gemm_helper.right_outer_begin();
         i < gemm_helper.right_outer_end(); ++i, ++x)
      result_size_vectors.get()[x] = other.size_vectors_.get()[i];
 
    // Compute the number of inner ranks
    const unsigned int k_rank =
        gemm_helper.left_inner_end() - gemm_helper.left_inner_begin();
 
    // Construct the result norm tensor
    Tensor<value_type> result_norms(
        gemm_helper.make_result_range<typename Tensor<T>::range_type>(
            tile_norms_.range(), other.tile_norms_.range()),
        0);
 
    if (k_rank > 0u) {
      // Compute size vector
      const vector_type k_sizes = recursive_outer_product(
          size_vectors_.get() + gemm_helper.left_inner_begin(), k_rank,
          [](const vector_type& size_vector) -> const vector_type& {
            return size_vector;
          });
 
      // TODO: Make this faster. It can be done without using temporaries
      // for the arguments, but requires a custom matrix multiply.
 
      Tensor<value_type> left(tile_norms_.range());
      const size_type mk = M * K;
      auto left_op = [](const value_type left, const value_type right) {
        return left * right;
      };
      for (size_type i = 0ul; i < mk; i += K)
        math::vector_op(left_op, K, left.data() + i, tile_norms_.data() + i,
                        k_sizes.data());
 
      Tensor<value_type> right(other.tile_norms_.range());
      for (integer i = 0ul, k = 0; k < K; i += N, ++k) {
        const value_type factor = k_sizes[k];
        auto right_op = [=](const value_type arg) { return arg * factor; };
        math::vector_op(right_op, N, right.data() + i,
                        other.tile_norms_.data() + i);
      }
 
      result_norms = left.gemm(right, abs_factor, gemm_helper);
 
      // Hard zero tiles that are below the zero threshold.
      result_norms.inplace_unary(
          [threshold, &zero_tile_count](value_type& value) {
            if (value < threshold) {
              value = value_type(0);
              ++zero_tile_count;
            }
          });
 
    } else {
      // This is an outer product, so the inputs can be used directly
      math::outer_fill(M, N, tile_norms_.data(), other.tile_norms_.data(),
                       result_norms.data(),
                       [threshold, &zero_tile_count, abs_factor](
                           const value_type left, const value_type right) {
                         value_type norm = left * right * abs_factor;
                         if (norm < threshold) {
                           norm = value_type(0);
                           ++zero_tile_count;
                         }
                         return norm;
                       });
    }
 
    return SparseShape_(result_norms, result_size_vectors, zero_tile_count);
  }
 
  template <typename Factor>
  SparseShape_ gemm(const SparseShape_& other, const Factor factor,
                    const math::GemmHelper& gemm_helper,
                    const Permutation& perm) const {
    return gemm(other, factor, gemm_helper).perm(perm);
  }
 
  template <typename Archive,
            typename std::enable_if<madness::archive::is_input_archive<
                Archive>::value>::type* = nullptr>
  void serialize(const Archive& ar) {
    ar& tile_norms_;
    const unsigned int dim = tile_norms_.range().rank();
    // allocate size_vectors_
    size_vectors_ = std::move(std::shared_ptr<vector_type>(
        new vector_type[dim], std::default_delete<vector_type[]>()));
    for (unsigned d = 0; d != dim; ++d) ar& size_vectors_.get()[d];
    ar& zero_tile_count_;
  }
 
  template <typename Archive,
            typename std::enable_if<madness::archive::is_output_archive<
                Archive>::value>::type* = nullptr>
  void serialize(const Archive& ar) const {
    ar& tile_norms_;
    const unsigned int dim = tile_norms_.range().rank();
    for (unsigned d = 0; d != dim; ++d) ar& size_vectors_.get()[d];
    ar& zero_tile_count_;
  }
 
 private:
  template <typename Factor>
  static value_type to_abs_factor(const Factor factor) {
    using std::abs;
    const auto cast_abs_factor = static_cast<value_type>(abs(factor));
    TA_ASSERT(std::isfinite(cast_abs_factor));
    return cast_abs_factor;
  }
 
};  // class SparseShape
 
// Static member initialization
template <typename T>
typename SparseShape<T>::value_type SparseShape<T>::threshold_ =
    std::numeric_limits<T>::epsilon();
 
 
template <typename T>
inline std::ostream& operator<<(std::ostream& os, const SparseShape<T>& shape) {
  os << "SparseShape<" << typeid(T).name() << ">:" << std::endl
     << shape.data() << std::endl;
  return os;
}
 
 
template <typename T>
bool is_replicated(World& world, const SparseShape<T>& shape) {
  const auto volume = shape.data().size();
  std::vector<T> data(shape.data().data(), shape.data().data() + volume);
  world.gop.max(data.data(), volume);
  for (size_t i = 0; i != data.size(); ++i) {
    if (data[i] != shape.data()[i]) return false;
  }
  return true;
}
 
#ifndef TILEDARRAY_HEADER_ONLY
 
extern template class SparseShape<float>;
 
#endif  // TILEDARRAY_HEADER_ONLY
 
}  // namespace TiledArray
 
#endif  // TILEDARRAY_SPASE_SHAPE_H__INCLUDED