gf_base.hpp

Go to the documentation of this file.

/*
 * Copyright (C) 1998-2018 ALPS Collaboration. See COPYRIGHT.TXT
 * All rights reserved. Use is subject to license terms. See LICENSE.TXT
 * For use in publications, see ACKNOWLEDGE.TXT
 */

#ifndef ALPSCORE_GF_H
#define ALPSCORE_GF_H


#include <tuple>
#include <vector>

#include <boost/preprocessor/cat.hpp>
#include <boost/preprocessor/repetition/repeat_from_to.hpp>

#include <alps/hdf5/archive.hpp>
#include <alps/hdf5/tensor.hpp>
#ifdef ALPS_HAVE_MPI
#include <alps/utilities/mpi.hpp>
#endif

#include <alps/gf/mesh.hpp>
#include <alps/numeric/tensors/tensor_base.hpp>
#include <alps/type_traits/index_sequence.hpp>
#include <alps/type_traits/tuple_traits.hpp>


namespace alps {
  namespace gf {
    namespace detail {
      template<class VTYPE, class Storage, class ...MESHES>
      class gf_base;

      template <typename Tuple, std::size_t...Is>
      std::ostream & output(std::ostream & os, const Tuple& t, index_sequence<Is...>) {
        using swallow = int[];
        (void)swallow{0, (void(os << std::get<Is>(t)), 0)...};
        return os;
      };

      template <typename ...Args>
      auto operator<<(std::ostream & os, const std::tuple<Args...>& t) ->
                                                DECLTYPE(detail::output(os, t, make_index_sequence<sizeof...(Args)>()))

      template<typename T> inline void print_no_complex(std::ostream &os, const T &z){
        os<<z;
      }
      template<> inline void print_no_complex(std::ostream &os, const std::complex<double> &z){
        //specialization for printing complex as ... real imag ...
        os<<z.real()<<" "<<z.imag();
      }
    }

    template<class VTYPE, class ...MESHES>
    using greenf = detail::gf_base<VTYPE, numerics::tensor<VTYPE, sizeof...(MESHES)>, MESHES...>;
    template<class VTYPE, class ...MESHES>
    using greenf_view = detail::gf_base<VTYPE, numerics::tensor_view<VTYPE, sizeof...(MESHES)>, MESHES...>;

    namespace detail {
      template<class VTYPE, class Storage, class ...MESHES>
      class gf_base {
        // types definition
      public:
        // Current GF types
        using value_type   = VTYPE;
        using storage_type = Storage;
        using mesh_types = std::tuple < MESHES... >;
        using index_types = std::tuple < typename MESHES::index_type... >;
      private:
        using gf_type   = gf_base < VTYPE, Storage, MESHES... >;
        using data_storage = numerics::tensor < VTYPE, sizeof...(MESHES) >;
        using data_view    = numerics::tensor_view < VTYPE, sizeof...(MESHES) >;
        template<typename St>
        using generic_gf   = gf_base < VTYPE, St, MESHES... >;
        template<typename RHS_VTYPE>
        using gf_op_type   = gf_base < decltype(RHS_VTYPE{} + VTYPE{}), numerics::tensor<decltype(RHS_VTYPE{} + VTYPE{}),
            sizeof...(MESHES)>, MESHES... >;

        // fields definition
      private:
        static constexpr int minor_version = 1;
        static constexpr int major_version = 0;
        static constexpr int N_ = sizeof...(MESHES);
        Storage data_;
        std::tuple < MESHES... > meshes_;
        bool empty_;

        // template hacks
        template <typename S, class Tup> struct subpack_impl;
        template<typename S, template<class...> class Tup, class... T>
        struct subpack_impl<S, Tup<T...> >
        {
          using type = gf_base<S, numerics::tensor_view < S, sizeof...(T)>,  T...>;
        };
        template<typename S, size_t I>
        using subpack = typename subpack_impl<S, decltype(tuple_tail < I, MESHES... >(meshes_) )>::type;

        template<typename RHS_VTYPE, typename LHS_VTYPE>
        struct convert {
          typedef typename std::conditional<std::is_integral<RHS_VTYPE>::value, VTYPE, RHS_VTYPE>::type type;
        };

        template<typename ...IndicesTypes>
        struct check_mesh {
          using type = typename std::conditional<std::is_convertible<std::tuple<IndicesTypes...>, index_types>::value,
              std::true_type, std::false_type>::type;
        };

      public:

        gf_base(const gf_type &g) = default;

        gf_base(gf_type &&g) = default;

        gf_base() : data_(std::array < size_t, N_ >{{0}}), empty_(true) {}

        gf_base(MESHES...meshes) : gf_base(std::make_tuple(meshes...)) {}
        gf_base(const mesh_types &meshes) : data_(get_sizes(meshes)), meshes_(meshes), empty_(false) {}
        gf_base(VTYPE* data, const mesh_types &meshes) : data_(data, get_sizes(meshes)), meshes_(meshes), empty_(false) {}
        gf_base(const data_storage &data, const mesh_types &meshes) : data_(data), meshes_(meshes), empty_(false) {}
        gf_base(data_storage &&data, const mesh_types &meshes) : data_(data), meshes_(meshes), empty_(false) {}
        gf_base(const data_storage &data, MESHES...meshes) : data_(data), meshes_(std::make_tuple(meshes...)), empty_(false) {}
        gf_base(data_storage && data, MESHES...meshes) : data_(std::move(data)), meshes_(std::make_tuple(meshes...)), empty_(false) {}

        template<typename St, typename = std::enable_if<!std::is_same<St, Storage>::value && std::is_same<St, data_view>::value > >
        gf_base(const gf_base<VTYPE, St, MESHES...> &g) : data_(g.data()), meshes_(g.meshes()), empty_(g.is_empty()) {}
        template<typename St, typename = std::enable_if<!std::is_same<St, Storage>::value && std::is_same<St, data_storage>::value > >
        gf_base(gf_base<VTYPE, St, MESHES...> &g) : data_(g.data()), meshes_(g.meshes()), empty_(g.is_empty()) {}

        template<typename RHS_VTYPE, typename St, typename = std::enable_if<std::is_same<data_storage, Storage>::value> >
        gf_base(const gf_base<RHS_VTYPE, St, MESHES...> &g) : data_(g.data()), meshes_(g.meshes()), empty_(g.is_empty()) {
          static_assert(std::is_convertible<RHS_VTYPE, VTYPE>::value, "Right-hand side data type is not convertible into left-hand side.");
        }
        template<typename RHS_VTYPE, typename St, typename = std::enable_if<std::is_same<data_storage, Storage>::value>>
        gf_base(gf_base<RHS_VTYPE, St, MESHES...> &&g) noexcept : data_(std::move(g.data())), meshes_(std::move(g.meshes())), empty_(g.is_empty()) {
          static_assert(std::is_convertible<RHS_VTYPE, VTYPE>::value, "Right-hand side data type is not convertible into left-hand side.");
        }

        template<typename St, typename...OLDMESHES, typename Index, typename ...Indices>
        gf_base(gf_base<VTYPE, numerics::detail::tensor_base < VTYPE, sizeof...(OLDMESHES), St >, OLDMESHES...> & g,
                std::tuple<OLDMESHES...>& oldmesh, const mesh_types &meshes, const Index ind, const Indices... idx) :
          data_(g.data()(ind(), idx()...)), meshes_(meshes), empty_(false) {}

        template<typename RHSTYPE ,typename St, typename...OLDMESHES, typename ...Indices>
        gf_base(const gf_base<RHSTYPE, numerics::detail::tensor_base < RHSTYPE, sizeof...(OLDMESHES), St >, OLDMESHES...> & g,
                const std::tuple<OLDMESHES...>& oldmesh, mesh_types &meshes, const Indices... idx) :
            data_(g.data()(idx()...)), meshes_(meshes), empty_(false) {}

        gf_type& operator=(const gf_type & rhs) {
          swap_meshes(rhs.meshes_, make_index_sequence<sizeof...(MESHES)>());
          data_ = rhs.data();
          empty_= rhs.empty_;
          return *this;
        }
        gf_type& operator=(gf_type && rhs) noexcept {
          swap_meshes(rhs.meshes_, make_index_sequence<sizeof...(MESHES)>());
          using std::swap;
          swap(data_, rhs.data_);
          swap(empty_, rhs.empty_);
          return *this;
        }

        void initialize() {
          data_.set_zero();
        }

        void check_meshes(const gf_type& rhs) const
        {
          if (meshes_ != rhs.meshes_) {
            throw std::invalid_argument("Green Functions have incompatible meshes");
          }
        }

        /*
         * Data access operations
         */
        template<class...Indices>
        typename std::enable_if < (sizeof...(Indices) == sizeof...(MESHES)), const VTYPE & >::type
        operator()(Indices...inds) const {
          // check that index types are the same as mesh indices
          static_assert(check_mesh<Indices...>::type::value, "Index type is inconsistent with mesh index type.");
          return value(std::forward < Indices >(inds)...);
        }

        template<class...Indices>
        typename std::enable_if < (sizeof...(Indices) == sizeof...(MESHES)), VTYPE & >::type
        operator()(Indices...inds) {
          // check that index types are the same as mesh indices
          static_assert(check_mesh<Indices...>::type::value, "Index type is inconsistent with mesh index type.");
          return value(std::forward < Indices >(inds)...);
        }

        template<class...Indices>
        auto operator()(typename std::enable_if<(sizeof...(Indices)+1 < N_),
            typename std::tuple_element<0,mesh_types>::type::index_type >::type ind, Indices...inds) ->
                                                                            subpack<VTYPE, sizeof...(Indices) + 1> {
          return subpack<VTYPE, sizeof...(Indices) + 1 >
                          (*this, meshes_, tuple_tail < sizeof...(Indices) + 1, MESHES...>(meshes_), ind, std::forward<Indices>(inds)...);
        }

        template<class...Indices>
        auto operator()(typename std::enable_if<(sizeof...(Indices)+1 < N_),
            typename std::tuple_element<0,mesh_types>::type::index_type >::type ind, Indices...inds) const ->
                                                                            subpack<const VTYPE, sizeof...(Indices) + 1> {
          auto t = tuple_tail < sizeof...(Indices) + 1, MESHES...>(meshes_);
          return subpack<const VTYPE, sizeof...(Indices) + 1>  (*this, meshes_, t, ind, std::forward<Indices>(inds)...);
        }


        /*
         * Basic GF arithmetics
         */

        template<typename RHS_VTYPE, typename RHS_STORAGE>
        gf_op_type<RHS_VTYPE> operator+(const gf_base<RHS_VTYPE, RHS_STORAGE, MESHES...> &rhs) const {
          throw_if_empty();
          gf_op_type<RHS_VTYPE> res(*this);
          return res += rhs;
        }

        template<typename RHS_GF>
        typename std::enable_if < std::is_convertible < RHS_GF, generic_gf<data_storage>>::value ||
                                  std::is_convertible < RHS_GF, generic_gf<data_view>>::value, gf_type & >::type
        operator+=(const RHS_GF &rhs) {
          throw_if_empty();
          data_ += rhs.data();
          return *this;
        }

        template<typename RHS_VTYPE, typename RHS_STORAGE>
        gf_op_type<RHS_VTYPE> operator-(const gf_base<RHS_VTYPE, RHS_STORAGE, MESHES...> &rhs) const {
          throw_if_empty();
          gf_op_type<RHS_VTYPE> res(*this);
          return res -= rhs;
        }

        template<typename RHS_GF>
        typename std::enable_if < std::is_convertible < RHS_GF, generic_gf<data_storage>>::value ||
                                  std::is_convertible < RHS_GF, generic_gf<data_view>>::value, gf_type & >::type
        operator-=(const RHS_GF &rhs) {
          throw_if_empty();
          data_ -= rhs.data();
          return *this;
        }

        template<typename RHS>
        typename std::enable_if < std::is_scalar < RHS >::value || std::is_same < VTYPE, RHS >::value, gf_type & >::type operator*=(RHS rhs) {
          throw_if_empty();
          data_ *= VTYPE(rhs);
          return *this;
        }

        template<typename RHS>
        gf_op_type<typename convert<RHS,VTYPE>::type> operator*(RHS rhs) const {
          throw_if_empty();
          gf_op_type<typename convert<RHS,VTYPE>::type> res(*this);
          return res *= rhs;
        }

        template<typename LHS>
        friend gf_op_type<typename convert<LHS,VTYPE>::type> operator*(LHS lhs, const gf_type & g) {
          return g * lhs;
        }

        template<typename RHS>
        typename std::enable_if < std::is_scalar < RHS >::value || std::is_same < VTYPE, RHS >::value, gf_type & >::type operator/=(RHS rhs) {
          throw_if_empty();
          data_ /= VTYPE(rhs);
          return *this;
        }

        template<typename RHS>
        gf_op_type<typename convert<RHS,VTYPE>::type> operator/(RHS rhs) const {
          throw_if_empty();
          gf_op_type<typename convert<RHS,VTYPE>::type> res(*this);
          return res /= rhs;
        }

        template<typename LHS>
        friend gf_op_type<typename convert<LHS,VTYPE>::type> operator/(LHS lhs, const gf_type & g) {
          return g / lhs;
        }

        gf_base < VTYPE, numerics::tensor<VTYPE, N_>, MESHES... > operator-() const {
          throw_if_empty();
          return (*this)*(VTYPE(-1.0));
        }

        template<typename RHS_GF>
        typename std::enable_if < std::is_convertible < RHS_GF, generic_gf<data_storage> >::value
                               || std::is_convertible < RHS_GF, generic_gf<data_view> >::value, bool >::type
        operator==(const RHS_GF &rhs) const {
          return (empty_ && rhs.is_empty()) || (data_.shape() == rhs.data().shape() && data_.storage() == rhs.data().storage() );
        }

        template<typename RHS_GF>
        bool operator!=(const RHS_GF &rhs) const {
          return !(*this==rhs);
        }

        double norm() const {
          throw_if_empty();
          return std::abs(*std::max_element(data_.data(), data_.num_elements() + data_.data(),
                                            [](VTYPE a, VTYPE b) {return std::abs(a) < std::abs(b);} ) );
        }

        // reshape green's function
        void reshape(MESHES...meshes) {
          std::array<size_t, N_> new_shape = get_sizes(std::make_tuple(meshes...));
          size_t new_size = std::accumulate(new_shape.begin(), new_shape.end(), size_t(1), std::multiplies<size_t>());
          if(std::is_same<data_view, Storage>::value) {
            if (data_.size() != new_size) {
              throw std::invalid_argument("The total size of resulting Green's function view should be the same as before.");
            }
          }
          data_.reshape(new_shape);
          meshes_ = std::make_tuple(meshes...);
          empty_ = false;
        }

        void save(alps::hdf5::archive &ar) const {
          save(ar, ar.get_context());
        }

        void load(alps::hdf5::archive &ar) {
          load(ar, ar.get_context());
        }

        void save(alps::hdf5::archive &ar, const std::string &path) const {
          throw_if_empty();
          save_version(ar, path);
          ar[path + "/data"] << data_;
          ar[path + "/mesh/N"] << int(N_);
          save_meshes(ar, path, make_index_sequence<sizeof...(MESHES)>());
        }

        void load(alps::hdf5::archive &ar, const std::string &path) {
          if (!check_version(ar, path)) throw std::runtime_error("Incompatible archive version");
          int ndim;
          ar[path + "/mesh/N"] >> ndim;
          if (ndim != N_) throw std::runtime_error("Wrong number of dimension reading GF, ndim=" + std::to_string(ndim)
                                                   + ", should be N=" + std::to_string(N_));
          load_meshes(ar, path, make_index_sequence<sizeof...(MESHES)>());
          data_ = numerics::tensor < VTYPE, N_ >(get_sizes(meshes_));
          ar[path + "/data"] >> data_;
          empty_ = false;
        }

        void save_version(alps::hdf5::archive &ar, const std::string &path) const {
          std::string vp = path + "/version/";
          ar[vp + "minor"] << int(minor_version);
          ar[vp + "major"] << int(major_version);
          ar[vp + "reference"] << "https://github.com/ALPSCore/H5GF/blob/master/H5GF.rst";
          ar[vp + "originator"] << "ALPSCore GF library, see http://www.alpscore.org";
        }

        bool check_version(alps::hdf5::archive &ar, const std::string &path) const {
          std::string vp = path + "/version/";
          int ver;
          ar[vp + "major"] >> ver;
          return (major_version == ver);
        }

        /*
         * Access to storage
         */
        const Storage &data() const { return data_; };
        Storage &data() { return data_; };

        bool is_empty() const {
          return empty_;
        }

        const mesh_types& meshes() const {return meshes_;}

        /*
         *  MPI routines
         */
#ifdef ALPS_HAVE_MPI

        void broadcast(const alps::mpi::communicator& comm, int root) {
          // check that root GF has been initilized
          if(comm.rank() == root) throw_if_empty();
          // broadcast grid meshes
          broadcast_mesh(comm, root);
          size_t root_sz=data_.size();
          alps::mpi::broadcast(comm, root_sz, root);
          // as long as all grids have been broadcasted we can define tensor object
          if(comm.rank() != root) data_ = numerics::tensor < VTYPE, N_ >(get_sizes(meshes_));
          alps::mpi::broadcast(comm, &data_.storage().data(0), root_sz, root);
        }
#endif

      private:

        /*
         *  MPI routines
         */
#ifdef ALPS_HAVE_MPI

        void broadcast_mesh(const alps::mpi::communicator& comm, int root) {
          // start from the zeroth mesh
          bcast_mesh(comm, root, std::integral_constant < int, 0 >());
        }
        template<int M>
        void bcast_mesh(const alps::mpi::communicator& comm, int root, std::integral_constant < int, M > i) {
          std::get<M>(meshes_).broadcast(comm, root);
          bcast_mesh(comm, root, std::integral_constant < int, M+1 >());
        }
        // Until we reach the last mesh object
        void bcast_mesh(const alps::mpi::communicator& comm, int root, std::integral_constant < int, N_ - 1> i) {
          std::get<N_ - 1>(meshes_).broadcast(comm, root);
        }
#endif

        inline void throw_if_empty() const {
#ifndef NDEBUG
          if (empty_) {
            throw std::runtime_error("gf is empty");
          }
#endif
        }

        template<class ... Indices>
        const VTYPE & value(Indices...inds) const {
          return data_(inds()...);
        }

        template<class ... Indices>
        VTYPE &value(Indices...inds) {
          return data_(inds()...);
        }

        std::array < size_t, N_ > get_sizes(MESHES...meshes) {
          return get_sizes(std::forward_as_tuple(meshes...));
        };

        std::array < size_t, N_ > get_sizes(const mesh_types &meshes) {
          std::array < size_t, N_ > sizes;
          sizes = fill_sizes(meshes, make_index_sequence<sizeof...(MESHES)>());
          return sizes;
        };

        template<size_t...Is>
        inline std::array < size_t, N_ > fill_sizes(const mesh_types &grids, index_sequence<Is...>) {
          return {{size_t(std::get<Is>(grids).extent())...}};
        };

        template<size_t...Is>
        void swap_meshes(const mesh_types & old_meshes, index_sequence<Is...>) {
          std::tie(std::get < Is >(meshes_) = std::get < Is >(old_meshes)...);
        }

        template<size_t...Is>
        void save_meshes(alps::hdf5::archive &ar, const std::string &path, index_sequence<Is...>) const {
          std::tie(ar[path + "/mesh/" + std::to_string(Is+1)] << std::get < Is >(meshes_)...);
        }


        template<size_t...Is>
        void load_meshes(alps::hdf5::archive &ar, const std::string &path, index_sequence<Is...>) {
          std::tie(ar[path + "/mesh/" + std::to_string(Is+1)] >> std::get < Is >(meshes_)...);
        }

        struct void_type {
          using type = void;
        };
        // Argument types.
        template<size_t i>
        struct args {
          using type
          = typename std::conditional < (i <= N_),
            std::tuple_element < i - 1, mesh_types >,
            void_type>::type::type;
        };

      public:

        #define MESH_FUNCTION(z, num, c) \
        template<typename = typename std::enable_if< (N_ >= num)> >\
        typename std::add_lvalue_reference<const typename args<num>::type>::type mesh##num() const {\
          return std::get<int(num-1)>(meshes_); \
        }

        #define MESH_TYPES(z, num, c) \
        /*template<typename = typename std::enable_if< (N_ >= num)> >*/\
        typedef typename args<num>::type mesh##num##_type;
        // = typename args<num>::type;

        /*
         * I guess 10 functions would be enough
         */
        BOOST_PP_REPEAT_FROM_TO (1, 11, MESH_FUNCTION, int)
        BOOST_PP_REPEAT_FROM_TO (1, 11, MESH_TYPES, int)

        friend std::ostream &operator<<(std::ostream &os, const gf_type & G ){
          using detail::operator<<;
          os<<G.meshes();
          for(int i=0;i<G.mesh1().extent();++i){
            os<<(G.mesh1().points()[i])<<" ";
            size_t points = G.data().size()/G.data().shape()[0];
            for(size_t j = 0; j< points; ++j) {
              detail::print_no_complex<value_type>(os, G.data().data()[j + G.data().index(i)]);
              os<<" ";
            }
            os<<std::endl;
          }
          return os;
        }
      };
    }
  }
}

#endif //ALPSCORE_GF_H