/[escript]/branches/numpy/escript/src/DataTypes.cpp
ViewVC logotype

Contents of /branches/numpy/escript/src/DataTypes.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2396 - (show annotations)
Thu Apr 23 23:58:29 2009 UTC (10 years ago) by jfenwick
File size: 19195 byte(s)
Branching to port escript to use numpy rather than numarray

1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2008 by University of Queensland
5 * Earth Systems Science Computational Center (ESSCC)
6 * http://www.uq.edu.au/esscc
7 *
8 * Primary Business: Queensland, Australia
9 * Licensed under the Open Software License version 3.0
10 * http://www.opensource.org/licenses/osl-3.0.php
11 *
12 *******************************************************/
13
14
15 #include "DataTypes.h"
16 #include <sstream>
17 #include <boost/python/extract.hpp>
18 #include <boost/python/tuple.hpp>
19 #include "DataException.h"
20
21 namespace {
22 using namespace boost::python;
23 using namespace escript;
24 using namespace escript::DataTypes;
25
26 /*
27 \brief
28 Calculate the slice range from the given python key object
29 Used by getSliceRegion - since it is not used anywhere else, I have elected not to declare it
30 in the header.
31 Returns the python slice object key as a pair of ints where the first
32 member is the start and the second member is the end. the presence of a possible
33 step attribute with value different from one will throw an exception
34
35 /param key - Input - key object specifying slice range.
36 */
37 std::pair<int,int>
38 getSliceRange(const boost::python::object& key,
39 const int shape)
40 {
41 /* default slice range is range of entire shape dimension */
42 int s0=0, s1=shape;;
43 extract<int> slice_int(key);
44 if (slice_int.check()) {
45 /* if the key is a single int set start=key and end=key */
46 /* in this case, we want to return a rank-1 dimension object from
47 this object, taken from a single index value for one of this
48 object's dimensions */
49 s0=slice_int();
50 s1=s0;
51 } else {
52 /* if key is a pair extract begin and end values */
53 extract<int> step(key.attr("step"));
54 if (step.check() && step()!=1) {
55 throw DataException("Error - Data does not support increments in slicing ");
56 } else {
57 extract<int> start(key.attr("start"));
58 if (start.check()) {
59 s0=start();
60 }
61 extract<int> stop(key.attr("stop"));
62 if (stop.check()) {
63 s1=stop();
64 }
65 }
66 }
67 if (s0 < 0)
68 throw DataException("Error - slice index out of range.");
69 if (s0 == s1 && s1 >= shape)
70 throw DataException("Error - slice index out of range.");
71 if (s0 != s1 && s1>shape)
72 throw DataException("Error - slice index out of range.");
73 if (s0 > s1)
74 throw DataException("Error - lower index must less or equal upper index.");
75 return std::pair<int,int>(s0,s1);
76 }
77 }
78
79
80 using namespace boost::python;
81
82 namespace escript
83 {
84 namespace DataTypes
85 {
86
87 int
88 noValues(const ShapeType& shape)
89 {
90 ShapeType::const_iterator i;
91 //
92 // An empty shape vector means rank 0 which contains 1 value
93 int noValues=1;
94 for (i=shape.begin();i!=shape.end();i++) {
95 noValues*=(*i);
96 }
97 return noValues;
98 }
99
100 int
101 noValues(const RegionLoopRangeType& region)
102 {
103 //
104 // An empty region vector means rank 0 which contains 1 value
105 int noValues=1;
106 unsigned int i;
107 for (i=0;i<region.size();i++) {
108 noValues*=region[i].second-region[i].first;
109 }
110 return noValues;
111 }
112
113 std::string
114 shapeToString(const DataTypes::ShapeType& shape)
115 {
116 std::stringstream temp;
117 temp << "(";
118 unsigned int i;
119 for (i=0;i<shape.size();i++) {
120 temp << shape[i];
121 if (i < shape.size()-1) {
122 temp << ",";
123 }
124 }
125 temp << ")";
126 return temp.str();
127 }
128
129
130
131
132
133 DataTypes::RegionType
134 getSliceRegion(const DataTypes::ShapeType& shape, const boost::python::object& key)
135 {
136 int slice_rank, i;
137 int this_rank=shape.size();
138 RegionType out(this_rank);
139 /* allow for case where key is singular eg: [1], this implies we
140 want to generate a rank-1 dimension object, as opposed to eg: [1,2]
141 which implies we want to take a rank dimensional object with one
142 dimension of size 1 */
143 extract<tuple> key_tuple(key);
144 if (key_tuple.check()) {
145 slice_rank=extract<int> (key.attr("__len__")());
146 /* ensure slice is correctly dimensioned */
147 if (slice_rank>this_rank) {
148 throw DataException("Error - rank of slices does not match rank of slicee");
149 } else {
150 /* calculate values for slice range */
151 for (i=0;i<slice_rank;i++) {
152 out[i]=getSliceRange(key[i],shape[i]);
153 }
154 }
155 } else {
156 slice_rank=1;
157 if (slice_rank>this_rank) {
158 throw DataException("Error - rank of slices does not match rank of slicee");
159 } else {
160 out[0]=getSliceRange(key,shape[0]);
161 }
162 }
163 for (i=slice_rank;i<this_rank;i++) {
164 out[i]=std::pair<int,int>(0,shape[i]);
165 }
166 return out;
167 }
168
169 DataTypes::ShapeType
170 getResultSliceShape(const RegionType& region)
171 {
172 int dimSize;
173 ShapeType result;
174 RegionType::const_iterator i;
175 for (i=region.begin();i!=region.end();i++) {
176 dimSize=((i->second)-(i->first));
177 if (dimSize!=0) {
178 result.push_back(dimSize);
179 }
180 }
181 return result;
182 }
183
184 DataTypes::RegionLoopRangeType
185 getSliceRegionLoopRange(const DataTypes::RegionType& region)
186 {
187 DataTypes::RegionLoopRangeType region_loop_range(region.size());
188 unsigned int i;
189 for (i=0;i<region.size();i++) {
190 if (region[i].first==region[i].second) {
191 region_loop_range[i].first=region[i].first;
192 region_loop_range[i].second=region[i].second+1;
193 } else {
194 region_loop_range[i].first=region[i].first;
195 region_loop_range[i].second=region[i].second;
196 }
197 }
198 return region_loop_range;
199 }
200
201
202 std::string
203 createShapeErrorMessage(const std::string& messagePrefix,
204 const DataTypes::ShapeType& other,
205 const DataTypes::ShapeType& thisShape)
206 {
207 std::stringstream temp;
208 temp << messagePrefix
209 << " This shape: " << shapeToString(thisShape)
210 << " Other shape: " << shapeToString(other);
211 return temp.str();
212 }
213
214
215 // Additional slice operations
216
217 inline
218 bool
219 checkOffset(ValueType::size_type offset, int size, int noval)
220 {
221 return (size >= (offset+noval));
222 }
223
224
225 void
226 copySlice(ValueType& left,
227 const ShapeType& leftShape,
228 ValueType::size_type thisOffset,
229 const ValueType& other,
230 const ShapeType& otherShape,
231 ValueType::size_type otherOffset,
232 const RegionLoopRangeType& region)
233 {
234 //
235 // Make sure views are not empty
236
237 EsysAssert(!left.size()==0,
238 "Error - left data is empty.");
239 EsysAssert(!other.size()==0,
240 "Error - other data is empty.");
241
242 //
243 // Check the view to be sliced from is compatible with the region to be sliced,
244 // and that the region to be sliced is compatible with this view:
245 EsysAssert(checkOffset(thisOffset,left.size(),noValues(leftShape)),
246 "Error - offset incompatible with this view.");
247 EsysAssert(otherOffset+noValues(leftShape)<=other.size(),
248 "Error - offset incompatible with other view.");
249
250 EsysAssert(getRank(otherShape)==region.size(),
251 "Error - slice not same rank as view to be sliced from.");
252
253 EsysAssert(noValues(leftShape)==noValues(getResultSliceShape(region)),
254 "Error - slice shape not compatible shape for this view.");
255
256 //
257 // copy the values in the specified region of the other view into this view
258
259 // the following loops cannot be parallelised due to the numCopy counter
260 int numCopy=0;
261
262 switch (region.size()) {
263 case 0:
264 /* this case should never be encountered,
265 as python will never pass us an empty region.
266 here for completeness only, allows slicing of a scalar */
267 // (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex()];
268
269 left[thisOffset+numCopy]=other[otherOffset];
270 numCopy++;
271 break;
272 case 1:
273 for (int i=region[0].first;i<region[0].second;i++) {
274 left[thisOffset+numCopy]=other[otherOffset+getRelIndex(otherShape,i)];
275 numCopy++;
276 }
277 break;
278 case 2:
279 for (int j=region[1].first;j<region[1].second;j++) {
280 for (int i=region[0].first;i<region[0].second;i++) {
281 /* (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i,j)];*/
282 left[thisOffset+numCopy]=other[otherOffset+getRelIndex(otherShape,i,j)];
283 numCopy++;
284 }
285 }
286 break;
287 case 3:
288 for (int k=region[2].first;k<region[2].second;k++) {
289 for (int j=region[1].first;j<region[1].second;j++) {
290 for (int i=region[0].first;i<region[0].second;i++) {
291 // (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i,j,k)];
292 left[thisOffset+numCopy]=other[otherOffset+getRelIndex(otherShape,i,j,k)];
293 numCopy++;
294 }
295 }
296 }
297 break;
298 case 4:
299 for (int l=region[3].first;l<region[3].second;l++) {
300 for (int k=region[2].first;k<region[2].second;k++) {
301 for (int j=region[1].first;j<region[1].second;j++) {
302 for (int i=region[0].first;i<region[0].second;i++) {
303 /* (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i,j,k,l)];*/
304 left[thisOffset+numCopy]=other[otherOffset+getRelIndex(otherShape,i,j,k,l)];
305 numCopy++;
306 }
307 }
308 }
309 }
310 break;
311 default:
312 std::stringstream mess;
313 mess << "Error - (copySlice) Invalid slice region rank: " << region.size();
314 throw DataException(mess.str());
315 }
316 }
317
318
319 void
320 copySliceFrom(ValueType& left,
321 const ShapeType& leftShape,
322 ValueType::size_type thisOffset,
323 const ValueType& other,
324 const ShapeType& otherShape,
325 ValueType::size_type otherOffset,
326 const RegionLoopRangeType& region)
327 {
328 //
329 // Make sure views are not empty
330
331 EsysAssert(left.size()!=0,
332 "Error - this view is empty.");
333 EsysAssert(other.size()!=0,
334 "Error - other view is empty.");
335
336 //
337 // Check this view is compatible with the region to be sliced,
338 // and that the region to be sliced is compatible with the other view:
339
340 EsysAssert(checkOffset(otherOffset,other.size(),noValues(otherShape)),
341 "Error - offset incompatible with other view.");
342 EsysAssert(thisOffset+noValues(otherShape)<=left.size(),
343 "Error - offset incompatible with this view.");
344
345 EsysAssert(getRank(leftShape)==region.size(),
346 "Error - slice not same rank as this view.");
347
348 EsysAssert(getRank(otherShape)==0 || noValues(otherShape)==noValues(getResultSliceShape(region)),
349 "Error - slice shape not compatible shape for other view.");
350
351 //
352 // copy the values in the other view into the specified region of this view
353
354 // allow for case where other view is a scalar
355 if (getRank(otherShape)==0) {
356
357 // the following loops cannot be parallelised due to the numCopy counter
358 int numCopy=0;
359
360 switch (region.size()) {
361 case 0:
362 /* this case should never be encountered,
363 as python will never pass us an empty region.
364 here for completeness only, allows slicing of a scalar */
365 //(*m_data)[thisOffset+relIndex()]=(*other.m_data)[otherOffset];
366 left[thisOffset]=other[otherOffset];
367 numCopy++;
368 break;
369 case 1:
370 for (int i=region[0].first;i<region[0].second;i++) {
371 left[thisOffset+getRelIndex(leftShape,i)]=other[otherOffset];
372 numCopy++;
373 }
374 break;
375 case 2:
376 for (int j=region[1].first;j<region[1].second;j++) {
377 for (int i=region[0].first;i<region[0].second;i++) {
378 left[thisOffset+getRelIndex(leftShape,i,j)]=other[otherOffset];
379 numCopy++;
380 }
381 }
382 break;
383 case 3:
384 for (int k=region[2].first;k<region[2].second;k++) {
385 for (int j=region[1].first;j<region[1].second;j++) {
386 for (int i=region[0].first;i<region[0].second;i++) {
387 left[thisOffset+getRelIndex(leftShape,i,j,k)]=other[otherOffset];
388 numCopy++;
389 }
390 }
391 }
392 break;
393 case 4:
394 for (int l=region[3].first;l<region[3].second;l++) {
395 for (int k=region[2].first;k<region[2].second;k++) {
396 for (int j=region[1].first;j<region[1].second;j++) {
397 for (int i=region[0].first;i<region[0].second;i++) {
398 left[thisOffset+getRelIndex(leftShape,i,j,k,l)]=other[otherOffset];
399 numCopy++;
400 }
401 }
402 }
403 }
404 break;
405 default:
406 std::stringstream mess;
407 mess << "Error - (copySliceFrom) Invalid slice region rank: " << region.size();
408 throw DataException(mess.str());
409 }
410
411 } else {
412
413 // the following loops cannot be parallelised due to the numCopy counter
414 int numCopy=0;
415
416 switch (region.size()) {
417 case 0:
418 /* this case should never be encountered,
419 as python will never pass us an empty region.
420 here for completeness only, allows slicing of a scalar */
421 //(*m_data)[thisOffset+relIndex()]=(*other.m_data)[otherOffset+numCopy];
422 left[thisOffset]=other[otherOffset+numCopy];
423 numCopy++;
424 break;
425 case 1:
426 for (int i=region[0].first;i<region[0].second;i++) {
427 left[thisOffset+getRelIndex(leftShape,i)]=other[otherOffset+numCopy];
428 numCopy++;
429 }
430 break;
431 case 2:
432 for (int j=region[1].first;j<region[1].second;j++) {
433 for (int i=region[0].first;i<region[0].second;i++) {
434 left[thisOffset+getRelIndex(leftShape,i,j)]=other[otherOffset+numCopy];
435 numCopy++;
436 }
437 }
438 break;
439 case 3:
440 for (int k=region[2].first;k<region[2].second;k++) {
441 for (int j=region[1].first;j<region[1].second;j++) {
442 for (int i=region[0].first;i<region[0].second;i++) {
443 left[thisOffset+getRelIndex(leftShape,i,j,k)]=other[otherOffset+numCopy];
444 numCopy++;
445 }
446 }
447 }
448 break;
449 case 4:
450 for (int l=region[3].first;l<region[3].second;l++) {
451 for (int k=region[2].first;k<region[2].second;k++) {
452 for (int j=region[1].first;j<region[1].second;j++) {
453 for (int i=region[0].first;i<region[0].second;i++) {
454 left[thisOffset+getRelIndex(leftShape,i,j,k,l)]=other[otherOffset+numCopy];
455 numCopy++;
456 }
457 }
458 }
459 }
460 break;
461 default:
462 std::stringstream mess;
463 mess << "Error - (copySliceFrom) Invalid slice region rank: " << region.size();
464 throw DataException(mess.str());
465 }
466
467 }
468
469 }
470
471 std::string
472 pointToString(const ValueType& data,const ShapeType& shape, int offset, const std::string& prefix)
473 {
474 using namespace std;
475 EsysAssert(data.size()>0,"Error - Data object is empty.");
476 stringstream temp;
477 string finalPrefix=prefix;
478 if (prefix.length() > 0) {
479 finalPrefix+=" ";
480 }
481 switch (getRank(shape)) {
482 case 0:
483 temp << finalPrefix << data[offset];
484 break;
485 case 1:
486 for (int i=0;i<shape[0];i++) {
487 temp << finalPrefix << "(" << i << ") " << data[i+offset];
488 if (i!=(shape[0]-1)) {
489 temp << endl;
490 }
491 }
492 break;
493 case 2:
494 for (int i=0;i<shape[0];i++) {
495 for (int j=0;j<shape[1];j++) {
496 temp << finalPrefix << "(" << i << "," << j << ") " << data[offset+getRelIndex(shape,i,j)];
497 if (!(i==(shape[0]-1) && j==(shape[1]-1))) {
498 temp << endl;
499 }
500 }
501 }
502 break;
503 case 3:
504 for (int i=0;i<shape[0];i++) {
505 for (int j=0;j<shape[1];j++) {
506 for (int k=0;k<shape[2];k++) {
507 temp << finalPrefix << "(" << i << "," << j << "," << k << ") " << data[offset+getRelIndex(shape,i,j,k)];
508 if (!(i==(shape[0]-1) && j==(shape[1]-1) && k==(shape[2]-1))) {
509 temp << endl;
510 }
511 }
512 }
513 }
514 break;
515 case 4:
516 for (int i=0;i<shape[0];i++) {
517 for (int j=0;j<shape[1];j++) {
518 for (int k=0;k<shape[2];k++) {
519 for (int l=0;l<shape[3];l++) {
520 temp << finalPrefix << "(" << i << "," << j << "," << k << "," << l << ") " << data[offset+getRelIndex(shape,i,j,k,l)];
521 if (!(i==(shape[0]-1) && j==(shape[1]-1) && k==(shape[2]-1) && l==(shape[3]-1))) {
522 temp << endl;
523 }
524 }
525 }
526 }
527 }
528 break;
529 default:
530 stringstream mess;
531 mess << "Error - (toString) Invalid rank: " << getRank(shape);
532 throw DataException(mess.str());
533 }
534 return temp.str();
535 }
536
537
538 void copyPoint(ValueType& dest, ValueType::size_type doffset, ValueType::size_type nvals, const ValueType& src, ValueType::size_type soffset)
539 {
540 EsysAssert((dest.size()>0&&src.size()>0&&checkOffset(doffset,dest.size(),nvals)),
541 "Error - Couldn't copy due to insufficient storage.");
542 // EsysAssert((checkShape(other.getShape())),
543 // createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",other.getShape(),m_shape));
544 if (checkOffset(doffset,dest.size(),nvals) && checkOffset(soffset,src.size(),nvals)) {
545 memcpy(&dest[doffset],&src[soffset],sizeof(double)*nvals);
546 } else {
547 throw DataException("Error - invalid offset specified.");
548 }
549
550
551
552 }
553
554 } // end namespace DataTypes
555 } // end namespace escript

  ViewVC Help
Powered by ViewVC 1.1.26