genericvectortools.hh

#ifndef DUNE_MATRIX_VECTOR_GENERICVECTORTOOLS_HH
#define DUNE_MATRIX_VECTOR_GENERICVECTORTOOLS_HH

/** \file
 *  \brief Various tools for working with istl vectors of arbitrary nesting depth
 */

#include <bitset>
#include <iostream>
#include <vector>

#include <dune/common/classname.hh>
#include <dune/common/concept.hh>
#include <dune/common/fvector.hh>
#include <dune/common/hybridutilities.hh>
#include <dune/common/typetraits.hh>

#include <dune/matrix-vector/algorithm.hh>
#include <dune/matrix-vector/matrixtraits.hh>

//! \brief Various tools for working with istl vectors of arbitrary nesting depth
namespace Dune {
namespace MatrixVector {
namespace Generic { // TODO change namespace name

// forward declaration for helper struct
namespace Impl {
template <class Vector, typename Enable = void>
struct ScalarSwitch;
} // end namespace Impl

//! Write vector to given stream
template <class Vector>
void writeBinary(std::ostream& s, const Vector& v) {
  Impl::ScalarSwitch<Vector>::writeBinary(s, v);
}

//! Read vector from given stream
template <class Vector>
void readBinary(std::istream& s, Vector& v) {
  Impl::ScalarSwitch<Vector>::readBinary(s, v);
}

//! Truncate vector by given bit set
template <class Vector, class BitVector>
void truncate(Vector& v, const BitVector& tr) {
  Impl::ScalarSwitch<Vector>::truncate(v, tr);
}

//! Resize vector recursively to size of given vector/matrix
// Note: We need the rvalue reference here, because we might obtain a
// temporary recursively, e.g. in std::bitset the operator[] yields a
// std::bitset::reference temporary due to its specific implementation.
template <class Vector, class SizeProvider>
void resize(Vector&& v, const SizeProvider& sizeProvider) {
  Impl::ScalarSwitch<SizeProvider>::resize(v, sizeProvider);
}
// Overload for setting from a std::bitset. We need this because it
// does not provide the range access the current implementation relies on.
template <class Vector, long unsigned int n>
void resize(Vector& v, const std::bitset<n>&) {
  FieldVector<double, n> sizeProvider;
  Impl::ScalarSwitch<decltype(sizeProvider)>::resize(v, sizeProvider);
}

namespace Impl {

/// helper concepts
namespace Concept {
struct HasResize {
  template <class C>
  auto require(C&& c) -> decltype(c.resize(0));
};

struct HasSize {
  // Note: Using Hybrid::size instead could be more general here. Make
  //       sure to update the getSize correspondingly if necessary.
  template <class C>
  auto require(C&& c) -> decltype(c.size());
};

struct HasN {
  template <class C>
  auto require(C&& c) -> decltype(c.N());
};

} // namespace Concept

//! Extract the size from the SizeProvider with these getSize overloads.
/// Note: We use PriorityTags to enforce trial order. Call getSize without
/// PriorityTag via the priority tag forward helper for standard behaviour.
// look for size method
template <class SizeProvider,
          typename = std::enable_if_t<models<Concept::HasSize, SizeProvider>()>>
auto getSize(const SizeProvider& sizeProvider, PriorityTag<20>) {
  return sizeProvider.size();
}
// look for N method (matrices should have this)
template <class SizeProvider,
          typename = std::enable_if_t<models<Concept::HasN, SizeProvider>()>>
auto getSize(const SizeProvider& sizeProvider, PriorityTag<10>) {
  return sizeProvider.N();
}
// fail if none of the above succeeds
template <class SizeProvider>
auto getSize(const SizeProvider&, PriorityTag<0>) {
  DUNE_THROW(NotImplemented, "Size extraction not implemented.");
  return -1;
}
// priority tag forward helper
template <class SizeProvider>
auto getSize(const SizeProvider& sizeProvider) {
  return getSize(sizeProvider, PriorityTag<42>());
}

/// compile-time helper traits
template <class M>
constexpr auto isMatrix() {
  return Std::bool_constant<MatrixTraits<M>::isMatrix>();
}

template <class V>
constexpr auto isVector() {
  // Note: At the time this comment is written, the only use of
  // isVector is in an overload after checking for isMatrix with
  // lower priority tag, so isNotMatrix should not be needed here.
  // Nevertheless, we keep it for clarity.
  constexpr bool isNotMatrix = not decltype(isMatrix<V>())::value;
  constexpr bool hasSize = models<Concept::HasSize, V>();
  return Std::bool_constant<hasSize and isNotMatrix>();
}

//! recursion helper for scalars nested in iterables (vectors)
template <class Vector>
struct ScalarSwitch<Vector,
                    typename std::enable_if_t<not IsNumber<Vector>::value>> {
  using This = ScalarSwitch<Vector>;

  static void writeBinary(std::ostream& s, const Vector& v) {
    Hybrid::forEach(v, [&s](auto&& vi) { Generic::writeBinary(s, vi); });
  }

  static void readBinary(std::istream& s, Vector& v) {
    Hybrid::forEach(v, [&s](auto&& vi) { Generic::readBinary(s, vi); });
  }

  template <class BitVector>
  static void truncate(Vector& v, const BitVector& tr) {
    sparseRangeFor(
        v, [&tr](auto&& vi, auto&& i) { Generic::truncate(vi, tr[i]); });
  }

  // Note: The scalarSwitch is with respect to the sizeProvider to avoid
  //       failures due to differing nesting depth. It might even be
  //       intentional to have a sizeProvider with less nesting than the
  //       vector to be resized.
  template <class Resizeable>
  static void resize(Resizeable& v, const Vector& sizeProvider) {
    static_assert(not IsNumber<Resizeable>::value, "SizeProvider seems to have nesting depth that is not applicable to given vector type.");

    using namespace Hybrid;
    auto s = getSize(sizeProvider);
    ifElse(models<Concept::HasResize, Resizeable>(),
           [&](auto&& id) { id(v).resize(s); },
           [&](auto) {
             if (size(v) != s) {
               std::stringstream msg;
               msg << "No resize method found when trying to resize " << className(v) << " with size " << size(v) 
                   << " from " << className(sizeProvider) << " with size " << s << ". "
                   << "Note that you cannot resize statically sized vectors." << std::endl;
               DUNE_THROW(RangeError, msg.str());
             }
           });

    // recursively resize (in this branch of ScalarSwitch, all size providers are nested)
    sparseRangeFor(sizeProvider, [&](auto&& row, auto&& i) {
      Generic::resize(v[i], This::subSizeProvider(row));
    });
  }

protected:
  //! Obtain the subSizeProvider from a sizeProvider row element with these overloads.
  /// Note: We use PriorityTags to enforce trial order w.r.t. failure fallback.
  /// Additionally we use it to circumvent ambiguous overloading.
  /// Call subSizeProvider without PriorityTag via the priority tag forward helper for standard behaviour.
  // ... from a sparse/dynamic matrix
  template <class SizeProviderRow, typename = std::enable_if_t<isMatrix<Vector>() and not isTupleOrDerived<Vector>()>>
  static decltype(auto) subSizeProvider(SizeProviderRow&& row, PriorityTag<21>) {
    return *row.begin();
  }
  // .. from a tuple like/static matrix
  template <class SizeProviderRow, typename = std::enable_if_t<isMatrix<Vector>() and isTupleOrDerived<Vector>()>>
  static decltype(auto) subSizeProvider(SizeProviderRow&& row, PriorityTag<20>) {
    return row[Indices::_0];
  }
  // ... from a vector
  template <class SizeProviderRow, typename = std::enable_if_t<isVector<Vector>()>>
  static decltype(auto) subSizeProvider(SizeProviderRow&& row, PriorityTag<10>) {
    return row;
  }
  // fail when none of the above can be applied to avoid unexpected resize behaviour
  template <class SizeProviderRow>
  static void subSizeProvider(SizeProviderRow&&, PriorityTag<0>) {
    static_assert(AlwaysFalse<Vector>::value, "MatrixVector::Generic::resize does not seem to work properly on given SizeProvider type. If it is a matrix, maybe you need to explicitly inject Dune::MatrixVector::MatrixTraits<SizeProvider>::isMatrix (see dune-matrix-vector/dune/matrix-vector/matrixtraits.hh).");
  }
  // priority tag forward helper
  template <class SizeProviderRow>
  static decltype(auto) subSizeProvider(SizeProviderRow&& row) {
    return This::subSizeProvider(std::forward<SizeProviderRow>(row), PriorityTag<42>());
  }
};

//! recursion anchors for scalars
template <class Scalar>
struct ScalarSwitch<Scalar,
                    typename std::enable_if_t<IsNumber<Scalar>::value>> {
  static void writeBinary(std::ostream& s, const Scalar& v) {
    s.write(reinterpret_cast<const char*>(&v), sizeof(Scalar));
  }

  static void readBinary(std::istream& s, Scalar& v) {
    for(auto&& vi: v)
      s.read(reinterpret_cast<char*>(&v), sizeof(Scalar));
  }

  template <class BitVector>
  static void truncate(Scalar& v, const BitVector& tr) {
    if (tr)
      v = 0;
  }

  template <class Resizeable>
  static void resize(Resizeable&, const Scalar&) {}
};

} // end namespace Impl

} // end namespace Generic
} // end namespace MatrixVector
} // end namespace Dune

#endif // DUNE_MATRIX_VECTOR_GENERICVECTORTOOLS_HH