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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1734 - (hide annotations)
Thu Aug 28 06:11:56 2008 UTC (11 years, 5 months ago) by jfenwick
Original Path: branches/arrayview_from_1695_trunk/escript/src/DataTypes.cpp
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 jfenwick 1698 /* $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 jfenwick 1714 #include <boost/python/extract.hpp>
17     #include <boost/python/tuple.hpp>
18     #include "DataException.h"
19 jfenwick 1698 #include "DataTypes.h"
20    
21 jfenwick 1714 using namespace boost::python;
22    
23 jfenwick 1698 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 jfenwick 1714
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 jfenwick 1715 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 jfenwick 1714
193 jfenwick 1724
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 jfenwick 1734
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 jfenwick 1698 } // end namespace DataTypes
471     } // end namespace escript

  ViewVC Help
Powered by ViewVC 1.1.26