/[escript]/branches/arrayview_from_1695_trunk/escript/src/DataArrayView.cpp
ViewVC logotype

Diff of /branches/arrayview_from_1695_trunk/escript/src/DataArrayView.cpp

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

trunk/esys2/escript/src/Data/DataArrayView.cpp revision 100 by jgs, Wed Dec 15 03:48:48 2004 UTC trunk/escript/src/DataArrayView.cpp revision 1388 by trankine, Fri Jan 11 07:45:58 2008 UTC
# Line 1  Line 1 
 // $Id$  
1    
2  /*  /* $Id$ */
  ******************************************************************************  
  *                                                                            *  
  *       COPYRIGHT  ACcESS 2004 -  All Rights Reserved                        *  
  *                                                                            *  
  * This software is the property of ACcESS. No part of this code              *  
  * may be copied in any form or by any means without the expressed written    *  
  * consent of ACcESS.  Copying, use or modification of this software          *  
  * by any unauthorised person is illegal unless that person has a software    *  
  * license agreement with ACcESS.                                             *  
  *                                                                            *  
  ******************************************************************************  
 */  
3    
4  #include "escript/Data/DataArrayView.h"  /*******************************************************
5  #include "escript/Data/DataException.h"   *
6     *           Copyright 2003-2007 by ACceSS MNRF
7     *       Copyright 2007 by University of Queensland
8     *
9     *                http://esscc.uq.edu.au
10     *        Primary Business: Queensland, Australia
11     *  Licensed under the Open Software License version 3.0
12     *     http://www.opensource.org/licenses/osl-3.0.php
13     *
14     *******************************************************/
15    
16    #include "DataArrayView.h"
17    #include "DataException.h"
18    
19  #include <sstream>  #include <sstream>
20  #include <vector>  
 #include <iostream>  
21  #include <boost/python/extract.hpp>  #include <boost/python/extract.hpp>
 #include <boost/python/object.hpp>  
22    
23  using namespace std;  using namespace std;
24  using namespace boost::python;  using namespace boost::python;
# Line 30  namespace escript { Line 27  namespace escript {
27    
28  DataArrayView::DataArrayView():  DataArrayView::DataArrayView():
29      m_shape(),      m_shape(),
30      m_offset(),      m_offset(0),
31      m_noValues()      m_noValues(0)
32  {  {
33      m_data=0;    m_data=0;
34  }  }
35    
36  DataArrayView::DataArrayView(ValueType& data,  DataArrayView::DataArrayView(ValueType& data,
# Line 45  DataArrayView::DataArrayView(ValueType& Line 42  DataArrayView::DataArrayView(ValueType&
42      m_noValues(noValues(viewShape))      m_noValues(noValues(viewShape))
43  {  {
44      //      //
45      // check the shape and throw an exception if an illegal shape is specified      // check the shape rank and size and throw an exception if a
46        // shape incompatible with the given data array is specified
47      if (viewShape.size() > m_maxRank) {      if (viewShape.size() > m_maxRank) {
48        stringstream temp;        stringstream temp;
49        temp << "Error- Couldn't construct DataArrayView, invalid rank."        temp << "Error- Couldn't construct DataArrayView, invalid rank."
50         << "Input data rank: " << viewShape.size()         << "Input data rank: " << viewShape.size() << " "
51         << " Maximum rank allowed for DataArrayView: "         << "Maximum rank allowed for DataArrayView: " << m_maxRank;
        << m_maxRank;  
52        throw DataException(temp.str());        throw DataException(temp.str());
53      }      }
54      if (m_data->size()<(noValues()+m_offset)) {      if (m_data->size() < (m_offset+noValues())) {
55        stringstream temp;        stringstream temp;
56        temp << "Error- Couldn't construct DataArrayView, insufficient storage."        temp << "Error- Couldn't construct DataArrayView, insufficient data values. "
57         << " Shape requires: " << noValues()+m_offset << " values."         << "Shape requires: " << noValues()+m_offset << " values."
58         << " Data values available: " << m_data->size();         << "Data values available: " << m_data->size();
59        throw DataException(temp.str());        throw DataException(temp.str());
60      }      }
61  }  }
# Line 74  DataArrayView::DataArrayView(const DataA Line 71  DataArrayView::DataArrayView(const DataA
71  bool  bool
72  DataArrayView::isEmpty() const  DataArrayView::isEmpty() const
73  {  {
74      return(m_data==0);    return (m_data==0);
75  }  }
76    
77  void  void
78  DataArrayView::copy(const boost::python::numeric::array& value)  DataArrayView::copy(const boost::python::numeric::array& value)
79  {  {
80      DataArrayView::ShapeType tempShape;          ShapeType tempShape;    
81      for (int i=0; i<value.getrank(); ++i) {      for (int i=0; i<value.getrank(); i++) {
82        tempShape.push_back(extract<int>(value.getshape()[i]));        tempShape.push_back(extract<int>(value.getshape()[i]));
83      }      }
84    
85      EsysAssert((!isEmpty()&&checkShape(tempShape)),      EsysAssert((!isEmpty()&&checkShape(tempShape)),
86             createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",tempShape));             createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",tempShape));
87    
88      if (value.getrank()==0) {      if (value.getrank()==0) {
89        (*this)()=extract<double>(value[value.getshape()]);        (*this)()=extract<double>(value[value.getshape()]);
90      } else if (value.getrank()==1) {      } else if (value.getrank()==1) {
91        for (ValueType::size_type i=0;i<tempShape[0];++i) {        for (ValueType::size_type i=0;i<tempShape[0];i++) {
92      (*this)(i)=extract<double>(value[i]);      (*this)(i)=extract<double>(value[i]);
93        }        }
94      } else if (value.getrank()==2) {      } else if (value.getrank()==2) {
95        for (ValueType::size_type i=0;i<tempShape[0];++i) {        for (ValueType::size_type i=0;i<tempShape[0];i++) {
96      for (ValueType::size_type j=0;j<tempShape[1];++j) {      for (ValueType::size_type j=0;j<tempShape[1];j++) {
97        (*this)(i,j)=extract<double>(value[i][j]);        (*this)(i,j)=extract<double>(value[i][j]);
98      }      }
99        }        }
100      } else if (value.getrank()==3) {      } else if (value.getrank()==3) {
101        for (ValueType::size_type i=0;i<tempShape[0];++i) {        for (ValueType::size_type i=0;i<tempShape[0];i++) {
102      for (ValueType::size_type j=0;j<tempShape[1];++j) {      for (ValueType::size_type j=0;j<tempShape[1];j++) {
103        for (ValueType::size_type k=0;k<tempShape[2];++k) {        for (ValueType::size_type k=0;k<tempShape[2];k++) {
104          (*this)(i,j,k)=extract<double>(value[i][j][k]);          (*this)(i,j,k)=extract<double>(value[i][j][k]);
105        }        }
106      }      }
107        }        }
108      } else if (value.getrank()==4) {      } else if (value.getrank()==4) {
109        for (ValueType::size_type i=0;i<tempShape[0];++i) {        for (ValueType::size_type i=0;i<tempShape[0];i++) {
110      for (ValueType::size_type j=0;j<tempShape[1];++j) {      for (ValueType::size_type j=0;j<tempShape[1];j++) {
111        for (ValueType::size_type k=0;k<tempShape[2];++k) {        for (ValueType::size_type k=0;k<tempShape[2];k++) {
112          for (ValueType::size_type m=0;m<tempShape[3];++m) {          for (ValueType::size_type l=0;l<tempShape[3];l++) {
113            (*this)(i,j,k,m)=extract<double>(value[i][j][k][m]);            (*this)(i,j,k,l)=extract<double>(value[i][j][k][l]);
114          }          }
115        }        }
116      }      }
117        }        }
118      }      }
119  }  }
120      
121  void  void
122  DataArrayView::copy(const DataArrayView& other)  DataArrayView::copy(const DataArrayView& other)
123  {  {
124      copy(m_offset,other);    copy(m_offset,other);
125  }  }
126    
127  void  void
128  DataArrayView::copy(ValueType::size_type offset,  DataArrayView::copy(ValueType::size_type offset,
129                      const DataArrayView& other)                      const DataArrayView& other)
130  {  {
131      EsysAssert((!isEmpty()&&checkShape(other.getShape())),      EsysAssert((!isEmpty()&&!other.isEmpty()&&checkOffset(offset)),
132               "Error - Couldn't copy due to insufficient storage.");
133        EsysAssert((checkShape(other.getShape())),
134             createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",other.getShape()));             createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",other.getShape()));
135      memcpy(&(*m_data)[offset],&(*other.m_data)[other.m_offset],sizeof(double)*noValues());      if (checkOffset(offset)) {
136          memcpy(&(*m_data)[offset],&(*other.m_data)[other.m_offset],sizeof(double)*noValues());
137        } else {
138          throw DataException("Error - invalid offset specified.");
139        }
140  }  }
141    
142  void  void
# Line 140  DataArrayView::copy(ValueType::size_type Line 144  DataArrayView::copy(ValueType::size_type
144                      const DataArrayView& other,                      const DataArrayView& other,
145                      ValueType::size_type otherOffset)                      ValueType::size_type otherOffset)
146  {  {
147      EsysAssert((!isEmpty()&&checkShape(other.getShape())),      EsysAssert((!isEmpty()&&!other.isEmpty()&&checkOffset(thisOffset)&&other.checkOffset(otherOffset)),
148               "Error - Couldn't copy due to insufficient storage.");
149        EsysAssert((checkShape(other.getShape())),
150             createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",other.getShape()));             createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",other.getShape()));
151      memcpy(&(*m_data)[thisOffset],&(*other.m_data)[otherOffset],sizeof(double)*noValues());      if (checkOffset(thisOffset)&&other.checkOffset(otherOffset)) {
152          memcpy(&(*m_data)[thisOffset],&(*other.m_data)[otherOffset],sizeof(double)*noValues());
153        } else {
154          throw DataException("Error - invalid offset specified.");
155        }
156  }  }
157    
158  void  void
159  DataArrayView::copy(ValueType::size_type offset,  DataArrayView::copy(ValueType::size_type offset,
160                      ValueType::value_type value)                      ValueType::value_type value)
161  {  {
162      //      EsysAssert((!isEmpty()&&checkOffset(offset)),
163      // fill the entire view with the single value             "Error - Couldn't copy due to insufficient storage.");
164      EsysAssert(!isEmpty(),"Error - View is empty");      if (checkOffset(offset)) {
165      ValueType temp(noValues(),value);        vector<double> temp(noValues(),value);
166      memcpy(&(*m_data)[offset],&temp[0],sizeof(double)*noValues());        memcpy(&(*m_data)[offset],&temp[0],sizeof(double)*noValues());
167        } else {
168          throw DataException("Error - invalid offset specified.");
169        }
170  }  }
171    
172  int  int
173  DataArrayView::getRank() const  DataArrayView::getRank() const
174  {  {
175      return m_shape.size();    return m_shape.size();
176  }  }
177    
178  const DataArrayView::ShapeType&  const
179    DataArrayView::ShapeType&
180  DataArrayView::getShape() const  DataArrayView::getShape() const
181  {  {
182      return m_shape;    return m_shape;
183  }  }
184    
185  int  int
186  DataArrayView::noValues(const ShapeType& shape)  DataArrayView::noValues(const ShapeType& shape)
187  {  {
188      ShapeType::const_iterator i;    ShapeType::const_iterator i;
189      //    //
190      // An empty shape vector means rank 0 which contains 1 value    // An empty shape vector means rank 0 which contains 1 value
191      int noValues=1;    int noValues=1;
192      for (i=shape.begin();i!=shape.end();++i) {    for (i=shape.begin();i!=shape.end();i++) {
193        noValues*=(*i);      noValues*=(*i);
194      }    }
195      return noValues;    return noValues;
196    }
197    
198    int
199    DataArrayView::noValues(const RegionLoopRangeType& region)
200    {
201      //
202      // An empty region vector means rank 0 which contains 1 value
203      int noValues=1;
204      for (int i=0;i<region.size();i++) {
205        noValues*=region[i].second-region[i].first;
206      }
207      return noValues;
208  }  }
209      
210  int  int
211  DataArrayView::noValues() const  DataArrayView::noValues() const
212  {  {
213      return m_noValues;    return m_noValues;
214  }  }
215    
216  bool  bool
217  DataArrayView::checkShape(const DataArrayView::ShapeType& other) const  DataArrayView::checkShape(const DataArrayView::ShapeType& other) const
218  {  {
219      return (m_shape==other);    return (m_shape==other);
220  }  }
221    
222  string  string
223  DataArrayView::createShapeErrorMessage(const string& messagePrefix,  DataArrayView::createShapeErrorMessage(const string& messagePrefix,
224                                         const DataArrayView::ShapeType& other) const                                         const DataArrayView::ShapeType& other) const
225  {  {
226      stringstream temp;    stringstream temp;
227      temp << messagePrefix    temp << messagePrefix
228       << " This shape: " << shapeToString(m_shape)         << " This shape: " << shapeToString(m_shape)
229       << " Other shape: " << shapeToString(other);         << " Other shape: " << shapeToString(other);
230      return temp.str();    return temp.str();
231  }  }
232    
233    DataArrayView::ValueType::size_type
234    DataArrayView::getOffset() const
235    {
236      return m_offset;
237    }
238    
239  void  void
240  DataArrayView::setOffset(ValueType::size_type offset)  DataArrayView::setOffset(ValueType::size_type offset)
241  {  {
242      EsysAssert((m_data->size()>=(noValues()+offset)), "Invalid offset");    EsysAssert((checkOffset(offset)), "Error - Invalid offset.");
243      if (checkOffset(offset)) {
244      m_offset=offset;      m_offset=offset;
245      } else {
246        throw DataException("Error - invalid offset specified.");
247      }
248    }
249    
250    void
251    DataArrayView::incrOffset()
252    {
253      EsysAssert((checkOffset(m_offset+noValues())), "Error - Cannot increment offset.");
254      if (checkOffset(m_offset+noValues())) {
255        m_offset=m_offset+noValues();
256      } else {
257        throw DataException("Error - Cannot increment offset.");
258      }
259    }
260    
261    bool
262    DataArrayView::checkOffset() const
263    {
264      return checkOffset(m_offset);
265    }
266    
267    bool
268    DataArrayView::checkOffset(ValueType::size_type offset) const
269    {
270      return (m_data->size() >= (offset+noValues()));
271    }
272    
273    DataArrayView::ValueType&
274    DataArrayView::getData() const
275    {
276      EsysAssert(!isEmpty(),"Error - View is empty");
277      return *m_data;
278    }
279    
280    DataArrayView::ValueType::reference
281    DataArrayView::getData(ValueType::size_type i) const
282    {
283      //
284      // don't do any index checking so allow one past the end of the
285      // vector providing the equivalent of end()
286      return (*m_data)[m_offset+i];
287  }  }
288    
289  DataArrayView::ShapeType  DataArrayView::ShapeType
290  DataArrayView::getResultSliceShape(const RegionType& region)  DataArrayView::getResultSliceShape(const RegionType& region)
291  {  {
292      RegionType::const_iterator i;    int dimSize;
293      DataArrayView::ShapeType result;    RegionType::const_iterator i;
294      for (i=region.begin();i!=region.end();++i) {    ShapeType result;
295        int dimSize=(i->second-i->first);    for (i=region.begin();i!=region.end();i++) {
296        if (dimSize!=0) {      dimSize=((i->second)-(i->first));
297      result.push_back(dimSize);      if (dimSize!=0) {
298          result.push_back(dimSize);
299        }
300      }
301      return result;
302    }
303    
304    DataArrayView::RegionType
305    DataArrayView::getSliceRegion(const boost::python::object& key) const
306    {
307      int slice_rank, i;
308      int this_rank=getRank();
309      RegionType out(this_rank);
310      /* allow for case where key is singular eg: [1], this implies we
311         want to generate a rank-1 dimension object, as opposed to eg: [1,2]
312         which implies we want to take a rank dimensional object with one
313         dimension of size 1 */
314      extract<tuple> key_tuple(key);
315      if (key_tuple.check()) {
316          slice_rank=extract<int> (key.attr("__len__")());
317          /* ensure slice is correctly dimensioned */
318          if (slice_rank>this_rank) {
319               throw DataException("Error - rank of slices does not match rank of slicee");
320          } else {
321             /* calculate values for slice range */
322             for (i=0;i<slice_rank;i++) {
323               out[i]=getSliceRange(key[i],getShape()[i]);
324             }
325          }
326      } else {
327          slice_rank=1;
328          if (slice_rank>this_rank) {
329               throw DataException("Error - rank of slices does not match rank of slicee");
330          } else {
331               out[0]=getSliceRange(key,getShape()[0]);
332          }
333      }
334      for (i=slice_rank;i<this_rank;i++) {
335        out[i]=std::pair<int,int>(0,getShape()[i]);
336      }
337      return out;
338    }
339    
340    std::pair<int,int>
341    getSliceRange(const boost::python::object& key,
342                  const int shape)
343    {
344      /* default slice range is range of entire shape dimension */
345      int s0=0, s1=shape;;
346      extract<int> slice_int(key);
347      if (slice_int.check()) {
348        /* if the key is a single int set start=key and end=key */
349        /* in this case, we want to return a rank-1 dimension object from
350           this object, taken from a single index value for one of this
351           object's dimensions */
352        s0=slice_int();
353        s1=s0;
354      } else {
355        /* if key is a pair extract begin and end values */
356        extract<int> step(key.attr("step"));
357        if (step.check() && step()!=1) {
358           throw DataException("Error - Data does not support increments in slicing ");
359        } else {
360          extract<int> start(key.attr("start"));
361          if (start.check()) {
362            s0=start();
363          }
364          extract<int> stop(key.attr("stop"));
365          if (stop.check()) {
366            s1=stop();
367          }
368        }
369      }
370      if (s0 < 0)
371         throw DataException("Error - slice index out of range.");
372      if (s0 == s1 && s1 >= shape)
373         throw DataException("Error - slice index out of range.");
374      if (s0 != s1 &&  s1>shape)
375         throw DataException("Error - slice index out of range.");
376      if (s0 > s1)
377         throw DataException("Error - lower index must less or equal upper index.");
378      return std::pair<int,int>(s0,s1);
379    }
380    
381    DataArrayView::RegionLoopRangeType
382    getSliceRegionLoopRange(const DataArrayView::RegionType& region)
383    {
384        DataArrayView::RegionLoopRangeType region_loop_range(region.size());
385        for (int i=0;i<region.size();i++) {
386          if (region[i].first==region[i].second) {
387            region_loop_range[i].first=region[i].first;
388            region_loop_range[i].second=region[i].second+1;
389          } else {
390            region_loop_range[i].first=region[i].first;
391            region_loop_range[i].second=region[i].second;
392        }        }
393      }      }
394      return result;      return region_loop_range;
395  }  }
396    
397  void  void
398  DataArrayView::copySlice(const DataArrayView& other,  DataArrayView::copySlice(const DataArrayView& other,
399                           const RegionType& region)                           const RegionLoopRangeType& region)
400  {  {
401      //    //
402      // version of slice that uses the default offset    // version of copySlice that uses the default offsets
403      copySlice(m_offset,other,other.m_offset,region);    copySlice(m_offset,other,other.m_offset,region);
404  }  }
405    
406  void  void
407  DataArrayView::copySlice(ValueType::size_type thisOffset,  DataArrayView::copySlice(ValueType::size_type thisOffset,
408                           const DataArrayView& other,                           const DataArrayView& other,
409                           ValueType::size_type otherOffset,                           ValueType::size_type otherOffset,
410                           const RegionType& region)                           const RegionLoopRangeType& region)
411  {  {
412    
413      //      //
414      // Assume region is okay, Other asserts will detect range errors      // Make sure views are not empty
415      EsysAssert((other.getRank()==region.size()),"Error - Invalid slice region.");  
416      EsysAssert((!isEmpty()&&checkShape(getResultSliceShape(region))),      EsysAssert(!isEmpty(),
417             createShapeErrorMessage("Error - Couldn't copy slice due to shape mismatch.",getResultSliceShape(region)));                 "Error - this view is empty.");
418        EsysAssert(!other.isEmpty(),
419                   "Error - other view is empty.");
420    
421      //      //
422      // take a local copy of the region and modify for the case      // Check the view to be sliced from is compatible with the region to be sliced,
423      // of 2:2 being equivalent to 2:3. The only difference in meaning being in      // and that the region to be sliced is compatible with this view:
424      // determining the output shape.  
425      RegionType localRegion=region;      EsysAssert(checkOffset(thisOffset),
426      for (int i=0;i<localRegion.size();++i) {                 "Error - offset incompatible with this view.");
427        if (localRegion[i].first==localRegion[i].second) {      EsysAssert(otherOffset+noValues()<=other.getData().size(),
428      ++(localRegion[i].second);                 "Error - offset incompatible with other view.");
429        }  
430      }      EsysAssert(other.getRank()==region.size(),
431                   "Error - slice not same rank as view to be sliced from.");
432    
433        EsysAssert(noValues()==noValues(getResultSliceShape(region)),
434                   "Error - slice shape not compatible shape for this view.");
435    
436        //
437        // copy the values in the specified region of the other view into this view
438    
439        // the following loops cannot be parallelised due to the numCopy counter
440      int numCopy=0;      int numCopy=0;
441      switch (localRegion.size()) {  
442        switch (region.size()) {
443      case 0:      case 0:
444        (*m_data)[thisOffset]=(*other.m_data)[otherOffset];        /* this case should never be encountered,
445             as python will never pass us an empty region.
446             here for completeness only, allows slicing of a scalar */
447          (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex()];
448          numCopy++;
449        break;        break;
450      case 1:      case 1:
451        for (int i=localRegion[0].first;i<localRegion[0].second;++i) {        for (int i=region[0].first;i<region[0].second;i++) {
452      (*m_data)[numCopy+thisOffset]=(*other.m_data)[i+otherOffset];      (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i)];
453      ++numCopy;          numCopy++;
454        }        }
455        break;        break;
456      case 2:      case 2:
457        for (int j=localRegion[1].first;j<localRegion[1].second;++j) {        for (int j=region[1].first;j<region[1].second;j++) {
458      for (int i=localRegion[0].first;i<localRegion[0].second;++i) {      for (int i=region[0].first;i<region[0].second;i++) {
459        (*m_data)[numCopy+thisOffset]=(*other.m_data)[other.relIndex(i,j)+otherOffset];        (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i,j)];
460        ++numCopy;        numCopy++;
461      }      }
462        }        }
463        break;        break;
464      case 3:      case 3:
465        for (int k=localRegion[2].first;k<localRegion[2].second;++k) {        for (int k=region[2].first;k<region[2].second;k++) {
466      for (int j=localRegion[1].first;j<localRegion[1].second;++j) {      for (int j=region[1].first;j<region[1].second;j++) {
467        for (int i=localRegion[0].first;i<localRegion[0].second;++i) {        for (int i=region[0].first;i<region[0].second;i++) {
468          (*m_data)[numCopy+thisOffset]=(*other.m_data)[other.relIndex(i,j,k)+otherOffset];          (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i,j,k)];
469          ++numCopy;          numCopy++;
470        }        }
471      }      }
472        }        }
473        break;        break;
474      case 4:      case 4:
475        for (int m=localRegion[2].first;m<localRegion[2].second;++m) {        for (int l=region[3].first;l<region[3].second;l++) {
476      for (int k=localRegion[2].first;k<localRegion[2].second;++k) {      for (int k=region[2].first;k<region[2].second;k++) {
477        for (int j=localRegion[1].first;j<localRegion[1].second;++j) {        for (int j=region[1].first;j<region[1].second;j++) {
478          for (int i=localRegion[0].first;i<localRegion[0].second;++i) {          for (int i=region[0].first;i<region[0].second;i++) {
479            (*m_data)[numCopy+thisOffset]=(*other.m_data)[other.relIndex(i,j,k,m)+otherOffset];            (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i,j,k,l)];
480            ++numCopy;            numCopy++;
481          }          }
482        }        }
483      }      }
# Line 298  DataArrayView::copySlice(ValueType::size Line 485  DataArrayView::copySlice(ValueType::size
485        break;        break;
486      default:      default:
487        stringstream mess;        stringstream mess;
488        mess << "Error - (copySlice) Invalid rank.";        mess << "Error - (copySlice) Invalid slice region rank: " << region.size();
489        throw DataException(mess.str());        throw DataException(mess.str());
490      }      }
491  }  }
492      
493  void  void
494  DataArrayView::copySliceFrom(const DataArrayView& other,  DataArrayView::copySliceFrom(const DataArrayView& other,
495                               const RegionType& region)                               const RegionLoopRangeType& region_loop_range)
496  {  {
497      //      //
498      // version of slice that uses the default offset      // version of copySliceFrom that uses the default offsets
499      copySliceFrom(m_offset,other,other.m_offset,region);      copySliceFrom(m_offset,other,other.m_offset,region_loop_range);
500  }  }
501    
502  void  void
503  DataArrayView::copySliceFrom(ValueType::size_type thisOffset,  DataArrayView::copySliceFrom(ValueType::size_type thisOffset,
504                               const DataArrayView& other,                               const DataArrayView& other,
505                               ValueType::size_type otherOffset,                               ValueType::size_type otherOffset,
506                               const RegionType& region)                               const RegionLoopRangeType& region)
507  {  {
508    
509      //      //
510      // Assume region is okay, Other asserts will detect range errors      // Make sure views are not empty
511      EsysAssert((getRank()==region.size()),"Error - Invalid slice region.");  
512      EsysAssert((!isEmpty()&&other.checkShape(getResultSliceShape(region))),      EsysAssert(!isEmpty(),
513             createShapeErrorMessage("Error - Couldn't copy slice due to shape mismatch.",                 "Error - this view is empty.");
514                         getResultSliceShape(region)));      EsysAssert(!other.isEmpty(),
515                   "Error - other view is empty.");
516    
517      //      //
518      // take a local copy of the region and modify for the case      // Check this view is compatible with the region to be sliced,
519      // of 2:2 being equivalent to 2:3. The only difference in meaning being in      // and that the region to be sliced is compatible with the other view:
520      // determining the output shape.  
521      RegionType localRegion=region;      EsysAssert(other.checkOffset(otherOffset),
522      for (int i=0;i<localRegion.size();++i) {                 "Error - offset incompatible with other view.");
523        if (localRegion[i].first==localRegion[i].second) {      EsysAssert(thisOffset+other.noValues()<=getData().size(),
524      ++(localRegion[i].second);                 "Error - offset incompatible with this view.");
525        }  
526      }      EsysAssert(getRank()==region.size(),
527      int numCopy=0;                 "Error - slice not same rank as this view.");
528      switch (localRegion.size()) {  
529      case 0:      EsysAssert(other.getRank()==0 || other.noValues()==noValues(getResultSliceShape(region)),
530        (*m_data)[thisOffset]=(*other.m_data)[otherOffset];                 "Error - slice shape not compatible shape for other view.");
531        break;  
532      case 1:      //
533        for (int i=localRegion[0].first;i<localRegion[0].second;++i) {      // copy the values in the other view into the specified region of this view
534      (*m_data)[i+thisOffset]=(*other.m_data)[numCopy+otherOffset];  
535      ++numCopy;      // allow for case where other view is a scalar
536        }      if (other.getRank()==0) {
537        break;  
538      case 2:          // the following loops cannot be parallelised due to the numCopy counter
539        for (int j=localRegion[1].first;j<localRegion[1].second;++j) {          int numCopy=0;
540      for (int i=localRegion[0].first;i<localRegion[0].second;++i) {  
541        (*m_data)[relIndex(i,j)+thisOffset]=(*other.m_data)[numCopy+otherOffset];          switch (region.size()) {
542        ++numCopy;          case 0:
543      }          /* this case should never be encountered,
544        }             as python will never pass us an empty region.
545        break;             here for completeness only, allows slicing of a scalar */
546      case 3:            (*m_data)[thisOffset+relIndex()]=(*other.m_data)[otherOffset];
547        for (int k=localRegion[2].first;k<localRegion[2].second;++k) {        numCopy++;
548      for (int j=localRegion[1].first;j<localRegion[1].second;++j) {            break;
549        for (int i=localRegion[0].first;i<localRegion[0].second;++i) {          case 1:
550          (*m_data)[relIndex(i,j,k)+thisOffset]=(*other.m_data)[numCopy+otherOffset];            for (int i=region[0].first;i<region[0].second;i++) {
551          ++numCopy;          (*m_data)[thisOffset+relIndex(i)]=(*other.m_data)[otherOffset];
552        }          numCopy++;
553      }            }
554        }            break;
555        break;          case 2:
556      case 4:            for (int j=region[1].first;j<region[1].second;j++) {
557        for (int m=localRegion[2].first;m<localRegion[2].second;++m) {          for (int i=region[0].first;i<region[0].second;i++) {
558      for (int k=localRegion[2].first;k<localRegion[2].second;++k) {            (*m_data)[thisOffset+relIndex(i,j)]=(*other.m_data)[otherOffset];
559        for (int j=localRegion[1].first;j<localRegion[1].second;++j) {            numCopy++;
         for (int i=localRegion[0].first;i<localRegion[0].second;++i) {  
           (*m_data)[relIndex(i,j,k,m)+thisOffset]=(*other.m_data)[numCopy+otherOffset];  
           ++numCopy;  
560          }          }
561        }            }
562      }            break;
563        }          case 3:
564        break;            for (int k=region[2].first;k<region[2].second;k++) {
565      default:          for (int j=region[1].first;j<region[1].second;j++) {
566        stringstream mess;            for (int i=region[0].first;i<region[0].second;i++) {
567        mess << "Error - (copySliceFrom) Invalid rank: " << localRegion.size() << " for region.";              (*m_data)[thisOffset+relIndex(i,j,k)]=(*other.m_data)[otherOffset];
568        throw DataException(mess.str());              numCopy++;
569      }            }
570  }          }
571                }
572  DataArrayView::ShapeType            break;
573  DataArrayView::determineResultShape(const DataArrayView& left,          case 4:
574                                      const DataArrayView& right)            for (int l=region[3].first;l<region[3].second;l++) {
575  {          for (int k=region[2].first;k<region[2].second;k++) {
576      DataArrayView::ShapeType temp;            for (int j=region[1].first;j<region[1].second;j++) {
577      for (int i=0; i<(left.getRank()-1); ++i) {              for (int i=region[0].first;i<region[0].second;i++) {
578        temp.push_back(left.getShape()[i]);                (*m_data)[thisOffset+relIndex(i,j,k,l)]=(*other.m_data)[otherOffset];
579      }                numCopy++;
580      for (int i=1; i<right.getRank(); ++i) {              }
581        temp.push_back(right.getShape()[i]);            }
582            }
583              }
584              break;
585            default:
586              stringstream mess;
587              mess << "Error - (copySliceFrom) Invalid slice region rank: " << region.size();
588              throw DataException(mess.str());
589            }
590    
591        } else {
592    
593            // the following loops cannot be parallelised due to the numCopy counter
594            int numCopy=0;
595    
596            switch (region.size()) {
597            case 0:
598            /* this case should never be encountered,
599               as python will never pass us an empty region.
600               here for completeness only, allows slicing of a scalar */
601              (*m_data)[thisOffset+relIndex()]=(*other.m_data)[otherOffset+numCopy];
602          numCopy++;
603              break;
604            case 1:
605              for (int i=region[0].first;i<region[0].second;i++) {
606            (*m_data)[thisOffset+relIndex(i)]=(*other.m_data)[otherOffset+numCopy];
607            numCopy++;
608              }
609              break;
610            case 2:
611              for (int j=region[1].first;j<region[1].second;j++) {
612            for (int i=region[0].first;i<region[0].second;i++) {
613              (*m_data)[thisOffset+relIndex(i,j)]=(*other.m_data)[otherOffset+numCopy];
614              numCopy++;
615            }
616              }
617              break;
618            case 3:
619              for (int k=region[2].first;k<region[2].second;k++) {
620            for (int j=region[1].first;j<region[1].second;j++) {
621              for (int i=region[0].first;i<region[0].second;i++) {
622                (*m_data)[thisOffset+relIndex(i,j,k)]=(*other.m_data)[otherOffset+numCopy];
623                numCopy++;
624              }
625            }
626              }
627              break;
628            case 4:
629              for (int l=region[3].first;l<region[3].second;l++) {
630            for (int k=region[2].first;k<region[2].second;k++) {
631              for (int j=region[1].first;j<region[1].second;j++) {
632                for (int i=region[0].first;i<region[0].second;i++) {
633                  (*m_data)[thisOffset+relIndex(i,j,k,l)]=(*other.m_data)[otherOffset+numCopy];
634                  numCopy++;
635                }
636              }
637            }
638              }
639              break;
640            default:
641              stringstream mess;
642              mess << "Error - (copySliceFrom) Invalid slice region rank: " << region.size();
643              throw DataException(mess.str());
644            }
645    
646      }      }
     return temp;  
 }  
647    
648    }
649    
650  string  string
651  DataArrayView::toString(const string& suffix) const  DataArrayView::toString(const string& suffix) const
652  {  {
653      EsysAssert(!isEmpty(),"Error - View is empty");      EsysAssert(!isEmpty(),"Error - View is empty.");
654      stringstream temp;      stringstream temp;
655      string finalSuffix=suffix;      string finalSuffix=suffix;
656      if (suffix.length() > 0) {      if (suffix.length() > 0) {
# Line 410  DataArrayView::toString(const string& su Line 661  DataArrayView::toString(const string& su
661        temp << finalSuffix << (*this)();        temp << finalSuffix << (*this)();
662        break;        break;
663      case 1:      case 1:
664        for (int i=0;i<getShape()[0];++i) {        for (int i=0;i<getShape()[0];i++) {
665      temp << "(" << i << ") " << finalSuffix << (*this)(i);      temp << finalSuffix << "(" << i << ") " << (*this)(i);
666      if (i!=(getShape()[0]-1)) {      if (i!=(getShape()[0]-1)) {
667        temp << endl;        temp << endl;
668      }      }
669        }        }
670        break;        break;
671      case 2:      case 2:
672        for (int i=0;i<getShape()[0];++i) {        for (int i=0;i<getShape()[0];i++) {
673      for (int j=0;j<getShape()[1];++j) {      for (int j=0;j<getShape()[1];j++) {
674        temp << "(" << i << "," << j << ") " << finalSuffix << (*this)(i,j);        temp << finalSuffix << "(" << i << "," << j << ") " << (*this)(i,j);
675        if (!(i==(getShape()[0]-1) && j==(getShape()[1]-1))) {        if (!(i==(getShape()[0]-1) && j==(getShape()[1]-1))) {
676          temp << endl;          temp << endl;
677        }        }
# Line 428  DataArrayView::toString(const string& su Line 679  DataArrayView::toString(const string& su
679        }        }
680        break;        break;
681      case 3:      case 3:
682        for (int i=0;i<getShape()[0];++i) {        for (int i=0;i<getShape()[0];i++) {
683      for (int j=0;j<getShape()[1];++j) {      for (int j=0;j<getShape()[1];j++) {
684        for (int k=0;k<getShape()[2];++k) {        for (int k=0;k<getShape()[2];k++) {
685          temp << "("          temp << finalSuffix << "(" << i << "," << j << "," << k << ") " << (*this)(i,j,k);
686           << i << ","          if (!(i==(getShape()[0]-1) && j==(getShape()[1]-1) && k==(getShape()[2]-1))) {
          << j << ","  
          << k << ") "  
          << finalSuffix << (*this)(i,j,k);  
         //  
         // don't put an endl after the last value  
         if (!(i==(getShape()[0]-1) &&  
           j==(getShape()[1]-1) &&  
           k==(getShape()[2]-1))) {  
687            temp << endl;            temp << endl;
688          }          }
689        }        }
# Line 448  DataArrayView::toString(const string& su Line 691  DataArrayView::toString(const string& su
691        }        }
692        break;        break;
693      case 4:      case 4:
694        // break;        for (int i=0;i<getShape()[0];i++) {
695        for (int j=0;j<getShape()[1];j++) {
696          for (int k=0;k<getShape()[2];k++) {
697            for (int l=0;l<getShape()[3];l++) {
698              temp << finalSuffix << "(" << i << "," << j << "," << k << "," << l << ") " << (*this)(i,j,k,l);
699              if (!(i==(getShape()[0]-1) && j==(getShape()[1]-1) && k==(getShape()[2]-1) && l==(getShape()[3]-1))) {
700                temp << endl;
701              }
702                }
703          }
704        }
705          }
706          break;
707      default:      default:
708        stringstream mess;        stringstream mess;
709        mess << "Error - (toString) Invalid rank: " << getRank();        mess << "Error - (toString) Invalid rank: " << getRank();
# Line 462  DataArrayView::shapeToString(const DataA Line 717  DataArrayView::shapeToString(const DataA
717  {  {
718      stringstream temp;      stringstream temp;
719      temp << "(";      temp << "(";
720      for (int i=0;i<shape.size();++i) {      for (int i=0;i<shape.size();i++) {
721        temp << shape[i];        temp << shape[i];
722        if (shape.size()-i != 1) {        if (i < shape.size()-1) {
723      temp << ",";      temp << ",";
724        }        }
725      }      }
# Line 477  DataArrayView::matMult(const DataArrayVi Line 732  DataArrayView::matMult(const DataArrayVi
732                         const DataArrayView& right,                         const DataArrayView& right,
733                         DataArrayView& result)                         DataArrayView& result)
734  {  {
735    
736      if (left.getRank()==0 || right.getRank()==0) {      if (left.getRank()==0 || right.getRank()==0) {
737        stringstream temp;        stringstream temp;
738        temp << "Error - (matMult) Invalid for rank 0 objects.";        temp << "Error - (matMult) Invalid for rank 0 objects.";
739        throw DataException(temp.str());        throw DataException(temp.str());
740      }      }
741    
742      if (left.getShape()[left.getRank()-1]!=right.getShape()[0]) {      if (left.getShape()[left.getRank()-1] != right.getShape()[0]) {
743        stringstream temp;        stringstream temp;
744        temp << "Error - (matMult) Dimension: " << left.getRank()        temp << "Error - (matMult) Dimension: " << left.getRank()
745         << ", size: " << left.getShape()[left.getRank()-1]         << ", size: " << left.getShape()[left.getRank()-1]
# Line 491  DataArrayView::matMult(const DataArrayVi Line 747  DataArrayView::matMult(const DataArrayVi
747         << " of RHS don't match.";         << " of RHS don't match.";
748        throw DataException(temp.str());        throw DataException(temp.str());
749      }      }
750      int outputRank=left.getRank()+right.getRank()-2;  
751            int outputRank = left.getRank()+right.getRank()-2;
752    
753      if (outputRank < 0) {      if (outputRank < 0) {
754        stringstream temp;        stringstream temp;
755        temp << "Error - (matMult) Left and right operands cannot be multiplied "        temp << "Error - (matMult) LHS and RHS cannot be multiplied "
756         << "as they have incompatable rank.";         << "as they have incompatable rank.";
757        throw DataException(temp.str());        throw DataException(temp.str());
758      }      }
# Line 508  DataArrayView::matMult(const DataArrayVi Line 765  DataArrayView::matMult(const DataArrayVi
765        throw DataException(temp.str());        throw DataException(temp.str());
766      }      }
767    
768      for (int i=0; i<(left.getRank()-1); ++i) {      for (int i=0; i<(left.getRank()-1); i++) {
769        if (left.getShape()[i]!=result.getShape()[i]) {        if (left.getShape()[i] != result.getShape()[i]) {
770      stringstream temp;      stringstream temp;
771      temp << "Error - (matMult) Dimension: " << i      temp << "Error - (matMult) Dimension: " << i
772           << " of LHS and output array don't match.";           << " of LHS and result array don't match.";
773      throw DataException(temp.str());      throw DataException(temp.str());
774        }        }
775      }      }
776    
777      for (int i=1; i<right.getRank(); ++i) {      for (int i=1; i<right.getRank(); i++) {
778        if (right.getShape()[i] != result.getShape()[i+left.getRank()-2]) {        if (right.getShape()[i] != result.getShape()[i+left.getRank()-2]) {
779      stringstream temp;      stringstream temp;
780      temp << "Error - (matMult) Dimension: " << i      temp << "Error - (matMult) Dimension: " << i
781           << ", size: " << right.getShape()[i]           << ", size: " << right.getShape()[i]
782           << " of RHS and dimension: " << i+left.getRank()-1           << " of RHS and dimension: " << i+left.getRank()-1
783           << ", size: " << result.getShape()[i+left.getRank()-1]           << ", size: " << result.getShape()[i+left.getRank()-1]
784           << " of output array don't match.";           << " of result array don't match.";
785      throw DataException(temp.str());      throw DataException(temp.str());
786        }        }
787      }      }
788    
789      switch (left.getRank()) {      switch (left.getRank()) {
790      case 1:  
       switch (right.getRank()) {  
791        case 1:        case 1:
792      result()=0;          switch (right.getRank()) {
793      for (int i=0;i<left.getShape()[0];++i) {            case 1:
794        result()+=left(i)*right(i);              result()=0;
795      }              for (int i=0;i<left.getShape()[0];i++) {
796      break;                result()+=left(i)*right(i);
797        case 2:              }
798          for (int i=0;i<result.getShape()[0];++i) {              break;
799        result(i)=0;            case 2:
800        for (int j=0;j<right.getShape()[0];++j) {              for (int i=0;i<result.getShape()[0];i++) {
801          result(i)+=left(j)*right(j,i);                result(i)=0;
802        }                for (int j=0;j<right.getShape()[0];j++) {
803      }                  result(i)+=left(j)*right(j,i);
804                  }
805                }
806                break;
807              default:
808                stringstream temp; temp << "Error - (matMult) Invalid rank. Programming error.";
809                throw DataException(temp.str());
810                break;
811            }
812          break;          break;
813        default:  
     stringstream temp;  
     temp << "Error - (matMult) Invalid rank. Programming error.";  
         throw DataException(temp.str());  
         break;  
       }  
       break;  
     case 2:  
       switch (right.getRank()) {  
       case 1:  
     result()=0;  
     for (int i=0;i<left.getShape()[0];++i) {  
       result(i)=0;  
       for (int j=0;j<left.getShape()[1];++j) {  
         result(i)+=left(i,j)*right(i);  
       }  
     }  
     break;  
814        case 2:        case 2:
815          for (int i=0;i<result.getShape()[0];++i) {          switch (right.getRank()) {
816        for (int j=0;j<result.getShape()[1];++j) {            case 1:
817          result(i,j)=0;              result()=0;
818          for (int jR=0;jR<right.getShape()[0];++jR) {                for (int i=0;i<left.getShape()[0];i++) {
819            result(i,j)+=left(i,jR)*right(jR,j);                 result(i)=0;
820          }                 for (int j=0;j<left.getShape()[1];j++) {
821        }                   result(i)+=left(i,j)*right(i);
822      }                 }
823                }
824            break;
825              case 2:
826                for (int i=0;i<result.getShape()[0];i++) {
827                  for (int j=0;j<result.getShape()[1];j++) {
828                    result(i,j)=0;
829                    for (int jR=0;jR<right.getShape()[0];jR++) {
830                      result(i,j)+=left(i,jR)*right(jR,j);
831                    }
832                  }
833                }
834                break;
835              default:
836                stringstream temp; temp << "Error - (matMult) Invalid rank. Programming error.";
837                throw DataException(temp.str());
838                break;
839            }
840          break;          break;
841    
842        default:        default:
843      stringstream temp;          stringstream temp; temp << "Error - (matMult) Not supported for rank: " << left.getRank();
     temp << "Error - (matMult) Invalid rank. Programming error.";  
844          throw DataException(temp.str());          throw DataException(temp.str());
845          break;          break;
       }  
       break;  
     default:  
       stringstream temp;  
       temp << "Error - (matMult) Not supported for rank: " << left.getRank();  
       throw DataException(temp.str());  
       break;  
846      }      }
847    
848  }  }
849    
850  DataArrayView::RegionType  DataArrayView::ShapeType
851  DataArrayView::getSliceRegion(const boost::python::object& key) const  DataArrayView::determineResultShape(const DataArrayView& left,
852                                        const DataArrayView& right)
853  {  {
854      int slice_len=0,i=0,rank=getRank(),out_len=0;      ShapeType temp;
855      extract<tuple> key_tuple(key);      for (int i=0; i<(left.getRank()-1); i++) {
856      if (key_tuple.check()) {        temp.push_back(left.getShape()[i]);
         slice_len=extract<int>(key.attr("__len__")());  
     } else {  
         slice_len=1;  
     }  
     if (slice_len>=rank) {  
          out_len=slice_len;  
     } else {  
          out_len=rank;  
857      }      }
858      DataArrayView::RegionType out(out_len);      for (int i=1; i<right.getRank(); i++) {
859      if (key_tuple.check()) {        temp.push_back(right.getShape()[i]);
          for (i=0;i<slice_len;i++) out[i]=getSliceRange(getShape()[i],key[i]);  
     } else {  
          out[0]=getSliceRange(getShape()[0],key);  
860      }      }
861      for (i=slice_len;i<rank;i++) out[i]=std::pair<int,int>(0,getShape()[i]);      return temp;
     // throw DataException("Error - number of slices is bigger then the rank.");  
     return out;  
 }  
   
 std::pair<int,int>  
 getSliceRange(const int s,  
               const boost::python::object& key)  
 {  
    int s1=s;  
    int s0=0;  
    /* if the key is an int we set start=end=key */  
    extract<int> slice_int(key);  
    if (slice_int.check()) {  
          s0=slice_int();  
          s1=s0;  
    } else {  
        extract<int> step(key.attr("step"));  
        if (step.check() && step()!=1) {  
             throw DataException("Error - Data does not support increments in slicing ");  
        } else {  
           extract<int> start(key.attr("start"));  
           if (start.check()) s0=start();  
           extract<int> stop(key.attr("stop"));  
           if (stop.check()) s1=stop();  
        }  
    }  
    return std::pair<int,int>(s0,s1);  
862  }  }
863    
864  bool operator==(const DataArrayView& left, const DataArrayView& right)  bool operator==(const DataArrayView& left, const DataArrayView& right)
# Line 647  bool operator==(const DataArrayView& lef Line 868  bool operator==(const DataArrayView& lef
868      // of the view is the same.      // of the view is the same.
869      bool result=(left.m_shape==right.m_shape);      bool result=(left.m_shape==right.m_shape);
870      if (result) {      if (result) {
871        for (int i=0;i<left.noValues();++i) {        for (int i=0;i<left.noValues();i++) {
872      result=result && (left.getData(i)==right.getData(i));      result=result && (left.getData(i)==right.getData(i));
873        }        }
874      }      }

Legend:
Removed from v.100  
changed lines
  Added in v.1388

  ViewVC Help
Powered by ViewVC 1.1.26