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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1734 - (show annotations)
Thu Aug 28 06:11:56 2008 UTC (11 years, 3 months ago) by jfenwick
File size: 16107 byte(s)
Added operations to Datatypes:
checkOffset
copySlice
copySliceFrom

Fixed some error reporting using EsysAssert.

Added two new c++ test suites:
DataTypesTest
DataMathsTest

Note that the test suite does not compile with dodebug=yes. There is an issue with linking one of the exception functions. I'm going to leave this 
until I have finished the rest of the work, perhaps Ken's scons changes will fix it.


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 //
228 // Make sure views are not empty
229
230 EsysAssert(!left.size()==0,
231 "Error - left data is empty.");
232 EsysAssert(!other.size()==0,
233 "Error - other data is empty.");
234
235 //
236 // Check the view to be sliced from is compatible with the region to be sliced,
237 // and that the region to be sliced is compatible with this view:
238
239 EsysAssert(checkOffset(thisOffset,left.size(),noValues(leftShape)),
240 "Error - offset incompatible with this view.");
241 EsysAssert(otherOffset+noValues(leftShape)<=other.size(),
242 "Error - offset incompatible with other view.");
243
244 EsysAssert(getRank(otherShape)==region.size(),
245 "Error - slice not same rank as view to be sliced from.");
246
247 EsysAssert(noValues(leftShape)==noValues(getResultSliceShape(region)),
248 "Error - slice shape not compatible shape for this view.");
249
250 //
251 // copy the values in the specified region of the other view into this view
252
253 // the following loops cannot be parallelised due to the numCopy counter
254 int numCopy=0;
255
256 switch (region.size()) {
257 case 0:
258 /* this case should never be encountered,
259 as python will never pass us an empty region.
260 here for completeness only, allows slicing of a scalar */
261 // (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex()];
262
263 left[thisOffset+numCopy]=other[otherOffset];
264 numCopy++;
265 break;
266 case 1:
267 for (int i=region[0].first;i<region[0].second;i++) {
268 left[thisOffset+numCopy]=other[otherOffset+getRelIndex(otherShape,i)];
269 numCopy++;
270 }
271 break;
272 case 2:
273 for (int j=region[1].first;j<region[1].second;j++) {
274 for (int i=region[0].first;i<region[0].second;i++) {
275 /* (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i,j)];*/
276 left[thisOffset+numCopy]=other[otherOffset+getRelIndex(otherShape,i,j)];
277 numCopy++;
278 }
279 }
280 break;
281 case 3:
282 for (int k=region[2].first;k<region[2].second;k++) {
283 for (int j=region[1].first;j<region[1].second;j++) {
284 for (int i=region[0].first;i<region[0].second;i++) {
285 // (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i,j,k)];
286 left[thisOffset+numCopy]=other[otherOffset+getRelIndex(otherShape,i,j,k)];
287 numCopy++;
288 }
289 }
290 }
291 break;
292 case 4:
293 for (int l=region[3].first;l<region[3].second;l++) {
294 for (int k=region[2].first;k<region[2].second;k++) {
295 for (int j=region[1].first;j<region[1].second;j++) {
296 for (int i=region[0].first;i<region[0].second;i++) {
297 /* (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i,j,k,l)];*/
298 left[thisOffset+numCopy]=other[otherOffset+getRelIndex(otherShape,i,j,k,l)];
299 numCopy++;
300 }
301 }
302 }
303 }
304 break;
305 default:
306 std::stringstream mess;
307 mess << "Error - (copySlice) Invalid slice region rank: " << region.size();
308 throw DataException(mess.str());
309 }
310 }
311
312
313 void
314 copySliceFrom(ValueType& left,
315 const ShapeType& leftShape,
316 ValueType::size_type thisOffset,
317 const ValueType& other,
318 const ShapeType& otherShape,
319 ValueType::size_type otherOffset,
320 const RegionLoopRangeType& region)
321 {
322
323 //
324 // Make sure views are not empty
325
326 EsysAssert(left.size()!=0,
327 "Error - this view is empty.");
328 EsysAssert(other.size()!=0,
329 "Error - other view is empty.");
330
331 //
332 // Check this view is compatible with the region to be sliced,
333 // and that the region to be sliced is compatible with the other view:
334
335 EsysAssert(checkOffset(otherOffset,other.size(),noValues(otherShape)),
336 "Error - offset incompatible with other view.");
337 EsysAssert(thisOffset+noValues(otherShape)<=left.size(),
338 "Error - offset incompatible with this view.");
339
340 EsysAssert(getRank(leftShape)==region.size(),
341 "Error - slice not same rank as this view.");
342
343 EsysAssert(getRank(otherShape)==0 || noValues(otherShape)==noValues(getResultSliceShape(region)),
344 "Error - slice shape not compatible shape for other view.");
345
346 //
347 // copy the values in the other view into the specified region of this view
348
349 // allow for case where other view is a scalar
350 if (getRank(otherShape)==0) {
351
352 // the following loops cannot be parallelised due to the numCopy counter
353 int numCopy=0;
354
355 switch (region.size()) {
356 case 0:
357 /* this case should never be encountered,
358 as python will never pass us an empty region.
359 here for completeness only, allows slicing of a scalar */
360 //(*m_data)[thisOffset+relIndex()]=(*other.m_data)[otherOffset];
361 left[thisOffset]=other[otherOffset];
362 numCopy++;
363 break;
364 case 1:
365 for (int i=region[0].first;i<region[0].second;i++) {
366 left[thisOffset+getRelIndex(leftShape,i)]=other[otherOffset];
367 numCopy++;
368 }
369 break;
370 case 2:
371 for (int j=region[1].first;j<region[1].second;j++) {
372 for (int i=region[0].first;i<region[0].second;i++) {
373 left[thisOffset+getRelIndex(leftShape,i,j)]=other[otherOffset];
374 numCopy++;
375 }
376 }
377 break;
378 case 3:
379 for (int k=region[2].first;k<region[2].second;k++) {
380 for (int j=region[1].first;j<region[1].second;j++) {
381 for (int i=region[0].first;i<region[0].second;i++) {
382 left[thisOffset+getRelIndex(leftShape,i,j,k)]=other[otherOffset];
383 numCopy++;
384 }
385 }
386 }
387 break;
388 case 4:
389 for (int l=region[3].first;l<region[3].second;l++) {
390 for (int k=region[2].first;k<region[2].second;k++) {
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[thisOffset+getRelIndex(leftShape,i,j,k,l)]=other[otherOffset];
394 numCopy++;
395 }
396 }
397 }
398 }
399 break;
400 default:
401 std::stringstream mess;
402 mess << "Error - (copySliceFrom) Invalid slice region rank: " << region.size();
403 throw DataException(mess.str());
404 }
405
406 } else {
407
408 // the following loops cannot be parallelised due to the numCopy counter
409 int numCopy=0;
410
411 switch (region.size()) {
412 case 0:
413 /* this case should never be encountered,
414 as python will never pass us an empty region.
415 here for completeness only, allows slicing of a scalar */
416 //(*m_data)[thisOffset+relIndex()]=(*other.m_data)[otherOffset+numCopy];
417 left[thisOffset]=other[otherOffset+numCopy];
418 numCopy++;
419 break;
420 case 1:
421 for (int i=region[0].first;i<region[0].second;i++) {
422 left[thisOffset+getRelIndex(leftShape,i)]=other[otherOffset+numCopy];
423 numCopy++;
424 }
425 break;
426 case 2:
427 for (int j=region[1].first;j<region[1].second;j++) {
428 for (int i=region[0].first;i<region[0].second;i++) {
429 left[thisOffset+getRelIndex(leftShape,i,j)]=other[otherOffset+numCopy];
430 numCopy++;
431 }
432 }
433 break;
434 case 3:
435 for (int k=region[2].first;k<region[2].second;k++) {
436 for (int j=region[1].first;j<region[1].second;j++) {
437 for (int i=region[0].first;i<region[0].second;i++) {
438 left[thisOffset+getRelIndex(leftShape,i,j,k)]=other[otherOffset+numCopy];
439 numCopy++;
440 }
441 }
442 }
443 break;
444 case 4:
445 for (int l=region[3].first;l<region[3].second;l++) {
446 for (int k=region[2].first;k<region[2].second;k++) {
447 for (int j=region[1].first;j<region[1].second;j++) {
448 for (int i=region[0].first;i<region[0].second;i++) {
449 left[thisOffset+getRelIndex(leftShape,i,j,k,l)]=other[otherOffset+numCopy];
450 numCopy++;
451 }
452 }
453 }
454 }
455 break;
456 default:
457 std::stringstream mess;
458 mess << "Error - (copySliceFrom) Invalid slice region rank: " << region.size();
459 throw DataException(mess.str());
460 }
461
462 }
463
464 }
465
466
467
468
469
470 } // end namespace DataTypes
471 } // end namespace escript

  ViewVC Help
Powered by ViewVC 1.1.26