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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26