/[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 1773 by jfenwick, Tue Sep 9 02:52:26 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          // 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       std::string
464       pointToString(const ValueType& data,const ShapeType& shape, int offset, const std::string& suffix)
465       {
466          using namespace std;
467          EsysAssert(data.size()>0,"Error - Data object is empty.");
468          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+offset];
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    
529    
530       void copyPoint(ValueType& dest, ValueType::size_type doffset, ValueType::size_type nvals, const ValueType& src, ValueType::size_type soffset)
531       {
532          EsysAssert((dest.size()>0&&src.size()>0&&checkOffset(doffset,dest.size(),nvals)),
533                     "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    
542    
543    
544       }
545    
546  }   // end namespace DataTypes  }   // end namespace DataTypes
547  }   // end namespace escript  }   // end namespace escript

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

  ViewVC Help
Powered by ViewVC 1.1.26