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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1751 - (hide 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 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     // 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 jfenwick 1747 std::string
464     pointToString(const ValueType& data,const ShapeType& shape, int offset, const std::string& suffix)
465     {
466     using namespace std;
467 jfenwick 1751 EsysAssert(data.size()>0,"Error - View is empty.");
468 jfenwick 1747 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 jfenwick 1734
529    
530 jfenwick 1747 void copyPoint(ValueType& dest, ValueType::size_type doffset, ValueType::size_type nvals, const ValueType& src, ValueType::size_type soffset)
531     {
532 jfenwick 1751 EsysAssert((dest.size()>0&&src.size()>0&&checkOffset(doffset,dest.size(),nvals)),
533 jfenwick 1747 "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 jfenwick 1734
542    
543 jfenwick 1747
544     }
545    
546 jfenwick 1698 } // end namespace DataTypes
547     } // end namespace escript

  ViewVC Help
Powered by ViewVC 1.1.26