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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1802 - (hide annotations)
Tue Sep 23 01:03:29 2008 UTC (11 years, 5 months ago) by jfenwick
File size: 19220 byte(s)
Added canTag methods to FunctionSpace and AbstractDomain (and its 
offspring).
This checks to see if the domain supports tags for the given type of 
function space.

Constructors for DataTagged now throw exceptions if you attempt to make 
a DataTagged with a FunctionSpace which does not support tags.
To allow the default constructor to work, NullDomain has a single 
functioncode which "supports" tagging.

Fixed a bug in DataTagged::toString and DataTypes::pointToString.

Added FunctionSpace::getListOfTagsSTL.

algorithm(DataTagged, BinaryFunction) in DataAlgorithm now only 
processes tags known to be in use.
This fixes mantis issue #0000186.

Added comment to Data.h intro warning about holding references if the 
underlying DataAbstract changes.

_python_ unit tests have been updated to test TaggedData with invalid 
FunctionSpaces and to give the correct answers to Lsup etc.


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 1802 temp << finalPrefix << data[offset];
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