/*******************************************************************************
*
* TRIQS: a Toolbox for Research in Interacting Quantum Systems
*
* Copyright (C) 2011-2013 by O. Parcollet
*
* TRIQS 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.
*
* TRIQS 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
* TRIQS. If not, see .
*
******************************************************************************/
#pragma once
#include "./common.hpp"
#include "../storages/shared_block.hpp"
#include "./assignment.hpp"
#include "./print.hpp"
#include "../indexmaps/cuboid/foreach.hpp"
#include "triqs/utility/exceptions.hpp"
#include "triqs/utility/typeid_name.hpp"
#include "triqs/utility/view_tools.hpp"
#include
#ifdef TRIQS_WITH_PYTHON_SUPPORT
#include "../python/numpy_extractor.hpp"
#include "../python/array_view_to_python.hpp"
#endif
namespace triqs { namespace arrays {
template AUTO_DECL get_shape (A const & x) RETURN(x.domain().lengths());
template size_t first_dim (A const & x) { return x.domain().lengths()[0];}
template size_t second_dim (A const & x) { return x.domain().lengths()[1];}
template size_t third_dim (A const & x) { return x.domain().lengths()[2];}
template size_t fourth_dim (A const & x) { return x.domain().lengths()[3];}
template size_t fifth_dim (A const & x) { return x.domain().lengths()[4];}
template size_t sixth_dim (A const & x) { return x.domain().lengths()[5];}
template size_t seventh_dim (A const & x) { return x.domain().lengths()[6];}
template size_t eighth_dim (A const & x) { return x.domain().lengths()[7];}
template size_t ninth_dim (A const & x) { return x.domain().lengths()[8];}
template class iterator_adapter;
template struct ISPViewType;
template
class indexmap_storage_pair : Tag::indexmap_storage_pair, TRIQS_CONCEPT_TAG_NAME(MutableCuboidArray) {
public :
using value_type=typename StorageType::value_type ;
static_assert(std::is_constructible::value, "array/array_view and const operate only on values");
static_assert(!std::is_const::value, "no const type");
using storage_type=StorageType ;
using indexmap_type=IndexMapType ;
using traversal_order_t = typename _get_traversal_order_t::type;
using traversal_order_arg = TraversalOrder;
static constexpr int rank = IndexMapType::domain_type::rank;
static constexpr bool is_const = IsConst;
protected:
indexmap_type indexmap_;
storage_type storage_;
// ------------------------------- constructors --------------------------------------------
indexmap_storage_pair() {}
indexmap_storage_pair(indexmap_type IM, storage_type ST) : indexmap_(std::move(IM)), storage_(std::move(ST)) {
#ifdef TRIQS_ARRAYS_CHECK_IM_STORAGE_COMPAT
if (ST.size() != IM.domain().number_of_elements())
TRIQS_RUNTIME_ERROR<<"index_storage_pair construction : storage and indices are not compatible";
#endif
}
/// The storage is allocated from the size of IM.
indexmap_storage_pair(const indexmap_type &IM) : indexmap_(IM), storage_() {
storage_ = StorageType(indexmap_.domain().number_of_elements());
}
public:
// Shallow copy
indexmap_storage_pair(indexmap_storage_pair const &X) = default;
indexmap_storage_pair(indexmap_storage_pair &&X) = default;
protected:
#ifdef TRIQS_WITH_PYTHON_SUPPORT
indexmap_storage_pair (PyObject * X, bool enforce_copy, const char * name ) {
try {
numpy_interface::numpy_extractor E(X, enforce_copy);
this->indexmap_ = E.indexmap(); this->storage_ = E.storage();
}
catch(numpy_interface::copy_exception s){// intercept only this one...
TRIQS_RUNTIME_ERROR<< " construction of a "<< name <<" from a numpy "
<<"\n T = "<< triqs::utility::typeid_name(value_type())
// lead to a nasty link pb ???
// linker search for IndexMapType::domain_type::rank in triqs.so
// and cannot resolve it ???
//<<"\n rank = "<< IndexMapType::domain_type::rank//this->rank
<<"\nfrom the python object \n"<< numpy_interface::object_to_string(X)
<<"\nThe error was :\n "<indexmap_,X.indexmap_); std::swap (this->storage_, X.storage_);
}
friend void swap( indexmap_storage_pair & A, indexmap_storage_pair & B) { A.swap_me(B);}
// ------------------------------- == --------------------------------------------
// at your own risk with floating value, but it is useful for int, string, etc....
// in particular for tests
friend bool operator==( indexmap_storage_pair const & A, indexmap_storage_pair const & B) {
if (A.shape() != B.shape()) return false;
auto ita = A.begin(); auto itb = B.begin();
for (;ita != A.end();++ita, ++itb) {if (!(*ita == *itb)) return false;}
return true;
}
friend bool operator!=( indexmap_storage_pair const & A, indexmap_storage_pair const & B) { return (!(A==B));}
public:
// ------------------------------- data access --------------------------------------------
indexmap_type const & indexmap() const {return indexmap_;}
storage_type const & storage() const {return storage_;}
storage_type & storage() {return storage_;}
#ifdef TRIQS_WITH_PYTHON_SUPPORT
PyObject *to_python() const {
if (is_empty()) TRIQS_RUNTIME_ERROR << "Error : trying to return an empty array/matrix/vector to python";
return numpy_interface::array_view_to_python(*this);
}
#endif
/// data_start is the starting point of the data of the object
/// this it NOT &storage()[0], which is the start of the underlying blokc
/// they are not equal for a view in general
value_type const * restrict data_start() const { return &storage_[indexmap_.start_shift()];}
value_type * restrict data_start() { return &storage_[indexmap_.start_shift()];}
using domain_type=typename indexmap_type::domain_type ;
domain_type const & domain() const { return indexmap_.domain();}
using shape_type=typename domain_type::index_value_type ;
shape_type const & shape() const { return domain().lengths();}
size_t shape(size_t i) const { return domain().lengths()[i];}
size_t num_elements() const { return domain().number_of_elements();}
//bool is_empty() const { return this->num_elements()==0;}
bool is_empty() const { return this->storage_.empty();}
// ------------------------------- operator () --------------------------------------------
using view_type = typename ISPViewType::type;
using weak_view_type = typename ISPViewType::type;
// Evaluation and slices
template
std14::enable_if_t<
(!clef::is_any_lazy::value) && (indexmaps::slicer::r_type::domain_type::rank==0) && (!IsConst)
, value_type &>
operator()(Args const & ... args) & { return storage_[indexmap_(args...)]; }
template
std14::enable_if_t<
(!clef::is_any_lazy::value) && (indexmaps::slicer::r_type::domain_type::rank==0)
, value_type const &>
operator()(Args const & ... args) const & { return storage_[indexmap_(args...)]; }
// && : return a & iif it is a non const view
template
std14::enable_if_t<
(!clef::is_any_lazy::value) && (indexmaps::slicer::r_type::domain_type::rank==0) && (!IsConst&&IsView)
, value_type &>
operator()(Args const & ... args) && {
// add here a security check in case it is a view, unique. For a regular type, move the result...
#ifdef TRIQS_ARRAYS_DEBUG
if (storage_.is_unique()) TRIQS_RUNTIME_ERROR <<"BUG : array : rvalue ref for an array...";
#endif
return storage_[indexmap_(args...)];
}
// && return a value if this is not a view (regular class) or it is a const_view
template
std14::enable_if_t<
(!clef::is_any_lazy::value) && (indexmaps::slicer::r_type::domain_type::rank==0) && (!(!IsConst&&IsView))
, value_type>
operator()(Args const & ... args) && { return storage_[indexmap_(args...)]; }
template struct result_of_call_as_view {
using IM2=typename indexmaps::slicer::r_type ;
//using V2=typename std::conditional::type, value_type>::type ;
using V2=value_type ;
static_assert(IM2::domain_type::rank !=0, "Internal error");
typedef typename ISPViewType < V2, IM2::domain_type::rank, typename IM2::traversal_order_in_template, ViewTag,
ForceBorrowed || StorageType::is_weak, is_const > ::type type;
};
template // non const version
typename boost::lazy_enable_if_c<
(!clef::is_any_lazy::value) && (indexmaps::slicer::r_type::domain_type::rank!=0) && (!IsConst)
, result_of_call_as_view
>::type // enable_if
operator()(Args const & ... args) & {
// simplify
return typename result_of_call_as_view::type ( indexmaps::slicer::invoke(indexmap_,args...), storage()); }
template // const version
typename boost::lazy_enable_if_c<
(!clef::is_any_lazy::value) && (indexmaps::slicer::r_type::domain_type::rank!=0)
, result_of_call_as_view
>::type // enable_if
operator()(Args const & ... args) const & {
return typename result_of_call_as_view::type ( indexmaps::slicer::invoke(indexmap_,args...), storage()); }
template // rvalue version : same value of weak as this
typename boost::lazy_enable_if_c<
(!clef::is_any_lazy::value) && (indexmaps::slicer::r_type::domain_type::rank!=0)
, result_of_call_as_view
>::type // enable_if
operator()(Args const & ... args) && {
//std::cerr << "slicing a temporary"<< this->storage().is_weak<< result_of_call_as_view::type::storage_type::is_weak << std::endl;
return typename result_of_call_as_view::type ( indexmaps::slicer::invoke(indexmap_,args...), std::move(storage()));
}
/// Equivalent to make_view
typename ISPViewType::type
operator()() const & { return *this; }
typename ISPViewType::type
operator()() & { return *this; }
typename ISPViewType::type
operator()() && { return *this; }
// Interaction with the CLEF library : calling with any clef expression as argument build a new clef expression
// NB : this is ok because indexmap_storage_pair has a shallow copy constructor ...
// Correction : no copy, just a ref...
// so A(i_) if A is an array will NOT copy the data....
template< typename... Args>
typename clef::_result_of::make_expr_call::type
operator()( Args&&... args ) const & {
static_assert(sizeof...(Args) <= indexmap_type::rank, "Incorrect number of variable in call");// not perfect : ellipsis ...
return make_expr_call(*this,std::forward(args)...);
}
template< typename... Args>
typename clef::_result_of::make_expr_call::type
operator()( Args&&... args ) & {
static_assert(sizeof...(Args) <= indexmap_type::rank, "Incorrect number of variable in call");// not perfect : ellipsis ...
return make_expr_call(*this,std::forward(args)...);
}
template< typename... Args>
typename clef::_result_of::make_expr_call::type
operator()( Args&&... args ) && {
static_assert(sizeof...(Args) <= indexmap_type::rank, "Incorrect number of variable in call");// not perfect : ellipsis ...
return make_expr_call(std::move(*this),std::forward(args)...);
}
// ------------------------------- clef auto assign --------------------------------------------
// For simple cases, it is assign_foreach. But when f(args...) is a function from a clef expression
// we make a chain call like cf clef vector adapter
template struct _worker {
indexmap_storage_pair &A;
Function const &f;
template void assign(T &x, RHS &&rhs) { x = std::forward(rhs); }
template void assign(T &x, clef::make_fun_impl &&rhs) {
triqs_clef_auto_assign(x, std::forward>(rhs));
}
template void operator()(Args const &... args) { this->assign(A(args...), f(args...)); }
};
template friend void triqs_clef_auto_assign(indexmap_storage_pair &x, Fnt f) {
foreach(x, _worker{x, f});
}
// template friend void triqs_clef_auto_assign (indexmap_storage_pair & x, Fnt f) { assign_foreach(x,f);}
// ------------------------------- Iterators --------------------------------------------
using const_iterator=iterator_adapter ;
using iterator=iterator_adapter ;
const_iterator begin() const {return const_iterator(indexmap(),storage(),false);}
const_iterator end() const {return const_iterator(indexmap(),storage(),true);}
const_iterator cbegin() const {return const_iterator(indexmap(),storage(),false);}
const_iterator cend() const {return const_iterator(indexmap(),storage(),true);}
iterator begin() {return iterator(indexmap(),storage(),false);}
iterator end() {return iterator(indexmap(),storage(),true);}
protected:
// ------------------------------- resize --------------------------------------------
//
void resize (domain_type const & d) {
indexmap_ = IndexMapType(d,indexmap_.get_memory_layout());
// build a new one with the lengths of IND BUT THE SAME layout !
// optimisation. Construct a storage only if the new index is not compatible (size mismatch).
if (storage_.size() != indexmap_.domain().number_of_elements())
storage_ = StorageType(indexmap_.domain().number_of_elements());
}
template
void resize_and_clone_data( Xtype const & X) { indexmap_ = X.indexmap(); storage_ = X.storage().clone(); }
// BOOST Serialization
friend class boost::serialization::access;
template
void serialize(Archive & ar, const unsigned int version) {
ar & TRIQS_MAKE_NVP("storage",this->storage_);
ar & TRIQS_MAKE_NVP("indexmap",this->indexmap_);
}
// pretty print of the array
friend std::ostream & operator << (std::ostream & out, const indexmap_storage_pair & A) {
if (A.storage().size()==0) out<<"empty ";
else pretty_print(out,A);
return out;
}
};// end class
}}//namespace triqs::arrays