/[escript]/branches/subworld2/escriptcore/src/DataTypes.cpp
ViewVC logotype

Contents of /branches/subworld2/escriptcore/src/DataTypes.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5504 - (show annotations)
Wed Mar 4 22:58:13 2015 UTC (4 years, 1 month ago) by jfenwick
File size: 21279 byte(s)
Again with a more up to date copy


1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2015 by University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development 2012-2013 by School of Earth Sciences
13 * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 *
15 *****************************************************************************/
16
17 #define ESNEEDPYTHON
18 #include "esysUtils/first.h"
19
20
21 #include "DataTypes.h"
22 #include <sstream>
23 #include <boost/python/extract.hpp>
24 #include <boost/python/tuple.hpp>
25 #include "DataException.h"
26
27 namespace {
28 using namespace boost::python;
29 using namespace escript;
30 using namespace escript::DataTypes;
31
32 /*
33 \brief
34 Calculate the slice range from the given python key object
35 Used by getSliceRegion - since it is not used anywhere else, I have elected not to declare it
36 in the header.
37 Returns the python slice object key as a pair of ints where the first
38 member is the start and the second member is the end. the presence of a possible
39 step attribute with value different from one will throw an exception
40
41 /param key - Input - key object specifying slice range.
42 */
43 std::pair<int,int>
44 getSliceRange(const boost::python::object& key,
45 const int shape)
46 {
47 /* default slice range is range of entire shape dimension */
48 int s0=0, s1=shape;;
49 extract<int> slice_int(key);
50 if (slice_int.check()) {
51 /* if the key is a single int set start=key and end=key */
52 /* in this case, we want to return a rank-1 dimension object from
53 this object, taken from a single index value for one of this
54 object's dimensions */
55 s0=slice_int();
56 s1=s0;
57 } else {
58 /* if key is a pair extract begin and end values */
59 extract<int> step(key.attr("step"));
60 if (step.check() && step()!=1) {
61 throw DataException("Error - Data does not support increments in slicing ");
62 } else {
63 extract<int> start(key.attr("start"));
64 if (start.check()) {
65 s0=start();
66 }
67 extract<int> stop(key.attr("stop"));
68 if (stop.check()) {
69 s1=stop();
70 }
71 }
72 }
73 if (s0 < 0)
74 throw DataException("Error - slice index out of range.");
75 if (s0 == s1 && s1 >= shape)
76 throw DataException("Error - slice index out of range.");
77 if (s0 != s1 && s1>shape)
78 throw DataException("Error - slice index out of range.");
79 if (s0 > s1)
80 throw DataException("Error - lower index must less or equal upper index.");
81 return std::pair<int,int>(s0,s1);
82 }
83 }
84
85
86 using namespace boost::python;
87
88 namespace escript
89 {
90 namespace DataTypes
91 {
92
93 int
94 noValues(const ShapeType& shape)
95 {
96 ShapeType::const_iterator i;
97 //
98 // An empty shape vector means rank 0 which contains 1 value
99 int noValues=1;
100 for (i=shape.begin();i!=shape.end();i++) {
101 noValues*=(*i);
102 }
103 return noValues;
104 }
105
106 int
107 noValues(const RegionLoopRangeType& region)
108 {
109 //
110 // An empty region vector means rank 0 which contains 1 value
111 int noValues=1;
112 unsigned int i;
113 for (i=0;i<region.size();i++) {
114 noValues*=region[i].second-region[i].first;
115 }
116 return noValues;
117 }
118
119 std::string
120 shapeToString(const DataTypes::ShapeType& shape)
121 {
122 std::stringstream temp;
123 temp << "(";
124 unsigned int i;
125 for (i=0;i<shape.size();i++) {
126 temp << shape[i];
127 if (i < shape.size()-1) {
128 temp << ",";
129 }
130 }
131 temp << ")";
132 return temp.str();
133 }
134
135
136
137
138
139 DataTypes::RegionType
140 getSliceRegion(const DataTypes::ShapeType& shape, const boost::python::object& key)
141 {
142 int slice_rank, i;
143 int this_rank=shape.size();
144 RegionType out(this_rank);
145 /* allow for case where key is singular eg: [1], this implies we
146 want to generate a rank-1 dimension object, as opposed to eg: [1,2]
147 which implies we want to take a rank dimensional object with one
148 dimension of size 1 */
149 extract<tuple> key_tuple(key);
150 if (key_tuple.check()) {
151 slice_rank=extract<int> (key.attr("__len__")());
152 /* ensure slice is correctly dimensioned */
153 if (slice_rank>this_rank) {
154 throw DataException("Error - rank of slices does not match rank of slicee");
155 } else {
156 /* calculate values for slice range */
157 for (i=0;i<slice_rank;i++) {
158 out[i]=getSliceRange(key[i],shape[i]);
159 }
160 }
161 } else {
162 slice_rank=1;
163 if (slice_rank>this_rank) {
164 throw DataException("Error - rank of slices does not match rank of slicee");
165 } else {
166 out[0]=getSliceRange(key,shape[0]);
167 }
168 }
169 for (i=slice_rank;i<this_rank;i++) {
170 out[i]=std::pair<int,int>(0,shape[i]);
171 }
172 return out;
173 }
174
175 DataTypes::ShapeType
176 getResultSliceShape(const RegionType& region)
177 {
178 int dimSize;
179 ShapeType result;
180 RegionType::const_iterator i;
181 for (i=region.begin();i!=region.end();i++) {
182 dimSize=((i->second)-(i->first));
183 if (dimSize!=0) {
184 result.push_back(dimSize);
185 }
186 }
187 return result;
188 }
189
190 DataTypes::RegionLoopRangeType
191 getSliceRegionLoopRange(const DataTypes::RegionType& region)
192 {
193 DataTypes::RegionLoopRangeType region_loop_range(region.size());
194 unsigned int i;
195 for (i=0;i<region.size();i++) {
196 if (region[i].first==region[i].second) {
197 region_loop_range[i].first=region[i].first;
198 region_loop_range[i].second=region[i].second+1;
199 } else {
200 region_loop_range[i].first=region[i].first;
201 region_loop_range[i].second=region[i].second;
202 }
203 }
204 return region_loop_range;
205 }
206
207
208 std::string
209 createShapeErrorMessage(const std::string& messagePrefix,
210 const DataTypes::ShapeType& other,
211 const DataTypes::ShapeType& thisShape)
212 {
213 std::stringstream temp;
214 temp << messagePrefix
215 << " This shape: " << shapeToString(thisShape)
216 << " Other shape: " << shapeToString(other);
217 return temp.str();
218 }
219
220
221 // Additional slice operations
222
223 inline
224 bool
225 checkOffset(ValueType::size_type offset, int size, int noval)
226 {
227 return (size >= (offset+noval));
228 }
229
230
231 void
232 copySlice(ValueType& left,
233 const ShapeType& leftShape,
234 ValueType::size_type thisOffset,
235 const ValueType& other,
236 const ShapeType& otherShape,
237 ValueType::size_type otherOffset,
238 const RegionLoopRangeType& region)
239 {
240 //
241 // Make sure views are not empty
242
243 EsysAssert(!left.size()==0,
244 "Error - left data is empty.");
245 EsysAssert(!other.size()==0,
246 "Error - other data is empty.");
247
248 //
249 // Check the view to be sliced from is compatible with the region to be sliced,
250 // and that the region to be sliced is compatible with this view:
251 EsysAssert(checkOffset(thisOffset,left.size(),noValues(leftShape)),
252 "Error - offset incompatible with this view.");
253 EsysAssert(otherOffset+noValues(leftShape)<=other.size(),
254 "Error - offset incompatible with other view.");
255
256 EsysAssert(getRank(otherShape)==region.size(),
257 "Error - slice not same rank as view to be sliced from.");
258
259 EsysAssert(noValues(leftShape)==noValues(getResultSliceShape(region)),
260 "Error - slice shape not compatible shape for this view.");
261
262 //
263 // copy the values in the specified region of the other view into this view
264
265 // the following loops cannot be parallelised due to the numCopy counter
266 int numCopy=0;
267
268 switch (region.size()) {
269 case 0:
270 /* this case should never be encountered,
271 as python will never pass us an empty region.
272 here for completeness only, allows slicing of a scalar */
273 // (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex()];
274
275 left[thisOffset+numCopy]=other[otherOffset];
276 numCopy++;
277 break;
278 case 1:
279 for (int i=region[0].first;i<region[0].second;i++) {
280 left[thisOffset+numCopy]=other[otherOffset+getRelIndex(otherShape,i)];
281 numCopy++;
282 }
283 break;
284 case 2:
285 for (int j=region[1].first;j<region[1].second;j++) {
286 for (int i=region[0].first;i<region[0].second;i++) {
287 /* (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i,j)];*/
288 left[thisOffset+numCopy]=other[otherOffset+getRelIndex(otherShape,i,j)];
289 numCopy++;
290 }
291 }
292 break;
293 case 3:
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)];
298 left[thisOffset+numCopy]=other[otherOffset+getRelIndex(otherShape,i,j,k)];
299 numCopy++;
300 }
301 }
302 }
303 break;
304 case 4:
305 for (int l=region[3].first;l<region[3].second;l++) {
306 for (int k=region[2].first;k<region[2].second;k++) {
307 for (int j=region[1].first;j<region[1].second;j++) {
308 for (int i=region[0].first;i<region[0].second;i++) {
309 /* (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i,j,k,l)];*/
310 left[thisOffset+numCopy]=other[otherOffset+getRelIndex(otherShape,i,j,k,l)];
311 numCopy++;
312 }
313 }
314 }
315 }
316 break;
317 default:
318 std::stringstream mess;
319 mess << "Error - (copySlice) Invalid slice region rank: " << region.size();
320 throw DataException(mess.str());
321 }
322 }
323
324
325 void
326 copySliceFrom(ValueType& left,
327 const ShapeType& leftShape,
328 ValueType::size_type thisOffset,
329 const ValueType& other,
330 const ShapeType& otherShape,
331 ValueType::size_type otherOffset,
332 const RegionLoopRangeType& region)
333 {
334 //
335 // Make sure views are not empty
336
337 EsysAssert(left.size()!=0,
338 "Error - this view is empty.");
339 EsysAssert(other.size()!=0,
340 "Error - other view is empty.");
341
342 //
343 // Check this view is compatible with the region to be sliced,
344 // and that the region to be sliced is compatible with the other view:
345
346 EsysAssert(checkOffset(otherOffset,other.size(),noValues(otherShape)),
347 "Error - offset incompatible with other view.");
348 EsysAssert(thisOffset+noValues(otherShape)<=left.size(),
349 "Error - offset incompatible with this view.");
350
351 EsysAssert(getRank(leftShape)==region.size(),
352 "Error - slice not same rank as this view.");
353
354 EsysAssert(getRank(otherShape)==0 || noValues(otherShape)==noValues(getResultSliceShape(region)),
355 "Error - slice shape not compatible shape for other view.");
356
357 //
358 // copy the values in the other view into the specified region of this view
359
360 // allow for case where other view is a scalar
361 if (getRank(otherShape)==0) {
362
363 // the following loops cannot be parallelised due to the numCopy counter
364 int numCopy=0;
365
366 switch (region.size()) {
367 case 0:
368 /* this case should never be encountered,
369 as python will never pass us an empty region.
370 here for completeness only, allows slicing of a scalar */
371 //(*m_data)[thisOffset+relIndex()]=(*other.m_data)[otherOffset];
372 left[thisOffset]=other[otherOffset];
373 numCopy++;
374 break;
375 case 1:
376 for (int i=region[0].first;i<region[0].second;i++) {
377 left[thisOffset+getRelIndex(leftShape,i)]=other[otherOffset];
378 numCopy++;
379 }
380 break;
381 case 2:
382 for (int j=region[1].first;j<region[1].second;j++) {
383 for (int i=region[0].first;i<region[0].second;i++) {
384 left[thisOffset+getRelIndex(leftShape,i,j)]=other[otherOffset];
385 numCopy++;
386 }
387 }
388 break;
389 case 3:
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)]=other[otherOffset];
394 numCopy++;
395 }
396 }
397 }
398 break;
399 case 4:
400 for (int l=region[3].first;l<region[3].second;l++) {
401 for (int k=region[2].first;k<region[2].second;k++) {
402 for (int j=region[1].first;j<region[1].second;j++) {
403 for (int i=region[0].first;i<region[0].second;i++) {
404 left[thisOffset+getRelIndex(leftShape,i,j,k,l)]=other[otherOffset];
405 numCopy++;
406 }
407 }
408 }
409 }
410 break;
411 default:
412 std::stringstream mess;
413 mess << "Error - (copySliceFrom) Invalid slice region rank: " << region.size();
414 throw DataException(mess.str());
415 }
416
417 } else {
418
419 // the following loops cannot be parallelised due to the numCopy counter
420 int numCopy=0;
421
422 switch (region.size()) {
423 case 0:
424 /* this case should never be encountered,
425 as python will never pass us an empty region.
426 here for completeness only, allows slicing of a scalar */
427 //(*m_data)[thisOffset+relIndex()]=(*other.m_data)[otherOffset+numCopy];
428 left[thisOffset]=other[otherOffset+numCopy];
429 numCopy++;
430 break;
431 case 1:
432 for (int i=region[0].first;i<region[0].second;i++) {
433 left[thisOffset+getRelIndex(leftShape,i)]=other[otherOffset+numCopy];
434 numCopy++;
435 }
436 break;
437 case 2:
438 for (int j=region[1].first;j<region[1].second;j++) {
439 for (int i=region[0].first;i<region[0].second;i++) {
440 left[thisOffset+getRelIndex(leftShape,i,j)]=other[otherOffset+numCopy];
441 numCopy++;
442 }
443 }
444 break;
445 case 3:
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)]=other[otherOffset+numCopy];
450 numCopy++;
451 }
452 }
453 }
454 break;
455 case 4:
456 for (int l=region[3].first;l<region[3].second;l++) {
457 for (int k=region[2].first;k<region[2].second;k++) {
458 for (int j=region[1].first;j<region[1].second;j++) {
459 for (int i=region[0].first;i<region[0].second;i++) {
460 left[thisOffset+getRelIndex(leftShape,i,j,k,l)]=other[otherOffset+numCopy];
461 numCopy++;
462 }
463 }
464 }
465 }
466 break;
467 default:
468 std::stringstream mess;
469 mess << "Error - (copySliceFrom) Invalid slice region rank: " << region.size();
470 throw DataException(mess.str());
471 }
472
473 }
474
475 }
476
477
478 void
479 pointToStream(std::ostream& os, const ValueType::ElementType* data,const ShapeType& shape, int offset, bool needsep, const std::string& sep)
480 {
481 using namespace std;
482 EsysAssert(data!=0, "Error - data is null");
483 // EsysAssert(data.size()>0,"Error - Data object is empty.");
484 switch (getRank(shape)) {
485 case 0:
486 if (needsep)
487 {
488 os << sep;
489 }
490 else
491 {
492 needsep=true;
493 }
494 os << data[offset];
495 break;
496 case 1:
497 for (int i=0;i<shape[0];i++) {
498 if (needsep)
499 {
500 os << sep;
501 }
502 else
503 {
504 needsep=true;
505 }
506 os << data[i+offset];
507 }
508 break;
509 case 2:
510 for (int i=0;i<shape[0];i++) {
511 for (int j=0;j<shape[1];j++) {
512 if (needsep)
513 {
514 os << sep;
515 }
516 else
517 {
518 needsep=true;
519 }
520 os << data[offset+getRelIndex(shape,i,j)];
521 }
522 }
523 break;
524 case 3:
525 for (int i=0;i<shape[0];i++) {
526 for (int j=0;j<shape[1];j++) {
527 for (int k=0;k<shape[2];k++) {
528 if (needsep)
529 {
530 os << sep;
531 }
532 else
533 {
534 needsep=true;
535 }
536 os << data[offset+getRelIndex(shape,i,j,k)];
537 }
538 }
539 }
540 break;
541 case 4:
542 for (int i=0;i<shape[0];i++) {
543 for (int j=0;j<shape[1];j++) {
544 for (int k=0;k<shape[2];k++) {
545 for (int l=0;l<shape[3];l++) {
546 if (needsep)
547 {
548 os << sep;
549 }
550 else
551 {
552 needsep=true;
553 }
554 os << data[offset+getRelIndex(shape,i,j,k,l)];
555 }
556 }
557 }
558 }
559 break;
560 default:
561 stringstream mess;
562 mess << "Error - (pointToStream) Invalid rank: " << getRank(shape);
563 throw DataException(mess.str());
564 }
565 }
566
567
568 std::string
569 pointToString(const ValueType& data,const ShapeType& shape, int offset, const std::string& prefix)
570 {
571 using namespace std;
572 EsysAssert(data.size()>0,"Error - Data object is empty.");
573 stringstream temp;
574 string finalPrefix=prefix;
575 if (prefix.length() > 0) {
576 finalPrefix+=" ";
577 }
578 switch (getRank(shape)) {
579 case 0:
580 temp << finalPrefix << data[offset];
581 break;
582 case 1:
583 for (int i=0;i<shape[0];i++) {
584 temp << finalPrefix << "(" << i << ") " << data[i+offset];
585 if (i!=(shape[0]-1)) {
586 temp << endl;
587 }
588 }
589 break;
590 case 2:
591 for (int i=0;i<shape[0];i++) {
592 for (int j=0;j<shape[1];j++) {
593 temp << finalPrefix << "(" << i << "," << j << ") " << data[offset+getRelIndex(shape,i,j)];
594 if (!(i==(shape[0]-1) && j==(shape[1]-1))) {
595 temp << endl;
596 }
597 }
598 }
599 break;
600 case 3:
601 for (int i=0;i<shape[0];i++) {
602 for (int j=0;j<shape[1];j++) {
603 for (int k=0;k<shape[2];k++) {
604 temp << finalPrefix << "(" << i << "," << j << "," << k << ") " << data[offset+getRelIndex(shape,i,j,k)];
605 if (!(i==(shape[0]-1) && j==(shape[1]-1) && k==(shape[2]-1))) {
606 temp << endl;
607 }
608 }
609 }
610 }
611 break;
612 case 4:
613 for (int i=0;i<shape[0];i++) {
614 for (int j=0;j<shape[1];j++) {
615 for (int k=0;k<shape[2];k++) {
616 for (int l=0;l<shape[3];l++) {
617 temp << finalPrefix << "(" << i << "," << j << "," << k << "," << l << ") " << data[offset+getRelIndex(shape,i,j,k,l)];
618 if (!(i==(shape[0]-1) && j==(shape[1]-1) && k==(shape[2]-1) && l==(shape[3]-1))) {
619 temp << endl;
620 }
621 }
622 }
623 }
624 }
625 break;
626 default:
627 stringstream mess;
628 mess << "Error - (toString) Invalid rank: " << getRank(shape);
629 throw DataException(mess.str());
630 }
631 return temp.str();
632 }
633
634
635 void copyPoint(ValueType& dest, ValueType::size_type doffset, ValueType::size_type nvals, const ValueType& src, ValueType::size_type soffset)
636 {
637 EsysAssert((dest.size()>0&&src.size()>0&&checkOffset(doffset,dest.size(),nvals)),
638 "Error - Couldn't copy due to insufficient storage.");
639 // EsysAssert((checkShape(other.getShape())),
640 // createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",other.getShape(),m_shape));
641 if (checkOffset(doffset,dest.size(),nvals) && checkOffset(soffset,src.size(),nvals)) {
642 memcpy(&dest[doffset],&src[soffset],sizeof(double)*nvals);
643 } else {
644 throw DataException("Error - invalid offset specified.");
645 }
646
647
648
649 }
650
651 } // end namespace DataTypes
652 } // end namespace escript

  ViewVC Help
Powered by ViewVC 1.1.26