/[escript]/trunk/escriptcore/src/Data.h
ViewVC logotype

Contents of /trunk/escriptcore/src/Data.h

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6424 - (show annotations)
Thu Nov 3 02:13:09 2016 UTC (2 years, 5 months ago) by jfenwick
File MIME type: text/plain
File size: 59666 byte(s)
Fixing copyWithMask to work with complex Data (issue #389)
1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2016 by The University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Apache License, version 2.0
9 * http://www.apache.org/licenses/LICENSE-2.0
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 /** \file Data.h */
18
19 #ifndef __ESCRIPT_DATA_H__
20 #define __ESCRIPT_DATA_H__
21
22 #include "system_dep.h"
23 #include "DataAbstract.h"
24 #include "DataException.h"
25 #include "DataTypes.h"
26 #include "EsysMPI.h"
27 #include "FunctionSpace.h"
28 #include "DataVectorOps.h"
29 #include <algorithm>
30 #include <string>
31 #include <sstream>
32
33 #include <boost/python/object.hpp>
34 #include <boost/python/tuple.hpp>
35 #include <boost/math/special_functions/bessel.hpp>
36
37 #ifndef ESCRIPT_MAX_DATA_RANK
38 #define ESCRIPT_MAX_DATA_RANK 4
39 #endif
40
41 namespace escript {
42
43 //
44 // Forward declaration for various implementations of Data.
45 class DataConstant;
46 class DataTagged;
47 class DataExpanded;
48 class DataLazy;
49
50 /**
51 \brief
52 Data represents a collection of datapoints.
53
54 Description:
55 Internally, the datapoints are actually stored by a DataAbstract object.
56 The specific instance of DataAbstract used may vary over the lifetime
57 of the Data object.
58 Some methods on this class return references (eg getShape()).
59 These references should not be used after an operation which changes the underlying DataAbstract object.
60 Doing so will lead to invalid memory access.
61 This should not affect any methods exposed via boost::python.
62 */
63 class Data {
64
65 public:
66
67 /**
68 Constructors.
69 */
70
71 /**
72 \brief
73 Default constructor.
74 Creates a DataEmpty object.
75 */
76 Data();
77
78 /**
79 \brief
80 Copy constructor.
81 WARNING: Only performs a shallow copy.
82 */
83 Data(const Data& inData);
84
85 /**
86 \brief
87 Constructor from another Data object. If "what" is different from the
88 function space of inData the inData are tried to be interpolated to what,
89 otherwise a shallow copy of inData is returned.
90 */
91 Data(const Data& inData,
92 const FunctionSpace& what);
93
94 /**
95 \brief Copy Data from an existing vector
96 */
97 Data(const DataTypes::RealVectorType& value,
98 const DataTypes::ShapeType& shape,
99 const FunctionSpace& what,
100 bool expanded);
101
102 /**
103 \brief
104 Constructor which creates a Data with points having the specified shape.
105
106 \param value - Input - Single real value applied to all Data.
107 \param dataPointShape - Input - The shape of each data point.
108 \param what - Input - A description of what this data represents.
109 \param expanded - Input - Flag, if true fill the entire container with
110 the given value. Otherwise a more efficient storage
111 mechanism will be used.
112 */
113 Data(DataTypes::real_t value,
114 const DataTypes::ShapeType& dataPointShape,
115 const FunctionSpace& what,
116 bool expanded);
117
118 /**
119 \brief
120 Constructor which creates a Data with points having the specified shape.
121
122 \param value - Input - Single complex value applied to all Data.
123 \param dataPointShape - Input - The shape of each data point.
124 \param what - Input - A description of what this data represents.
125 \param expanded - Input - Flag, if true fill the entire container with
126 the given value. Otherwise a more efficient storage
127 mechanism will be used.
128 */
129 explicit
130 Data(DataTypes::cplx_t value,
131 const DataTypes::ShapeType& dataPointShape,
132 const FunctionSpace& what,
133 bool expanded);
134
135 /**
136 \brief
137 Constructor which performs a deep copy of a region from another Data object.
138
139 \param inData - Input - Input Data object.
140 \param region - Input - Region to copy.
141 */
142 Data(const Data& inData,
143 const DataTypes::RegionType& region);
144
145
146 /**
147 \brief
148 Constructor which copies data from a wrapped array.
149
150 \param w - Input - Input data.
151 \param what - Input - A description of what this data represents.
152 \param expanded - Input - Flag, if true fill the entire container with
153 the value. Otherwise a more efficient storage
154 mechanism will be used.
155 */
156 Data(const WrappedArray& w, const FunctionSpace& what,
157 bool expanded);
158
159
160 /**
161 \brief
162 Constructor which creates a DataConstant.
163 Copies data from any object that can be treated like a python array/sequence.
164 All other parameters are copied from other.
165
166 \param value - Input - Input data.
167 \param other - Input - contains all other parameters.
168 */
169 Data(const boost::python::object& value,
170 const Data& other);
171
172 /**
173 This constructor subsumes a number of previous python ones.
174
175 Data(const boost::python::object& value,
176 const FunctionSpace& what=FunctionSpace(),
177 bool expanded=false);
178
179 Data(DataTypes::real_t value,
180 const boost::python::tuple& shape=boost::python::make_tuple(),
181 const FunctionSpace& what=FunctionSpace(),
182 bool expanded=false);
183
184 and a new
185
186 Data(cplx_t value,
187 const boost::python::tuple& shape=boost::python::make_tuple(),
188 const FunctionSpace& what=FunctionSpace(),
189 bool expanded=false);
190
191 */
192 Data(boost::python::object value,
193 boost::python::object par1=boost::python::object(),
194 boost::python::object par2=boost::python::object(),
195 boost::python::object par3=boost::python::object());
196
197
198
199
200
201 /**
202 \brief Create a Data using an existing DataAbstract. Warning: The new object assumes ownership of the pointer!
203 Once you have passed the pointer, do not delete it.
204 */
205 explicit Data(DataAbstract* underlyingdata);
206
207 /**
208 \brief Create a Data based on the supplied DataAbstract
209 */
210 explicit Data(DataAbstract_ptr underlyingdata);
211
212 /**
213 \brief
214 Destructor
215 */
216 ~Data();
217
218 /**
219 \brief Make this object a deep copy of "other".
220 */
221 void
222 copy(const Data& other);
223
224 /**
225 \brief Return a pointer to a deep copy of this object.
226 */
227 Data
228 copySelf() const;
229
230
231 /**
232 \brief produce a delayed evaluation version of this Data.
233 */
234 Data
235 delay();
236
237 /**
238 \brief convert the current data into lazy data.
239 */
240 void
241 delaySelf();
242
243
244 /**
245 Member access methods.
246 */
247
248 /**
249 \brief
250 switches on update protection
251
252 */
253 void
254 setProtection();
255
256 /**
257 \brief
258 Returns true, if the data object is protected against update
259
260 */
261 bool
262 isProtected() const;
263
264
265 /**
266 \brief
267 Return the value of a data point as a python tuple.
268 */
269 const boost::python::object
270 getValueOfDataPointAsTuple(int dataPointNo);
271
272 /**
273 \brief
274 sets the values of a data-point from a python object on this process
275 */
276 void
277 setValueOfDataPointToPyObject(int dataPointNo, const boost::python::object& py_object);
278
279 /**
280 \brief
281 sets the values of a data-point from a array-like object on this process
282 */
283 void
284 setValueOfDataPointToArray(int dataPointNo, const boost::python::object&);
285
286 /**
287 \brief
288 sets the values of a data-point on this process
289 */
290 void
291 setValueOfDataPoint(int dataPointNo, const DataTypes::real_t);
292
293 void
294 setValueOfDataPointC(int dataPointNo, const DataTypes::cplx_t);
295
296
297 /**
298 \brief Return a data point across all processors as a python tuple.
299 */
300 const boost::python::object
301 getValueOfGlobalDataPointAsTuple(int procNo, int dataPointNo);
302
303
304 /**
305 \brief Set the value of a global data point
306 */
307 void
308 setTupleForGlobalDataPoint(int id, int proc, boost::python::object);
309
310 /**
311 \brief
312 Return the tag number associated with the given data-point.
313
314 */
315 int
316 getTagNumber(int dpno);
317
318
319 /**
320 \brief
321 Write the data as a string. For large amounts of data, a summary is printed.
322 */
323 std::string
324 toString() const;
325
326 /**
327 \brief
328 Whatever the current Data type make this into a DataExpanded.
329 */
330 void
331 expand();
332
333 /**
334 \brief
335 If possible convert this Data to DataTagged. This will only allow
336 Constant data to be converted to tagged. An attempt to convert
337 Expanded data to tagged will throw an exception.
338 */
339 void
340 tag();
341
342 /**
343 \brief If this data is lazy, then convert it to ready data.
344 What type of ready data depends on the expression. For example, Constant+Tagged==Tagged.
345 */
346 void
347 resolve();
348
349 /**
350 \brief returns return true if data contains NaN.
351 \warning This is dependent on the ability to reliably detect NaNs on your compiler.
352 See the nancheck function in LocalOps for details.
353 */
354 bool
355 hasNaN();
356
357 /**
358 \brief replaces all NaN values with value
359 */
360 void
361 replaceNaN(DataTypes::real_t value);
362
363 /**
364 \brief replaces all NaN values with value
365 */
366 void
367 replaceNaN(DataTypes::cplx_t value);
368
369 /**
370 \brief replaces all NaN values with value
371 */
372 void
373 replaceNaNPython(boost::python::object obj);
374
375
376
377
378 /**
379 \brief Ensures data is ready for write access.
380 This means that the data will be resolved if lazy and will be copied if shared with another Data object.
381 \warning This method should only be called in single threaded sections of code. (It modifies m_data).
382 Do not create any Data objects from this one between calling requireWrite and getSampleDataRW.
383 Doing so might introduce additional sharing.
384 */
385 void
386 requireWrite();
387
388 /**
389 \brief
390 Return true if this Data is expanded.
391 \note To determine if a sample will contain separate values for each datapoint. Use actsExpanded instead.
392 */
393 bool
394 isExpanded() const;
395
396 /**
397 \brief
398 Return true if this Data is expanded or resolves to expanded.
399 That is, if it has a separate value for each datapoint in the sample.
400 */
401 bool
402 actsExpanded() const;
403
404
405 /**
406 \brief
407 Return true if this Data is tagged.
408 */
409 bool
410 isTagged() const;
411
412 /**
413 \brief
414 Return true if this Data is constant.
415 */
416 bool
417 isConstant() const;
418
419 /**
420 \brief Return true if this Data is lazy.
421 */
422 bool
423 isLazy() const;
424
425 /**
426 \brief Return true if this data is ready.
427 */
428 bool
429 isReady() const;
430
431 /**
432 \brief
433 Return true if this Data holds an instance of DataEmpty. This is _not_ the same as asking if the object
434 contains datapoints.
435 */
436 bool
437 isEmpty() const;
438
439 /**
440 \brief
441 True if components of this data are stored as complex
442 */
443 bool
444 isComplex() const;
445
446 /**
447 \brief
448 Return the function space.
449 */
450 inline
451 const FunctionSpace&
452 getFunctionSpace() const
453 {
454 return m_data->getFunctionSpace();
455 }
456
457 /**
458 \brief
459 Return the domain.
460 */
461 inline
462 // const AbstractDomain&
463 const_Domain_ptr
464 getDomain() const
465 {
466 return getFunctionSpace().getDomain();
467 }
468
469
470 /**
471 \brief
472 Return the domain.
473 TODO: For internal use only. This should be removed.
474 */
475 inline
476 // const AbstractDomain&
477 Domain_ptr
478 getDomainPython() const
479 {
480 return getFunctionSpace().getDomainPython();
481 }
482
483 /**
484 \brief
485 Return the rank of the point data.
486 */
487 inline
488 unsigned int
489 getDataPointRank() const
490 {
491 return m_data->getRank();
492 }
493
494 /**
495 \brief
496 Return the number of data points
497 */
498 inline
499 int
500 getNumDataPoints() const
501 {
502 return getNumSamples() * getNumDataPointsPerSample();
503 }
504 /**
505 \brief
506 Return the number of samples.
507 */
508 inline
509 int
510 getNumSamples() const
511 {
512 return m_data->getNumSamples();
513 }
514
515 /**
516 \brief
517 Return the number of data points per sample.
518 */
519 inline
520 int
521 getNumDataPointsPerSample() const
522 {
523 return m_data->getNumDPPSample();
524 }
525
526 /**
527 \brief
528 Returns true if the number of data points per sample and the number of
529 samples match the respective argument. DataEmpty always returns true.
530 */
531 inline
532 bool numSamplesEqual(int numDataPointsPerSample, int numSamples) const
533 {
534 return (isEmpty() ||
535 (numDataPointsPerSample==getNumDataPointsPerSample() && numSamples==getNumSamples()));
536 }
537
538 /**
539 \brief
540 Returns true if the shape matches the vector (dimensions[0],...,
541 dimensions[rank-1]). DataEmpty always returns true.
542 */
543 inline
544 bool isDataPointShapeEqual(int rank, const int* dimensions) const
545 {
546 if (isEmpty())
547 return true;
548 const DataTypes::ShapeType givenShape(&dimensions[0],&dimensions[rank]);
549 return (getDataPointShape()==givenShape);
550 }
551
552 /**
553 \brief
554 Return the number of values in the shape for this object.
555 */
556 int
557 getNoValues() const
558 {
559 return m_data->getNoValues();
560 }
561
562
563 /**
564 \brief
565 dumps the object into a netCDF file
566 */
567 void
568 dump(const std::string fileName) const;
569
570 /**
571 \brief returns the values of the object as a list of tuples (one for each datapoint).
572
573 \param scalarastuple If true, scalar data will produce single valued tuples [(1,) (2,) ...]
574 If false, the result is a list of scalars [1, 2, ...]
575 */
576 const boost::python::object
577 toListOfTuples(bool scalarastuple=true);
578
579
580 /**
581 \brief
582 Return the sample data for the given sample no.
583 Please do not use this unless you NEED to access samples individually
584 \param sampleNo - Input - the given sample no.
585 \return pointer to the sample data.
586 */
587 const DataTypes::real_t*
588 getSampleDataRO(DataTypes::RealVectorType::size_type sampleNo, DataTypes::real_t dummy=0) const;
589
590 const DataTypes::cplx_t*
591 getSampleDataRO(DataTypes::CplxVectorType::size_type sampleNo, DataTypes::cplx_t dummy) const;
592
593
594 /**
595 \brief
596 Return the sample data for the given sample no.
597 Please do not use this unless you NEED to access samples individually
598 \param sampleNo - Input - the given sample no.
599 \return pointer to the sample data.
600 */
601 DataTypes::real_t*
602 getSampleDataRW(DataTypes::RealVectorType::size_type sampleNo, DataTypes::real_t dummy=0);
603
604 DataTypes::cplx_t*
605 getSampleDataRW(DataTypes::RealVectorType::size_type sampleNo, DataTypes::cplx_t dummy);
606
607
608
609 /**
610 \brief
611 Return a pointer to the beginning of the underlying data
612 \warning please avoid using this method since it by-passes possible lazy improvements. May be removed without notice.
613 \return pointer to the data.
614 */
615 const DataTypes::real_t*
616 getDataRO(DataTypes::real_t dummy=0) const;
617
618 const DataTypes::cplx_t*
619 getDataRO(DataTypes::cplx_t dummy) const;
620
621
622
623 /**
624 \brief
625 Return the sample data for the given tag. If an attempt is made to
626 access data that isn't tagged an exception will be thrown.
627 \param tag - Input - the tag key.
628 */
629 inline
630 DataTypes::real_t*
631 getSampleDataByTag(int tag, DataTypes::real_t dummy=0)
632 {
633 return m_data->getSampleDataByTag(tag, dummy);
634 }
635
636 inline
637 DataTypes::cplx_t*
638 getSampleDataByTag(int tag, DataTypes::cplx_t dummy)
639 {
640 return m_data->getSampleDataByTag(tag, dummy);
641 }
642
643
644 /**
645 \brief
646 Return a reference into the DataVector which points to the specified data point.
647 \param sampleNo - Input -
648 \param dataPointNo - Input -
649 */
650 DataTypes::RealVectorType::const_reference
651 getDataPointRO(int sampleNo, int dataPointNo);
652
653 /**
654 \brief
655 Return a reference into the DataVector which points to the specified data point.
656 \param sampleNo - Input -
657 \param dataPointNo - Input -
658 */
659 DataTypes::RealVectorType::reference
660 getDataPointRW(int sampleNo, int dataPointNo);
661
662
663
664 /**
665 \brief
666 Return the offset for the given sample and point within the sample
667 */
668 inline
669 DataTypes::RealVectorType::size_type
670 getDataOffset(int sampleNo,
671 int dataPointNo)
672 {
673 return m_data->getPointOffset(sampleNo,dataPointNo);
674 }
675
676 /**
677 \brief
678 Return a reference to the data point shape.
679 */
680 inline
681 const DataTypes::ShapeType&
682 getDataPointShape() const
683 {
684 return m_data->getShape();
685 }
686
687 /**
688 \brief
689 Return the data point shape as a tuple of integers.
690 */
691 const boost::python::tuple
692 getShapeTuple() const;
693
694 /**
695 \brief
696 Return the size of the data point. It is the product of the
697 data point shape dimensions.
698 */
699 int
700 getDataPointSize() const;
701
702 /**
703 \brief
704 Return the number of doubles stored for this Data.
705 */
706 DataTypes::RealVectorType::size_type
707 getLength() const;
708
709 /**
710 \brief Return true if this object contains no samples.
711 This is not the same as isEmpty()
712 */
713 bool
714 hasNoSamples() const
715 {
716 return m_data->getNumSamples()==0;
717 }
718
719 /**
720 \brief
721 Assign the given value to the tag assocciated with name. Implicitly converts this
722 object to type DataTagged. Throws an exception if this object
723 cannot be converted to a DataTagged object or name cannot be mapped onto a tag key.
724 \param name - Input - name of tag.
725 \param value - Input - Value to associate with given key.
726 */
727 void
728 setTaggedValueByName(std::string name,
729 const boost::python::object& value);
730
731 /**
732 \brief
733 Assign the given value to the tag. Implicitly converts this
734 object to type DataTagged if it is constant.
735
736 \param tagKey - Input - Integer key.
737 \param value - Input - Value to associate with given key.
738 ==>*
739 */
740 void
741 setTaggedValue(int tagKey,
742 const boost::python::object& value);
743
744 /**
745 \brief
746 Assign the given value to the tag. Implicitly converts this
747 object to type DataTagged if it is constant.
748
749 \param tagKey - Input - Integer key.
750 \param pointshape - Input - The shape of the value parameter
751 \param value - Input - Value to associate with given key.
752 \param dataOffset - Input - Offset of the begining of the point within the value parameter
753 */
754 void
755 setTaggedValueFromCPP(int tagKey,
756 const DataTypes::ShapeType& pointshape,
757 const DataTypes::RealVectorType& value,
758 int dataOffset=0);
759
760
761 void
762 setTaggedValueFromCPP(int tagKey,
763 const DataTypes::ShapeType& pointshape,
764 const DataTypes::CplxVectorType& value,
765 int dataOffset=0);
766
767 /**
768 \brief
769 Copy other Data object into this Data object where mask is positive.
770 */
771 void
772 copyWithMask(const Data& other,
773 const Data& mask);
774
775 /**
776 Data object operation methods and operators.
777 */
778
779 /**
780 \brief
781 set all values to zero
782 *
783 */
784 void
785 setToZero();
786
787 /**
788 \brief
789 Interpolates this onto the given functionspace and returns
790 the result as a Data object.
791 *
792 */
793 Data
794 interpolate(const FunctionSpace& functionspace) const;
795
796 Data
797 interpolateFromTable3D(const WrappedArray& table, DataTypes::real_t Amin, DataTypes::real_t Astep,
798 DataTypes::real_t undef, Data& B, DataTypes::real_t Bmin, DataTypes::real_t Bstep, Data& C,
799 DataTypes::real_t Cmin, DataTypes::real_t Cstep, bool check_boundaries);
800
801 Data
802 interpolateFromTable2D(const WrappedArray& table, DataTypes::real_t Amin, DataTypes::real_t Astep,
803 DataTypes::real_t undef, Data& B, DataTypes::real_t Bmin, DataTypes::real_t Bstep,bool check_boundaries);
804
805 Data
806 interpolateFromTable1D(const WrappedArray& table, DataTypes::real_t Amin, DataTypes::real_t Astep,
807 DataTypes::real_t undef,bool check_boundaries);
808
809
810 Data
811 interpolateFromTable3DP(boost::python::object table, DataTypes::real_t Amin, DataTypes::real_t Astep,
812 Data& B, DataTypes::real_t Bmin, DataTypes::real_t Bstep, Data& C, DataTypes::real_t Cmin, DataTypes::real_t Cstep, DataTypes::real_t undef,bool check_boundaries);
813
814
815 Data
816 interpolateFromTable2DP(boost::python::object table, DataTypes::real_t Amin, DataTypes::real_t Astep,
817 Data& B, DataTypes::real_t Bmin, DataTypes::real_t Bstep, DataTypes::real_t undef,bool check_boundaries);
818
819 Data
820 interpolateFromTable1DP(boost::python::object table, DataTypes::real_t Amin, DataTypes::real_t Astep,
821 DataTypes::real_t undef,bool check_boundaries);
822
823 Data
824 nonuniforminterp(boost::python::object in, boost::python::object out, bool check_boundaries);
825
826 Data
827 nonuniformslope(boost::python::object in, boost::python::object out, bool check_boundaries);
828
829 /**
830 \brief
831 Calculates the gradient of the data at the data points of functionspace.
832 If functionspace is not present the function space of Function(getDomain()) is used.
833 *
834 */
835 Data
836 gradOn(const FunctionSpace& functionspace) const;
837
838 Data
839 grad() const;
840
841 /**
842 \brief
843 Calculate the integral over the function space domain as a python tuple.
844 */
845 boost::python::object
846 integrateToTuple_const() const;
847
848
849 /**
850 \brief
851 Calculate the integral over the function space domain as a python tuple.
852 */
853 boost::python::object
854 integrateToTuple();
855
856
857
858 /**
859 \brief
860 Returns 1./ Data object
861 *
862 */
863 Data
864 oneOver() const;
865 /**
866 \brief
867 Return a Data with a 1 for +ive values and a 0 for 0 or -ive values.
868 *
869 */
870 Data
871 wherePositive() const;
872
873 /**
874 \brief
875 Return a Data with a 1 for -ive values and a 0 for +ive or 0 values.
876 *
877 */
878 Data
879 whereNegative() const;
880
881 /**
882 \brief
883 Return a Data with a 1 for +ive or 0 values and a 0 for -ive values.
884 *
885 */
886 Data
887 whereNonNegative() const;
888
889 /**
890 \brief
891 Return a Data with a 1 for -ive or 0 values and a 0 for +ive values.
892 *
893 */
894 Data
895 whereNonPositive() const;
896
897 /**
898 \brief
899 Return a Data with a 1 for 0 values and a 0 for +ive or -ive values.
900 *
901 */
902 Data
903 whereZero(DataTypes::real_t tol=0.0) const;
904
905 /**
906 \brief
907 Return a Data with a 0 for 0 values and a 1 for +ive or -ive values.
908 *
909 */
910 Data
911 whereNonZero(DataTypes::real_t tol=0.0) const;
912
913 /**
914 \brief
915 Return the maximum absolute value of this Data object.
916
917 The method is not const because lazy data needs to be expanded before Lsup can be computed.
918 The _const form can be used when the Data object is const, however this will only work for
919 Data which is not Lazy.
920
921 For Data which contain no samples (or tagged Data for which no tags in use have a value)
922 zero is returned.
923 */
924 DataTypes::real_t
925 Lsup();
926
927 DataTypes::real_t
928 Lsup_const() const;
929
930
931 /**
932 \brief
933 Return the maximum value of this Data object.
934
935 The method is not const because lazy data needs to be expanded before sup can be computed.
936 The _const form can be used when the Data object is const, however this will only work for
937 Data which is not Lazy.
938
939 For Data which contain no samples (or tagged Data for which no tags in use have a value)
940 a large negative value is returned.
941 */
942 DataTypes::real_t
943 sup();
944
945 DataTypes::real_t
946 sup_const() const;
947
948
949 /**
950 \brief
951 Return the minimum value of this Data object.
952
953 The method is not const because lazy data needs to be expanded before inf can be computed.
954 The _const form can be used when the Data object is const, however this will only work for
955 Data which is not Lazy.
956
957 For Data which contain no samples (or tagged Data for which no tags in use have a value)
958 a large positive value is returned.
959 */
960 DataTypes::real_t
961 inf();
962
963 DataTypes::real_t
964 inf_const() const;
965
966
967
968 /**
969 \brief
970 Return the absolute value of each data point of this Data object.
971 *
972 */
973 Data
974 abs() const;
975
976 /**
977 \brief
978 Return the phase/arg/angular-part of complex values.
979 *
980 */
981 Data
982 phase() const;
983
984
985 /**
986 \brief
987 Return the maximum value of each data point of this Data object.
988 *
989 */
990 Data
991 maxval() const;
992
993 /**
994 \brief
995 Return the minimum value of each data point of this Data object.
996 *
997 */
998 Data
999 minval() const;
1000
1001 /**
1002 \brief
1003 Return the (sample number, data-point number) of the data point with
1004 the minimum component value in this Data object.
1005 \note If you are working in python, please consider using Locator
1006 instead of manually manipulating process and point IDs.
1007 */
1008 const boost::python::tuple
1009 minGlobalDataPoint() const;
1010
1011 /**
1012 \brief
1013 Return the (sample number, data-point number) of the data point with
1014 the minimum component value in this Data object.
1015 \note If you are working in python, please consider using Locator
1016 instead of manually manipulating process and point IDs.
1017 */
1018 const boost::python::tuple
1019 maxGlobalDataPoint() const;
1020
1021
1022
1023 /**
1024 \brief
1025 Return the sign of each data point of this Data object.
1026 -1 for negative values, zero for zero values, 1 for positive values.
1027 *
1028 */
1029 Data
1030 sign() const;
1031
1032 /**
1033 \brief
1034 Return the symmetric part of a matrix which is half the matrix plus its transpose.
1035 *
1036 */
1037 Data
1038 symmetric() const;
1039
1040 /**
1041 \brief
1042 Return the antisymmetric part of a matrix which is half the matrix minus its transpose.
1043 *
1044 */
1045 Data
1046 antisymmetric() const;
1047
1048
1049 /**
1050 \brief
1051 Return the hermitian part of a matrix which is half the matrix plus its adjoint.
1052 *
1053 */
1054 Data
1055 hermitian() const;
1056
1057 /**
1058 \brief
1059 Return the anti-hermitian part of a matrix which is half the matrix minus its hermitian.
1060 *
1061 */
1062 Data
1063 antihermitian() const;
1064
1065 /**
1066 \brief
1067 Return the trace of a matrix
1068 *
1069 */
1070 Data
1071 trace(int axis_offset) const;
1072
1073 /**
1074 \brief
1075 Transpose each data point of this Data object around the given axis.
1076 *
1077 */
1078 Data
1079 transpose(int axis_offset) const;
1080
1081 /**
1082 \brief
1083 Return the eigenvalues of the symmetric part at each data point of this Data object in increasing values.
1084 Currently this function is restricted to rank 2, square shape, and dimension 3.
1085 *
1086 */
1087 Data
1088 eigenvalues() const;
1089
1090 /**
1091 \brief
1092 Return the eigenvalues and corresponding eigenvcetors of the symmetric part at each data point of this Data object.
1093 the eigenvalues are ordered in increasing size where eigenvalues with relative difference less than
1094 tol are treated as equal. The eigenvectors are orthogonal, normalized and the sclaed such that the
1095 first non-zero entry is positive.
1096 Currently this function is restricted to rank 2, square shape, and dimension 3
1097 *
1098 */
1099 const boost::python::tuple
1100 eigenvalues_and_eigenvectors(const DataTypes::real_t tol=1.e-12) const;
1101
1102 /**
1103 \brief
1104 swaps the components axis0 and axis1
1105 *
1106 */
1107 Data
1108 swapaxes(const int axis0, const int axis1) const;
1109
1110 /**
1111 \brief
1112 Return the error function erf of each data point of this Data object.
1113 *
1114 */
1115 Data
1116 erf() const;
1117
1118
1119 /**
1120 \brief
1121 For complex values return the conjugate values.
1122 For non-complex data return a copy
1123 */
1124 Data
1125 conjugate() const;
1126
1127 Data
1128 real() const;
1129
1130 Data
1131 imag() const;
1132
1133 /**
1134 \brief
1135 Return the sin of each data point of this Data object.
1136 *
1137 */
1138 Data
1139 sin() const;
1140
1141 /**
1142 \brief
1143 Return the cos of each data point of this Data object.
1144 *
1145 */
1146 Data
1147 cos() const;
1148
1149 /**
1150 \brief
1151 Bessel worker function.
1152 *
1153 */
1154 Data
1155 bessel(int order, DataTypes::real_t (*besselfunc) (int,DataTypes::real_t) );
1156
1157
1158 /**
1159 \brief
1160 Return the Bessel function of the first kind for each data point of this Data object.
1161 *
1162 */
1163 Data
1164 besselFirstKind(int order);
1165
1166 /**
1167 \brief
1168 Return the Bessel function of the second kind for each data point of this Data object.
1169 *
1170 */
1171 Data
1172 besselSecondKind(int order);
1173
1174
1175 /**
1176 \brief
1177 Return the tan of each data point of this Data object.
1178 *
1179 */
1180 Data
1181 tan() const;
1182
1183 /**
1184 \brief
1185 Return the asin of each data point of this Data object.
1186 *
1187 */
1188 Data
1189 asin() const;
1190
1191 /**
1192 \brief
1193 Return the acos of each data point of this Data object.
1194 *
1195 */
1196 Data
1197 acos() const;
1198
1199 /**
1200 \brief
1201 Return the atan of each data point of this Data object.
1202 *
1203 */
1204 Data
1205 atan() const;
1206
1207 /**
1208 \brief
1209 Return the sinh of each data point of this Data object.
1210 *
1211 */
1212 Data
1213 sinh() const;
1214
1215 /**
1216 \brief
1217 Return the cosh of each data point of this Data object.
1218 *
1219 */
1220 Data
1221 cosh() const;
1222
1223 /**
1224 \brief
1225 Return the tanh of each data point of this Data object.
1226 *
1227 */
1228 Data
1229 tanh() const;
1230
1231 /**
1232 \brief
1233 Return the asinh of each data point of this Data object.
1234 *
1235 */
1236 Data
1237 asinh() const;
1238
1239 /**
1240 \brief
1241 Return the acosh of each data point of this Data object.
1242 *
1243 */
1244 Data
1245 acosh() const;
1246
1247 /**
1248 \brief
1249 Return the atanh of each data point of this Data object.
1250 *
1251 */
1252 Data
1253 atanh() const;
1254
1255 /**
1256 \brief
1257 Return the log to base 10 of each data point of this Data object.
1258 *
1259 */
1260 Data
1261 log10() const;
1262
1263 /**
1264 \brief
1265 Return the natural log of each data point of this Data object.
1266 *
1267 */
1268 Data
1269 log() const;
1270
1271 /**
1272 \brief
1273 Return the exponential function of each data point of this Data object.
1274 *
1275 */
1276 Data
1277 exp() const;
1278
1279 /**
1280 \brief
1281 Return the square root of each data point of this Data object.
1282 *
1283 */
1284 Data
1285 sqrt() const;
1286
1287 /**
1288 \brief
1289 Return the negation of each data point of this Data object.
1290 *
1291 */
1292 Data
1293 neg() const;
1294
1295 /**
1296 \brief
1297 Return the identity of each data point of this Data object.
1298 Simply returns this object unmodified.
1299 *
1300 */
1301 Data
1302 pos() const;
1303
1304 /**
1305 \brief
1306 Return the given power of each data point of this Data object.
1307
1308 \param right Input - the power to raise the object to.
1309 *
1310 */
1311 Data
1312 powD(const Data& right) const;
1313
1314 /**
1315 \brief
1316 Return the given power of each data point of this boost python object.
1317
1318 \param right Input - the power to raise the object to.
1319 *
1320 */
1321 Data
1322 powO(const boost::python::object& right) const;
1323
1324 /**
1325 \brief
1326 Return the given power of each data point of this boost python object.
1327
1328 \param left Input - the bases
1329 *
1330 */
1331
1332 Data
1333 rpowO(const boost::python::object& left) const;
1334
1335 /**
1336 \brief
1337 Overloaded operator +=
1338 \param right - Input - The right hand side.
1339 *
1340 */
1341 Data& operator+=(const Data& right);
1342 Data& operator+=(const boost::python::object& right);
1343
1344 Data& operator=(const Data& other);
1345
1346 /**
1347 \brief
1348 Overloaded operator -=
1349 \param right - Input - The right hand side.
1350 *
1351 */
1352 Data& operator-=(const Data& right);
1353 Data& operator-=(const boost::python::object& right);
1354
1355 /**
1356 \brief
1357 Overloaded operator *=
1358 \param right - Input - The right hand side.
1359 *
1360 */
1361 Data& operator*=(const Data& right);
1362 Data& operator*=(const boost::python::object& right);
1363
1364 /**
1365 \brief
1366 Overloaded operator /=
1367 \param right - Input - The right hand side.
1368 *
1369 */
1370 Data& operator/=(const Data& right);
1371 Data& operator/=(const boost::python::object& right);
1372
1373 /**
1374 \brief
1375 Newer style division operator for python
1376 */
1377 Data truedivD(const Data& right);
1378
1379 /**
1380 \brief
1381 Newer style division operator for python
1382 */
1383 Data truedivO(const boost::python::object& right);
1384
1385 /**
1386 \brief
1387 Newer style division operator for python
1388 */
1389 Data rtruedivO(const boost::python::object& left);
1390
1391 /**
1392 \brief
1393 wrapper for python add operation
1394 */
1395 boost::python::object __add__(const boost::python::object& right);
1396
1397
1398 /**
1399 \brief
1400 wrapper for python subtract operation
1401 */
1402 boost::python::object __sub__(const boost::python::object& right);
1403
1404 /**
1405 \brief
1406 wrapper for python reverse subtract operation
1407 */
1408 boost::python::object __rsub__(const boost::python::object& right);
1409
1410 /**
1411 \brief
1412 wrapper for python multiply operation
1413 */
1414 boost::python::object __mul__(const boost::python::object& right);
1415
1416 /**
1417 \brief
1418 wrapper for python divide operation
1419 */
1420 boost::python::object __div__(const boost::python::object& right);
1421
1422 /**
1423 \brief
1424 wrapper for python reverse divide operation
1425 */
1426 boost::python::object __rdiv__(const boost::python::object& right);
1427
1428 /**
1429 \brief return inverse of matricies.
1430 */
1431 Data
1432 matrixInverse() const;
1433
1434 /**
1435 \brief
1436 Returns true if this can be interpolated to functionspace.
1437 */
1438 bool
1439 probeInterpolation(const FunctionSpace& functionspace) const;
1440
1441 /**
1442 Data object slicing methods.
1443 */
1444
1445 /**
1446 \brief
1447 Returns a slice from this Data object.
1448
1449 /description
1450 Implements the [] get operator in python.
1451 Calls getSlice.
1452
1453 \param key - Input - python slice tuple specifying
1454 slice to return.
1455 */
1456 Data
1457 getItem(const boost::python::object& key) const;
1458
1459 /**
1460 \brief
1461 Copies slice from value into this Data object.
1462
1463 Implements the [] set operator in python.
1464 Calls setSlice.
1465
1466 \param key - Input - python slice tuple specifying
1467 slice to copy from value.
1468 \param value - Input - Data object to copy from.
1469 */
1470 void
1471 setItemD(const boost::python::object& key,
1472 const Data& value);
1473
1474 void
1475 setItemO(const boost::python::object& key,
1476 const boost::python::object& value);
1477
1478 // These following public methods should be treated as private.
1479
1480 /**
1481 \brief
1482 Perform the given unary operation on every element of every data point in
1483 this Data object.
1484 */
1485 template <class UnaryFunction>
1486 inline
1487 void
1488 unaryOp2(UnaryFunction operation);
1489
1490 /**
1491 \brief
1492 Return a Data object containing the specified slice of
1493 this Data object.
1494 \param region - Input - Region to copy.
1495 *
1496 */
1497 Data
1498 getSlice(const DataTypes::RegionType& region) const;
1499
1500 /**
1501 \brief
1502 Copy the specified slice from the given value into this
1503 Data object.
1504 \param value - Input - Data to copy from.
1505 \param region - Input - Region to copy.
1506 *
1507 */
1508 void
1509 setSlice(const Data& value,
1510 const DataTypes::RegionType& region);
1511
1512 /**
1513 \brief
1514 print the data values to stdout. Used for debugging
1515 */
1516 void
1517 print(void);
1518
1519 /**
1520 \brief
1521 return the MPI rank number of the local data
1522 MPI_COMM_WORLD is assumed and the result of MPI_Comm_size()
1523 is returned
1524 */
1525 int
1526 get_MPIRank(void) const;
1527
1528 /**
1529 \brief
1530 return the MPI rank number of the local data
1531 MPI_COMM_WORLD is assumed and the result of MPI_Comm_rank()
1532 is returned
1533 */
1534 int
1535 get_MPISize(void) const;
1536
1537 /**
1538 \brief
1539 return the MPI rank number of the local data
1540 MPI_COMM_WORLD is assumed and returned.
1541 */
1542 MPI_Comm
1543 get_MPIComm(void) const;
1544
1545 /**
1546 \brief
1547 return the object produced by the factory, which is a DataConstant or DataExpanded
1548 TODO Ownership of this object should be explained in doco.
1549 */
1550 DataAbstract*
1551 borrowData(void) const;
1552
1553 DataAbstract_ptr
1554 borrowDataPtr(void) const;
1555
1556 DataReady_ptr
1557 borrowReadyPtr(void) const;
1558
1559
1560
1561 /**
1562 \brief
1563 Return a pointer to the beginning of the datapoint at the specified offset.
1564 TODO Eventually these should be inlined.
1565 \param i - position(offset) in the underlying datastructure
1566 */
1567
1568 DataTypes::RealVectorType::const_reference
1569 getDataAtOffsetRO(DataTypes::RealVectorType::size_type i, DataTypes::real_t dummy=0);
1570
1571 DataTypes::RealVectorType::reference
1572 getDataAtOffsetRW(DataTypes::RealVectorType::size_type i, DataTypes::real_t dummy=0);
1573
1574 DataTypes::CplxVectorType::const_reference
1575 getDataAtOffsetRO(DataTypes::CplxVectorType::size_type i, DataTypes::cplx_t dummy);
1576
1577 DataTypes::CplxVectorType::reference
1578 getDataAtOffsetRW(DataTypes::CplxVectorType::size_type i, DataTypes::cplx_t dummy);
1579
1580
1581 /**
1582 \brief Ensures that the Data is expanded and returns its underlying vector
1583 Does not check for exclusive write so do that before calling if sharing
1584 Is a posibility.
1585 \warning For domain implementors only. Using this function will
1586 avoid using optimisations like lazy evaluation. It is intended
1587 to allow quick initialisation of Data by domain; not as a bypass around
1588 escript's other mechanisms.
1589 */
1590 DataTypes::RealVectorType&
1591 getExpandedVectorReference(DataTypes::real_t dummy=0);
1592
1593 DataTypes::CplxVectorType&
1594 getExpandedVectorReference(DataTypes::cplx_t dummy);
1595
1596
1597 /**
1598 * \brief For tagged Data returns the number of tags with values.
1599 * For non-tagged data will return 0 (even Data which has been expanded from tagged).
1600 */
1601 size_t
1602 getNumberOfTaggedValues() const;
1603
1604 /*
1605 * \brief make the data complex
1606 */
1607 void complicate();
1608
1609 protected:
1610
1611 private:
1612 void init_from_data_and_fs(const Data& inData,
1613 const FunctionSpace& functionspace);
1614
1615 template <typename S>
1616 void
1617 maskWorker(Data& other2, Data& mask2, S sentinel);
1618
1619 template <class BinaryOp>
1620 DataTypes::real_t
1621 #ifdef ESYS_MPI
1622 lazyAlgWorker(DataTypes::real_t init, MPI_Op mpiop_type);
1623 #else
1624 lazyAlgWorker(DataTypes::real_t init);
1625 #endif
1626
1627 DataTypes::real_t
1628 LsupWorker() const;
1629
1630 DataTypes::real_t
1631 supWorker() const;
1632
1633 DataTypes::real_t
1634 infWorker() const;
1635
1636 boost::python::object
1637 integrateWorker() const;
1638
1639 void
1640 calc_minGlobalDataPoint(int& ProcNo, int& DataPointNo) const;
1641
1642 void
1643 calc_maxGlobalDataPoint(int& ProcNo, int& DataPointNo) const;
1644
1645 // For internal use in Data.cpp only!
1646 // other uses should call the main entry points and allow laziness
1647 Data
1648 minval_nonlazy() const;
1649
1650 // For internal use in Data.cpp only!
1651 Data
1652 maxval_nonlazy() const;
1653
1654
1655 /**
1656 \brief
1657 Check *this and the right operand are compatible. Throws
1658 an exception if they aren't.
1659 \param right - Input - The right hand side.
1660 */
1661 inline
1662 void
1663 operandCheck(const Data& right) const
1664 {
1665 return m_data->operandCheck(*(right.m_data.get()));
1666 }
1667
1668 /**
1669 \brief
1670 Perform the specified reduction algorithm on every element of every data point in
1671 this Data object according to the given function and return the single value result.
1672 */
1673 template <class BinaryFunction>
1674 inline
1675 DataTypes::real_t
1676 reduction(BinaryFunction operation,
1677 DataTypes::real_t initial_value) const;
1678
1679 /**
1680 \brief
1681 Reduce each data-point in this Data object using the given operation. Return a Data
1682 object with the same number of data-points, but with each data-point containing only
1683 one value - the result of the reduction operation on the corresponding data-point in
1684 this Data object
1685 */
1686 template <class BinaryFunction>
1687 inline
1688 Data
1689 dp_algorithm(BinaryFunction operation,
1690 DataTypes::real_t initial_value) const;
1691
1692 /**
1693 \brief
1694 Convert the data type of the RHS to match this.
1695 \param right - Input - data type to match.
1696 */
1697 void
1698 typeMatchLeft(Data& right) const;
1699
1700 /**
1701 \brief
1702 Convert the data type of this to match the RHS.
1703 \param right - Input - data type to match.
1704 */
1705 void
1706 typeMatchRight(const Data& right);
1707
1708 /**
1709 \brief
1710 Construct a Data object of the appropriate type.
1711 */
1712
1713 void
1714 initialise(const DataTypes::RealVectorType& value,
1715 const DataTypes::ShapeType& shape,
1716 const FunctionSpace& what,
1717 bool expanded);
1718
1719 void
1720 initialise(const DataTypes::CplxVectorType& value,
1721 const DataTypes::ShapeType& shape,
1722 const FunctionSpace& what,
1723 bool expanded);
1724
1725 void
1726 initialise(const WrappedArray& value,
1727 const FunctionSpace& what,
1728 bool expanded);
1729
1730 void
1731 initialise(const DataTypes::real_t value,
1732 const DataTypes::ShapeType& shape,
1733 const FunctionSpace& what,
1734 bool expanded);
1735
1736 void
1737 initialise(const DataTypes::cplx_t value,
1738 const DataTypes::ShapeType& shape,
1739 const FunctionSpace& what,
1740 bool expanded);
1741 //
1742 // flag to protect the data object against any update
1743 bool m_protected;
1744 bool m_lazy;
1745
1746 //
1747 // pointer to the actual data object
1748 // boost::shared_ptr<DataAbstract> m_data;
1749 DataAbstract_ptr m_data;
1750
1751 // If possible please use getReadyPtr instead.
1752 // But see warning below.
1753 const DataReady*
1754 getReady() const
1755 {
1756 const DataReady* dr=dynamic_cast<const DataReady*>(m_data.get());
1757 ESYS_ASSERT(dr!=0, "error casting to DataReady.");
1758 return dr;
1759 }
1760
1761 DataReady*
1762 getReady()
1763 {
1764 DataReady* dr=dynamic_cast<DataReady*>(m_data.get());
1765 ESYS_ASSERT(dr!=0, "error casting to DataReady.");
1766 return dr;
1767 }
1768
1769
1770 // Be wary of using this for local operations since it (temporarily) increases reference count.
1771 // If you are just using this to call a method on DataReady instead of DataAbstract consider using
1772 // getReady() instead
1773 DataReady_ptr
1774 getReadyPtr()
1775 {
1776 DataReady_ptr dr=REFCOUNTNS::dynamic_pointer_cast<DataReady>(m_data);
1777 ESYS_ASSERT(dr.get()!=0, "error casting to DataReady.");
1778 return dr;
1779 }
1780
1781 const_DataReady_ptr
1782 getReadyPtr() const
1783 {
1784 const_DataReady_ptr dr=REFCOUNTNS::dynamic_pointer_cast<const DataReady>(m_data);
1785 ESYS_ASSERT(dr.get()!=0, "error casting to DataReady.");
1786 return dr;
1787 }
1788
1789 // In the isShared() method below:
1790 // A problem would occur if m_data (the address pointed to) were being modified
1791 // while the call m_data->is_shared is being executed.
1792 //
1793 // Q: So why do I think this code can be thread safe/correct?
1794 // A: We need to make some assumptions.
1795 // 1. We assume it is acceptable to return true under some conditions when we aren't shared.
1796 // 2. We assume that no constructions or assignments which will share previously unshared
1797 // will occur while this call is executing. This is consistent with the way Data:: and C are written.
1798 //
1799 // This means that the only transition we need to consider, is when a previously shared object is
1800 // not shared anymore. ie. the other objects have been destroyed or a deep copy has been made.
1801 // In those cases the m_shared flag changes to false after m_data has completed changing.
1802 // For any threads executing before the flag switches they will assume the object is still shared.
1803 bool isShared() const
1804 {
1805 #ifdef SLOWSHARECHECK
1806 return m_data->isShared(); // single threadsafe check for this
1807 #else
1808 return !m_data.unique();
1809 #endif
1810 }
1811
1812 void forceResolve()
1813 {
1814 if (isLazy())
1815 {
1816 #ifdef _OPENMP
1817 if (omp_in_parallel())
1818 { // Yes this is throwing an exception out of an omp thread which is forbidden.
1819 throw DataException("Please do not call forceResolve() in a parallel region.");
1820 }
1821 #endif
1822 resolve();
1823 }
1824 }
1825
1826 /**
1827 \brief if another object is sharing out member data make a copy to work with instead.
1828 This code should only be called from single threaded sections of code.
1829 */
1830 void exclusiveWrite()
1831 {
1832 #ifdef _OPENMP
1833 if (omp_in_parallel())
1834 {
1835 throw DataException("Programming error. Please do not run exclusiveWrite() in multi-threaded sections.");
1836 }
1837 #endif
1838 forceResolve();
1839 if (isShared())
1840 {
1841 DataAbstract* t=m_data->deepCopy();
1842 set_m_data(DataAbstract_ptr(t));
1843 }
1844 #ifdef EXWRITECHK
1845 m_data->exclusivewritecalled=true;
1846 #endif
1847 }
1848
1849 /**
1850 \brief checks if caller can have exclusive write to the object
1851 */
1852 void checkExclusiveWrite()
1853 {
1854 if (isLazy() || isShared())
1855 {
1856 std::ostringstream oss;
1857 oss << "Programming error. ExclusiveWrite required - please call requireWrite() isLazy=" << isLazy() << " isShared()=" << isShared();
1858 throw DataException(oss.str());
1859 }
1860 }
1861
1862 /**
1863 \brief Modify the data abstract hosted by this Data object
1864 For internal use only.
1865 Passing a pointer to null is permitted (do this in the destructor)
1866 \warning Only to be called in single threaded code or inside a single/critical section. This method needs to be atomic.
1867 */
1868 void set_m_data(DataAbstract_ptr p);
1869
1870
1871 void TensorSelfUpdateBinaryOperation(const Data& right, escript::ES_optype operation);
1872
1873 friend class DataAbstract; // To allow calls to updateShareStatus
1874 friend class TestDomain; // so its getX will work quickly
1875 #ifdef IKNOWWHATIMDOING
1876 friend Data applyBinaryCFunction(boost::python::object cfunc, boost::python::tuple shape, escript::Data& d, escript::Data& e);
1877 #endif
1878 friend Data condEval(escript::Data& mask, escript::Data& trueval, escript::Data& falseval);
1879 friend Data randomData(const boost::python::tuple& shape, const FunctionSpace& what, long seed, const boost::python::tuple& filter);
1880
1881 };
1882
1883
1884 #ifdef IKNOWWHATIMDOING
1885 Data
1886 applyBinaryCFunction(boost::python::object func, boost::python::tuple shape, escript::Data& d, escript::Data& e);
1887 #endif
1888
1889 Data
1890 condEval(escript::Data& mask, escript::Data& trueval, escript::Data& falseval);
1891
1892
1893
1894 /**
1895 \brief Create a new Expanded Data object filled with pseudo-random data.
1896 */
1897 Data randomData(const boost::python::tuple& shape,
1898 const FunctionSpace& what,
1899 long seed, const boost::python::tuple& filter);
1900
1901
1902 } // end namespace escript
1903
1904
1905 // No, this is not supposed to be at the top of the file
1906 // DataAbstact needs to be declared first, then DataReady needs to be fully declared
1907 // so that I can dynamic cast between them below.
1908 #include "DataReady.h"
1909 #include "DataLazy.h"
1910 #include "DataExpanded.h"
1911 #include "DataConstant.h"
1912 #include "DataTagged.h"
1913
1914 namespace escript
1915 {
1916
1917
1918
1919 inline
1920 DataTypes::real_t*
1921 Data::getSampleDataRW(DataTypes::RealVectorType::size_type sampleNo, DataTypes::real_t dummy)
1922 {
1923 if (isLazy())
1924 {
1925 throw DataException("Error, attempt to acquire RW access to lazy data. Please call requireWrite() first.");
1926 }
1927 #ifdef EXWRITECHK
1928 if (!getReady()->exclusivewritecalled)
1929 {
1930 throw DataException("Error, call to Data::getSampleDataRW without a preceeding call to requireWrite/exclusiveWrite.");
1931 }
1932 #endif
1933 return getReady()->getSampleDataRW(sampleNo, dummy);
1934 }
1935
1936 inline
1937 DataTypes::cplx_t*
1938 Data::getSampleDataRW(DataTypes::CplxVectorType::size_type sampleNo, DataTypes::cplx_t dummy)
1939 {
1940 if (isLazy())
1941 {
1942 throw DataException("Error, attempt to acquire RW access to lazy data. Please call requireWrite() first.");
1943 }
1944 #ifdef EXWRITECHK
1945 if (!getReady()->exclusivewritecalled)
1946 {
1947 throw DataException("Error, call to Data::getSampleDataRW without a preceeding call to requireWrite/exclusiveWrite.");
1948 }
1949 #endif
1950 return getReady()->getSampleDataRW(sampleNo, dummy);
1951 }
1952
1953
1954 inline
1955 const DataTypes::real_t*
1956 Data::getSampleDataRO(DataTypes::RealVectorType::size_type sampleNo,DataTypes::real_t dummy) const
1957 {
1958 DataLazy* l=dynamic_cast<DataLazy*>(m_data.get());
1959 if (l!=0)
1960 {
1961 size_t offset=0;
1962 const DataTypes::RealVectorType* res=l->resolveSample(sampleNo,offset);
1963 return &((*res)[offset]);
1964 }
1965 return getReady()->getSampleDataRO(sampleNo, dummy);
1966 }
1967
1968 inline
1969 const DataTypes::cplx_t*
1970 Data::getSampleDataRO(DataTypes::RealVectorType::size_type sampleNo, DataTypes::cplx_t dummy) const
1971 {
1972 DataLazy* l=dynamic_cast<DataLazy*>(m_data.get());
1973 if (l!=0)
1974 {
1975 throw DataException("Programming error: complex lazy objects are not supported.");
1976 }
1977 return getReady()->getSampleDataRO(sampleNo, dummy);
1978 }
1979
1980
1981 inline
1982 const DataTypes::real_t*
1983 Data::getDataRO(DataTypes::real_t dummy) const
1984 {
1985 if (isLazy())
1986 {
1987 throw DataException("Programmer error - getDataRO must not be called on Lazy Data.");
1988 }
1989 if (getNumSamples()==0)
1990 {
1991 return 0;
1992 }
1993 else
1994 {
1995 return &(getReady()->getTypedVectorRO(0)[0]);
1996 }
1997 }
1998
1999 inline
2000 const DataTypes::cplx_t*
2001 Data::getDataRO(DataTypes::cplx_t dummy) const
2002 {
2003 if (isLazy())
2004 {
2005 throw DataException("Programmer error - getDataRO must not be called on Lazy Data.");
2006 }
2007 if (getNumSamples()==0)
2008 {
2009 return 0;
2010 }
2011 else
2012 {
2013 return &(getReady()->getTypedVectorRO(dummy)[0]);
2014 }
2015 }
2016
2017
2018 /**
2019 Binary Data object operators.
2020 */
2021 inline DataTypes::real_t rpow(DataTypes::real_t x,DataTypes::real_t y)
2022 {
2023 return pow(y,x);
2024 }
2025
2026 /**
2027 \brief
2028 Operator+
2029 Takes two Data objects.
2030 */
2031 Data operator+(const Data& left, const Data& right);
2032
2033 /**
2034 \brief
2035 Operator-
2036 Takes two Data objects.
2037 */
2038 Data operator-(const Data& left, const Data& right);
2039
2040 /**
2041 \brief
2042 Operator*
2043 Takes two Data objects.
2044 */
2045 Data operator*(const Data& left, const Data& right);
2046
2047 /**
2048 \brief
2049 Operator/
2050 Takes two Data objects.
2051 */
2052 Data operator/(const Data& left, const Data& right);
2053
2054 /**
2055 \brief
2056 Operator+
2057 Takes LHS Data object and RHS python::object.
2058 python::object must be convertable to Data type.
2059 */
2060 Data operator+(const Data& left, const boost::python::object& right);
2061
2062 /**
2063 \brief
2064 Operator-
2065 Takes LHS Data object and RHS python::object.
2066 python::object must be convertable to Data type.
2067 */
2068 Data operator-(const Data& left, const boost::python::object& right);
2069
2070 /**
2071 \brief
2072 Operator*
2073 Takes LHS Data object and RHS python::object.
2074 python::object must be convertable to Data type.
2075 */
2076 Data operator*(const Data& left, const boost::python::object& right);
2077
2078 /**
2079 \brief
2080 Operator/
2081 Takes LHS Data object and RHS python::object.
2082 python::object must be convertable to Data type.
2083 */
2084 Data operator/(const Data& left, const boost::python::object& right);
2085
2086 /**
2087 \brief
2088 Operator+
2089 Takes LHS python::object and RHS Data object.
2090 python::object must be convertable to Data type.
2091 */
2092 Data operator+(const boost::python::object& left, const Data& right);
2093
2094 /**
2095 \brief
2096 Operator-
2097 Takes LHS python::object and RHS Data object.
2098 python::object must be convertable to Data type.
2099 */
2100 Data operator-(const boost::python::object& left, const Data& right);
2101
2102 /**
2103 \brief
2104 Operator*
2105 Takes LHS python::object and RHS Data object.
2106 python::object must be convertable to Data type.
2107 */
2108 Data operator*(const boost::python::object& left, const Data& right);
2109
2110 /**
2111 \brief
2112 Operator/
2113 Takes LHS python::object and RHS Data object.
2114 python::object must be convertable to Data type.
2115 */
2116 Data operator/(const boost::python::object& left, const Data& right);
2117
2118
2119
2120 /**
2121 \brief
2122 Output operator
2123 */
2124 std::ostream& operator<<(std::ostream& o, const Data& data);
2125
2126 /**
2127 \brief
2128 Compute a tensor product of two Data objects
2129 \param arg_0 - Input - Data object
2130 \param arg_1 - Input - Data object
2131 \param axis_offset - Input - axis offset
2132 \param transpose - Input - 0: transpose neither, 1: transpose arg0, 2: transpose arg1
2133 */
2134 Data
2135 C_GeneralTensorProduct(Data& arg_0,
2136 Data& arg_1,
2137 int axis_offset=0,
2138 int transpose=0);
2139
2140 /**
2141 \brief
2142 Operator/
2143 Takes RHS Data object.
2144 */
2145 inline
2146 Data
2147 Data::truedivD(const Data& right)
2148 {
2149 return *this / right;
2150 }
2151
2152 /**
2153 \brief
2154 Operator/
2155 Takes RHS python::object.
2156 */
2157 inline
2158 Data
2159 Data::truedivO(const boost::python::object& right)
2160 {
2161 Data tmp(right, getFunctionSpace(), false);
2162 return truedivD(tmp);
2163 }
2164
2165 /**
2166 \brief
2167 Operator/
2168 Takes LHS python::object.
2169 */
2170 inline
2171 Data
2172 Data::rtruedivO(const boost::python::object& left)
2173 {
2174 Data tmp(left, getFunctionSpace(), false);
2175 return tmp.truedivD(*this);
2176 }
2177
2178
2179
2180 /**
2181 \brief
2182 Perform the given Data object reduction algorithm on this and return the result.
2183 Given operation combines each element of each data point, thus argument
2184 object (*this) is a rank n Data object, and returned object is a scalar.
2185 Calls escript::algorithm.
2186 */
2187 template <class BinaryFunction>
2188 inline
2189 DataTypes::real_t
2190 Data::reduction(BinaryFunction operation, DataTypes::real_t initial_value) const
2191 {
2192 if (isExpanded()) {
2193 DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());
2194 ESYS_ASSERT(leftC!=0, "Programming error - casting to DataExpanded.");
2195
2196 DataExpanded& data=*leftC;
2197 int i,j;
2198 int numDPPSample=data.getNumDPPSample();
2199 int numSamples=data.getNumSamples();
2200 DataTypes::real_t global_current_value=initial_value;
2201 DataTypes::real_t local_current_value;
2202 const auto& vec=data.getTypedVectorRO(typename BinaryFunction::first_argument_type(0));
2203 const DataTypes::ShapeType& shape=data.getShape();
2204 // calculate the reduction operation value for each data point
2205 // reducing the result for each data-point into the current_value variables
2206 #pragma omp parallel private(local_current_value)
2207 {
2208 local_current_value=initial_value;
2209 #pragma omp for private(i,j) schedule(static)
2210 for (i=0;i<numSamples;i++) {
2211 for (j=0;j<numDPPSample;j++) {
2212 local_current_value=operation(local_current_value,escript::reductionOpVector(vec,shape,data.getPointOffset(i,j),operation,initial_value));
2213
2214 }
2215 }
2216 #pragma omp critical
2217 global_current_value=operation(global_current_value,local_current_value);
2218 }
2219 return global_current_value;
2220 } else if (isTagged()) {
2221 DataTagged* leftC=dynamic_cast<DataTagged*>(m_data.get());
2222 ESYS_ASSERT(leftC!=0, "Programming error - casting to DataTagged.");
2223
2224 DataTagged& data=*leftC;
2225 DataTypes::real_t current_value=initial_value;
2226
2227 const auto& vec=data.getTypedVectorRO(typename BinaryFunction::first_argument_type(0));
2228 const DataTypes::ShapeType& shape=data.getShape();
2229 const DataTagged::DataMapType& lookup=data.getTagLookup();
2230 const std::list<int> used=data.getFunctionSpace().getListOfTagsSTL();
2231 for (std::list<int>::const_iterator i=used.begin();i!=used.end();++i)
2232 {
2233 int tag=*i;
2234 if (tag==0) // check for the default tag
2235 {
2236 current_value=operation(current_value,escript::reductionOpVector(vec,shape,data.getDefaultOffset(),operation,initial_value));
2237 }
2238 else
2239 {
2240 DataTagged::DataMapType::const_iterator it=lookup.find(tag);
2241 if (it!=lookup.end())
2242 {
2243 current_value=operation(current_value,escript::reductionOpVector(vec,shape,it->second,operation,initial_value));
2244 }
2245 }
2246 }
2247 return current_value;
2248 } else if (isConstant()) {
2249 DataConstant* leftC=dynamic_cast<DataConstant*>(m_data.get());
2250 ESYS_ASSERT(leftC!=0, "Programming error - casting to DataConstant.");
2251 return escript::reductionOpVector(leftC->getTypedVectorRO(typename BinaryFunction::first_argument_type(0)),leftC->getShape(),0,operation,initial_value);
2252 } else if (isEmpty()) {
2253 throw DataException("Error - Operations (algorithm) not permitted on instances of DataEmpty.");
2254 } else if (isLazy()) {
2255 throw DataException("Error - Operations not permitted on instances of DataLazy.");
2256 } else {
2257 throw DataException("Error - Data encapsulates an unknown type.");
2258 }
2259 }
2260
2261 /**
2262 \brief
2263 Perform the given data point reduction algorithm on data and return the result.
2264 Given operation combines each element within each data point into a scalar,
2265 thus argument object is a rank n Data object, and returned object is a
2266 rank 0 Data object.
2267 Calls escript::dp_algorithm.
2268 */
2269 template <class BinaryFunction>
2270 inline
2271 Data
2272 Data::dp_algorithm(BinaryFunction operation, DataTypes::real_t initial_value) const
2273 {
2274 if (isEmpty()) {
2275 throw DataException("Error - Operations (dp_algorithm) not permitted on instances of DataEmpty.");
2276 }
2277 else if (isExpanded()) {
2278 Data result(0,DataTypes::ShapeType(),getFunctionSpace(),isExpanded());
2279 DataExpanded* dataE=dynamic_cast<DataExpanded*>(m_data.get());
2280 DataExpanded* resultE=dynamic_cast<DataExpanded*>(result.m_data.get());
2281 ESYS_ASSERT(dataE!=0, "Programming error - casting data to DataExpanded.");
2282 ESYS_ASSERT(resultE!=0, "Programming error - casting result to DataExpanded.");
2283
2284
2285
2286 int i,j;
2287 int numSamples=dataE->getNumSamples();
2288 int numDPPSample=dataE->getNumDPPSample();
2289 // DataArrayView dataView=data.getPointDataView();
2290 // DataArrayView resultView=result.getPointDataView();
2291 const auto& dataVec=dataE->getTypedVectorRO(initial_value);
2292 const DataTypes::ShapeType& shape=dataE->getShape();
2293 auto& resultVec=resultE->getTypedVectorRW(initial_value);
2294 // perform the operation on each data-point and assign
2295 // this to the corresponding element in result
2296 #pragma omp parallel for private(i,j) schedule(static)
2297 for (i=0;i<numSamples;i++) {
2298 for (j=0;j<numDPPSample;j++) {
2299 resultVec[resultE->getPointOffset(i,j)] =
2300 escript::reductionOpVector(dataVec, shape, dataE->getPointOffset(i,j),operation,initial_value);
2301
2302 }
2303 }
2304 //escript::dp_algorithm(*dataE,*resultE,operation,initial_value);
2305 return result;
2306 }
2307 else if (isTagged()) {
2308 DataTagged* dataT=dynamic_cast<DataTagged*>(m_data.get());
2309 ESYS_ASSERT(dataT!=0, "Programming error - casting data to DataTagged.");
2310 DataTypes::RealVectorType defval(1);
2311 defval[0]=0;
2312 DataTagged* resultT=new DataTagged(getFunctionSpace(), DataTypes::scalarShape, defval, dataT);
2313
2314
2315 const DataTypes::ShapeType& shape=dataT->getShape();
2316 const auto& vec=dataT->getTypedVectorRO(initial_value);
2317 const DataTagged::DataMapType& lookup=dataT->getTagLookup();
2318 for (DataTagged::DataMapType::const_iterator i=lookup.begin(); i!=lookup.end(); i++) {
2319 resultT->getDataByTagRW(i->first,0) =
2320 escript::reductionOpVector(vec,shape,dataT->getOffsetForTag(i->first),operation,initial_value);
2321 }
2322 resultT->getTypedVectorRW(initial_value)[resultT->getDefaultOffset()] = escript::reductionOpVector(dataT->getTypedVectorRO(initial_value),dataT->getShape(),dataT->getDefaultOffset(),operation,initial_value);
2323
2324
2325
2326
2327 //escript::dp_algorithm(*dataT,*resultT,operation,initial_value);
2328 return Data(resultT); // note: the Data object now owns the resultT pointer
2329 }
2330 else if (isConstant()) {
2331 Data result(0,DataTypes::ShapeType(),getFunctionSpace(),isExpanded());
2332 DataConstant* dataC=dynamic_cast<DataConstant*>(m_data.get());
2333 DataConstant* resultC=dynamic_cast<DataConstant*>(result.m_data.get());
2334 ESYS_ASSERT(dataC!=0, "Programming error - casting data to DataConstant.");
2335 ESYS_ASSERT(resultC!=0, "Programming error - casting result to DataConstant.");
2336
2337 DataConstant& data=*dataC;
2338 resultC->getTypedVectorRW(initial_value)[0] =
2339 escript::reductionOpVector(data.getTypedVectorRO(initial_value),data.getShape(),0,operation,initial_value);
2340
2341 //escript::dp_algorithm(*dataC,*resultC,operation,initial_value);
2342 return result;
2343 } else if (isLazy()) {
2344 throw DataException("Error - Operations not permitted on instances of DataLazy.");
2345 } else {
2346 throw DataException("Error - Data encapsulates an unknown type.");
2347 }
2348 }
2349
2350
2351 /**
2352 \brief
2353 Compute a tensor operation with two Data objects
2354 \param arg_0 - Input - Data object
2355 \param arg_1 - Input - Data object
2356 \param operation - Input - Binary op functor
2357 */
2358 Data
2359 C_TensorBinaryOperation(Data const &arg_0,
2360 Data const &arg_1,
2361 ES_optype operation);
2362
2363
2364 Data
2365 C_TensorUnaryOperation(Data const &arg_0,
2366 escript::ES_optype operation,
2367 DataTypes::real_t tol=0);
2368
2369 } // namespace escript
2370
2371 #endif // __ESCRIPT_DATA_H__
2372

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26