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