1
0
mirror of https://github.com/paboyle/Grid.git synced 2025-04-11 14:40:46 +01:00

599 lines
22 KiB
C++

/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./lib/serialisation/BaseIO.h
Copyright (C) 2015
Author: Antonin Portelli <antonin.portelli@me.com>
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
Author: Guido Cossu <guido.cossu@ed.ac.uk>
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 2 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, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#ifndef GRID_SERIALISATION_ABSTRACT_READER_H
#define GRID_SERIALISATION_ABSTRACT_READER_H
#include <type_traits>
#include <Grid/tensors/Tensors.h>
#include <Grid/serialisation/VectorUtils.h>
#include <Grid/Eigen/unsupported/CXX11/Tensor>
namespace Grid {
// Abstract writer/reader classes ////////////////////////////////////////////
// static polymorphism implemented using CRTP idiom
class Serializable;
// which types are supported scalar types for Eigen::Tensor
template<typename T> struct is_eigen_tensor_scalar : std::integral_constant<bool,
std::is_arithmetic<T>::value || is_complex<T>::value> {};
// which types are grid tensors
template <typename T> struct is_grid_tensor : public std::false_type {
static constexpr unsigned int rank = 0;
static constexpr unsigned int dim = 1;
};
template <typename T> struct is_grid_tensor<iScalar<T>> : public std::true_type {};
template <typename T, int N> struct is_grid_tensor<iVector<T, N>> : public std::true_type {};
template <typename T, int N> struct is_grid_tensor<iMatrix<T, N>> : public std::true_type {};
// Rank and dimension of grid tensors, i.e. compositions of iScalar, iVector and iMatrix
template <typename T> struct grid_tensor_att {
static constexpr unsigned int depth = 0; // How many levels of Grid Tensor there are
static constexpr unsigned int rank = 0; // The rank of the grid tensor (i.e. how many indices used)
static constexpr unsigned int rank_non_trivial = 0; // As per rank, but excludes those of dimension 1
static constexpr unsigned int count = 1; // total number of elements (i.e. product of dimensions)
typedef T scalar_type; // Type of the underlying scalar
static constexpr size_t scalar_size = sizeof(T); // Size of the underlying scalar in bytes
static constexpr size_t size = scalar_size * count; // total size of elements in bytes
// e.g. iScalar<iVector<Complex,1>>
// depth = 2
// rank = 1
// rank_non_trivial = 0
// count = 1
// e.g. iVector<iMatrix<Complex,3>,4>
// depth = 2
// rank = 3
// rank_non_trivial = 3
// count = 36
// e.g. iScalar<iVector<iMatrix<Complex,4>,3>>
// depth = 3
// rank = 3
// rank_non_trivial = 3
// count = 48
};
template <typename T> struct grid_tensor_att<iScalar<T>> {
static constexpr unsigned int depth = 1 + grid_tensor_att<T>::depth;
static constexpr unsigned int rank = 0 + grid_tensor_att<T>::rank;
static constexpr unsigned int rank_non_trivial = 0 + grid_tensor_att<T>::rank_non_trivial;
static constexpr unsigned int count = 1 * grid_tensor_att<T>::count;
typedef typename grid_tensor_att<T>::scalar_type scalar_type;
static constexpr size_t scalar_size = grid_tensor_att<T>::scalar_size;
static constexpr size_t size = scalar_size * count;
};
template <typename T, int N> struct grid_tensor_att<iVector<T, N>> {
static constexpr unsigned int depth = 1 + grid_tensor_att<T>::depth;
static constexpr unsigned int rank = 1 + grid_tensor_att<T>::rank;
static constexpr unsigned int rank_non_trivial = (N>1 ? 1 : 0) + grid_tensor_att<T>::rank_non_trivial;
static constexpr unsigned int count = N * grid_tensor_att<T>::count;
typedef typename grid_tensor_att<T>::scalar_type scalar_type;
static constexpr size_t scalar_size = grid_tensor_att<T>::scalar_size;
static constexpr size_t size = scalar_size * count;
};
template <typename T, int N> struct grid_tensor_att<iMatrix<T, N>> {
static constexpr unsigned int depth = 1 + grid_tensor_att<T>::depth;
static constexpr unsigned int rank = 2 + grid_tensor_att<T>::rank;
static constexpr unsigned int rank_non_trivial = (N>1 ? 2 : 0) + grid_tensor_att<T>::rank_non_trivial;
static constexpr unsigned int count = N * N * grid_tensor_att<T>::count;
typedef typename grid_tensor_att<T>::scalar_type scalar_type;
static constexpr size_t scalar_size = grid_tensor_att<T>::scalar_size;
static constexpr size_t size = scalar_size * count;
};
// Static abstract writer
template <typename T>
class Writer
{
public:
Writer(void);
virtual ~Writer(void) = default;
void push(const std::string &s);
void pop(void);
template <typename U>
typename std::enable_if<std::is_base_of<Serializable, U>::value, void>::type
write(const std::string& s, const U &output);
template <typename U>
typename std::enable_if<!std::is_base_of<Serializable, U>::value
&& !std::is_base_of<Eigen::TensorBase<U, Eigen::ReadOnlyAccessors>, U>::value, void>::type
write(const std::string& s, const U &output);
template <typename U>
void write(const std::string &s, const iScalar<U> &output);
template <typename U, int N>
void write(const std::string &s, const iVector<U, N> &output);
template <typename U, int N>
void write(const std::string &s, const iMatrix<U, N> &output);
template <typename ETensor>
typename std::enable_if<std::is_base_of<Eigen::TensorBase<ETensor, Eigen::ReadOnlyAccessors>, ETensor>::value && (std::is_arithmetic<typename ETensor::Scalar>::value || Grid::is_complex<typename ETensor::Scalar>::value), void>::type
write(const std::string &s, const ETensor &output);
template <typename ETensor/*, typename U, int N*/>
typename std::enable_if<std::is_base_of<Eigen::TensorBase<ETensor, Eigen::ReadOnlyAccessors>, ETensor>::value
&& is_grid_tensor<typename ETensor::Scalar>::value
//&& !(std::is_arithmetic<typename ETensor::Scalar>::value || Grid::is_complex<typename ETensor::Scalar>::value)
/*&& ( std::is_base_of<typename ETensor::Scalar, iScalar<U> >::value
|| std::is_base_of<typename ETensor::Scalar, iVector<U, N>>::value
|| std::is_base_of<typename ETensor::Scalar, iMatrix<U, N>>::value )*/, void>::type
write(const std::string &s, const ETensor &output);
/*template <typename ETensor, typename U, int N>
typename std::enable_if<std::is_base_of<Eigen::TensorBase<ETensor, Eigen::ReadOnlyAccessors>, ETensor>::value
&& std::is_base_of<typename ETensor::Scalar, iVector<U, N>>::value, void>::type
write(const std::string &s, const ETensor &output);
template <typename ETensor, typename U, int N>
typename std::enable_if<std::is_base_of<Eigen::TensorBase<ETensor, Eigen::ReadOnlyAccessors>, ETensor>::value
&& std::is_base_of<typename ETensor::Scalar, iMatrix<U, N>>::value, void>::type
write(const std::string &s, const ETensor &output);*/
void scientificFormat(const bool set);
bool isScientific(void);
void setPrecision(const unsigned int prec);
unsigned int getPrecision(void);
private:
T *upcast;
bool scientific_{false};
unsigned int prec_{0};
};
// Static abstract reader
template <typename T>
class Reader
{
public:
Reader(void);
virtual ~Reader(void) = default;
bool push(const std::string &s);
void pop(void);
template <typename U>
typename std::enable_if<std::is_base_of<Serializable, U>::value, void>::type
read(const std::string& s, U &output);
template <typename U>
typename std::enable_if<!std::is_base_of<Serializable, U>::value, void>::type
read(const std::string& s, U &output);
template <typename U>
void read(const std::string &s, iScalar<U> &output);
template <typename U, int N>
void read(const std::string &s, iVector<U, N> &output);
template <typename U, int N>
void read(const std::string &s, iMatrix<U, N> &output);
protected:
template <typename U>
void fromString(U &output, const std::string &s);
private:
T *upcast;
};
// What is the vtype
template<typename T> struct isReader {
static const bool value = false;
};
template<typename T> struct isWriter {
static const bool value = false;
};
// Writer template implementation
template <typename T>
Writer<T>::Writer(void)
{
upcast = static_cast<T *>(this);
}
template <typename T>
void Writer<T>::push(const std::string &s)
{
upcast->push(s);
}
template <typename T>
void Writer<T>::pop(void)
{
upcast->pop();
}
template <typename T>
template <typename U>
typename std::enable_if<std::is_base_of<Serializable, U>::value, void>::type
Writer<T>::write(const std::string &s, const U &output)
{
U::write(*this, s, output);
}
template <typename T>
template <typename U>
typename std::enable_if<!std::is_base_of<Serializable, U>::value
&& !std::is_base_of<Eigen::TensorBase<U, Eigen::ReadOnlyAccessors>, U>::value, void>::type
Writer<T>::write(const std::string &s, const U &output)
{
upcast->writeDefault(s, output);
}
template <typename T>
template <typename U>
void Writer<T>::write(const std::string &s, const iScalar<U> &output)
{
upcast->writeDefault(s, tensorToVec(output));
}
template <typename T>
template <typename U, int N>
void Writer<T>::write(const std::string &s, const iVector<U, N> &output)
{
upcast->writeDefault(s, tensorToVec(output));
}
template <typename T>
template <typename U, int N>
void Writer<T>::write(const std::string &s, const iMatrix<U, N> &output)
{
upcast->writeDefault(s, tensorToVec(output));
}
// Eigen::Tensors of arithmetic/complex base type
template <typename T>
template <typename ETensor>
typename std::enable_if<std::is_base_of<Eigen::TensorBase<ETensor, Eigen::ReadOnlyAccessors>, ETensor>::value && (std::is_arithmetic<typename ETensor::Scalar>::value || Grid::is_complex<typename ETensor::Scalar>::value), void>::type
Writer<T>::write(const std::string &s, const ETensor &output)
{
std::cout << "Eigen::Tensors of arithmetic/complex base type" << std::endl;
const typename ETensor::Index NumElements{output.size()};
assert( NumElements > 0 );
if( NumElements == 1 )
upcast->writeDefault(s, * output.data());
else {
// We're not interested in trivial dimensions, i.e. dimensions = 1
unsigned int TrivialDimCount{0};
std::vector<size_t> ReducedDims;
for(auto i = 0; i < output.NumDimensions; i++ ) {
auto dim = output.dimension(i);
if( dim <= 1 ) {
TrivialDimCount++;
assert( dim == 1 ); // Not expecting dimension to be <= 0
} else {
size_t s = static_cast<size_t>(dim);
assert( s == dim ); // check we didn't lose anything in the conversion
ReducedDims.push_back(s);
}
}
assert( output.NumDimensions > TrivialDimCount > 0 ); // NB: NumElements > 1 implies this is not a scalar, so some dims should be left
// Create a single, flat vector to hold all the data
std::vector<typename ETensor::Scalar> flat(NumElements);
// Now copy all the data to my flat vector
// Regardless of the Eigen::Tensor storage order, the copy will be Row Major
std::array<typename ETensor::Index, ETensor::NumIndices> MyIndex;
for( int i = 0 ; i < output.NumDimensions ; i++ ) MyIndex[i] = 0;
for( typename ETensor::Index n = 0; n < NumElements; n++ ) {
flat[n] = output( MyIndex );
// Now increment the index
for( int i = output.NumDimensions - 1; i >= 0 && ++MyIndex[i] == output.dimension(i); i-- )
MyIndex[i] = 0;
}
upcast->template writeMultiDim<typename ETensor::Scalar>(s, ReducedDims, flat);
}
}
// Eigen::Tensors of iScalar<U>
template <typename T>
template <typename ETensor/*, typename U, int N*/>
typename std::enable_if<std::is_base_of<Eigen::TensorBase<ETensor, Eigen::ReadOnlyAccessors>, ETensor>::value
&& is_grid_tensor<typename ETensor::Scalar>::value
//&& !(std::is_arithmetic<typename ETensor::Scalar>::value || Grid::is_complex<typename ETensor::Scalar>::value)
/*&& ( std::is_base_of<typename ETensor::Scalar, iScalar<U> >::value
|| std::is_base_of<typename ETensor::Scalar, iVector<U, N>>::value
|| std::is_base_of<typename ETensor::Scalar, iMatrix<U, N>>::value )*/, void>::type
Writer<T>::write(const std::string &s, const ETensor &output)
{
std::cout << "Eigen::Tensors of iScalar<U>" << std::endl;
const typename ETensor::Index NumElements{output.size()};
assert( NumElements > 0 );
if( NumElements == 1 )
upcast->writeDefault(s, tensorToVec(* output.data()));
else {
// We're not interested in trivial dimensions, i.e. dimensions = 1
unsigned int TrivialDimCount{0};
std::vector<size_t> ReducedDims;
for(auto i = 0; i < output.NumDimensions; i++ ) {
auto dim = output.dimension(i);
if( dim <= 1 ) {
TrivialDimCount++;
assert( dim == 1 ); // Not expecting dimension to be <= 0
} else {
size_t s = static_cast<size_t>(dim);
assert( s == dim ); // check we didn't lose anything in the conversion
ReducedDims.push_back(s);
}
}
assert( output.NumDimensions > TrivialDimCount > 0 ); // NB: NumElements > 1 implies this is not a scalar, so some dims should be left
// Now add the extra dimensions, based on object zero
typename TensorToVec<typename ETensor::Scalar>::type ttv = tensorToVec(* output.data());
Flatten<typename TensorToVec<typename ETensor::Scalar>::type> f(ttv);
const std::vector<size_t> & ExtraDims{f.getDim()};
size_t ExtraCount{1};
for( auto i : ExtraDims ) {
assert( i > 0 );
ExtraCount *= i;
ReducedDims.push_back(i);
}
typedef typename ETensor::Scalar::scalar_type Scalar;
assert( sizeof( typename ETensor::Scalar ) == ExtraCount * sizeof( Scalar ) );
// Create a single, flat vector to hold all the data
const typename ETensor::Index TotalNumElements = NumElements * ExtraCount;
std::vector<Scalar> flat(TotalNumElements);
// Now copy all the data to my flat vector
// Regardless of the Eigen::Tensor storage order, the copy will be Row Major
std::array<typename ETensor::Index, ETensor::NumIndices> MyIndex;
for( int i = 0 ; i < output.NumDimensions ; i++ ) MyIndex[i] = 0;
for( typename ETensor::Index n = 0; n < TotalNumElements; ) {
const Scalar * p = reinterpret_cast<const Scalar *>( &output( MyIndex ));
for( auto j = 0; j < ExtraCount ; j++ )
flat[n++] = * p++;
// Now increment the index
for( int i = output.NumDimensions - 1; i >= 0 && ++MyIndex[i] == output.dimension(i); i-- )
MyIndex[i] = 0;
}
upcast->template writeMultiDim<Scalar>(s, ReducedDims, flat);
}
}
// Eigen::Tensors of iVector<U, N>
/*template <typename T>
template <typename ETensor, typename U, int N>
typename std::enable_if<std::is_base_of<Eigen::TensorBase<ETensor, Eigen::ReadOnlyAccessors>, ETensor>::value
&& std::is_base_of<typename ETensor::Scalar, iVector<U, N>>::value, void>::type
Writer<T>::write(const std::string &s, const ETensor &output)
{
//upcast->writeDefault(s, tensorToVec(output));
std::cout << "I really should add code to write Eigen::Tensor (iVector) ..." << std::endl;
}
// Eigen::Tensors of iMatrix<U, N>
template <typename T>
template <typename ETensor, typename U, int N>
typename std::enable_if<std::is_base_of<Eigen::TensorBase<ETensor, Eigen::ReadOnlyAccessors>, ETensor>::value
&& std::is_base_of<typename ETensor::Scalar, iMatrix<U, N>>::value, void>::type
Writer<T>::write(const std::string &s, const ETensor &output)
{
//upcast->writeDefault(s, tensorToVec(output));
std::cout << "I really should add code to write Eigen::Tensor (iMatrix) ..." << std::endl;
}*/
template <typename T>
void Writer<T>::scientificFormat(const bool set)
{
scientific_ = set;
}
template <typename T>
bool Writer<T>::isScientific(void)
{
return scientific_;
}
template <typename T>
void Writer<T>::setPrecision(const unsigned int prec)
{
prec_ = prec;
}
template <typename T>
unsigned int Writer<T>::getPrecision(void)
{
return prec_;
}
// Reader template implementation
template <typename T>
Reader<T>::Reader(void)
{
upcast = static_cast<T *>(this);
}
template <typename T>
bool Reader<T>::push(const std::string &s)
{
return upcast->push(s);
}
template <typename T>
void Reader<T>::pop(void)
{
upcast->pop();
}
template <typename T>
template <typename U>
typename std::enable_if<std::is_base_of<Serializable, U>::value, void>::type
Reader<T>::read(const std::string &s, U &output)
{
U::read(*this, s, output);
}
template <typename T>
template <typename U>
typename std::enable_if<!std::is_base_of<Serializable, U>::value, void>::type
Reader<T>::read(const std::string &s, U &output)
{
upcast->readDefault(s, output);
}
template <typename T>
template <typename U>
void Reader<T>::read(const std::string &s, iScalar<U> &output)
{
typename TensorToVec<iScalar<U>>::type v;
upcast->readDefault(s, v);
vecToTensor(output, v);
}
template <typename T>
template <typename U, int N>
void Reader<T>::read(const std::string &s, iVector<U, N> &output)
{
typename TensorToVec<iVector<U, N>>::type v;
upcast->readDefault(s, v);
vecToTensor(output, v);
}
template <typename T>
template <typename U, int N>
void Reader<T>::read(const std::string &s, iMatrix<U, N> &output)
{
typename TensorToVec<iMatrix<U, N>>::type v;
upcast->readDefault(s, v);
vecToTensor(output, v);
}
template <typename T>
template <typename U>
void Reader<T>::fromString(U &output, const std::string &s)
{
std::istringstream is(s);
is.exceptions(std::ios::failbit);
try
{
is >> std::boolalpha >> output;
}
catch(std::ios_base::failure &e)
{
std::cerr << "numerical conversion failure on '" << s << "' ";
std::cerr << "(typeid: " << typeid(U).name() << ")" << std::endl;
abort();
}
}
// serializable base class ///////////////////////////////////////////////////
class Serializable
{
public:
template <typename T>
static inline void write(Writer<T> &WR,const std::string &s,
const Serializable &obj)
{}
template <typename T>
static inline void read(Reader<T> &RD,const std::string &s,
Serializable &obj)
{}
friend inline std::ostream & operator<<(std::ostream &os,
const Serializable &obj)
{
return os;
}
template <typename T>
static inline typename std::enable_if<!std::is_base_of<Eigen::TensorBase<T, Eigen::ReadOnlyAccessors>, T>::value, bool>::type
CompareMember(const T &lhs, const T &rhs) {
return lhs == rhs;
}
template <typename T>
static inline typename std::enable_if<std::is_base_of<Eigen::TensorBase<T, Eigen::ReadOnlyAccessors>, T>::value, bool>::type
CompareMember(const T &lhs, const T &rhs) {
Eigen::Tensor<bool, 0> bResult = (lhs == rhs).all();
return bResult(0);
}
template <typename T>
static inline typename std::enable_if<std::is_base_of<Eigen::TensorBase<T, Eigen::ReadOnlyAccessors>, T>::value, bool>::type
CompareMember(const std::vector<T> &lhs, const std::vector<T> &rhs) {
const auto NumElements{lhs.size()};
bool bResult = ( NumElements == rhs.size() );
for( auto i = 0 ; i < NumElements && bResult ; i++ ) {
Eigen::Tensor<bool, 0> b = (lhs[i] == rhs[i]).all();
bResult = b(0);
}
return bResult;
}
template <typename T>
static inline typename std::enable_if<!std::is_base_of<Eigen::TensorBase<T, Eigen::ReadOnlyAccessors>, T>::value, void>::type
WriteMember(std::ostream &os, const T &object) {
os << object;
}
template <typename T>
static inline typename std::enable_if<std::is_base_of<Eigen::TensorBase<T, Eigen::ReadOnlyAccessors>, T>::value, void>::type
WriteMember(std::ostream &os, const T &object) {
os << "Eigen::Tensor";
}
};
// Generic writer interface //////////////////////////////////////////////////
template <typename T>
inline void push(Writer<T> &w, const std::string &s) {
w.push(s);
}
template <typename T>
inline void push(Writer<T> &w, const char *s)
{
w.push(std::string(s));
}
template <typename T>
inline void pop(Writer<T> &w)
{
w.pop();
}
template <typename T, typename U>
inline void write(Writer<T> &w, const std::string& s, const U &output)
{
w.write(s, output);
}
// Generic reader interface //////////////////////////////////////////////////
template <typename T>
inline bool push(Reader<T> &r, const std::string &s)
{
return r.push(s);
}
template <typename T>
inline bool push(Reader<T> &r, const char *s)
{
return r.push(std::string(s));
}
template <typename T>
inline void pop(Reader<T> &r)
{
r.pop();
}
template <typename T, typename U>
inline void read(Reader<T> &r, const std::string &s, U &output)
{
r.read(s, output);
}
}
#endif