/[escript]/trunk/escript/src/DataTypes.cpp
ViewVC logotype

Contents of /trunk/escript/src/DataTypes.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2005 - (show annotations)
Mon Nov 10 01:21:39 2008 UTC (10 years, 10 months ago) by jfenwick
File size: 19195 byte(s)
Bringing all changes across from schroedinger.
(Note this does not mean development is done, just that it will happen
on the trunk for now).
If anyone notices any problems please contact me.


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