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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1781 - (hide annotations)
Thu Sep 11 05:03:14 2008 UTC (11 years, 5 months ago) by jfenwick
Original Path: branches/arrayview_from_1695_trunk/escript/src/DataTypes.cpp
File size: 19215 byte(s)
Branch commit

Merged changes from trunk version 1695 up to and including version 1779.


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 1779 namespace {
22 jfenwick 1714 using namespace boost::python;
23 jfenwick 1779 using namespace escript;
24     using namespace escript::DataTypes;
25 jfenwick 1714
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 jfenwick 1779 }
78 jfenwick 1714
79    
80 jfenwick 1779 using namespace boost::python;
81 jfenwick 1714
82 jfenwick 1779 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 jfenwick 1714 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 jfenwick 1715 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 jfenwick 1714
201 jfenwick 1724
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 jfenwick 1734
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 jfenwick 1747 std::string
472 jfenwick 1779 pointToString(const ValueType& data,const ShapeType& shape, int offset, const std::string& prefix)
473 jfenwick 1747 {
474     using namespace std;
475 jfenwick 1773 EsysAssert(data.size()>0,"Error - Data object is empty.");
476 jfenwick 1747 stringstream temp;
477 jfenwick 1779 string finalPrefix=prefix;
478     if (prefix.length() > 0) {
479     finalPrefix+=" ";
480 jfenwick 1747 }
481     switch (getRank(shape)) {
482     case 0:
483 jfenwick 1779 temp << finalPrefix << data[0];
484 jfenwick 1747 break;
485     case 1:
486     for (int i=0;i<shape[0];i++) {
487 jfenwick 1779 temp << finalPrefix << "(" << i << ") " << data[i+offset];
488 jfenwick 1747 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 jfenwick 1779 temp << finalPrefix << "(" << i << "," << j << ") " << data[offset+getRelIndex(shape,i,j)];
497 jfenwick 1747 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 jfenwick 1779 temp << finalPrefix << "(" << i << "," << j << "," << k << ") " << data[offset+getRelIndex(shape,i,j,k)];
508 jfenwick 1747 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 jfenwick 1779 temp << finalPrefix << "(" << i << "," << j << "," << k << "," << l << ") " << data[offset+getRelIndex(shape,i,j,k,l)];
521 jfenwick 1747 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 jfenwick 1734
537    
538 jfenwick 1747 void copyPoint(ValueType& dest, ValueType::size_type doffset, ValueType::size_type nvals, const ValueType& src, ValueType::size_type soffset)
539     {
540 jfenwick 1751 EsysAssert((dest.size()>0&&src.size()>0&&checkOffset(doffset,dest.size(),nvals)),
541 jfenwick 1747 "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 jfenwick 1734
550    
551 jfenwick 1747
552     }
553    
554 jfenwick 1698 } // end namespace DataTypes
555 jfenwick 1781 } // end namespace escript

  ViewVC Help
Powered by ViewVC 1.1.26