escript  Revision_
DataVector.h
Go to the documentation of this file.
1 
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2020 by The University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Apache License, version 2.0
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development 2012-2013 by School of Earth Sciences
13 * Development from 2014-2017 by Centre for Geoscience Computing (GeoComp)
14 * Development from 2019 by School of Earth and Environmental Sciences
15 **
16 *****************************************************************************/
17 
18 #ifndef __ESCRIPT_DATAVECTOR_H__
19 #define __ESCRIPT_DATAVECTOR_H__
20 
21 #include "system_dep.h"
22 #include "DataTypes.h"
23 #include "Assert.h"
24 #include "DataVectorAlt.h"
25 #include "DataVectorTaipan.h"
26 
27 #ifdef ESYS_HAVE_BOOST_NUMPY
28 #include <boost/python.hpp>
29 #include <boost/python/numpy.hpp>
30 #endif
31 
32 // ensure that nobody else tries to instantiate the complex version
34 
35 
36 namespace escript {
37 
38 // Functions in DataTypes:: which manipulate DataVectors
39 namespace DataTypes
40 {
41 
42  // This is the main version we had
43  //typedef DataVectorTaipan DataVector;
46 
62  void
63  pointToStream(std::ostream& os, const RealVectorType::ElementType* data,const ShapeType& shape, int offset, bool needsep=true, const std::string& sep=",");
64 
80  void
81  pointToStream(std::ostream& os, const CplxVectorType::ElementType* data,const ShapeType& shape, int offset, bool needsep=true, const std::string& sep=",");
82 
83 
84 #ifdef ESYS_HAVE_BOOST_NUMPY
85 
86  void
87  pointToNumpyArrayOld(boost::python::numpy::ndarray& dataArray, const RealVectorType::ElementType* data, const ShapeType& shape, long offset, long numsamples, long dpps, long numdata);
88 
89 
90  void
91  pointToNumpyArrayOld(boost::python::numpy::ndarray& dataArray, const CplxVectorType::ElementType* data, const ShapeType& shape, long offset, long numsamples, long dpps, long numdata);
92 
93 
105  void
106  pointToNumpyArray(boost::python::numpy::ndarray& dataArray, const RealVectorType::ElementType* data, const ShapeType& shape, int offset, int d, int index);
107 
108 
109  void
110  pointToNumpyArray(boost::python::numpy::ndarray& dataArray, const CplxVectorType::ElementType* data, const ShapeType& shape, int offset, int d, int index);
111 #endif
112 
121  std::string
122  pointToString(const RealVectorType& data,const ShapeType& shape, int offset, const std::string& prefix);
123 
124 
125  std::string
126  pointToString(const CplxVectorType& data,const ShapeType& shape, int offset, const std::string& prefix);
127 
137  void copyPoint(RealVectorType& dest, vec_size_type doffset, vec_size_type nvals, const RealVectorType& src, vec_size_type soffset);
138 
148  void copyPoint(CplxVectorType& dest, vec_size_type doffset, vec_size_type nvals, const CplxVectorType& src, vec_size_type soffset);
149 
156 
171  template <class VEC>
173  void
174  copySlice(VEC& left,
175  const ShapeType& leftShape,
176  typename VEC::size_type leftOffset,
177  const VEC& other,
178  const ShapeType& otherShape,
179  typename VEC::size_type otherOffset,
180  const RegionLoopRangeType& region)
181  {
182  //
183  // Make sure vectors are not empty
184 
185  ESYS_ASSERT(!left.size()==!1, "left data is empty."); //Note: !1=0, but clang returns an error if the rhs is 0 here
186  ESYS_ASSERT(!other.size()==!1, "other data is empty.");
187 
188  //
189  // Check the vector to be sliced from is compatible with the region to be sliced,
190  // and that the region to be sliced is compatible with this vector:
191  ESYS_ASSERT(checkOffset(leftOffset,left.size(),noValues(leftShape)),
192  "offset incompatible with this vector.");
193  ESYS_ASSERT(otherOffset+noValues(leftShape)<=other.size(),
194  "offset incompatible with other vector.");
195 
196  ESYS_ASSERT(getRank(otherShape)==region.size(),
197  "slice not same rank as vector to be sliced from.");
198 
199  ESYS_ASSERT(noValues(leftShape)==noValues(getResultSliceShape(region)),
200  "slice shape not compatible shape for this vector.");
201 
202  //
203  // copy the values in the specified region of the other vector into this vector
204 
205  // the following loops cannot be parallelised due to the numCopy counter
206  int numCopy=0;
207 
208  switch (region.size()) {
209  case 0:
210  /* this case should never be encountered,
211  as python will never pass us an empty region.
212  here for completeness only, allows slicing of a scalar */
213 // (*m_data)[leftOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex()];
214 
215  left[leftOffset+numCopy]=other[otherOffset];
216  numCopy++;
217  break;
218  case 1:
219  for (int i=region[0].first;i<region[0].second;i++) {
220  left[leftOffset+numCopy]=other[otherOffset+getRelIndex(otherShape,i)];
221  numCopy++;
222  }
223  break;
224  case 2:
225  for (int j=region[1].first;j<region[1].second;j++) {
226  for (int i=region[0].first;i<region[0].second;i++) {
227 /* (*m_data)[leftOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i,j)];*/
228  left[leftOffset+numCopy]=other[otherOffset+getRelIndex(otherShape,i,j)];
229  numCopy++;
230  }
231  }
232  break;
233  case 3:
234  for (int k=region[2].first;k<region[2].second;k++) {
235  for (int j=region[1].first;j<region[1].second;j++) {
236  for (int i=region[0].first;i<region[0].second;i++) {
237 // (*m_data)[leftOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i,j,k)];
238  left[leftOffset+numCopy]=other[otherOffset+getRelIndex(otherShape,i,j,k)];
239  numCopy++;
240  }
241  }
242  }
243  break;
244  case 4:
245  for (int l=region[3].first;l<region[3].second;l++) {
246  for (int k=region[2].first;k<region[2].second;k++) {
247  for (int j=region[1].first;j<region[1].second;j++) {
248  for (int i=region[0].first;i<region[0].second;i++) {
249 /* (*m_data)[leftOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i,j,k,l)];*/
250  left[leftOffset+numCopy]=other[otherOffset+getRelIndex(otherShape,i,j,k,l)];
251  numCopy++;
252  }
253  }
254  }
255  }
256  break;
257  default:
258  std::stringstream mess;
259  mess << "Error - (copySlice) Invalid slice region rank: " << region.size();
260  throw DataException(mess.str());
261  }
262  }
263 
278  template<typename VEC>
280  void
281  copySliceFrom(VEC& left,
282  const ShapeType& leftShape,
283  typename VEC::size_type leftOffset,
284  const VEC& other,
285  const ShapeType& otherShape,
286  typename VEC::size_type otherOffset,
287  const RegionLoopRangeType& region)
288  {
289  //
290  // Make sure vectors are not empty
291 
292  ESYS_ASSERT(left.size()!=0, "this vector is empty.");
293  ESYS_ASSERT(other.size()!=0, "other vector is empty.");
294 
295  //
296  // Check this vector is compatible with the region to be sliced,
297  // and that the region to be sliced is compatible with the other vector:
298 
299  ESYS_ASSERT(checkOffset(otherOffset,other.size(),noValues(otherShape)),
300  "offset incompatible with other vector.");
301  ESYS_ASSERT(leftOffset+noValues(otherShape)<=left.size(),
302  "offset incompatible with this vector.");
303 
304  ESYS_ASSERT(getRank(leftShape)==region.size(),
305  "slice not same rank as this vector.");
306 
307  ESYS_ASSERT(getRank(otherShape)==0 || noValues(otherShape)==noValues(getResultSliceShape(region)),
308  "slice shape not compatible shape for other vector.");
309 
310  //
311  // copy the values in the other vector into the specified region of this vector
312 
313  // allow for case where other vector is a scalar
314  if (getRank(otherShape)==0) {
315 
316  // the following loops cannot be parallelised due to the numCopy counter
317  int numCopy=0;
318 
319  switch (region.size()) {
320  case 0:
321  /* this case should never be encountered,
322  as python will never pass us an empty region.
323  here for completeness only, allows slicing of a scalar */
324  //(*m_data)[leftOffset+relIndex()]=(*other.m_data)[otherOffset];
325  left[leftOffset]=other[otherOffset];
326  numCopy++;
327  break;
328  case 1:
329  for (int i=region[0].first;i<region[0].second;i++) {
330  left[leftOffset+getRelIndex(leftShape,i)]=other[otherOffset];
331  numCopy++;
332  }
333  break;
334  case 2:
335  for (int j=region[1].first;j<region[1].second;j++) {
336  for (int i=region[0].first;i<region[0].second;i++) {
337  left[leftOffset+getRelIndex(leftShape,i,j)]=other[otherOffset];
338  numCopy++;
339  }
340  }
341  break;
342  case 3:
343  for (int k=region[2].first;k<region[2].second;k++) {
344  for (int j=region[1].first;j<region[1].second;j++) {
345  for (int i=region[0].first;i<region[0].second;i++) {
346  left[leftOffset+getRelIndex(leftShape,i,j,k)]=other[otherOffset];
347  numCopy++;
348  }
349  }
350  }
351  break;
352  case 4:
353  for (int l=region[3].first;l<region[3].second;l++) {
354  for (int k=region[2].first;k<region[2].second;k++) {
355  for (int j=region[1].first;j<region[1].second;j++) {
356  for (int i=region[0].first;i<region[0].second;i++) {
357  left[leftOffset+getRelIndex(leftShape,i,j,k,l)]=other[otherOffset];
358  numCopy++;
359  }
360  }
361  }
362  }
363  break;
364  default:
365  std::stringstream mess;
366  mess << "Error - (copySliceFrom) Invalid slice region rank: " << region.size();
367  throw DataException(mess.str());
368  }
369 
370  } else {
371 
372  // the following loops cannot be parallelised due to the numCopy counter
373  int numCopy=0;
374 
375  switch (region.size()) {
376  case 0:
377  /* this case should never be encountered,
378  as python will never pass us an empty region.
379  here for completeness only, allows slicing of a scalar */
380  //(*m_data)[leftOffset+relIndex()]=(*other.m_data)[otherOffset+numCopy];
381  left[leftOffset]=other[otherOffset+numCopy];
382  numCopy++;
383  break;
384  case 1:
385  for (int i=region[0].first;i<region[0].second;i++) {
386  left[leftOffset+getRelIndex(leftShape,i)]=other[otherOffset+numCopy];
387  numCopy++;
388  }
389  break;
390  case 2:
391  for (int j=region[1].first;j<region[1].second;j++) {
392  for (int i=region[0].first;i<region[0].second;i++) {
393  left[leftOffset+getRelIndex(leftShape,i,j)]=other[otherOffset+numCopy];
394  numCopy++;
395  }
396  }
397  break;
398  case 3:
399  for (int k=region[2].first;k<region[2].second;k++) {
400  for (int j=region[1].first;j<region[1].second;j++) {
401  for (int i=region[0].first;i<region[0].second;i++) {
402  left[leftOffset+getRelIndex(leftShape,i,j,k)]=other[otherOffset+numCopy];
403  numCopy++;
404  }
405  }
406  }
407  break;
408  case 4:
409  for (int l=region[3].first;l<region[3].second;l++) {
410  for (int k=region[2].first;k<region[2].second;k++) {
411  for (int j=region[1].first;j<region[1].second;j++) {
412  for (int i=region[0].first;i<region[0].second;i++) {
413  left[leftOffset+getRelIndex(leftShape,i,j,k,l)]=other[otherOffset+numCopy];
414  numCopy++;
415  }
416  }
417  }
418  }
419  break;
420  default:
421  std::stringstream mess;
422  mess << "Error - (copySliceFrom) Invalid slice region rank: " << region.size();
423  throw DataException(mess.str());
424  }
425 
426  }
427 
428  }
429 }
430 
431 } // end of namespace
432 
433 #endif // __ESCRIPT_DATAVECTOR_H__
434 
ESCRIPT_DLL_API
#define ESCRIPT_DLL_API
Definition: escriptcore/src/system_dep.h:29
escript::DataTypes::real_t
double real_t
type of all real-valued scalars in escript
Definition: DataTypes.h:76
escript::DataTypes::DataVectorAlt::ElementType
T ElementType
Definition: DataVectorAlt.h:81
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:173
escript::DataTypes::getRank
int getRank(const DataTypes::ShapeType &shape)
Return the rank (number of dimensions) of the given shape.
Definition: DataTypes.h:244
escript::DataTypes::checkOffset
bool checkOffset(vec_size_type offset, int size, int noval)
Definition: DataTypes.h:349
escript::DataTypes::DataVectorAlt
Definition: DataVectorAlt.h:62
Assert.h
system_dep.h
escript::DataTypes::vec_size_type
long vec_size_type
Definition: DataTypes.h:73
escript::DataException
Definition: DataException.h:39
escript::DataTypes::ShapeType
std::vector< int > ShapeType
The shape of a single datapoint.
Definition: DataTypes.h:68
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:280
escript::DataTypes::noValues
int noValues(const ShapeType &shape)
Calculate the number of values in a datapoint with the given shape.
Definition: DataTypes.cpp:90
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:172
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:70
DataVectorTaipan.h
escript::DataTypes::RealVectorType
escript::DataTypes::DataVectorAlt< real_t > RealVectorType
Vector to store underlying data.
Definition: DataVector.h:43
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:259
DataTypes.h
escript::DataTypes::cplx_t
std::complex< real_t > cplx_t
complex data type
Definition: DataTypes.h:79
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:78
DataException.h
DataVectorAlt.h
escript::DataTypes::CplxVectorType
escript::DataTypes::DataVectorAlt< cplx_t > CplxVectorType
Definition: DataVector.h:44