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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26