DataVector.h
Go to the documentation of this file.
/*****************************************************************************
*
* Copyright (c) 2003-2020 by The University of Queensland
* http://www.uq.edu.au
*
* Primary Business: Queensland, Australia
* Licensed under the Apache License, version 2.0
* http://www.apache.org/licenses/LICENSE-2.0
*
* Development until 2012 by Earth Systems Science Computational Center (ESSCC)
* Development 2012-2013 by School of Earth Sciences
* Development from 2014-2017 by Centre for Geoscience Computing (GeoComp)
* Development from 2019 by School of Earth and Environmental Sciences
**
*****************************************************************************/
#ifndef __ESCRIPT_DATAVECTOR_H__
#define __ESCRIPT_DATAVECTOR_H__
#include "system_dep.h"
#include "DataTypes.h"
#include "Assert.h"
#include "DataVectorAlt.h"
#include "DataVectorTaipan.h"
#ifdef ESYS_HAVE_BOOST_NUMPY
#include <boost/python.hpp>
#include <boost/python/numpy.hpp>
#endif
// ensure that nobody else tries to instantiate the complex version
extern template class escript::DataTypes::DataVectorAlt<escript::DataTypes::cplx_t>;
namespace escript {
// Functions in DataTypes:: which manipulate DataVectors
namespace DataTypes
{
// This is the main version we had
//typedef DataVectorTaipan DataVector;
typedef escript::DataTypes::DataVectorAlt<real_t> RealVectorType;
typedef escript::DataTypes::DataVectorAlt<cplx_t> CplxVectorType;
void
pointToStream(std::ostream& os, const RealVectorType::ElementType* data,const ShapeType& shape, int offset, bool needsep=true, const std::string& sep=",");
void
pointToStream(std::ostream& os, const CplxVectorType::ElementType* data,const ShapeType& shape, int offset, bool needsep=true, const std::string& sep=",");
#ifdef ESYS_HAVE_BOOST_NUMPY
void
pointToNumpyArrayOld(boost::python::numpy::ndarray& dataArray, const RealVectorType::ElementType* data, const ShapeType& shape, long offset, long numsamples, long dpps, long numdata);
void
pointToNumpyArrayOld(boost::python::numpy::ndarray& dataArray, const CplxVectorType::ElementType* data, const ShapeType& shape, long offset, long numsamples, long dpps, long numdata);
void
pointToNumpyArray(boost::python::numpy::ndarray& dataArray, const RealVectorType::ElementType* data, const ShapeType& shape, int offset, int d, int index);
void
pointToNumpyArray(boost::python::numpy::ndarray& dataArray, const CplxVectorType::ElementType* data, const ShapeType& shape, int offset, int d, int index);
#endif
std::string
pointToString(const RealVectorType& data,const ShapeType& shape, int offset, const std::string& prefix);
std::string
pointToString(const CplxVectorType& data,const ShapeType& shape, int offset, const std::string& prefix);
void copyPoint(RealVectorType& dest, vec_size_type doffset, vec_size_type nvals, const RealVectorType& src, vec_size_type soffset);
void copyPoint(CplxVectorType& dest, vec_size_type doffset, vec_size_type nvals, const CplxVectorType& src, vec_size_type soffset);
void fillComplexFromReal(const RealVectorType& r, CplxVectorType& c);
template <class VEC>
ESCRIPT_DLL_API
void
copySlice(VEC& left,
const ShapeType& leftShape,
typename VEC::size_type leftOffset,
const VEC& other,
const ShapeType& otherShape,
typename VEC::size_type otherOffset,
const RegionLoopRangeType& region)
{
//
// Make sure vectors are not empty
ESYS_ASSERT(!left.size()==!1, "left data is empty."); //Note: !1=0, but clang returns an error if the rhs is 0 here
ESYS_ASSERT(!other.size()==!1, "other data is empty.");
//
// Check the vector to be sliced from is compatible with the region to be sliced,
// and that the region to be sliced is compatible with this vector:
ESYS_ASSERT(checkOffset(leftOffset,left.size(),noValues(leftShape)),
"offset incompatible with this vector.");
ESYS_ASSERT(otherOffset+noValues(leftShape)<=other.size(),
"offset incompatible with other vector.");
ESYS_ASSERT(getRank(otherShape)==region.size(),
"slice not same rank as vector to be sliced from.");
ESYS_ASSERT(noValues(leftShape)==noValues(getResultSliceShape(region)),
"slice shape not compatible shape for this vector.");
//
// copy the values in the specified region of the other vector into this vector
// the following loops cannot be parallelised due to the numCopy counter
int numCopy=0;
switch (region.size()) {
case 0:
/* this case should never be encountered,
as python will never pass us an empty region.
here for completeness only, allows slicing of a scalar */
// (*m_data)[leftOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex()];
left[leftOffset+numCopy]=other[otherOffset];
numCopy++;
break;
case 1:
for (int i=region[0].first;i<region[0].second;i++) {
left[leftOffset+numCopy]=other[otherOffset+getRelIndex(otherShape,i)];
numCopy++;
}
break;
case 2:
for (int j=region[1].first;j<region[1].second;j++) {
for (int i=region[0].first;i<region[0].second;i++) {
/* (*m_data)[leftOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i,j)];*/
left[leftOffset+numCopy]=other[otherOffset+getRelIndex(otherShape,i,j)];
numCopy++;
}
}
break;
case 3:
for (int k=region[2].first;k<region[2].second;k++) {
for (int j=region[1].first;j<region[1].second;j++) {
for (int i=region[0].first;i<region[0].second;i++) {
// (*m_data)[leftOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i,j,k)];
left[leftOffset+numCopy]=other[otherOffset+getRelIndex(otherShape,i,j,k)];
numCopy++;
}
}
}
break;
case 4:
for (int l=region[3].first;l<region[3].second;l++) {
for (int k=region[2].first;k<region[2].second;k++) {
for (int j=region[1].first;j<region[1].second;j++) {
for (int i=region[0].first;i<region[0].second;i++) {
/* (*m_data)[leftOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i,j,k,l)];*/
left[leftOffset+numCopy]=other[otherOffset+getRelIndex(otherShape,i,j,k,l)];
numCopy++;
}
}
}
}
break;
default:
std::stringstream mess;
mess << "Error - (copySlice) Invalid slice region rank: " << region.size();
throw DataException(mess.str());
}
}
template<typename VEC>
ESCRIPT_DLL_API
void
copySliceFrom(VEC& left,
const ShapeType& leftShape,
typename VEC::size_type leftOffset,
const VEC& other,
const ShapeType& otherShape,
typename VEC::size_type otherOffset,
const RegionLoopRangeType& region)
{
//
// Make sure vectors are not empty
ESYS_ASSERT(left.size()!=0, "this vector is empty.");
ESYS_ASSERT(other.size()!=0, "other vector is empty.");
//
// Check this vector is compatible with the region to be sliced,
// and that the region to be sliced is compatible with the other vector:
ESYS_ASSERT(checkOffset(otherOffset,other.size(),noValues(otherShape)),
"offset incompatible with other vector.");
ESYS_ASSERT(leftOffset+noValues(otherShape)<=left.size(),
"offset incompatible with this vector.");
ESYS_ASSERT(getRank(leftShape)==region.size(),
"slice not same rank as this vector.");
ESYS_ASSERT(getRank(otherShape)==0 || noValues(otherShape)==noValues(getResultSliceShape(region)),
"slice shape not compatible shape for other vector.");
//
// copy the values in the other vector into the specified region of this vector
// allow for case where other vector is a scalar
if (getRank(otherShape)==0) {
// the following loops cannot be parallelised due to the numCopy counter
int numCopy=0;
switch (region.size()) {
case 0:
/* this case should never be encountered,
as python will never pass us an empty region.
here for completeness only, allows slicing of a scalar */
//(*m_data)[leftOffset+relIndex()]=(*other.m_data)[otherOffset];
left[leftOffset]=other[otherOffset];
numCopy++;
break;
case 1:
for (int i=region[0].first;i<region[0].second;i++) {
left[leftOffset+getRelIndex(leftShape,i)]=other[otherOffset];
numCopy++;
}
break;
case 2:
for (int j=region[1].first;j<region[1].second;j++) {
for (int i=region[0].first;i<region[0].second;i++) {
left[leftOffset+getRelIndex(leftShape,i,j)]=other[otherOffset];
numCopy++;
}
}
break;
case 3:
for (int k=region[2].first;k<region[2].second;k++) {
for (int j=region[1].first;j<region[1].second;j++) {
for (int i=region[0].first;i<region[0].second;i++) {
left[leftOffset+getRelIndex(leftShape,i,j,k)]=other[otherOffset];
numCopy++;
}
}
}
break;
case 4:
for (int l=region[3].first;l<region[3].second;l++) {
for (int k=region[2].first;k<region[2].second;k++) {
for (int j=region[1].first;j<region[1].second;j++) {
for (int i=region[0].first;i<region[0].second;i++) {
left[leftOffset+getRelIndex(leftShape,i,j,k,l)]=other[otherOffset];
numCopy++;
}
}
}
}
break;
default:
std::stringstream mess;
mess << "Error - (copySliceFrom) Invalid slice region rank: " << region.size();
throw DataException(mess.str());
}
} else {
// the following loops cannot be parallelised due to the numCopy counter
int numCopy=0;
switch (region.size()) {
case 0:
/* this case should never be encountered,
as python will never pass us an empty region.
here for completeness only, allows slicing of a scalar */
//(*m_data)[leftOffset+relIndex()]=(*other.m_data)[otherOffset+numCopy];
left[leftOffset]=other[otherOffset+numCopy];
numCopy++;
break;
case 1:
for (int i=region[0].first;i<region[0].second;i++) {
left[leftOffset+getRelIndex(leftShape,i)]=other[otherOffset+numCopy];
numCopy++;
}
break;
case 2:
for (int j=region[1].first;j<region[1].second;j++) {
for (int i=region[0].first;i<region[0].second;i++) {
left[leftOffset+getRelIndex(leftShape,i,j)]=other[otherOffset+numCopy];
numCopy++;
}
}
break;
case 3:
for (int k=region[2].first;k<region[2].second;k++) {
for (int j=region[1].first;j<region[1].second;j++) {
for (int i=region[0].first;i<region[0].second;i++) {
left[leftOffset+getRelIndex(leftShape,i,j,k)]=other[otherOffset+numCopy];
numCopy++;
}
}
}
break;
case 4:
for (int l=region[3].first;l<region[3].second;l++) {
for (int k=region[2].first;k<region[2].second;k++) {
for (int j=region[1].first;j<region[1].second;j++) {
for (int i=region[0].first;i<region[0].second;i++) {
left[leftOffset+getRelIndex(leftShape,i,j,k,l)]=other[otherOffset+numCopy];
numCopy++;
}
}
}
}
break;
default:
std::stringstream mess;
mess << "Error - (copySliceFrom) Invalid slice region rank: " << region.size();
throw DataException(mess.str());
}
}
}
}
} // end of namespace
#endif // __ESCRIPT_DATAVECTOR_H__
ESCRIPT_DLL_API
#define ESCRIPT_DLL_API
Definition: escriptcore/src/system_dep.h:30
escript::DataTypes::real_t
double real_t
type of all real-valued scalars in escript
Definition: DataTypes.h:79
escript::DataTypes::DataVectorAlt::ElementType
T ElementType
Definition: DataVectorAlt.h:85
escript::DataTypes::copySlice
void copySlice(VEC &left, const ShapeType &leftShape, typename VEC::size_type leftOffset, const VEC &other, const ShapeType &otherShape, typename VEC::size_type otherOffset, const RegionLoopRangeType &region)
Copy a data slice specified by the given region and offset from the "other" vector into the "left" ve...
Definition: DataVector.h:174
escript::DataTypes::getRank
int getRank(const DataTypes::ShapeType &shape)
Return the rank (number of dimensions) of the given shape.
Definition: DataTypes.h:247
escript::DataTypes::checkOffset
bool checkOffset(vec_size_type offset, int size, int noval)
Definition: DataTypes.h:352
escript::DataTypes::DataVectorAlt
Definition: DataVectorAlt.h:65
Assert.h
system_dep.h
escript::DataTypes::vec_size_type
long vec_size_type
Definition: DataTypes.h:76
escript::DataException
Definition: DataException.h:41
escript::DataTypes::ShapeType
std::vector< int > ShapeType
The shape of a single datapoint.
Definition: DataTypes.h:71
escript::DataTypes::fillComplexFromReal
void fillComplexFromReal(const RealVectorType &r, CplxVectorType &c)
copy data from a real vector to a complex vector The complex vector will be resized as needed and any...
escript::DataTypes::copySliceFrom
void copySliceFrom(VEC &left, const ShapeType &leftShape, typename VEC::size_type leftOffset, const VEC &other, const ShapeType &otherShape, typename VEC::size_type otherOffset, const RegionLoopRangeType &region)
Copy data into a slice specified by the given region and offset in the left vector from the other vec...
Definition: DataVector.h:281
escript::DataTypes::noValues
int noValues(const ShapeType &shape)
Calculate the number of values in a datapoint with the given shape.
Definition: DataTypes.cpp:91
escript::DataTypes::pointToStream
void pointToStream(std::ostream &os, const RealVectorType::ElementType *data, const ShapeType &shape, int offset, bool needsep=true, const std::string &sep=",")
Display a single value (with the specified shape) from the data.
escript::DataTypes::getResultSliceShape
DataTypes::ShapeType getResultSliceShape(const RegionType &region)
Determine the shape of the specified slice region.
Definition: DataTypes.cpp:173
escript::DataTypes::copyPoint
void copyPoint(RealVectorType &dest, vec_size_type doffset, vec_size_type nvals, const RealVectorType &src, vec_size_type soffset)
Copy a point from one vector to another. Note: This version does not check to see if shapes are the s...
Taipan.h
WrappedArray.h
escript::DataTypes::RegionLoopRangeType
std::vector< std::pair< int, int > > RegionLoopRangeType
Definition: DataTypes.h:73
DataVectorTaipan.h
escript::DataTypes::RealVectorType
escript::DataTypes::DataVectorAlt< real_t > RealVectorType
Vector to store underlying data.
Definition: DataVector.h:44
escript
Definition: AbstractContinuousDomain.cpp:23
DataVector.h
escript::DataTypes::pointToString
std::string pointToString(const RealVectorType &data, const ShapeType &shape, int offset, const std::string &prefix)
Display a single value (with the specified shape) from the data.
escript::DataTypes::getRelIndex
vec_size_type getRelIndex(const DataTypes::ShapeType &shape, vec_size_type i)
Compute the offset (in 1D vector) of a given subscript with a shape.
Definition: DataTypes.h:262
DataTypes.h
escript::DataTypes::cplx_t
std::complex< real_t > cplx_t
complex data type
Definition: DataTypes.h:82
ESYS_ASSERT
#define ESYS_ASSERT(a, b)
EsysAssert is a MACRO that will throw an exception if the boolean condition specified is false.
Definition: Assert.h:79
DataException.h
DataVectorAlt.h
escript::DataTypes::CplxVectorType
escript::DataTypes::DataVectorAlt< cplx_t > CplxVectorType
Definition: DataVector.h:45