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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1698 by jfenwick, Tue Aug 12 01:13:16 2008 UTC revision 1747 by jfenwick, Wed Sep 3 04:56:50 2008 UTC
# Line 13  Line 13 
13  *******************************************************/  *******************************************************/
14    
15  #include <sstream>  #include <sstream>
16    #include <boost/python/extract.hpp>
17    #include <boost/python/tuple.hpp>
18    #include "DataException.h"
19  #include "DataTypes.h"  #include "DataTypes.h"
20    
21    using namespace boost::python;
22    
23  namespace escript  namespace escript
24  {  {
25  namespace DataTypes  namespace DataTypes
# Line 62  namespace DataTypes Line 67  namespace DataTypes
67        return temp.str();        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       std::string
467       pointToString(const ValueType& data,const ShapeType& shape, int offset, const std::string& suffix)
468       {
469          using namespace std;
470          EsysAssert(!data.isEmpty(),"Error - View is empty.");
471          stringstream temp;
472          string finalSuffix=suffix;
473          if (suffix.length() > 0) {
474             finalSuffix+=" ";
475          }
476          switch (getRank(shape)) {
477          case 0:
478             temp << finalSuffix << data[0];
479             break;
480          case 1:
481             for (int i=0;i<shape[0];i++) {
482                temp << finalSuffix << "(" << i << ") " << data[i];
483                if (i!=(shape[0]-1)) {
484                   temp << endl;
485                }
486             }
487             break;
488          case 2:
489             for (int i=0;i<shape[0];i++) {
490                for (int j=0;j<shape[1];j++) {
491                   temp << finalSuffix << "(" << i << "," << j << ") " << data[offset+getRelIndex(shape,i,j)];
492                   if (!(i==(shape[0]-1) && j==(shape[1]-1))) {
493                      temp << endl;
494                   }
495                }
496             }
497             break;
498          case 3:
499             for (int i=0;i<shape[0];i++) {
500                for (int j=0;j<shape[1];j++) {
501                   for (int k=0;k<shape[2];k++) {
502                      temp << finalSuffix << "(" << i << "," << j << "," << k << ") " << data[offset+getRelIndex(shape,i,j,k)];
503                      if (!(i==(shape[0]-1) && j==(shape[1]-1) && k==(shape[2]-1))) {
504                         temp << endl;
505                      }
506                   }
507                }
508             }
509             break;
510          case 4:
511             for (int i=0;i<shape[0];i++) {
512                for (int j=0;j<shape[1];j++) {
513                   for (int k=0;k<shape[2];k++) {
514                      for (int l=0;l<shape[3];l++) {
515                         temp << finalSuffix << "(" << i << "," << j << "," << k << "," << l << ") " << data[offset+getRelIndex(shape,i,j,k,l)];
516                         if (!(i==(shape[0]-1) && j==(shape[1]-1) && k==(shape[2]-1) && l==(shape[3]-1))) {
517                            temp << endl;
518                         }
519                      }
520                   }
521                }
522             }
523             break;
524          default:
525             stringstream mess;
526             mess << "Error - (toString) Invalid rank: " << getRank(shape);
527             throw DataException(mess.str());
528          }
529          return temp.str();
530       }
531    
532    
533       void copyPoint(ValueType& dest, ValueType::size_type doffset, ValueType::size_type nvals, const ValueType& src, ValueType::size_type soffset)
534       {
535          EsysAssert((!dest.isEmpty()&&!other.isEmpty()&&checkOffset(doffset,dest.size(),nvals)),
536                     "Error - Couldn't copy due to insufficient storage.");
537    //       EsysAssert((checkShape(other.getShape())),
538    //                  createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",other.getShape(),m_shape));
539          if (checkOffset(doffset,dest.size(),nvals) && checkOffset(soffset,src.size(),nvals)) {
540             memcpy(&dest[doffset],&src[soffset],sizeof(double)*nvals);
541          } else {
542             throw DataException("Error - invalid offset specified.");
543          }
544    
545    
546    
547       }
548    
549  }   // end namespace DataTypes  }   // end namespace DataTypes
550  }   // end namespace escript  }   // end namespace escript

Legend:
Removed from v.1698  
changed lines
  Added in v.1747

  ViewVC Help
Powered by ViewVC 1.1.26