00001 #ifndef NPSTAT_ARRAYND_HH_
00002 #define NPSTAT_ARRAYND_HH_
00003
00014 #include <cassert>
00015
00016 #include "Alignment/Geners/interface/ClassId.hh"
00017
00018 #include "JetMETCorrections/InterpolationTables/interface/SimpleFunctors.h"
00019 #include "JetMETCorrections/InterpolationTables/interface/ArrayRange.h"
00020 #include "JetMETCorrections/InterpolationTables/interface/AbsArrayProjector.h"
00021 #include "JetMETCorrections/InterpolationTables/interface/AbsVisitor.h"
00022 #include "JetMETCorrections/InterpolationTables/interface/PreciseType.h"
00023 #include "JetMETCorrections/InterpolationTables/interface/ProperDblFromCmpl.h"
00024
00025 namespace npstat {
00047 template <typename Numeric, unsigned StackLen=1U, unsigned StackDim=10U>
00048 class ArrayND
00049 {
00050 template <typename Num2, unsigned Len2, unsigned Dim2>
00051 friend class ArrayND;
00052
00053 public:
00054 typedef Numeric value_type;
00055 typedef typename ProperDblFromCmpl<Numeric>::type proper_double;
00056
00071 ArrayND();
00072
00074
00082 explicit ArrayND(const ArrayShape& shape);
00083 ArrayND(const unsigned* shape, unsigned dim);
00085
00087 ArrayND(const ArrayND&);
00088
00097 template <typename Num2, unsigned Len2, unsigned Dim2>
00098 ArrayND(const ArrayND<Num2, Len2, Dim2>&);
00099
00104 template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
00105 ArrayND(const ArrayND<Num2, Len2, Dim2>&, Functor f);
00106
00108 template <typename Num2, unsigned Len2, unsigned Dim2>
00109 ArrayND(const ArrayND<Num2, Len2, Dim2>& from,
00110 const ArrayRange& fromRange);
00111
00113 template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
00114 ArrayND(const ArrayND<Num2, Len2, Dim2>& from,
00115 const ArrayRange& fromRange, Functor f);
00116
00128 template <typename Num2, unsigned Len2, unsigned Dim2>
00129 ArrayND(const ArrayND<Num2, Len2, Dim2>& slicedArray,
00130 const unsigned *indices, unsigned nIndices);
00131
00133 template <typename Num1, unsigned Len1, unsigned Dim1,
00134 typename Num2, unsigned Len2, unsigned Dim2>
00135 ArrayND(const ArrayND<Num1, Len1, Dim1>& a1,
00136 const ArrayND<Num2, Len2, Dim2>& a2);
00137
00139
00143 explicit ArrayND(unsigned n0);
00144 ArrayND(unsigned n0, unsigned n1);
00145 ArrayND(unsigned n0, unsigned n1, unsigned n2);
00146 ArrayND(unsigned n0, unsigned n1, unsigned n2, unsigned n3);
00147 ArrayND(unsigned n0, unsigned n1, unsigned n2, unsigned n3,
00148 unsigned n4);
00149 ArrayND(unsigned n0, unsigned n1, unsigned n2, unsigned n3,
00150 unsigned n4, unsigned n5);
00151 ArrayND(unsigned n0, unsigned n1, unsigned n2, unsigned n3,
00152 unsigned n4, unsigned n5, unsigned n6);
00153 ArrayND(unsigned n0, unsigned n1, unsigned n2, unsigned n3,
00154 unsigned n4, unsigned n5, unsigned n6, unsigned n7);
00155 ArrayND(unsigned n0, unsigned n1, unsigned n2, unsigned n3,
00156 unsigned n4, unsigned n5, unsigned n6, unsigned n7,
00157 unsigned n8);
00158 ArrayND(unsigned n0, unsigned n1, unsigned n2, unsigned n3,
00159 unsigned n4, unsigned n5, unsigned n6, unsigned n7,
00160 unsigned n8, unsigned n9);
00162
00164 ~ArrayND();
00165
00174 ArrayND& operator=(const ArrayND&);
00175
00177 template <typename Num2, unsigned Len2, unsigned Dim2>
00178 ArrayND& operator=(const ArrayND<Num2,Len2,Dim2>&);
00179
00181 template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
00182 ArrayND& assign(const ArrayND<Num2, Len2, Dim2>&, Functor f);
00183
00189 ArrayND& uninitialize();
00190
00192
00197 Numeric& value(const unsigned *index, unsigned indexLen);
00198 const Numeric& value(const unsigned *index, unsigned indexLen) const;
00200
00202
00206 Numeric& valueAt(const unsigned *index, unsigned indexLen);
00207 const Numeric& valueAt(const unsigned *index, unsigned indexLen) const;
00209
00211
00212 Numeric& linearValue(unsigned long index);
00213 const Numeric& linearValue(unsigned long index) const;
00215
00217
00218 Numeric& linearValueAt(unsigned long index);
00219 const Numeric& linearValueAt(unsigned long index) const;
00221
00223 void convertLinearIndex(unsigned long l, unsigned* index,
00224 unsigned indexLen) const;
00225
00227 unsigned long linearIndex(const unsigned* idx, unsigned idxLen) const;
00228
00229
00231 inline unsigned long length() const {return len_;}
00232
00234 inline const Numeric* data() const {return data_;}
00235
00237 inline bool isShapeKnown() const {return shapeIsKnown_;}
00238
00240 inline unsigned rank() const {return dim_;}
00241
00243 ArrayShape shape() const;
00244
00246 inline const unsigned *shapeData() const {return shape_;}
00247
00249 ArrayRange fullRange() const;
00250
00252 unsigned span(unsigned dim) const;
00253
00255 unsigned maximumSpan() const;
00256
00258 unsigned minimumSpan() const;
00259
00261 inline const unsigned long* strides() const {return strides_;}
00262
00264 bool isZero() const;
00265
00271 bool isDensity() const;
00272
00274 template <typename Num2>
00275 ArrayND& setData(const Num2* data, unsigned long dataLength);
00276
00278 template <unsigned Len2, unsigned Dim2>
00279 bool operator==(const ArrayND<Numeric,Len2,Dim2>&) const;
00280
00282 template <unsigned Len2, unsigned Dim2>
00283 bool operator!=(const ArrayND<Numeric,Len2,Dim2>&) const;
00284
00286 template <unsigned Len2, unsigned Dim2>
00287 double maxAbsDifference(const ArrayND<Numeric,Len2,Dim2>&) const;
00288
00290 ArrayND operator+() const;
00291
00293 ArrayND operator-() const;
00294
00296 template <unsigned Len2, unsigned Dim2>
00297 ArrayND operator+(const ArrayND<Numeric,Len2,Dim2>& r) const;
00298
00300 template <unsigned Len2, unsigned Dim2>
00301 ArrayND operator-(const ArrayND<Numeric,Len2,Dim2>& r) const;
00302
00304 template <typename Num2>
00305 ArrayND operator*(const Num2& r) const;
00306
00308 template <typename Num2>
00309 ArrayND operator/(const Num2& r) const;
00310
00312
00316 template <typename Num2>
00317 ArrayND& operator*=(const Num2& r);
00318
00319 template <typename Num2>
00320 ArrayND& operator/=(const Num2& r);
00321
00322 template <typename Num2, unsigned Len2, unsigned Dim2>
00323 ArrayND& operator+=(const ArrayND<Num2,Len2,Dim2>& r);
00324
00325 template <typename Num2, unsigned Len2, unsigned Dim2>
00326 ArrayND& operator-=(const ArrayND<Num2,Len2,Dim2>& r);
00328
00330 template <typename Num3, typename Num2, unsigned Len2, unsigned Dim2>
00331 ArrayND& addmul(const ArrayND<Num2,Len2,Dim2>& r, const Num3& c);
00332
00334 template <typename Num2, unsigned Len2, unsigned Dim2>
00335 ArrayND outer(const ArrayND<Num2,Len2,Dim2>& r) const;
00336
00341 ArrayND contract(unsigned pos1, unsigned pos2) const;
00342
00348 template <typename Num2, unsigned Len2, unsigned Dim2>
00349 ArrayND dot(const ArrayND<Num2,Len2,Dim2>& r) const;
00350
00365 template <typename Num2, unsigned Len2, unsigned Dim2>
00366 ArrayND marginalize(const ArrayND<Num2,Len2,Dim2>& prior,
00367 const unsigned* indexMap, unsigned mapLen) const;
00368
00370 ArrayND transpose(unsigned pos1, unsigned pos2) const;
00371
00373 ArrayND transpose() const;
00374
00381 template <typename Num2>
00382 Num2 sum() const;
00383
00388 template <typename Num2>
00389 Num2 sumsq() const;
00390
00399 template <typename Num2>
00400 ArrayND derivative(double scale=1.0) const;
00401
00406 template <typename Num2>
00407 ArrayND cdfArray(double scale=1.0) const;
00408
00413 template <typename Num2>
00414 Num2 cdfValue(const unsigned *index, unsigned indexLen) const;
00415
00427 template <typename Num2>
00428 void convertToLastDimCdf(ArrayND* sumSlice, bool useTrapezoids);
00429
00431 Numeric min() const;
00432
00434 Numeric min(unsigned *index, unsigned indexLen) const;
00435
00437 Numeric max() const;
00438
00440 Numeric max(unsigned *index, unsigned indexLen) const;
00441
00443
00451 Numeric& closest(const double *x, unsigned xDim);
00452 const Numeric& closest(const double *x, unsigned xDim) const;
00454
00463 Numeric interpolate1(const double *x, unsigned xDim) const;
00464
00471 Numeric interpolate3(const double *x, unsigned xDim) const;
00472
00482 template <class Functor>
00483 ArrayND& apply(Functor f);
00484
00492 template <class Functor>
00493 ArrayND& scanInPlace(Functor f);
00494
00496 ArrayND& constFill(Numeric c);
00497
00499 ArrayND& clear();
00500
00507 ArrayND& linearFill(const double* coeff, unsigned coeffLen, double c);
00508
00515 template <class Functor>
00516 ArrayND& functorFill(Functor f);
00517
00524 ArrayND& makeUnit();
00525
00527 ArrayND& makeNonNegative();
00528
00538 unsigned makeCopulaSteps(double tolerance, unsigned maxIterations);
00539
00544 template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
00545 void jointScan(ArrayND<Num2, Len2, Dim2>& other, Functor binaryFunct);
00546
00549 template <typename Num2, unsigned Len2, unsigned Dim2>
00550 inline ArrayND& inPlaceMul(const ArrayND<Num2,Len2,Dim2>& r)
00551 {
00552 jointScan(const_cast<ArrayND<Num2,Len2,Dim2>&>(r),
00553 multeq_left<Numeric,Num2>());
00554 return *this;
00555 }
00556
00574 template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
00575 void jointSubrangeScan(ArrayND<Num2, Len2, Dim2>& other,
00576 const unsigned* thisCorner,
00577 const unsigned* range,
00578 const unsigned* otherCorner,
00579 unsigned arrLen,
00580 Functor binaryFunct);
00581
00587 template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
00588 void dualCircularScan(ArrayND<Num2, Len2, Dim2>& other,
00589 const unsigned* thisCorner,
00590 const unsigned* range,
00591 const unsigned* otherCorner,
00592 unsigned arrLen,
00593 Functor binaryFunct);
00594
00599 template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
00600 void flatCircularScan(ArrayND<Num2, Len2, Dim2>& other,
00601 const unsigned* thisCorner,
00602 const unsigned* range,
00603 const unsigned* otherCorner,
00604 unsigned arrLen,
00605 Functor binaryFunct);
00606
00611 template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
00612 void circularFlatScan(ArrayND<Num2, Len2, Dim2>& other,
00613 const unsigned* thisCorner,
00614 const unsigned* range,
00615 const unsigned* otherCorner,
00616 unsigned arrLen,
00617 Functor binaryFunct);
00618
00625 template <typename Num2, typename Integer>
00626 void processSubrange(AbsArrayProjector<Numeric,Num2>& f,
00627 const BoxND<Integer>& subrange) const;
00628
00636 template <typename Num2, unsigned Len2, unsigned Dim2>
00637 void exportSubrange(const unsigned* fromCorner, unsigned lenCorner,
00638 ArrayND<Num2, Len2, Dim2>* dest) const;
00639
00641 template <typename Num2, unsigned Len2, unsigned Dim2>
00642 void importSubrange(const unsigned* fromCorner, unsigned lenCorner,
00643 const ArrayND<Num2, Len2, Dim2>& from);
00644
00650 template <typename Num2, unsigned Len2, unsigned Dim2>
00651 bool isClose(const ArrayND<Num2,Len2,Dim2>& r, double eps) const;
00652
00654 bool isCompatible(const ArrayShape& shape) const;
00655
00660 template <typename Num2, unsigned Len2, unsigned Dim2>
00661 bool isShapeCompatible(const ArrayND<Num2,Len2,Dim2>& r) const;
00662
00672 template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
00673 void jointSliceScan(ArrayND<Num2,Len2,Dim2>& slice,
00674 const unsigned *fixedIndices,
00675 const unsigned *fixedIndexValues,
00676 unsigned nFixedIndices,
00677 Functor binaryFunct);
00678
00680 template <typename Num2, unsigned Len2, unsigned Dim2>
00681 inline void exportSlice(ArrayND<Num2,Len2,Dim2>* slice,
00682 const unsigned *fixedIndices,
00683 const unsigned *fixedIndexValues,
00684 unsigned nFixedIndices) const
00685 {
00686 assert(slice);
00687 (const_cast<ArrayND*>(this))->jointSliceScan(
00688 *slice, fixedIndices, fixedIndexValues, nFixedIndices,
00689 scast_assign_right<Numeric,Num2>());
00690 }
00691
00693 template <typename Num2, unsigned Len2, unsigned Dim2>
00694 inline void importSlice(const ArrayND<Num2,Len2,Dim2>& slice,
00695 const unsigned *fixedIndices,
00696 const unsigned *fixedIndexValues,
00697 unsigned nFixedIndices)
00698 {
00699 jointSliceScan(const_cast<ArrayND<Num2,Len2,Dim2>&>(slice),
00700 fixedIndices, fixedIndexValues, nFixedIndices,
00701 scast_assign_left<Numeric,Num2>());
00702 }
00703
00711 template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
00712 void applySlice(ArrayND<Num2,Len2,Dim2>& slice,
00713 const unsigned *fixedIndices, unsigned nFixedIndices,
00714 Functor binaryFunct);
00715
00720 template <typename Num2, unsigned Len2, unsigned Dim2>
00721 inline ArrayND& multiplyBySlice(const ArrayND<Num2,Len2,Dim2>& slice,
00722 const unsigned *fixedIndices,
00723 unsigned nFixedIndices)
00724 {
00725 applySlice(const_cast<ArrayND<Num2,Len2,Dim2>&>(slice),
00726 fixedIndices, nFixedIndices,
00727 multeq_left<Numeric,Num2>());
00728 return *this;
00729 }
00730
00732
00739 template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
00740 void project(ArrayND<Num2,Len2,Dim2>* projection,
00741 AbsArrayProjector<Numeric,Num3>& projector,
00742 const unsigned *projectedIndices,
00743 unsigned nProjectedIndices) const;
00744
00745 template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
00746 void project(ArrayND<Num2,Len2,Dim2>* projection,
00747 AbsVisitor<Numeric,Num3>& projector,
00748 const unsigned *projectedIndices,
00749 unsigned nProjectedIndices) const;
00751
00753
00758 template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
00759 void addToProjection(ArrayND<Num2,Len2,Dim2>* projection,
00760 AbsArrayProjector<Numeric,Num3>& projector,
00761 const unsigned *projectedIndices,
00762 unsigned nProjectedIndices) const;
00763
00764 template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
00765 void subtractFromProjection(ArrayND<Num2,Len2,Dim2>* projection,
00766 AbsArrayProjector<Numeric,Num3>& projector,
00767 const unsigned *projectedIndices,
00768 unsigned nProjectedIndices) const;
00769
00770 template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
00771 void addToProjection(ArrayND<Num2,Len2,Dim2>* projection,
00772 AbsVisitor<Numeric,Num3>& projector,
00773 const unsigned *projectedIndices,
00774 unsigned nProjectedIndices) const;
00775
00776 template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
00777 void subtractFromProjection(ArrayND<Num2,Len2,Dim2>* projection,
00778 AbsVisitor<Numeric,Num3>& projector,
00779 const unsigned *projectedIndices,
00780 unsigned nProjectedIndices) const;
00782
00788 template <typename Num2, unsigned Len2, unsigned Dim2>
00789 void rotate(const unsigned* shifts, unsigned lenShifts,
00790 ArrayND<Num2, Len2, Dim2>* rotated) const;
00791
00797 template <typename Num2, unsigned Len2, unsigned Dim2>
00798 void multiMirror(ArrayND<Num2, Len2, Dim2>* out) const;
00799
00801
00805 Numeric& operator()();
00806 const Numeric& operator()() const;
00807
00808 Numeric& operator()(unsigned i0);
00809 const Numeric& operator()(unsigned i0) const;
00810
00811 Numeric& operator()(unsigned i0, unsigned i1);
00812 const Numeric& operator()(unsigned i0, unsigned i1) const;
00813
00814 Numeric& operator()(unsigned i0, unsigned i1, unsigned i2);
00815 const Numeric& operator()(unsigned i0, unsigned i1, unsigned i2) const;
00816
00817 Numeric& operator()(unsigned i0, unsigned i1,
00818 unsigned i2, unsigned i3);
00819 const Numeric& operator()(unsigned i0, unsigned i1,
00820 unsigned i2, unsigned i3) const;
00821
00822 Numeric& operator()(unsigned i0, unsigned i1,
00823 unsigned i2, unsigned i3, unsigned i4);
00824 const Numeric& operator()(unsigned i0, unsigned i1,
00825 unsigned i2, unsigned i3, unsigned i4) const;
00826
00827 Numeric& operator()(unsigned i0, unsigned i1, unsigned i2,
00828 unsigned i3, unsigned i4, unsigned i5);
00829 const Numeric& operator()(unsigned i0, unsigned i1, unsigned i2,
00830 unsigned i3, unsigned i4, unsigned i5) const;
00831
00832 Numeric& operator()(unsigned i0, unsigned i1, unsigned i2,
00833 unsigned i3, unsigned i4, unsigned i5,
00834 unsigned i6);
00835 const Numeric& operator()(unsigned i0, unsigned i1, unsigned i2,
00836 unsigned i3, unsigned i4, unsigned i5,
00837 unsigned i6) const;
00838
00839 Numeric& operator()(unsigned i0, unsigned i1, unsigned i2,
00840 unsigned i3, unsigned i4, unsigned i5,
00841 unsigned i6, unsigned i7);
00842 const Numeric& operator()(unsigned i0, unsigned i1, unsigned i2,
00843 unsigned i3, unsigned i4, unsigned i5,
00844 unsigned i6, unsigned i7) const;
00845
00846 Numeric& operator()(unsigned i0, unsigned i1, unsigned i2,
00847 unsigned i3, unsigned i4, unsigned i5,
00848 unsigned i6, unsigned i7, unsigned i8);
00849 const Numeric& operator()(unsigned i0, unsigned i1, unsigned i2,
00850 unsigned i3, unsigned i4, unsigned i5,
00851 unsigned i6, unsigned i7, unsigned i8) const;
00852
00853 Numeric& operator()(unsigned i0, unsigned i1, unsigned i2,
00854 unsigned i3, unsigned i4, unsigned i5,
00855 unsigned i6, unsigned i7, unsigned i8,
00856 unsigned i9);
00857 const Numeric& operator()(unsigned i0, unsigned i1, unsigned i2,
00858 unsigned i3, unsigned i4, unsigned i5,
00859 unsigned i6, unsigned i7, unsigned i8,
00860 unsigned i9) const;
00862
00864
00868 Numeric& at();
00869 const Numeric& at() const;
00870
00871 Numeric& at(unsigned i0);
00872 const Numeric& at(unsigned i0) const;
00873
00874 Numeric& at(unsigned i0, unsigned i1);
00875 const Numeric& at(unsigned i0, unsigned i1) const;
00876
00877 Numeric& at(unsigned i0, unsigned i1, unsigned i2);
00878 const Numeric& at(unsigned i0, unsigned i1, unsigned i2) const;
00879
00880 Numeric& at(unsigned i0, unsigned i1,
00881 unsigned i2, unsigned i3);
00882 const Numeric& at(unsigned i0, unsigned i1,
00883 unsigned i2, unsigned i3) const;
00884
00885 Numeric& at(unsigned i0, unsigned i1,
00886 unsigned i2, unsigned i3, unsigned i4);
00887 const Numeric& at(unsigned i0, unsigned i1,
00888 unsigned i2, unsigned i3, unsigned i4) const;
00889
00890 Numeric& at(unsigned i0, unsigned i1, unsigned i2,
00891 unsigned i3, unsigned i4, unsigned i5);
00892 const Numeric& at(unsigned i0, unsigned i1, unsigned i2,
00893 unsigned i3, unsigned i4, unsigned i5) const;
00894
00895 Numeric& at(unsigned i0, unsigned i1, unsigned i2,
00896 unsigned i3, unsigned i4, unsigned i5,
00897 unsigned i6);
00898 const Numeric& at(unsigned i0, unsigned i1, unsigned i2,
00899 unsigned i3, unsigned i4, unsigned i5,
00900 unsigned i6) const;
00901
00902 Numeric& at(unsigned i0, unsigned i1, unsigned i2,
00903 unsigned i3, unsigned i4, unsigned i5,
00904 unsigned i6, unsigned i7);
00905 const Numeric& at(unsigned i0, unsigned i1, unsigned i2,
00906 unsigned i3, unsigned i4, unsigned i5,
00907 unsigned i6, unsigned i7) const;
00908
00909 Numeric& at(unsigned i0, unsigned i1, unsigned i2,
00910 unsigned i3, unsigned i4, unsigned i5,
00911 unsigned i6, unsigned i7, unsigned i8);
00912 const Numeric& at(unsigned i0, unsigned i1, unsigned i2,
00913 unsigned i3, unsigned i4, unsigned i5,
00914 unsigned i6, unsigned i7, unsigned i8) const;
00915
00916 Numeric& at(unsigned i0, unsigned i1, unsigned i2,
00917 unsigned i3, unsigned i4, unsigned i5,
00918 unsigned i6, unsigned i7, unsigned i8,
00919 unsigned i9);
00920 const Numeric& at(unsigned i0, unsigned i1, unsigned i2,
00921 unsigned i3, unsigned i4, unsigned i5,
00922 unsigned i6, unsigned i7, unsigned i8,
00923 unsigned i9) const;
00925
00927
00931 Numeric& cl();
00932 const Numeric& cl() const;
00933
00934 Numeric& cl(double x0);
00935 const Numeric& cl(double x0) const;
00936
00937 Numeric& cl(double x0, double x1);
00938 const Numeric& cl(double x0, double x1) const;
00939
00940 Numeric& cl(double x0, double x1, double x2);
00941 const Numeric& cl(double x0, double x1, double x2) const;
00942
00943 Numeric& cl(double x0, double x1,
00944 double x2, double x3);
00945 const Numeric& cl(double x0, double x1,
00946 double x2, double x3) const;
00947
00948 Numeric& cl(double x0, double x1,
00949 double x2, double x3, double x4);
00950 const Numeric& cl(double x0, double x1,
00951 double x2, double x3, double x4) const;
00952
00953 Numeric& cl(double x0, double x1, double x2,
00954 double x3, double x4, double x5);
00955 const Numeric& cl(double x0, double x1, double x2,
00956 double x3, double x4, double x5) const;
00957
00958 Numeric& cl(double x0, double x1, double x2,
00959 double x3, double x4, double x5,
00960 double x6);
00961 const Numeric& cl(double x0, double x1, double x2,
00962 double x3, double x4, double x5,
00963 double x6) const;
00964
00965 Numeric& cl(double x0, double x1, double x2,
00966 double x3, double x4, double x5,
00967 double x6, double x7);
00968 const Numeric& cl(double x0, double x1, double x2,
00969 double x3, double x4, double x5,
00970 double x6, double x7) const;
00971
00972 Numeric& cl(double x0, double x1, double x2,
00973 double x3, double x4, double x5,
00974 double x6, double x7, double x8);
00975 const Numeric& cl(double x0, double x1, double x2,
00976 double x3, double x4, double x5,
00977 double x6, double x7, double x8) const;
00978
00979 Numeric& cl(double x0, double x1, double x2,
00980 double x3, double x4, double x5,
00981 double x6, double x7, double x8,
00982 double x9);
00983 const Numeric& cl(double x0, double x1, double x2,
00984 double x3, double x4, double x5,
00985 double x6, double x7, double x8,
00986 double x9) const;
00988
00990
00991 inline gs::ClassId classId() const {return gs::ClassId(*this);}
00992 bool write(std::ostream& of) const;
00994
00995 static const char* classname();
00996 static inline unsigned version() {return 1;}
00997 static void restore(const gs::ClassId& id, std::istream& in,
00998 ArrayND* array);
00999 private:
01000 Numeric localData_[StackLen];
01001 Numeric* data_;
01002
01003 unsigned long localStrides_[StackDim];
01004 unsigned long *strides_;
01005
01006 unsigned localShape_[StackDim];
01007 unsigned *shape_;
01008
01009 unsigned long len_;
01010 unsigned dim_;
01011
01012 bool shapeIsKnown_;
01013
01014
01015 void buildFromShapePtr(const unsigned*, unsigned);
01016
01017
01018 void buildStrides();
01019
01020
01021 void linearFillLoop(unsigned level, double s0,
01022 unsigned long idx, double shift,
01023 const double* coeffs);
01024
01025
01026 template <class Functor>
01027 void functorFillLoop(unsigned level, unsigned long idx,
01028 Functor f, unsigned* farg);
01029
01030
01031 Numeric interpolateLoop(unsigned level, const double *x,
01032 const Numeric* base) const;
01033
01034
01035 template <typename Num1, unsigned Len1, unsigned Dim1,
01036 typename Num2, unsigned Len2, unsigned Dim2>
01037 void outerProductLoop(unsigned level, unsigned long idx0,
01038 unsigned long idx1, unsigned long idx2,
01039 const ArrayND<Num1, Len1, Dim1>& a1,
01040 const ArrayND<Num2, Len2, Dim2>& a2);
01041
01042
01043 void contractLoop(unsigned thisLevel, unsigned resLevel,
01044 unsigned pos1, unsigned pos2,
01045 unsigned long idxThis, unsigned long idxRes,
01046 ArrayND& result) const;
01047
01048
01049 void transposeLoop(unsigned level, unsigned pos1, unsigned pos2,
01050 unsigned long idxThis, unsigned long idxRes,
01051 ArrayND& result) const;
01052
01053
01054 template <typename Num2, unsigned Len2, unsigned Dim2>
01055 void dotProductLoop(unsigned level, unsigned long idx0,
01056 unsigned long idx1, unsigned long idx2,
01057 const ArrayND<Num2, Len2, Dim2>& r,
01058 ArrayND& result) const;
01059
01060
01061 template <typename Num2, unsigned Len2, unsigned Dim2>
01062 Numeric marginalizeInnerLoop(unsigned long idx,
01063 unsigned levelPr, unsigned long idxPr,
01064 const ArrayND<Num2,Len2,Dim2>& prior,
01065 const unsigned* indexMap) const;
01066 template <typename Num2, unsigned Len2, unsigned Dim2>
01067 void marginalizeLoop(unsigned level, unsigned long idx,
01068 unsigned levelRes, unsigned long idxRes,
01069 const ArrayND<Num2,Len2,Dim2>& prior,
01070 const unsigned* indexMap, ArrayND& res) const;
01071
01072
01073
01074 template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
01075 void copyRangeLoopFunct(unsigned level, unsigned long idx0,
01076 unsigned long idx1,
01077 const ArrayND<Num2, Len2, Dim2>& r,
01078 const ArrayRange& range, Functor f);
01079
01080
01081
01082
01083 template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
01084 void commonSubrangeLoop(unsigned level, unsigned long idx0,
01085 unsigned long idx1,
01086 const unsigned* thisCorner,
01087 const unsigned* range,
01088 const unsigned* otherCorner,
01089 ArrayND<Num2, Len2, Dim2>& other,
01090 Functor binaryFunct);
01091
01092
01093
01094 template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
01095 void dualCircularLoop(unsigned level, unsigned long idx0,
01096 unsigned long idx1,
01097 const unsigned* thisCorner,
01098 const unsigned* range,
01099 const unsigned* otherCorner,
01100 ArrayND<Num2, Len2, Dim2>& other,
01101 Functor binaryFunct);
01102
01103
01104
01105
01106
01107 template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
01108 void flatCircularLoop(unsigned level, unsigned long idx0,
01109 unsigned long idx1,
01110 const unsigned* thisCorner,
01111 const unsigned* range,
01112 const unsigned* otherCorner,
01113 ArrayND<Num2, Len2, Dim2>& other,
01114 Functor binaryFunct);
01115
01116
01117
01118
01119 template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
01120 void circularFlatLoop(unsigned level, unsigned long idx0,
01121 unsigned long idx1,
01122 const unsigned* thisCorner,
01123 const unsigned* range,
01124 const unsigned* otherCorner,
01125 ArrayND<Num2, Len2, Dim2>& other,
01126 Functor binaryFunct);
01127
01128
01129 template <typename Num2, unsigned Len2, unsigned Dim2>
01130 unsigned long verifySliceCompatibility(
01131 const ArrayND<Num2,Len2,Dim2>& slice,
01132 const unsigned *fixedIndices,
01133 const unsigned *fixedIndexValues,
01134 unsigned nFixedIndices) const;
01135
01136
01137 template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
01138 void jointSliceLoop(unsigned level, unsigned long idx0,
01139 unsigned level1, unsigned long idx1,
01140 ArrayND<Num2,Len2,Dim2>& slice,
01141 const unsigned *fixedIndices,
01142 const unsigned *fixedIndexValues,
01143 unsigned nFixedIndices, Functor binaryFunctor);
01144
01145
01146 template <typename Num2, class Functor>
01147 void scaleBySliceInnerLoop(unsigned level, unsigned long idx0,
01148 Num2& scale,
01149 const unsigned* projectedIndices,
01150 unsigned nProjectedIndices,
01151 Functor binaryFunct);
01152
01153 template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
01154 void scaleBySliceLoop(unsigned level, unsigned long idx0,
01155 unsigned level1, unsigned long idx1,
01156 ArrayND<Num2,Len2,Dim2>& slice,
01157 const unsigned *fixedIndices,
01158 unsigned nFixedIndices,
01159 Functor binaryFunct);
01160
01161
01162 template <typename Num2>
01163 void projectInnerLoop(unsigned level, unsigned long idx0,
01164 unsigned* currentIndex,
01165 AbsArrayProjector<Numeric,Num2>& projector,
01166 const unsigned* projectedIndices,
01167 unsigned nProjectedIndices) const;
01168
01169 template <typename Num2, unsigned Len2, unsigned Dim2,
01170 typename Num3, class Op>
01171 void projectLoop(unsigned level, unsigned long idx0,
01172 unsigned level1, unsigned long idx1,
01173 unsigned* currentIndex,
01174 ArrayND<Num2,Len2,Dim2>* projection,
01175 AbsArrayProjector<Numeric,Num3>& projector,
01176 const unsigned* projectedIndices,
01177 unsigned nProjectedIndices, Op fcn) const;
01178
01179
01180
01181
01182
01183
01184
01185
01186
01187
01188
01189 template <typename Num2, unsigned Len2, unsigned Dim2,
01190 typename Num3, class Op>
01191 void projectLoop2(unsigned level, unsigned long idx0,
01192 unsigned level1, unsigned long idx1,
01193 ArrayND<Num2,Len2,Dim2>* projection,
01194 AbsVisitor<Numeric,Num3>& projector,
01195 const unsigned* projectedIndices,
01196 unsigned nProjectedIndices, Op fcn) const;
01197
01198 template <typename Num2>
01199 void projectInnerLoop2(unsigned level, unsigned long idx0,
01200 AbsVisitor<Numeric,Num2>& projector,
01201 const unsigned* projectedIndices,
01202 unsigned nProjectedIndices) const;
01203
01204 template <typename Num2, typename Integer>
01205 void processSubrangeLoop(unsigned level, unsigned long idx0,
01206 unsigned* currentIndex,
01207 AbsArrayProjector<Numeric,Num2>& f,
01208 const BoxND<Integer>& subrange) const;
01209
01210
01211 template <typename Accumulator>
01212 Accumulator sumBelowLoop(unsigned level, unsigned long idx0,
01213 const unsigned* limit) const;
01214
01215
01216 template <typename Accumulator>
01217 void convertToLastDimCdfLoop(ArrayND* sumSlice, unsigned level,
01218 unsigned long idx0,
01219 unsigned long idxSlice,
01220 bool useTrapezoids);
01221
01222
01223
01224 unsigned coordToIndex(double coord, unsigned idim) const;
01225
01226
01227 template <typename Num2, unsigned Len2, unsigned Dim2>
01228 void verifyProjectionCompatibility(
01229 const ArrayND<Num2,Len2,Dim2>& projection,
01230 const unsigned *projectedIndices,
01231 unsigned nProjectedIndices) const;
01232
01233 };
01234 }
01235
01236 #include <cmath>
01237 #include <climits>
01238 #include <algorithm>
01239 #include <sstream>
01240 #include "JetMETCorrections/InterpolationTables/interface/NpstatException.h"
01241
01242 #include "Alignment/Geners/interface/GenericIO.hh"
01243 #include "Alignment/Geners/interface/IOIsUnsigned.hh"
01244
01245 #include "JetMETCorrections/InterpolationTables/interface/allocators.h"
01246
01247 #include "JetMETCorrections/InterpolationTables/interface/interpolate.h"
01248 #include "JetMETCorrections/InterpolationTables/interface/absDifference.h"
01249 #include "JetMETCorrections/InterpolationTables/interface/ComplexComparesFalse.h"
01250 #include "JetMETCorrections/InterpolationTables/interface/ComplexComparesAbs.h"
01251
01252 #define me_macro_check_loop_prerequisites(METHOD, INNERLOOP) \
01253 template<typename Numeric, unsigned Len, unsigned Dim> \
01254 template <typename Num2, unsigned Len2, unsigned Dim2, class Functor> \
01255 void ArrayND<Numeric,Len,Dim>:: METHOD ( \
01256 ArrayND<Num2, Len2, Dim2>& other, \
01257 const unsigned* thisCorner, \
01258 const unsigned* range, \
01259 const unsigned* otherCorner, \
01260 const unsigned arrLen, \
01261 Functor binaryFunct) \
01262 { \
01263 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument( \
01264 "Initialize npstat::ArrayND before calling method \"" \
01265 #METHOD "\""); \
01266 if (!other.shapeIsKnown_) throw npstat::NpstatInvalidArgument( \
01267 "In npstat::ArrayND::" #METHOD ": uninitialized argument array");\
01268 if (dim_ != other.dim_) throw npstat::NpstatInvalidArgument( \
01269 "In npstat::ArrayND::" #METHOD ": incompatible argument array rank");\
01270 if (arrLen != dim_) throw npstat::NpstatInvalidArgument( \
01271 "In npstat::ArrayND::" #METHOD ": incompatible index length"); \
01272 if (dim_) \
01273 { \
01274 assert(thisCorner); \
01275 assert(range); \
01276 assert(otherCorner); \
01277 INNERLOOP (0U, 0UL, 0UL, thisCorner, range, \
01278 otherCorner, other, binaryFunct); \
01279 } \
01280 else \
01281 binaryFunct(localData_[0], other.localData_[0]); \
01282 }
01283
01284 namespace npstat {
01285 template<typename Numeric, unsigned Len, unsigned Dim>
01286 template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
01287 void ArrayND<Numeric,Len,Dim>::commonSubrangeLoop(
01288 unsigned level, unsigned long idx0,
01289 unsigned long idx1,
01290 const unsigned* thisCorner,
01291 const unsigned* range,
01292 const unsigned* otherCorner,
01293 ArrayND<Num2, Len2, Dim2>& r,
01294 Functor binaryFunct)
01295 {
01296 const unsigned imax = range[level];
01297
01298 if (level == dim_ - 1)
01299 {
01300 Numeric* left = data_ + (idx0 + thisCorner[level]);
01301 Numeric* const lMax = data_ + (idx0 + shape_[level]);
01302 Num2* right = r.data_ + (idx1 + otherCorner[level]);
01303 Num2* const rMax = r.data_ + (idx1 + r.shape_[level]);
01304
01305 for (unsigned i=0; i<imax && left<lMax && right<rMax; ++i)
01306 binaryFunct(*left++, *right++);
01307 }
01308 else
01309 {
01310 const unsigned long leftStride = strides_[level];
01311 const unsigned long leftMax = idx0 + shape_[level]*leftStride;
01312 idx0 += thisCorner[level]*leftStride;
01313 const unsigned long rightStride = r.strides_[level];
01314 const unsigned long rightMax = idx1 + r.shape_[level]*rightStride;
01315 idx1 += otherCorner[level]*rightStride;
01316
01317 for (unsigned i=0; i<imax && idx0 < leftMax && idx1 < rightMax;
01318 ++i, idx0 += leftStride, idx1 += rightStride)
01319 commonSubrangeLoop(level+1, idx0, idx1, thisCorner, range,
01320 otherCorner, r, binaryFunct);
01321 }
01322 }
01323
01324 template<typename Numeric, unsigned Len, unsigned Dim>
01325 template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
01326 void ArrayND<Numeric,Len,Dim>::dualCircularLoop(
01327 unsigned level, unsigned long idx0,
01328 unsigned long idx1,
01329 const unsigned* thisCorner,
01330 const unsigned* range,
01331 const unsigned* otherCorner,
01332 ArrayND<Num2, Len2, Dim2>& r,
01333 Functor binaryFunct)
01334 {
01335 const unsigned imax = range[level];
01336 const unsigned leftShift = thisCorner[level];
01337 const unsigned leftPeriod = shape_[level];
01338 const unsigned rightShift = otherCorner[level];
01339 const unsigned rightPeriod = r.shape_[level];
01340
01341 if (level == dim_ - 1)
01342 {
01343 Numeric* left = data_ + idx0;
01344 Num2* right = r.data_ + idx1;
01345 for (unsigned i=0; i<imax; ++i)
01346 binaryFunct(left[(i+leftShift)%leftPeriod],
01347 right[(i+rightShift)%rightPeriod]);
01348 }
01349 else
01350 {
01351 const unsigned long leftStride = strides_[level];
01352 const unsigned long rightStride = r.strides_[level];
01353 for (unsigned i=0; i<imax; ++i)
01354 dualCircularLoop(
01355 level+1, idx0+((i+leftShift)%leftPeriod)*leftStride,
01356 idx1+((i+rightShift)%rightPeriod)*rightStride,
01357 thisCorner, range, otherCorner, r, binaryFunct);
01358 }
01359 }
01360
01361 template<typename Numeric, unsigned Len, unsigned Dim>
01362 template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
01363 void ArrayND<Numeric,Len,Dim>::flatCircularLoop(
01364 unsigned level, unsigned long idx0,
01365 unsigned long idx1,
01366 const unsigned* thisCorner,
01367 const unsigned* range,
01368 const unsigned* otherCorner,
01369 ArrayND<Num2, Len2, Dim2>& r,
01370 Functor binaryFunct)
01371 {
01372 const unsigned imax = range[level];
01373 const unsigned rightShift = otherCorner[level];
01374 const unsigned rightPeriod = r.shape_[level];
01375
01376 if (level == dim_ - 1)
01377 {
01378 Numeric* left = data_ + (idx0 + thisCorner[level]);
01379 Numeric* const lMax = data_ + (idx0 + shape_[level]);
01380 Num2* right = r.data_ + idx1;
01381
01382 for (unsigned i=0; i<imax && left<lMax; ++i)
01383 binaryFunct(*left++, right[(i+rightShift)%rightPeriod]);
01384 }
01385 else
01386 {
01387 const unsigned long leftStride = strides_[level];
01388 const unsigned long leftMax = idx0 + shape_[level]*leftStride;
01389 idx0 += thisCorner[level]*leftStride;
01390 const unsigned long rightStride = r.strides_[level];
01391
01392 for (unsigned i=0; i<imax && idx0 < leftMax; ++i, idx0+=leftStride)
01393 flatCircularLoop(
01394 level+1, idx0,
01395 idx1+((i+rightShift)%rightPeriod)*rightStride,
01396 thisCorner, range, otherCorner, r, binaryFunct);
01397 }
01398 }
01399
01400 template<typename Numeric, unsigned Len, unsigned Dim>
01401 template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
01402 void ArrayND<Numeric,Len,Dim>::circularFlatLoop(
01403 unsigned level, unsigned long idx0,
01404 unsigned long idx1,
01405 const unsigned* thisCorner,
01406 const unsigned* range,
01407 const unsigned* otherCorner,
01408 ArrayND<Num2, Len2, Dim2>& r,
01409 Functor binaryFunct)
01410 {
01411 const unsigned imax = range[level];
01412 const unsigned leftShift = thisCorner[level];
01413 const unsigned leftPeriod = shape_[level];
01414
01415 if (level == dim_ - 1)
01416 {
01417 Numeric* left = data_ + idx0;
01418 Num2* right = r.data_ + (idx1 + otherCorner[level]);
01419 Num2* const rMax = r.data_ + (idx1 + r.shape_[level]);
01420
01421 for (unsigned i=0; i<imax && right<rMax; ++i)
01422 binaryFunct(left[(i+leftShift)%leftPeriod], *right++);
01423 }
01424 else
01425 {
01426 const unsigned long leftStride = strides_[level];
01427 const unsigned long rightStride = r.strides_[level];
01428 const unsigned long rightMax = idx1 + r.shape_[level]*rightStride;
01429 idx1 += otherCorner[level]*rightStride;
01430
01431 for (unsigned i=0; i<imax && idx1<rightMax; ++i, idx1+=rightStride)
01432 circularFlatLoop(
01433 level+1, idx0+((i+leftShift)%leftPeriod)*leftStride,
01434 idx1, thisCorner, range, otherCorner, r, binaryFunct);
01435 }
01436 }
01437
01438 me_macro_check_loop_prerequisites(jointSubrangeScan, commonSubrangeLoop)
01439 me_macro_check_loop_prerequisites(dualCircularScan, dualCircularLoop)
01440 me_macro_check_loop_prerequisites(flatCircularScan, flatCircularLoop)
01441 me_macro_check_loop_prerequisites(circularFlatScan, circularFlatLoop)
01442
01443 template <typename Numeric, unsigned StackLen, unsigned StackDim>
01444 template <typename Num2, unsigned Len2, unsigned Dim2>
01445 Numeric ArrayND<Numeric,StackLen,StackDim>::marginalizeInnerLoop(
01446 unsigned long idx, const unsigned levelPr, unsigned long idxPr,
01447 const ArrayND<Num2,Len2,Dim2>& prior,
01448 const unsigned* indexMap) const
01449 {
01450 Numeric sum = Numeric();
01451 const unsigned long myStride = strides_[indexMap[levelPr]];
01452 const unsigned imax = prior.shape_[levelPr];
01453 assert(imax == shape_[indexMap[levelPr]]);
01454 if (levelPr == prior.dim_ - 1)
01455 {
01456 for (unsigned i=0; i<imax; ++i)
01457 sum += data_[idx+i*myStride]*prior.data_[idxPr++];
01458 }
01459 else
01460 {
01461 const unsigned long priorStride = prior.strides_[levelPr];
01462 for (unsigned i=0; i<imax; ++i)
01463 {
01464 sum += marginalizeInnerLoop(idx, levelPr+1U, idxPr,
01465 prior, indexMap);
01466 idx += myStride;
01467 idxPr += priorStride;
01468 }
01469 }
01470 return sum;
01471 }
01472
01473 template <typename Numeric, unsigned StackLen, unsigned StackDim>
01474 template <typename Num2, unsigned Len2, unsigned Dim2>
01475 void ArrayND<Numeric,StackLen,StackDim>::marginalizeLoop(
01476 const unsigned level, unsigned long idx,
01477 const unsigned levelRes, unsigned long idxRes,
01478 const ArrayND<Num2,Len2,Dim2>& prior,
01479 const unsigned* indexMap, ArrayND& result) const
01480 {
01481 if (level == dim_)
01482 {
01483 const Numeric res = marginalizeInnerLoop(
01484 idx, 0U, 0UL, prior, indexMap);
01485 if (result.dim_)
01486 result.data_[idxRes] = res;
01487 else
01488 result.localData_[0] = res;
01489 }
01490 else
01491 {
01492
01493 bool mapped = false;
01494 for (unsigned i=0; i<prior.dim_; ++i)
01495 if (level == indexMap[i])
01496 {
01497 mapped = true;
01498 break;
01499 }
01500 if (mapped)
01501 marginalizeLoop(level+1U, idx, levelRes, idxRes,
01502 prior, indexMap, result);
01503 else
01504 {
01505 const unsigned imax = shape_[level];
01506 const unsigned long myStride = strides_[level];
01507 const unsigned long resStride = result.strides_[levelRes];
01508 for (unsigned i=0; i<imax; ++i)
01509 {
01510 marginalizeLoop(level+1U, idx, levelRes+1U, idxRes,
01511 prior, indexMap, result);
01512 idx += myStride;
01513 idxRes += resStride;
01514 }
01515 }
01516 }
01517 }
01518
01519 template<typename Numeric, unsigned StackLen, unsigned StackDim>
01520 template <typename Num2, unsigned Len2, unsigned Dim2>
01521 ArrayND<Numeric,StackLen,StackDim>
01522 ArrayND<Numeric,StackLen,StackDim>::marginalize(
01523 const ArrayND<Num2,Len2,Dim2>& prior,
01524 const unsigned* indexMap, const unsigned mapLen) const
01525 {
01526 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
01527 "Initialize npstat::ArrayND before calling method \"marginalize\"");
01528 if (!(prior.dim_ && prior.dim_ <= dim_)) throw npstat::NpstatInvalidArgument(
01529 "In npstat::ArrayND::marginalize: incompatible argument array rank");
01530 const unsigned resultDim = dim_ - prior.dim_;
01531
01532
01533 if (mapLen != prior.dim_) throw npstat::NpstatInvalidArgument(
01534 "In npstat::ArrayND::marginalize: incompatible index map length");
01535 assert(indexMap);
01536 for (unsigned i=0; i<mapLen; ++i)
01537 {
01538 const unsigned thisInd = indexMap[i];
01539 if (shape_[thisInd] != prior.shape_[i]) throw npstat::NpstatInvalidArgument(
01540 "In npstat::ArrayND::marginalize: "
01541 "incompatible argument array dimensions");
01542 if (thisInd >= dim_) throw npstat::NpstatOutOfRange(
01543 "In npstat::ArrayND::marginalize: index map entry out of range");
01544 for (unsigned j=0; j<i; ++j)
01545 if (indexMap[j] == thisInd) throw npstat::NpstatInvalidArgument(
01546 "In npstat::ArrayND::marginalize: "
01547 "duplicate entry in the index map");
01548 }
01549
01550
01551 ArrayShape newShape;
01552 newShape.reserve(resultDim);
01553 for (unsigned i=0; i<dim_; ++i)
01554 {
01555 bool mapped = false;
01556 for (unsigned j=0; j<mapLen; ++j)
01557 if (indexMap[j] == i)
01558 {
01559 mapped = true;
01560 break;
01561 }
01562 if (!mapped)
01563 newShape.push_back(shape_[i]);
01564 }
01565
01566 ArrayND result(newShape);
01567 assert(result.dim_ == resultDim);
01568 bool calculated = false;
01569 if (resultDim == 0)
01570 {
01571 calculated = true;
01572 for (unsigned i=0; i<dim_; ++i)
01573 if (indexMap[i] != i)
01574 {
01575 calculated = false;
01576 break;
01577 }
01578 if (calculated)
01579 {
01580 Numeric sum = Numeric();
01581 for (unsigned long i=0; i<len_; ++i)
01582 sum += data_[i]*prior.data_[i];
01583 result.localData_[0] = sum;
01584 }
01585 }
01586
01587 if (!calculated)
01588 marginalizeLoop(0U, 0UL, 0U, 0UL, prior, indexMap, result);
01589
01590 return result;
01591 }
01592
01593 template<typename Numeric, unsigned Len, unsigned Dim>
01594 void ArrayND<Numeric,Len,Dim>::buildFromShapePtr(
01595 const unsigned* sizes, const unsigned dim)
01596 {
01597 dim_ = dim;
01598 if (dim_)
01599 {
01600 assert(sizes);
01601 for (unsigned i=0; i<dim_; ++i)
01602 if (sizes[i] == 0)
01603 throw npstat::NpstatInvalidArgument(
01604 "In npstat::ArrayND::buildFromShapePtr: "
01605 "detected span of zero");
01606
01607
01608 shape_ = makeBuffer(dim_, localShape_, Dim);
01609 for (unsigned i=0; i<dim_; ++i)
01610 {
01611 shape_[i] = sizes[i];
01612 len_ *= shape_[i];
01613 }
01614
01615
01616 buildStrides();
01617
01618
01619 data_ = makeBuffer(len_, localData_, Len);
01620 }
01621 else
01622 {
01623 localData_[0] = Numeric();
01624 data_ = localData_;
01625 }
01626 }
01627
01628 template<typename Numeric, unsigned StackLen, unsigned StackDim>
01629 template <typename Num2, unsigned Len2, unsigned Dim2>
01630 ArrayND<Numeric,StackLen,StackDim>::ArrayND(
01631 const ArrayND<Num2, Len2, Dim2>& slicedArray,
01632 const unsigned *fixedIndices, const unsigned nFixedIndices)
01633 : data_(0), strides_(0), shape_(0),
01634 len_(1UL), dim_(slicedArray.dim_ - nFixedIndices),
01635 shapeIsKnown_(true)
01636 {
01637 if (nFixedIndices)
01638 {
01639 assert(fixedIndices);
01640 if (nFixedIndices > slicedArray.dim_) throw npstat::NpstatInvalidArgument(
01641 "In npstat::ArrayND slicing constructor: too many fixed indices");
01642 if (!slicedArray.shapeIsKnown_) throw npstat::NpstatInvalidArgument(
01643 "In npstat::ArrayND slicing constructor: "
01644 "uninitialized argument array");
01645
01646
01647 for (unsigned j=0; j<nFixedIndices; ++j)
01648 if (fixedIndices[j] >= slicedArray.dim_)
01649 throw npstat::NpstatOutOfRange("In npstat::ArrayND slicing "
01650 "constructor: fixed index out of range");
01651
01652
01653 shape_ = makeBuffer(dim_, localShape_, StackDim);
01654 unsigned idim = 0;
01655 for (unsigned i=0; i<slicedArray.dim_; ++i)
01656 {
01657 bool fixed = false;
01658 for (unsigned j=0; j<nFixedIndices; ++j)
01659 if (fixedIndices[j] == i)
01660 {
01661 fixed = true;
01662 break;
01663 }
01664 if (!fixed)
01665 {
01666 assert(idim < dim_);
01667 shape_[idim++] = slicedArray.shape_[i];
01668 }
01669 }
01670 assert(idim == dim_);
01671
01672 if (dim_)
01673 {
01674
01675 for (unsigned i=0; i<dim_; ++i)
01676 len_ *= shape_[i];
01677
01678
01679 buildStrides();
01680
01681
01682 data_ = makeBuffer(len_, localData_, StackLen);
01683 }
01684 else
01685 {
01686 localData_[0] = Numeric();
01687 data_ = localData_;
01688 }
01689 }
01690 else
01691 {
01692 new (this) ArrayND(slicedArray);
01693 }
01694 }
01695
01696 template<typename Numeric, unsigned StackLen, unsigned StackDim>
01697 template <typename Num2, unsigned Len2, unsigned Dim2>
01698 unsigned long ArrayND<Numeric,StackLen,StackDim>::verifySliceCompatibility(
01699 const ArrayND<Num2,Len2,Dim2>& slice,
01700 const unsigned *fixedIndices,
01701 const unsigned *fixedIndexValues,
01702 const unsigned nFixedIndices) const
01703 {
01704 if (!(nFixedIndices && nFixedIndices <= dim_))
01705 throw npstat::NpstatInvalidArgument(
01706 "In npstat::ArrayND::verifySliceCompatibility: "
01707 "invalid number of fixed indices");
01708 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
01709 "Initialize npstat::ArrayND before calling "
01710 "method \"verifySliceCompatibility\"");
01711 if (!slice.shapeIsKnown_) throw npstat::NpstatInvalidArgument(
01712 "In npstat::ArrayND::verifySliceCompatibility: "
01713 "uninitialized argument array");
01714 if (slice.dim_ != dim_ - nFixedIndices) throw npstat::NpstatInvalidArgument(
01715 "In npstat::ArrayND::verifySliceCompatibility: "
01716 "incompatible argument array rank");
01717 assert(fixedIndices);
01718 assert(fixedIndexValues);
01719
01720 for (unsigned j=0; j<nFixedIndices; ++j)
01721 if (fixedIndices[j] >= dim_) throw npstat::NpstatOutOfRange(
01722 "In npstat::ArrayND::verifySliceCompatibility: "
01723 "fixed index out of range");
01724
01725
01726 unsigned long idx = 0UL;
01727 unsigned sliceDim = 0U;
01728 for (unsigned i=0; i<dim_; ++i)
01729 {
01730 bool fixed = false;
01731 for (unsigned j=0; j<nFixedIndices; ++j)
01732 if (fixedIndices[j] == i)
01733 {
01734 fixed = true;
01735 if (fixedIndexValues[j] >= shape_[i])
01736 throw npstat::NpstatOutOfRange(
01737 "In npstat::ArrayND::verifySliceCompatibility: "
01738 "fixed index value out of range");
01739 idx += fixedIndexValues[j]*strides_[i];
01740 break;
01741 }
01742 if (!fixed)
01743 {
01744 if (shape_[i] != slice.shape_[sliceDim])
01745 throw npstat::NpstatInvalidArgument(
01746 "In npstat::ArrayND::verifySliceCompatibility: "
01747 "incompatible argument array dimensions");
01748 ++sliceDim;
01749 }
01750 }
01751 assert(sliceDim == slice.dim_);
01752 return idx;
01753 }
01754
01755 template<typename Numeric, unsigned StackLen, unsigned StackDim>
01756 template <typename Num2, unsigned Len2, unsigned Dim2, class Op>
01757 void ArrayND<Numeric,StackLen,StackDim>::jointSliceLoop(
01758 const unsigned level, const unsigned long idx0,
01759 const unsigned level1, const unsigned long idx1,
01760 ArrayND<Num2,Len2,Dim2>& slice,
01761 const unsigned *fixedIndices,
01762 const unsigned *fixedIndexValues,
01763 const unsigned nFixedIndices,
01764 Op fcn)
01765 {
01766 bool fixed = false;
01767 for (unsigned j=0; j<nFixedIndices; ++j)
01768 if (fixedIndices[j] == level)
01769 {
01770 fixed = true;
01771 break;
01772 }
01773 if (fixed)
01774 {
01775 jointSliceLoop(level+1, idx0, level1, idx1,
01776 slice, fixedIndices, fixedIndexValues,
01777 nFixedIndices, fcn);
01778 }
01779 else
01780 {
01781 const unsigned imax = shape_[level];
01782 assert(imax == slice.shape_[level1]);
01783 const unsigned long stride = strides_[level];
01784
01785 if (level1 == slice.dim_ - 1)
01786 {
01787 Num2* to = slice.data_ + idx1;
01788 for (unsigned i = 0; i<imax; ++i)
01789 fcn(data_[idx0 + i*stride], to[i]);
01790 }
01791 else
01792 {
01793 const unsigned long stride2 = slice.strides_[level1];
01794 for (unsigned i = 0; i<imax; ++i)
01795 jointSliceLoop(level+1, idx0+i*stride,
01796 level1+1, idx1+i*stride2,
01797 slice, fixedIndices, fixedIndexValues,
01798 nFixedIndices, fcn);
01799 }
01800 }
01801 }
01802
01803 template<typename Numeric, unsigned StackLen, unsigned StackDim>
01804 template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
01805 void ArrayND<Numeric,StackLen,StackDim>::jointSliceScan(
01806 ArrayND<Num2,Len2,Dim2>& slice,
01807 const unsigned *fixedIndices,
01808 const unsigned *fixedIndexValues,
01809 const unsigned nFixedIndices,
01810 Functor binaryFunct)
01811 {
01812 const unsigned long idx = verifySliceCompatibility(
01813 slice, fixedIndices, fixedIndexValues, nFixedIndices);
01814 if (slice.dim_)
01815 jointSliceLoop(0U, idx, 0U, 0UL, slice, fixedIndices,
01816 fixedIndexValues, nFixedIndices, binaryFunct);
01817 else
01818 binaryFunct(data_[idx], slice.localData_[0]);
01819 }
01820
01821 template<typename Numeric, unsigned StackLen, unsigned StackDim>
01822 template <typename Num2>
01823 void ArrayND<Numeric,StackLen,StackDim>::projectInnerLoop(
01824 const unsigned level, unsigned long idx0,
01825 unsigned* currentIndex,
01826 AbsArrayProjector<Numeric,Num2>& projector,
01827 const unsigned *projectedIndices,
01828 const unsigned nProjectedIndices) const
01829 {
01830
01831 const unsigned idx = projectedIndices[level];
01832 const unsigned imax = shape_[idx];
01833 const unsigned long stride = strides_[idx];
01834 const bool last = (level == nProjectedIndices - 1);
01835
01836 for (unsigned i = 0; i<imax; ++i)
01837 {
01838 currentIndex[idx] = i;
01839 if (last)
01840 projector.process(currentIndex, dim_, idx0, data_[idx0]);
01841 else
01842 projectInnerLoop(level+1, idx0, currentIndex, projector,
01843 projectedIndices, nProjectedIndices);
01844 idx0 += stride;
01845 }
01846 }
01847
01848 template<typename Numeric, unsigned StackLen, unsigned StackDim>
01849 template <typename Num2, unsigned Len2, unsigned Dim2,
01850 typename Num3, class Op>
01851 void ArrayND<Numeric,StackLen,StackDim>::projectLoop(
01852 const unsigned level, const unsigned long idx0,
01853 const unsigned level1, const unsigned long idx1,
01854 unsigned* currentIndex,
01855 ArrayND<Num2,Len2,Dim2>* projection,
01856 AbsArrayProjector<Numeric,Num3>& projector,
01857 const unsigned *projectedIndices,
01858 const unsigned nProjectedIndices, Op fcn) const
01859 {
01860
01861
01862
01863
01864
01865
01866
01867 if (level == dim_)
01868 {
01869 assert(level1 == projection->dim_);
01870 projector.clear();
01871 projectInnerLoop(0U, idx0, currentIndex, projector,
01872 projectedIndices, nProjectedIndices);
01873 if (projection->dim_)
01874 fcn(projection->data_[idx1], projector.result());
01875 else
01876 fcn(projection->localData_[0], projector.result());
01877 }
01878 else
01879 {
01880 bool iterated = false;
01881 for (unsigned j=0; j<nProjectedIndices; ++j)
01882 if (projectedIndices[j] == level)
01883 {
01884 iterated = true;
01885 break;
01886 }
01887 if (iterated)
01888 {
01889
01890 projectLoop(level+1, idx0, level1, idx1,
01891 currentIndex, projection, projector,
01892 projectedIndices, nProjectedIndices, fcn);
01893 }
01894 else
01895 {
01896 const unsigned imax = shape_[level];
01897 const unsigned long stride = strides_[level];
01898
01899
01900 const unsigned long stride2 = projection->strides_[level1];
01901 for (unsigned i = 0; i<imax; ++i)
01902 {
01903 currentIndex[level] = i;
01904 projectLoop(level+1, idx0+i*stride,
01905 level1+1, idx1+i*stride2,
01906 currentIndex, projection, projector,
01907 projectedIndices, nProjectedIndices, fcn);
01908 }
01909 }
01910 }
01911 }
01912
01913 template<typename Numeric, unsigned StackLen, unsigned StackDim>
01914 template <typename Num2, unsigned Len2, unsigned Dim2>
01915 void ArrayND<Numeric,StackLen,StackDim>::verifyProjectionCompatibility(
01916 const ArrayND<Num2,Len2,Dim2>& projection,
01917 const unsigned *projectedIndices,
01918 const unsigned nProjectedIndices) const
01919 {
01920 if (!(nProjectedIndices && nProjectedIndices <= dim_))
01921 throw npstat::NpstatInvalidArgument(
01922 "In npstat::ArrayND::verifyProjectionCompatibility: "
01923 "invalid number of projected indices");
01924 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
01925 "Initialize npstat::ArrayND before calling "
01926 "method \"verifyProjectionCompatibility\"");
01927 if (!projection.shapeIsKnown_) throw npstat::NpstatInvalidArgument(
01928 "In npstat::ArrayND::verifyProjectionCompatibility: "
01929 "uninitialized argument array");
01930 if (projection.dim_ != dim_ - nProjectedIndices)
01931 throw npstat::NpstatInvalidArgument(
01932 "In npstat::ArrayND::verifyProjectionCompatibility: "
01933 "incompatible argument array rank");
01934 assert(projectedIndices);
01935
01936 for (unsigned j=0; j<nProjectedIndices; ++j)
01937 if (projectedIndices[j] >= dim_) throw npstat::NpstatOutOfRange(
01938 "In npstat::ArrayND::verifyProjectionCompatibility: "
01939 "projected index out of range");
01940
01941
01942 unsigned sliceDim = 0U;
01943 for (unsigned i=0; i<dim_; ++i)
01944 {
01945 bool fixed = false;
01946 for (unsigned j=0; j<nProjectedIndices; ++j)
01947 if (projectedIndices[j] == i)
01948 {
01949 fixed = true;
01950 break;
01951 }
01952 if (!fixed)
01953 {
01954 if (shape_[i] != projection.shape_[sliceDim])
01955 throw npstat::NpstatInvalidArgument(
01956 "In npstat::ArrayND::verifyProjectionCompatibility: "
01957 "incompatible argument array dimensions");
01958 ++sliceDim;
01959 }
01960 }
01961 assert(sliceDim == projection.dim_);
01962 }
01963
01964 template<typename Numeric, unsigned StackLen, unsigned StackDim>
01965 template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
01966 void ArrayND<Numeric,StackLen,StackDim>::project(
01967 ArrayND<Num2,Len2,Dim2>* projection,
01968 AbsArrayProjector<Numeric,Num3>& projector,
01969 const unsigned *projectedIndices,
01970 const unsigned nProjectedIndices) const
01971 {
01972 assert(projection);
01973 verifyProjectionCompatibility(*projection, projectedIndices,
01974 nProjectedIndices);
01975 unsigned ibuf[StackDim];
01976 unsigned* buf = makeBuffer(dim_, ibuf, StackDim);
01977 for (unsigned i=0; i<dim_; ++i)
01978 buf[i] = 0U;
01979 projectLoop(0U, 0UL, 0U, 0UL, buf, projection,
01980 projector, projectedIndices, nProjectedIndices,
01981 scast_assign_left<Num2,Num3>());
01982 destroyBuffer(buf, ibuf);
01983 }
01984
01985 template<typename Numeric, unsigned StackLen, unsigned StackDim>
01986 template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
01987 void ArrayND<Numeric,StackLen,StackDim>::addToProjection(
01988 ArrayND<Num2,Len2,Dim2>* projection,
01989 AbsArrayProjector<Numeric,Num3>& projector,
01990 const unsigned *projectedIndices,
01991 const unsigned nProjectedIndices) const
01992 {
01993 assert(projection);
01994 verifyProjectionCompatibility(*projection, projectedIndices,
01995 nProjectedIndices);
01996 unsigned ibuf[StackDim];
01997 unsigned* buf = makeBuffer(dim_, ibuf, StackDim);
01998 for (unsigned i=0; i<dim_; ++i)
01999 buf[i] = 0U;
02000 projectLoop(0U, 0UL, 0U, 0UL, buf, projection,
02001 projector, projectedIndices, nProjectedIndices,
02002 scast_pluseq_left<Num2,Num3>());
02003 destroyBuffer(buf, ibuf);
02004 }
02005
02006 template<typename Numeric, unsigned StackLen, unsigned StackDim>
02007 template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
02008 void ArrayND<Numeric,StackLen,StackDim>::subtractFromProjection(
02009 ArrayND<Num2,Len2,Dim2>* projection,
02010 AbsArrayProjector<Numeric,Num3>& projector,
02011 const unsigned *projectedIndices,
02012 const unsigned nProjectedIndices) const
02013 {
02014 assert(projection);
02015 verifyProjectionCompatibility(*projection, projectedIndices,
02016 nProjectedIndices);
02017 unsigned ibuf[StackDim];
02018 unsigned* buf = makeBuffer(dim_, ibuf, StackDim);
02019 for (unsigned i=0; i<dim_; ++i)
02020 buf[i] = 0U;
02021 projectLoop(0U, 0UL, 0U, 0UL, buf, projection,
02022 projector, projectedIndices, nProjectedIndices,
02023 scast_minuseq_left<Num2,Num3>());
02024 destroyBuffer(buf, ibuf);
02025 }
02026
02027 template<typename Numeric, unsigned StackLen, unsigned StackDim>
02028 template <typename Num2>
02029 void ArrayND<Numeric,StackLen,StackDim>::projectInnerLoop2(
02030 const unsigned level, const unsigned long idx0,
02031 AbsVisitor<Numeric,Num2>& projector,
02032 const unsigned *projectedIndices,
02033 const unsigned nProjectedIndices) const
02034 {
02035 const unsigned idx = projectedIndices[level];
02036 const unsigned imax = shape_[idx];
02037 const unsigned long stride = strides_[idx];
02038 const bool last = (level == nProjectedIndices - 1);
02039
02040 for (unsigned i = 0; i<imax; ++i)
02041 {
02042 if (last)
02043 projector.process(data_[idx0+i*stride]);
02044 else
02045 projectInnerLoop2(level+1, idx0+i*stride, projector,
02046 projectedIndices, nProjectedIndices);
02047 }
02048 }
02049
02050 template<typename Numeric, unsigned StackLen, unsigned StackDim>
02051 template <typename Num2, unsigned Len2, unsigned Dim2,
02052 typename Num3, class Op>
02053 void ArrayND<Numeric,StackLen,StackDim>::projectLoop2(
02054 const unsigned level, const unsigned long idx0,
02055 const unsigned level1, const unsigned long idx1,
02056 ArrayND<Num2,Len2,Dim2>* projection,
02057 AbsVisitor<Numeric,Num3>& projector,
02058 const unsigned *projectedIndices,
02059 const unsigned nProjectedIndices, Op fcn) const
02060 {
02061 if (level == dim_)
02062 {
02063 assert(level1 == projection->dim_);
02064 projector.clear();
02065 projectInnerLoop2(0U, idx0, projector,
02066 projectedIndices, nProjectedIndices);
02067 if (projection->dim_)
02068 fcn(projection->data_[idx1], projector.result());
02069 else
02070 fcn(projection->localData_[0], projector.result());
02071 }
02072 else
02073 {
02074 bool fixed = false;
02075 for (unsigned j=0; j<nProjectedIndices; ++j)
02076 if (projectedIndices[j] == level)
02077 {
02078 fixed = true;
02079 break;
02080 }
02081 if (fixed)
02082 {
02083 projectLoop2(level+1, idx0, level1, idx1,
02084 projection, projector,
02085 projectedIndices, nProjectedIndices, fcn);
02086 }
02087 else
02088 {
02089 const unsigned imax = shape_[level];
02090 const unsigned long stride = strides_[level];
02091 const unsigned long stride2 = projection->strides_[level1];
02092 for (unsigned i = 0; i<imax; ++i)
02093 projectLoop2(level+1, idx0+i*stride,
02094 level1+1, idx1+i*stride2,
02095 projection, projector,
02096 projectedIndices, nProjectedIndices, fcn);
02097 }
02098 }
02099 }
02100
02101 template<typename Numeric, unsigned StackLen, unsigned StackDim>
02102 template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
02103 void ArrayND<Numeric,StackLen,StackDim>::project(
02104 ArrayND<Num2,Len2,Dim2>* projection,
02105 AbsVisitor<Numeric,Num3>& projector,
02106 const unsigned *projectedIndices,
02107 const unsigned nProjectedIndices) const
02108 {
02109 assert(projection);
02110 verifyProjectionCompatibility(*projection, projectedIndices,
02111 nProjectedIndices);
02112 projectLoop2(0U, 0UL, 0U, 0UL, projection,
02113 projector, projectedIndices, nProjectedIndices,
02114 scast_assign_left<Num2,Num3>());
02115 }
02116
02117 template<typename Numeric, unsigned StackLen, unsigned StackDim>
02118 template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
02119 void ArrayND<Numeric,StackLen,StackDim>::addToProjection(
02120 ArrayND<Num2,Len2,Dim2>* projection,
02121 AbsVisitor<Numeric,Num3>& projector,
02122 const unsigned *projectedIndices,
02123 const unsigned nProjectedIndices) const
02124 {
02125 assert(projection);
02126 verifyProjectionCompatibility(*projection, projectedIndices,
02127 nProjectedIndices);
02128 projectLoop2(0U, 0UL, 0U, 0UL, projection,
02129 projector, projectedIndices, nProjectedIndices,
02130 scast_pluseq_left<Num2,Num3>());
02131 }
02132
02133 template<typename Numeric, unsigned StackLen, unsigned StackDim>
02134 template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
02135 void ArrayND<Numeric,StackLen,StackDim>::subtractFromProjection(
02136 ArrayND<Num2,Len2,Dim2>* projection,
02137 AbsVisitor<Numeric,Num3>& projector,
02138 const unsigned *projectedIndices,
02139 const unsigned nProjectedIndices) const
02140 {
02141 assert(projection);
02142 verifyProjectionCompatibility(*projection, projectedIndices,
02143 nProjectedIndices);
02144 projectLoop2(0U, 0UL, 0U, 0UL, projection,
02145 projector, projectedIndices, nProjectedIndices,
02146 scast_minuseq_left<Num2,Num3>());
02147 }
02148
02149 template<typename Numeric, unsigned StackLen, unsigned StackDim>
02150 template <typename Num2, class Functor>
02151 void ArrayND<Numeric,StackLen,StackDim>::scaleBySliceInnerLoop(
02152 const unsigned level, const unsigned long idx0,
02153 Num2& scale, const unsigned *projectedIndices,
02154 const unsigned nProjectedIndices, Functor binaryFunct)
02155 {
02156 const unsigned idx = projectedIndices[level];
02157 const unsigned imax = shape_[idx];
02158 const unsigned long stride = strides_[idx];
02159
02160 if (level == nProjectedIndices - 1)
02161 {
02162 Numeric* data = data_ + idx0;
02163 for (unsigned i = 0; i<imax; ++i)
02164 binaryFunct(data[i*stride], scale);
02165 }
02166 else
02167 for (unsigned i = 0; i<imax; ++i)
02168 scaleBySliceInnerLoop(level+1, idx0+i*stride, scale,
02169 projectedIndices, nProjectedIndices,
02170 binaryFunct);
02171 }
02172
02173 template<typename Numeric, unsigned StackLen, unsigned StackDim>
02174 template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
02175 void ArrayND<Numeric,StackLen,StackDim>::scaleBySliceLoop(
02176 unsigned level, unsigned long idx0,
02177 unsigned level1, unsigned long idx1,
02178 ArrayND<Num2,Len2,Dim2>& slice,
02179 const unsigned *projectedIndices,
02180 const unsigned nProjectedIndices,
02181 Functor binaryFunct)
02182 {
02183 if (level == dim_)
02184 {
02185 assert(level1 == slice.dim_);
02186 Num2& scaleFactor = slice.dim_ ? slice.data_[idx1] :
02187 slice.localData_[0];
02188 scaleBySliceInnerLoop(0U, idx0, scaleFactor, projectedIndices,
02189 nProjectedIndices, binaryFunct);
02190 }
02191 else
02192 {
02193 bool fixed = false;
02194 for (unsigned j=0; j<nProjectedIndices; ++j)
02195 if (projectedIndices[j] == level)
02196 {
02197 fixed = true;
02198 break;
02199 }
02200 if (fixed)
02201 {
02202 scaleBySliceLoop(level+1, idx0, level1, idx1, slice,
02203 projectedIndices, nProjectedIndices,
02204 binaryFunct);
02205 }
02206 else
02207 {
02208 const unsigned imax = shape_[level];
02209 const unsigned long stride = strides_[level];
02210 const unsigned long stride2 = slice.strides_[level1];
02211 for (unsigned i = 0; i<imax; ++i)
02212 scaleBySliceLoop(level+1, idx0+i*stride, level1+1,
02213 idx1+i*stride2, slice, projectedIndices,
02214 nProjectedIndices, binaryFunct);
02215 }
02216 }
02217 }
02218
02219 template<typename Numeric, unsigned StackLen, unsigned StackDim>
02220 template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
02221 void ArrayND<Numeric,StackLen,StackDim>::jointScan(
02222 ArrayND<Num2, Len2, Dim2>& r, Functor binaryFunct)
02223 {
02224 if (!isShapeCompatible(r)) throw npstat::NpstatInvalidArgument(
02225 "In npstat::ArrayND::jointScan: incompatible argument array shape");
02226 if (dim_)
02227 for (unsigned long i=0; i<len_; ++i)
02228 binaryFunct(data_[i], r.data_[i]);
02229 else
02230 binaryFunct(localData_[0], r.localData_[0]);
02231 }
02232
02233 template<typename Numeric, unsigned StackLen, unsigned StackDim>
02234 template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
02235 void ArrayND<Numeric,StackLen,StackDim>::applySlice(
02236 ArrayND<Num2,Len2,Dim2>& slice,
02237 const unsigned *fixedIndices, const unsigned nFixedIndices,
02238 Functor binaryFunct)
02239 {
02240 if (nFixedIndices)
02241 {
02242 verifyProjectionCompatibility(slice, fixedIndices, nFixedIndices);
02243 if (slice.dim_ == 0U)
02244 for (unsigned long i=0; i<len_; ++i)
02245 binaryFunct(data_[i], slice.localData_[0]);
02246 else
02247 scaleBySliceLoop(0U, 0UL, 0U, 0UL, slice,
02248 fixedIndices, nFixedIndices, binaryFunct);
02249 }
02250 else
02251 jointScan(slice, binaryFunct);
02252 }
02253
02254 template<typename Numeric, unsigned StackLen, unsigned StackDim>
02255 inline bool ArrayND<Numeric,StackLen,StackDim>::isCompatible(
02256 const ArrayShape& shape) const
02257 {
02258 if (!shapeIsKnown_)
02259 return false;
02260 if (dim_ != shape.size())
02261 return false;
02262 if (dim_)
02263 {
02264 for (unsigned i=0; i<dim_; ++i)
02265 if (shape_[i] != shape[i])
02266 return false;
02267 }
02268 else
02269 assert(len_ == 1UL);
02270 return true;
02271 }
02272
02273 template<typename Numeric, unsigned StackLen, unsigned StackDim>
02274 template <typename Num2, unsigned Len2, unsigned Dim2>
02275 inline bool ArrayND<Numeric,StackLen,StackDim>::isShapeCompatible(
02276 const ArrayND<Num2,Len2,Dim2>& r) const
02277 {
02278 if (!shapeIsKnown_)
02279 return false;
02280 if (!r.shapeIsKnown_)
02281 return false;
02282 if (dim_ != r.dim_)
02283 return false;
02284 if (len_ != r.len_)
02285 return false;
02286 if (dim_)
02287 {
02288 assert(shape_);
02289 assert(r.shape_);
02290 for (unsigned i=0; i<dim_; ++i)
02291 if (shape_[i] != r.shape_[i])
02292 return false;
02293 }
02294 else
02295 assert(len_ == 1UL);
02296 return true;
02297 }
02298
02299 template<typename Numeric, unsigned Len, unsigned Dim>
02300 template<typename Num2, typename Integer>
02301 void ArrayND<Numeric,Len,Dim>::processSubrangeLoop(
02302 const unsigned level, unsigned long idx0,
02303 unsigned* currentIndex,
02304 AbsArrayProjector<Numeric,Num2>& f,
02305 const BoxND<Integer>& subrange) const
02306 {
02307
02308 const Interval<Integer>& levelRange(subrange[level]);
02309 long long int iminl = static_cast<long long int>(levelRange.min());
02310 if (iminl < 0LL) iminl = 0LL;
02311 long long int imaxl = static_cast<long long int>(levelRange.max());
02312 if (imaxl < 0LL) imaxl = 0LL;
02313
02314
02315 const unsigned imin = static_cast<unsigned>(iminl);
02316 unsigned imax = static_cast<unsigned>(imaxl);
02317 if (imax > shape_[level])
02318 imax = shape_[level];
02319
02320 if (level == dim_ - 1)
02321 {
02322 idx0 += imin;
02323 for (unsigned i=imin; i<imax; ++i, ++idx0)
02324 {
02325 currentIndex[level] = i;
02326 f.process(currentIndex, dim_, idx0, data_[idx0]);
02327 }
02328 }
02329 else
02330 {
02331 const unsigned long stride = strides_[level];
02332 idx0 += imin*stride;
02333 for (unsigned i=imin; i<imax; ++i)
02334 {
02335 currentIndex[level] = i;
02336 processSubrangeLoop(level+1U, idx0, currentIndex, f, subrange);
02337 idx0 += stride;
02338 }
02339 }
02340 }
02341
02342 template<typename Numeric, unsigned Len, unsigned StackDim>
02343 template <typename Num2, typename Integer>
02344 void ArrayND<Numeric,Len,StackDim>::processSubrange(
02345 AbsArrayProjector<Numeric,Num2>& f,
02346 const BoxND<Integer>& subrange) const
02347 {
02348 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
02349 "Initialize npstat::ArrayND before calling method \"processSubrange\"");
02350 if (!dim_) throw npstat::NpstatInvalidArgument(
02351 "npstat::ArrayND::processSubrange method "
02352 "can not be used with array of 0 rank");
02353 if (dim_ != subrange.dim()) throw npstat::NpstatInvalidArgument(
02354 "In npstat::ArrayND::processSubrange: incompatible subrange rank");
02355 unsigned ibuf[StackDim];
02356 unsigned* buf = makeBuffer(dim_, ibuf, StackDim);
02357 for (unsigned i=0; i<dim_; ++i)
02358 buf[i] = 0U;
02359 processSubrangeLoop(0U, 0UL, buf, f, subrange);
02360 destroyBuffer(buf, ibuf);
02361 }
02362
02363 template<typename Numeric, unsigned Len, unsigned Dim>
02364 template<typename Num2>
02365 inline ArrayND<Numeric,Len,Dim>& ArrayND<Numeric,Len,Dim>::setData(
02366 const Num2* data, const unsigned long dataLength)
02367 {
02368 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
02369 "Initialize npstat::ArrayND before calling method \"setData\"");
02370 if (dataLength != len_) throw npstat::NpstatInvalidArgument(
02371 "In npstat::ArrayND::setData: incompatible input data length");
02372 copyBuffer(data_, data, dataLength);
02373 return *this;
02374 }
02375
02376 template<typename Numeric, unsigned Len, unsigned Dim>
02377 void ArrayND<Numeric,Len,Dim>::buildStrides()
02378 {
02379 assert(dim_);
02380 if (strides_ == 0)
02381 strides_ = makeBuffer(dim_, localStrides_, Dim);
02382 strides_[dim_ - 1] = 1UL;
02383 for (unsigned j=dim_ - 1; j>0; --j)
02384 strides_[j - 1] = strides_[j]*shape_[j];
02385 }
02386
02387 template<typename Numeric, unsigned StackLen, unsigned StackDim>
02388 inline ArrayND<Numeric,StackLen,StackDim>::ArrayND()
02389 : data_(0), strides_(0), shape_(0),
02390 len_(0UL), dim_(0U), shapeIsKnown_(false)
02391 {
02392 localData_[0] = Numeric();
02393 data_ = localData_;
02394 }
02395
02396 template<typename Numeric, unsigned Len, unsigned Dim>
02397 ArrayND<Numeric,Len,Dim>::ArrayND(const ArrayND& r)
02398 : data_(0), strides_(0), shape_(0),
02399 len_(r.len_), dim_(r.dim_), shapeIsKnown_(r.shapeIsKnown_)
02400 {
02401 if (dim_)
02402 {
02403
02404 shape_ = makeBuffer(dim_, localShape_, Dim);
02405 copyBuffer(shape_, r.shape_, dim_);
02406
02407
02408 strides_ = makeBuffer(dim_, localStrides_, Dim);
02409 copyBuffer(strides_, r.strides_, dim_);
02410
02411
02412 data_ = makeBuffer(len_, localData_, Len);
02413 copyBuffer(data_, r.data_, len_);
02414 }
02415 else
02416 {
02417 assert(len_ == 1UL);
02418 localData_[0] = r.localData_[0];
02419 data_ = localData_;
02420 }
02421 }
02422
02423 template<typename Numeric, unsigned Len, unsigned Dim>
02424 template <typename Num2, unsigned Len2, unsigned Dim2>
02425 ArrayND<Numeric,Len,Dim>::ArrayND(const ArrayND<Num2, Len2, Dim2>& r)
02426 : data_(0), strides_(0), shape_(0),
02427 len_(r.len_), dim_(r.dim_), shapeIsKnown_(r.shapeIsKnown_)
02428 {
02429 if (dim_)
02430 {
02431
02432 shape_ = makeBuffer(dim_, localShape_, Dim);
02433 copyBuffer(shape_, r.shape_, dim_);
02434
02435
02436 strides_ = makeBuffer(dim_, localStrides_, Dim);
02437 copyBuffer(strides_, r.strides_, dim_);
02438
02439
02440 data_ = makeBuffer(len_, localData_, Len);
02441 copyBuffer(data_, r.data_, len_);
02442 }
02443 else
02444 {
02445 assert(len_ == 1UL);
02446 localData_[0] = static_cast<Numeric>(r.localData_[0]);
02447 data_ = localData_;
02448 }
02449 }
02450
02451 template<typename Numeric, unsigned Len, unsigned Dim>
02452 template<typename Num2, unsigned Len2, unsigned Dim2, class Functor>
02453 ArrayND<Numeric,Len,Dim>::ArrayND(const ArrayND<Num2, Len2, Dim2>& r,
02454 Functor f)
02455 : data_(0), strides_(0), shape_(0),
02456 len_(r.len_), dim_(r.dim_), shapeIsKnown_(r.shapeIsKnown_)
02457 {
02458 if (dim_)
02459 {
02460
02461 shape_ = makeBuffer(dim_, localShape_, Dim);
02462 copyBuffer(shape_, r.shape_, dim_);
02463
02464
02465 strides_ = makeBuffer(dim_, localStrides_, Dim);
02466 copyBuffer(strides_, r.strides_, dim_);
02467
02468
02469 data_ = makeBuffer(len_, localData_, Len);
02470 for (unsigned long i=0; i<len_; ++i)
02471 data_[i] = static_cast<Numeric>(f(r.data_[i]));
02472 }
02473 else
02474 {
02475 assert(len_ == 1UL);
02476 localData_[0] = static_cast<Numeric>(f(r.localData_[0]));
02477 data_ = localData_;
02478 }
02479 }
02480
02481 template<typename Numeric, unsigned Len, unsigned Dim>
02482 template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
02483 void ArrayND<Numeric,Len,Dim>::copyRangeLoopFunct(
02484 const unsigned level, unsigned long idx0,
02485 unsigned long idx1,
02486 const ArrayND<Num2, Len2, Dim2>& r,
02487 const ArrayRange& range, Functor f)
02488 {
02489 const unsigned imax = shape_[level];
02490 if (level == dim_ - 1)
02491 {
02492 Numeric* to = data_ + idx0;
02493 const Num2* from = r.data_ + (idx1 + range[level].min());
02494 for (unsigned i=0; i<imax; ++i)
02495 *to++ = static_cast<Numeric>(f(*from++));
02496 }
02497 else
02498 {
02499 const unsigned long fromstride = r.strides_[level];
02500 const unsigned long tostride = strides_[level];
02501 idx1 += range[level].min()*fromstride;
02502 for (unsigned i=0; i<imax; ++i)
02503 {
02504 copyRangeLoopFunct(level+1, idx0, idx1, r, range, f);
02505 idx0 += tostride;
02506 idx1 += fromstride;
02507 }
02508 }
02509 }
02510
02511 template<typename Numeric, unsigned Len, unsigned Dim>
02512 template <typename Num2, unsigned Len2, unsigned Dim2>
02513 ArrayND<Numeric,Len,Dim>::ArrayND(
02514 const ArrayND<Num2, Len2, Dim2>& r, const ArrayRange& range)
02515 : data_(0), strides_(0), shape_(0),
02516 len_(r.len_), dim_(r.dim_), shapeIsKnown_(r.shapeIsKnown_)
02517 {
02518 if (!range.isCompatible(r.shape_, r.dim_))
02519 throw npstat::NpstatInvalidArgument(
02520 "In npstat::ArrayND subrange constructor: invalid subrange");
02521 if (dim_)
02522 {
02523 len_ = range.rangeSize();
02524 if (!len_)
02525 throw npstat::NpstatInvalidArgument(
02526 "In npstat::ArrayND subrange constructor: empty subrange");
02527
02528
02529 shape_ = makeBuffer(dim_, localShape_, Dim);
02530 range.rangeLength(shape_, dim_);
02531
02532
02533 buildStrides();
02534
02535
02536 data_ = makeBuffer(len_, localData_, Len);
02537
02538
02539 if (dim_ > CHAR_BIT*sizeof(unsigned long))
02540 throw npstat::NpstatInvalidArgument(
02541 "In npstat::ArrayND subrange constructor: "
02542 "input array rank is too large");
02543 unsigned lolim[CHAR_BIT*sizeof(unsigned long)];
02544 range.lowerLimits(lolim, dim_);
02545 unsigned toBuf[CHAR_BIT*sizeof(unsigned long)];
02546 clearBuffer(toBuf, dim_);
02547 (const_cast<ArrayND<Num2, Len2, Dim2>&>(r)).commonSubrangeLoop(
02548 0U, 0UL, 0UL, lolim, shape_, toBuf, *this,
02549 scast_assign_right<Num2,Numeric>());
02550 }
02551 else
02552 {
02553 assert(len_ == 1UL);
02554 localData_[0] = static_cast<Numeric>(r.localData_[0]);
02555 data_ = localData_;
02556 }
02557 }
02558
02559 template<typename Numeric, unsigned Len, unsigned Dim>
02560 template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
02561 ArrayND<Numeric,Len,Dim>::ArrayND(
02562 const ArrayND<Num2, Len2, Dim2>& r, const ArrayRange& range,
02563 Functor f)
02564 : data_(0), strides_(0), shape_(0),
02565 len_(r.len_), dim_(r.dim_), shapeIsKnown_(r.shapeIsKnown_)
02566 {
02567 if (!range.isCompatible(r.shape_, r.dim_))
02568 throw npstat::NpstatInvalidArgument(
02569 "In npstat::ArrayND transforming subrange constructor: "
02570 "incompatible subrange");
02571 if (dim_)
02572 {
02573 len_ = range.rangeSize();
02574 if (!len_)
02575 throw npstat::NpstatInvalidArgument(
02576 "In npstat::ArrayND transforming subrange constructor: "
02577 "empty subrange");
02578
02579
02580 shape_ = makeBuffer(dim_, localShape_, Dim);
02581 for (unsigned i=0; i<dim_; ++i)
02582 shape_[i] = range[i].length();
02583
02584
02585 buildStrides();
02586
02587
02588 data_ = makeBuffer(len_, localData_, Len);
02589
02590
02591 copyRangeLoopFunct(0U, 0UL, 0UL, r, range, f);
02592 }
02593 else
02594 {
02595 assert(len_ == 1UL);
02596 localData_[0] = static_cast<Numeric>(f(r.localData_[0]));
02597 data_ = localData_;
02598 }
02599 }
02600
02601 template<typename Numeric, unsigned Len, unsigned Dim>
02602 ArrayND<Numeric,Len,Dim>::ArrayND(const ArrayShape& sh)
02603 : data_(0), strides_(0), shape_(0), len_(1UL), shapeIsKnown_(true)
02604 {
02605 const unsigned sz = sh.size();
02606 buildFromShapePtr(sz ? &sh[0] : 0, sz);
02607 }
02608
02609 template<typename Numeric, unsigned Len, unsigned Dim>
02610 ArrayND<Numeric,Len,Dim>::ArrayND(const unsigned* sizes,
02611 const unsigned dim)
02612 : data_(0), strides_(0), shape_(0), len_(1UL), shapeIsKnown_(true)
02613 {
02614 buildFromShapePtr(sizes, dim);
02615 }
02616
02617 template<typename Numeric, unsigned Len, unsigned Dim>
02618 ArrayND<Numeric,Len,Dim>::ArrayND(const unsigned n0)
02619 : data_(0), strides_(0), shape_(0), len_(1UL), shapeIsKnown_(true)
02620 {
02621 const unsigned dim = 1U;
02622 unsigned sizes[dim];
02623 sizes[0] = n0;
02624 buildFromShapePtr(sizes, dim);
02625 }
02626
02627 template<typename Numeric, unsigned Len, unsigned Dim>
02628 ArrayND<Numeric,Len,Dim>::ArrayND(const unsigned n0,
02629 const unsigned n1)
02630 : data_(0), strides_(0), shape_(0), len_(1UL), shapeIsKnown_(true)
02631 {
02632 const unsigned dim = 2U;
02633 unsigned sizes[dim];
02634 sizes[0] = n0;
02635 sizes[1] = n1;
02636 buildFromShapePtr(sizes, dim);
02637 }
02638
02639 template<typename Numeric, unsigned Len, unsigned Dim>
02640 ArrayND<Numeric,Len,Dim>::ArrayND(const unsigned n0,
02641 const unsigned n1,
02642 const unsigned n2)
02643 : data_(0), strides_(0), shape_(0), len_(1UL), shapeIsKnown_(true)
02644 {
02645 const unsigned dim = 3U;
02646 unsigned sizes[dim];
02647 sizes[0] = n0;
02648 sizes[1] = n1;
02649 sizes[2] = n2;
02650 buildFromShapePtr(sizes, dim);
02651 }
02652
02653 template<typename Numeric, unsigned Len, unsigned Dim>
02654 ArrayND<Numeric,Len,Dim>::ArrayND(const unsigned n0,
02655 const unsigned n1,
02656 const unsigned n2,
02657 const unsigned n3)
02658 : data_(0), strides_(0), shape_(0), len_(1UL), shapeIsKnown_(true)
02659 {
02660 const unsigned dim = 4U;
02661 unsigned sizes[dim];
02662 sizes[0] = n0;
02663 sizes[1] = n1;
02664 sizes[2] = n2;
02665 sizes[3] = n3;
02666 buildFromShapePtr(sizes, dim);
02667 }
02668
02669 template<typename Numeric, unsigned Len, unsigned Dim>
02670 ArrayND<Numeric,Len,Dim>::ArrayND(const unsigned n0,
02671 const unsigned n1,
02672 const unsigned n2,
02673 const unsigned n3,
02674 const unsigned n4)
02675 : data_(0), strides_(0), shape_(0), len_(1UL), shapeIsKnown_(true)
02676 {
02677 const unsigned dim = 5U;
02678 unsigned sizes[dim];
02679 sizes[0] = n0;
02680 sizes[1] = n1;
02681 sizes[2] = n2;
02682 sizes[3] = n3;
02683 sizes[4] = n4;
02684 buildFromShapePtr(sizes, dim);
02685 }
02686
02687 template<typename Numeric, unsigned Len, unsigned Dim>
02688 ArrayND<Numeric,Len,Dim>::ArrayND(const unsigned n0,
02689 const unsigned n1,
02690 const unsigned n2,
02691 const unsigned n3,
02692 const unsigned n4,
02693 const unsigned n5)
02694 : data_(0), strides_(0), shape_(0), len_(1UL), shapeIsKnown_(true)
02695 {
02696 const unsigned dim = 6U;
02697 unsigned sizes[dim];
02698 sizes[0] = n0;
02699 sizes[1] = n1;
02700 sizes[2] = n2;
02701 sizes[3] = n3;
02702 sizes[4] = n4;
02703 sizes[5] = n5;
02704 buildFromShapePtr(sizes, dim);
02705 }
02706
02707 template<typename Numeric, unsigned Len, unsigned Dim>
02708 ArrayND<Numeric,Len,Dim>::ArrayND(const unsigned n0,
02709 const unsigned n1,
02710 const unsigned n2,
02711 const unsigned n3,
02712 const unsigned n4,
02713 const unsigned n5,
02714 const unsigned n6)
02715 : data_(0), strides_(0), shape_(0), len_(1UL), shapeIsKnown_(true)
02716 {
02717 const unsigned dim = 7U;
02718 unsigned sizes[dim];
02719 sizes[0] = n0;
02720 sizes[1] = n1;
02721 sizes[2] = n2;
02722 sizes[3] = n3;
02723 sizes[4] = n4;
02724 sizes[5] = n5;
02725 sizes[6] = n6;
02726 buildFromShapePtr(sizes, dim);
02727 }
02728
02729 template<typename Numeric, unsigned Len, unsigned Dim>
02730 ArrayND<Numeric,Len,Dim>::ArrayND(const unsigned n0,
02731 const unsigned n1,
02732 const unsigned n2,
02733 const unsigned n3,
02734 const unsigned n4,
02735 const unsigned n5,
02736 const unsigned n6,
02737 const unsigned n7)
02738 : data_(0), strides_(0), shape_(0), len_(1UL), shapeIsKnown_(true)
02739 {
02740 const unsigned dim = 8U;
02741 unsigned sizes[dim];
02742 sizes[0] = n0;
02743 sizes[1] = n1;
02744 sizes[2] = n2;
02745 sizes[3] = n3;
02746 sizes[4] = n4;
02747 sizes[5] = n5;
02748 sizes[6] = n6;
02749 sizes[7] = n7;
02750 buildFromShapePtr(sizes, dim);
02751 }
02752
02753 template<typename Numeric, unsigned Len, unsigned Dim>
02754 ArrayND<Numeric,Len,Dim>::ArrayND(const unsigned n0,
02755 const unsigned n1,
02756 const unsigned n2,
02757 const unsigned n3,
02758 const unsigned n4,
02759 const unsigned n5,
02760 const unsigned n6,
02761 const unsigned n7,
02762 const unsigned n8)
02763 : data_(0), strides_(0), shape_(0), len_(1UL), shapeIsKnown_(true)
02764 {
02765 const unsigned dim = 9U;
02766 unsigned sizes[dim];
02767 sizes[0] = n0;
02768 sizes[1] = n1;
02769 sizes[2] = n2;
02770 sizes[3] = n3;
02771 sizes[4] = n4;
02772 sizes[5] = n5;
02773 sizes[6] = n6;
02774 sizes[7] = n7;
02775 sizes[8] = n8;
02776 buildFromShapePtr(sizes, dim);
02777 }
02778
02779 template<typename Numeric, unsigned Len, unsigned Dim>
02780 ArrayND<Numeric,Len,Dim>::ArrayND(const unsigned n0,
02781 const unsigned n1,
02782 const unsigned n2,
02783 const unsigned n3,
02784 const unsigned n4,
02785 const unsigned n5,
02786 const unsigned n6,
02787 const unsigned n7,
02788 const unsigned n8,
02789 const unsigned n9)
02790 : data_(0), strides_(0), shape_(0), len_(1UL), shapeIsKnown_(true)
02791 {
02792 const unsigned dim = 10U;
02793 unsigned sizes[dim];
02794 sizes[0] = n0;
02795 sizes[1] = n1;
02796 sizes[2] = n2;
02797 sizes[3] = n3;
02798 sizes[4] = n4;
02799 sizes[5] = n5;
02800 sizes[6] = n6;
02801 sizes[7] = n7;
02802 sizes[8] = n8;
02803 sizes[9] = n9;
02804 buildFromShapePtr(sizes, dim);
02805 }
02806
02807 template<typename Numeric, unsigned Len, unsigned Dim>
02808 template<typename Num1, unsigned Len1, unsigned Dim1,
02809 typename Num2, unsigned Len2, unsigned Dim2>
02810 void ArrayND<Numeric,Len,Dim>::outerProductLoop(
02811 const unsigned level, unsigned long idx0,
02812 unsigned long idx1, unsigned long idx2,
02813 const ArrayND<Num1, Len1, Dim1>& a1,
02814 const ArrayND<Num2, Len2, Dim2>& a2)
02815 {
02816 const unsigned imax = shape_[level];
02817 if (level == dim_ - 1)
02818 {
02819 for (unsigned i=0; i<imax; ++i)
02820 data_[idx0 + i] = a1.data_[idx1]*a2.data_[idx2 + i];
02821 }
02822 else
02823 {
02824 for (unsigned i=0; i<imax; ++i)
02825 {
02826 outerProductLoop(level+1, idx0, idx1, idx2, a1, a2);
02827 idx0 += strides_[level];
02828 if (level < a1.dim_)
02829 idx1 += a1.strides_[level];
02830 else
02831 idx2 += a2.strides_[level - a1.dim_];
02832 }
02833 }
02834 }
02835
02836 template<typename Numeric, unsigned Len, unsigned Dim>
02837 template<typename Num1, unsigned Len1, unsigned Dim1,
02838 typename Num2, unsigned Len2, unsigned Dim2>
02839 ArrayND<Numeric,Len,Dim>::ArrayND(const ArrayND<Num1, Len1, Dim1>& a1,
02840 const ArrayND<Num2, Len2, Dim2>& a2)
02841 : data_(0), strides_(0), shape_(0),
02842 len_(1UL), dim_(a1.dim_ + a2.dim_), shapeIsKnown_(true)
02843 {
02844 if (!(a1.shapeIsKnown_ && a2.shapeIsKnown_))
02845 throw npstat::NpstatInvalidArgument(
02846 "In npstat::ArrayND outer product constructor: "
02847 "uninitialized argument array");
02848 if (dim_)
02849 {
02850 shape_ = makeBuffer(dim_, localShape_, Dim);
02851 copyBuffer(shape_, a1.shape_, a1.dim_);
02852 copyBuffer(shape_+a1.dim_, a2.shape_, a2.dim_);
02853
02854 for (unsigned i=0; i<dim_; ++i)
02855 {
02856 assert(shape_[i]);
02857 len_ *= shape_[i];
02858 }
02859
02860
02861 buildStrides();
02862
02863
02864 data_ = makeBuffer(len_, localData_, Len);
02865
02866
02867 if (a1.dim_ == 0)
02868 {
02869 for (unsigned long i=0; i<len_; ++i)
02870 data_[i] = a1.localData_[0] * a2.data_[i];
02871 }
02872 else if (a2.dim_ == 0)
02873 {
02874 for (unsigned long i=0; i<len_; ++i)
02875 data_[i] = a1.data_[i] * a2.localData_[0];
02876 }
02877 else
02878 outerProductLoop(0U, 0UL, 0UL, 0UL, a1, a2);
02879 }
02880 else
02881 {
02882 localData_[0] = a1.localData_[0] * a2.localData_[0];
02883 data_ = localData_;
02884 }
02885 }
02886
02887 template<typename Numeric, unsigned Len, unsigned Dim>
02888 inline ArrayND<Numeric,Len,Dim>::~ArrayND()
02889 {
02890 destroyBuffer(data_, localData_);
02891 destroyBuffer(strides_, localStrides_);
02892 destroyBuffer(shape_, localShape_);
02893 }
02894
02895 template<typename Numeric, unsigned Len, unsigned Dim>
02896 ArrayND<Numeric,Len,Dim>&
02897 ArrayND<Numeric,Len,Dim>::operator=(const ArrayND& r)
02898 {
02899 if (this == &r)
02900 return *this;
02901 if (shapeIsKnown_)
02902 {
02903 if (!r.shapeIsKnown_) throw npstat::NpstatInvalidArgument(
02904 "In npstat::ArrayND assignment operator: "
02905 "uninitialized argument array");
02906 if (!isShapeCompatible(r)) throw npstat::NpstatInvalidArgument(
02907 "In npstat::ArrayND assignment operator: "
02908 "incompatible argument array shape");
02909 if (dim_)
02910 copyBuffer(data_, r.data_, len_);
02911 else
02912 localData_[0] = r.localData_[0];
02913 }
02914 else
02915 {
02916
02917
02918 if (r.shapeIsKnown_)
02919 new (this) ArrayND(r);
02920 }
02921 return *this;
02922 }
02923
02924 template<typename Numeric, unsigned Len, unsigned Dim>
02925 template <typename Num2, unsigned Len2, unsigned Dim2>
02926 ArrayND<Numeric,Len,Dim>&
02927 ArrayND<Numeric,Len,Dim>::operator=(const ArrayND<Num2,Len2,Dim2>& r)
02928 {
02929 if ((void*)this == (void*)(&r))
02930 return *this;
02931 if (shapeIsKnown_)
02932 {
02933 if (!r.shapeIsKnown_) throw npstat::NpstatInvalidArgument(
02934 "In npstat::ArrayND assignment operator: "
02935 "uninitialized argument array");
02936 if (!isShapeCompatible(r)) throw npstat::NpstatInvalidArgument(
02937 "In npstat::ArrayND assignment operator: "
02938 "incompatible argument array shape");
02939 if (dim_)
02940 copyBuffer(data_, r.data_, len_);
02941 else
02942 localData_[0] = static_cast<Numeric>(r.localData_[0]);
02943 }
02944 else
02945 {
02946
02947
02948 if (r.shapeIsKnown_)
02949 new (this) ArrayND(r);
02950 }
02951 return *this;
02952 }
02953
02954 template<typename Numeric, unsigned Len, unsigned Dim>
02955 template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
02956 ArrayND<Numeric,Len,Dim>&
02957 ArrayND<Numeric,Len,Dim>::assign(const ArrayND<Num2,Len2,Dim2>& r,
02958 Functor f)
02959 {
02960 if (shapeIsKnown_)
02961 {
02962 if (!r.shapeIsKnown_) throw npstat::NpstatInvalidArgument(
02963 "In npstat::ArrayND::assign: uninitialized argument array");
02964 if (!isShapeCompatible(r)) throw npstat::NpstatInvalidArgument(
02965 "In npstat::ArrayND::assign: incompatible argument array shape");
02966 if (dim_)
02967 for (unsigned long i=0; i<len_; ++i)
02968 data_[i] = static_cast<Numeric>(f(r.data_[i]));
02969 else
02970 localData_[0] = static_cast<Numeric>(f(r.localData_[0]));
02971 }
02972 else
02973 {
02974
02975
02976 if (r.shapeIsKnown_)
02977 new (this) ArrayND(r, f);
02978 }
02979 return *this;
02980 }
02981
02982 template<typename Numeric, unsigned Len, unsigned Dim>
02983 inline ArrayShape ArrayND<Numeric,Len,Dim>::shape() const
02984 {
02985 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
02986 "Initialize npstat::ArrayND before calling method \"shape\"");
02987 return ArrayShape(shape_, shape_+dim_);
02988 }
02989
02990 template<typename Numeric, unsigned Len, unsigned Dim>
02991 inline ArrayRange ArrayND<Numeric,Len,Dim>::fullRange() const
02992 {
02993 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
02994 "Initialize npstat::ArrayND before calling method \"fullRange\"");
02995 ArrayRange range;
02996 if (dim_)
02997 {
02998 range.reserve(dim_);
02999 for (unsigned i=0; i<dim_; ++i)
03000 range.push_back(Interval<unsigned>(0U, shape_[i]));
03001 }
03002 return range;
03003 }
03004
03005 template<typename Numeric, unsigned Len, unsigned Dim>
03006 bool ArrayND<Numeric,Len,Dim>::isDensity() const
03007 {
03008 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
03009 "Initialize npstat::ArrayND before calling method \"isDensity\"");
03010 const Numeric zero = Numeric();
03011 bool hasPositive = false;
03012 if (dim_)
03013 for (unsigned long i=0; i<len_; ++i)
03014 {
03015
03016
03017
03018
03019
03020 if (data_[i] == zero)
03021 continue;
03022 if (ComplexComparesFalse<Numeric>::less(zero, data_[i]))
03023 hasPositive = true;
03024 else
03025 return false;
03026 }
03027 else
03028 hasPositive = ComplexComparesFalse<Numeric>::less(
03029 zero, localData_[0]);
03030 return hasPositive;
03031 }
03032
03033 template<typename Numeric, unsigned Len, unsigned Dim>
03034 bool ArrayND<Numeric,Len,Dim>::isZero() const
03035 {
03036 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
03037 "Initialize npstat::ArrayND before calling method \"isZero\"");
03038 const Numeric zero = Numeric();
03039 if (dim_)
03040 {
03041 for (unsigned long i=0; i<len_; ++i)
03042 if (data_[i] != zero)
03043 return false;
03044 }
03045 else
03046 if (localData_[0] != zero)
03047 return false;
03048 return true;
03049 }
03050
03051 template<typename Numeric, unsigned Len, unsigned Dim>
03052 void ArrayND<Numeric,Len,Dim>::convertLinearIndex(
03053 unsigned long l, unsigned* idx, const unsigned idxLen) const
03054 {
03055 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
03056 "Initialize npstat::ArrayND before calling "
03057 "method \"convertLinearIndex\"");
03058 if (!dim_) throw npstat::NpstatInvalidArgument(
03059 "npstat::ArrayND::convertLinearIndex method "
03060 "can not be used with array of 0 rank");
03061 if (idxLen != dim_) throw npstat::NpstatInvalidArgument(
03062 "In npstat::ArrayND::convertLinearIndex: incompatible index length");
03063 if (l >= len_) throw npstat::NpstatOutOfRange(
03064 "In npstat::ArrayND::convertLinearIndex: linear index out of range");
03065 assert(idx);
03066
03067 for (unsigned i=0; i<dim_; ++i)
03068 {
03069 idx[i] = l / strides_[i];
03070 l -= (idx[i] * strides_[i]);
03071 }
03072 }
03073
03074 template<typename Numeric, unsigned Len, unsigned Dim>
03075 unsigned long ArrayND<Numeric,Len,Dim>::linearIndex(
03076 const unsigned* index, unsigned idxLen) const
03077 {
03078 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
03079 "Initialize npstat::ArrayND before calling method \"linearIndex\"");
03080 if (!dim_) throw npstat::NpstatInvalidArgument(
03081 "npstat::ArrayND::linearIndex method "
03082 "can not be used with array of 0 rank");
03083 if (idxLen != dim_) throw npstat::NpstatInvalidArgument(
03084 "In npstat::ArrayND::linearIndex: incompatible index length");
03085 assert(index);
03086
03087 unsigned long idx = 0UL;
03088 for (unsigned i=0; i<dim_; ++i)
03089 {
03090 if (index[i] >= shape_[i])
03091 throw npstat::NpstatOutOfRange(
03092 "In npstat::ArrayND::linearIndex: index out of range");
03093 idx += index[i]*strides_[i];
03094 }
03095 return idx;
03096 }
03097
03098 template<typename Numeric, unsigned Len, unsigned Dim>
03099 inline Numeric& ArrayND<Numeric,Len,Dim>::value(
03100 const unsigned *index, const unsigned dim)
03101 {
03102 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
03103 "Initialize npstat::ArrayND before calling method \"value\"");
03104 if (dim != dim_) throw npstat::NpstatInvalidArgument(
03105 "In npstat::ArrayND::value: incompatible index length");
03106 if (dim)
03107 {
03108 assert(index);
03109 unsigned long idx = 0UL;
03110 for (unsigned i=0; i<dim_; ++i)
03111 idx += index[i]*strides_[i];
03112 return data_[idx];
03113 }
03114 else
03115 return localData_[0];
03116 }
03117
03118 template<typename Numeric, unsigned Len, unsigned Dim>
03119 inline const Numeric& ArrayND<Numeric,Len,Dim>::value(
03120 const unsigned *index, const unsigned dim) const
03121 {
03122 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
03123 "Initialize npstat::ArrayND before calling method \"value\"");
03124 if (dim != dim_) throw npstat::NpstatInvalidArgument(
03125 "In npstat::ArrayND::value: incompatible index length");
03126 if (dim)
03127 {
03128 assert(index);
03129 unsigned long idx = 0UL;
03130 for (unsigned i=0; i<dim_; ++i)
03131 idx += index[i]*strides_[i];
03132 return data_[idx];
03133 }
03134 else
03135 return localData_[0];
03136 }
03137
03138 template<typename Numeric, unsigned Len, unsigned Dim>
03139 inline Numeric& ArrayND<Numeric,Len,Dim>::linearValue(
03140 const unsigned long index)
03141 {
03142 return data_[index];
03143 }
03144
03145 template<typename Numeric, unsigned Len, unsigned Dim>
03146 inline const Numeric& ArrayND<Numeric,Len,Dim>::linearValue(
03147 const unsigned long index) const
03148 {
03149 return data_[index];
03150 }
03151
03152 template<typename Numeric, unsigned Len, unsigned Dim>
03153 inline Numeric& ArrayND<Numeric,Len,Dim>::linearValueAt(
03154 const unsigned long index)
03155 {
03156 if (index >= len_)
03157 throw npstat::NpstatOutOfRange(
03158 "In npstat::ArrayND::linearValueAt: linear index out of range");
03159 return data_[index];
03160 }
03161
03162 template<typename Numeric, unsigned Len, unsigned Dim>
03163 inline const Numeric& ArrayND<Numeric,Len,Dim>::linearValueAt(
03164 const unsigned long index) const
03165 {
03166 if (index >= len_)
03167 throw npstat::NpstatOutOfRange(
03168 "In npstat::ArrayND::linearValueAt: linear index out of range");
03169 return data_[index];
03170 }
03171
03172 template<typename Numeric, unsigned Len, unsigned Dim>
03173 inline unsigned ArrayND<Numeric,Len,Dim>::coordToIndex(
03174 const double x, const unsigned idim) const
03175 {
03176 if (x <= 0.0)
03177 return 0;
03178 else if (x >= static_cast<double>(shape_[idim] - 1))
03179 return shape_[idim] - 1;
03180 else
03181 return static_cast<unsigned>(std::floor(x + 0.5));
03182 }
03183
03184 template<typename Numeric, unsigned Len, unsigned Dim>
03185 inline const Numeric& ArrayND<Numeric,Len,Dim>::closest(
03186 const double *x, const unsigned dim) const
03187 {
03188 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
03189 "Initialize npstat::ArrayND before calling method \"closest\"");
03190 if (dim != dim_) throw npstat::NpstatInvalidArgument(
03191 "In npstat::ArrayND::closest: incompatible data length");
03192 if (dim)
03193 {
03194 assert(x);
03195 unsigned long idx = 0UL;
03196 for (unsigned i=0; i<dim_; ++i)
03197 idx += coordToIndex(x[i], i)*strides_[i];
03198 return data_[idx];
03199 }
03200 else
03201 return localData_[0];
03202 }
03203
03204 template<typename Numeric, unsigned Len, unsigned Dim>
03205 inline Numeric& ArrayND<Numeric,Len,Dim>::closest(
03206 const double *x, const unsigned dim)
03207 {
03208 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
03209 "Initialize npstat::ArrayND before calling method \"closest\"");
03210 if (dim != dim_) throw npstat::NpstatInvalidArgument(
03211 "In npstat::ArrayND::closest: incompatible data length");
03212 if (dim)
03213 {
03214 assert(x);
03215 unsigned long idx = 0UL;
03216 for (unsigned i=0; i<dim_; ++i)
03217 idx += coordToIndex(x[i], i)*strides_[i];
03218 return data_[idx];
03219 }
03220 else
03221 return localData_[0];
03222 }
03223
03224 template<typename Numeric, unsigned Len, unsigned Dim>
03225 inline const Numeric& ArrayND<Numeric,Len,Dim>::valueAt(
03226 const unsigned *index, const unsigned dim) const
03227 {
03228 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
03229 "Initialize npstat::ArrayND before calling method \"valueAt\"");
03230 if (dim != dim_) throw npstat::NpstatInvalidArgument(
03231 "In npstat::ArrayND::valueAt: incompatible index length");
03232 if (dim)
03233 {
03234 assert(index);
03235 unsigned long idx = 0UL;
03236 for (unsigned i=0; i<dim_; ++i)
03237 {
03238 if (index[i] >= shape_[i]) throw npstat::NpstatOutOfRange(
03239 "In npstat::ArrayND::valueAt: index out of range");
03240 idx += index[i]*strides_[i];
03241 }
03242 return data_[idx];
03243 }
03244 else
03245 return localData_[0];
03246 }
03247
03248 template<typename Numeric, unsigned Len, unsigned Dim>
03249 inline Numeric& ArrayND<Numeric,Len,Dim>::valueAt(
03250 const unsigned *index, const unsigned dim)
03251 {
03252 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
03253 "Initialize npstat::ArrayND before calling method \"valueAt\"");
03254 if (dim != dim_) throw npstat::NpstatInvalidArgument(
03255 "In npstat::ArrayND::valueAt: incompatible index length");
03256 if (dim)
03257 {
03258 assert(index);
03259 unsigned long idx = 0UL;
03260 for (unsigned i=0; i<dim_; ++i)
03261 {
03262 if (index[i] >= shape_[i]) throw npstat::NpstatOutOfRange(
03263 "In npstat::ArrayND::valueAt: index out of range");
03264 idx += index[i]*strides_[i];
03265 }
03266 return data_[idx];
03267 }
03268 else
03269 return localData_[0];
03270 }
03271
03272 template<typename Numeric, unsigned Len, unsigned Dim>
03273 inline Numeric& ArrayND<Numeric,Len,Dim>::operator()()
03274 {
03275 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
03276 "Initialize npstat::ArrayND before calling method \"operator()\"");
03277 if (dim_) throw npstat::NpstatInvalidArgument(
03278 "In npstat::ArrayND::operator(): wrong # of args (not rank 0 array)");
03279 return localData_[0];
03280 }
03281
03282 template<typename Numeric, unsigned Len, unsigned Dim>
03283 inline const Numeric& ArrayND<Numeric,Len,Dim>::operator()() const
03284 {
03285 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
03286 "Initialize npstat::ArrayND before calling method \"operator()\"");
03287 if (dim_) throw npstat::NpstatInvalidArgument(
03288 "In npstat::ArrayND::operator(): wrong # of args (not rank 0 array)");
03289 return localData_[0];
03290 }
03291
03292 template<typename Numeric, unsigned Len, unsigned Dim>
03293 inline Numeric& ArrayND<Numeric,Len,Dim>::operator()(
03294 const unsigned i)
03295 {
03296 if (1U != dim_) throw npstat::NpstatInvalidArgument(
03297 "In npstat::ArrayND::operator(): wrong # of args (not rank 1 array)");
03298 return data_[i];
03299 }
03300
03301 template<typename Numeric, unsigned Len, unsigned Dim>
03302 inline const Numeric& ArrayND<Numeric,Len,Dim>::operator()(
03303 const unsigned i) const
03304 {
03305 if (1U != dim_) throw npstat::NpstatInvalidArgument(
03306 "In npstat::ArrayND::operator(): wrong # of args (not rank 1 array)");
03307 return data_[i];
03308 }
03309
03310 template<typename Numeric, unsigned Len, unsigned Dim>
03311 const Numeric& ArrayND<Numeric,Len,Dim>::at() const
03312 {
03313 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
03314 "Initialize npstat::ArrayND before calling method \"at\"");
03315 if (dim_) throw npstat::NpstatInvalidArgument(
03316 "In npstat::ArrayND::at: wrong # of args (not rank 0 array)");
03317 return localData_[0];
03318 }
03319
03320 template<typename Numeric, unsigned Len, unsigned Dim>
03321 Numeric& ArrayND<Numeric,Len,Dim>::at()
03322 {
03323 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
03324 "Initialize npstat::ArrayND before calling method \"at\"");
03325 if (dim_) throw npstat::NpstatInvalidArgument(
03326 "In npstat::ArrayND::at: wrong # of args (not rank 0 array)");
03327 return localData_[0];
03328 }
03329
03330 template<typename Numeric, unsigned Len, unsigned Dim>
03331 const Numeric& ArrayND<Numeric,Len,Dim>::at(
03332 const unsigned i0) const
03333 {
03334 if (1U != dim_) throw npstat::NpstatInvalidArgument(
03335 "In npstat::ArrayND::at: wrong # of args (not rank 1 array)");
03336 if (i0 >= shape_[0]) throw npstat::NpstatOutOfRange(
03337 "In npstat::ArrayND::at: index 0 out of range (rank 1)");
03338 return data_[i0];
03339 }
03340
03341 template<typename Numeric, unsigned Len, unsigned Dim>
03342 Numeric& ArrayND<Numeric,Len,Dim>::at(
03343 const unsigned i0)
03344 {
03345 if (1U != dim_) throw npstat::NpstatInvalidArgument(
03346 "In npstat::ArrayND::at: wrong # of args (not rank 1 array)");
03347 if (i0 >= shape_[0]) throw npstat::NpstatOutOfRange(
03348 "In npstat::ArrayND::at: index 0 out of range (rank 1)");
03349 return data_[i0];
03350 }
03351
03352 template<typename Numeric, unsigned Len, unsigned Dim>
03353 inline Numeric& ArrayND<Numeric,Len,Dim>::operator()(
03354 const unsigned i0,
03355 const unsigned i1)
03356 {
03357 if (2U != dim_) throw npstat::NpstatInvalidArgument(
03358 "In npstat::ArrayND::operator(): wrong # of args (not rank 2 array)");
03359 return data_[i0*strides_[0] + i1];
03360 }
03361
03362 template<typename Numeric, unsigned Len, unsigned Dim>
03363 inline const Numeric& ArrayND<Numeric,Len,Dim>::operator()(
03364 const unsigned i0,
03365 const unsigned i1) const
03366 {
03367 if (2U != dim_) throw npstat::NpstatInvalidArgument(
03368 "In npstat::ArrayND::operator(): wrong # of args (not rank 2 array)");
03369 return data_[i0*strides_[0] + i1];
03370 }
03371
03372 template<typename Numeric, unsigned Len, unsigned Dim>
03373 const Numeric& ArrayND<Numeric,Len,Dim>::at(
03374 const unsigned i0,
03375 const unsigned i1) const
03376 {
03377 if (2U != dim_) throw npstat::NpstatInvalidArgument(
03378 "In npstat::ArrayND::at: wrong # of args (not rank 2 array)");
03379 if (i0 >= shape_[0]) throw npstat::NpstatOutOfRange(
03380 "In npstat::ArrayND::at: index 0 out of range (rank 2)");
03381 if (i1 >= shape_[1]) throw npstat::NpstatOutOfRange(
03382 "In npstat::ArrayND::at: index 1 out of range (rank 2)");
03383 return data_[i0*strides_[0] + i1];
03384 }
03385
03386 template<typename Numeric, unsigned Len, unsigned Dim>
03387 Numeric& ArrayND<Numeric,Len,Dim>::at(
03388 const unsigned i0,
03389 const unsigned i1)
03390 {
03391 if (2U != dim_) throw npstat::NpstatInvalidArgument(
03392 "In npstat::ArrayND::at: wrong # of args (not rank 2 array)");
03393 if (i0 >= shape_[0]) throw npstat::NpstatOutOfRange(
03394 "In npstat::ArrayND::at: index 0 out of range (rank 2)");
03395 if (i1 >= shape_[1]) throw npstat::NpstatOutOfRange(
03396 "In npstat::ArrayND::at: index 1 out of range (rank 2)");
03397 return data_[i0*strides_[0] + i1];
03398 }
03399
03400 template<typename Numeric, unsigned Len, unsigned Dim>
03401 inline const Numeric& ArrayND<Numeric,Len,Dim>::operator()(
03402 const unsigned i0,
03403 const unsigned i1,
03404 const unsigned i2) const
03405 {
03406 if (3U != dim_) throw npstat::NpstatInvalidArgument(
03407 "In npstat::ArrayND::operator(): wrong # of args (not rank 3 array)");
03408 return data_[i0*strides_[0] + i1*strides_[1] + i2];
03409 }
03410
03411 template<typename Numeric, unsigned Len, unsigned Dim>
03412 inline const Numeric& ArrayND<Numeric,Len,Dim>::operator()(
03413 const unsigned i0,
03414 const unsigned i1,
03415 const unsigned i2,
03416 const unsigned i3) const
03417 {
03418 if (4U != dim_) throw npstat::NpstatInvalidArgument(
03419 "In npstat::ArrayND::operator(): wrong # of args (not rank 4 array)");
03420 return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] + i3];
03421 }
03422
03423 template<typename Numeric, unsigned Len, unsigned Dim>
03424 inline const Numeric& ArrayND<Numeric,Len,Dim>::operator()(
03425 const unsigned i0,
03426 const unsigned i1,
03427 const unsigned i2,
03428 const unsigned i3,
03429 const unsigned i4) const
03430 {
03431 if (5U != dim_) throw npstat::NpstatInvalidArgument(
03432 "In npstat::ArrayND::operator(): wrong # of args (not rank 5 array)");
03433 return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] +
03434 i3*strides_[3] + i4];
03435 }
03436
03437 template<typename Numeric, unsigned Len, unsigned Dim>
03438 inline const Numeric& ArrayND<Numeric,Len,Dim>::operator()(
03439 const unsigned i0,
03440 const unsigned i1,
03441 const unsigned i2,
03442 const unsigned i3,
03443 const unsigned i4,
03444 const unsigned i5) const
03445 {
03446 if (6U != dim_) throw npstat::NpstatInvalidArgument(
03447 "In npstat::ArrayND::operator(): wrong # of args (not rank 6 array)");
03448 return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] +
03449 i3*strides_[3] + i4*strides_[4] + i5];
03450 }
03451
03452 template<typename Numeric, unsigned Len, unsigned Dim>
03453 inline const Numeric& ArrayND<Numeric,Len,Dim>::operator()(
03454 const unsigned i0,
03455 const unsigned i1,
03456 const unsigned i2,
03457 const unsigned i3,
03458 const unsigned i4,
03459 const unsigned i5,
03460 const unsigned i6) const
03461 {
03462 if (7U != dim_) throw npstat::NpstatInvalidArgument(
03463 "In npstat::ArrayND::operator(): wrong # of args (not rank 7 array)");
03464 return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] +
03465 i3*strides_[3] + i4*strides_[4] + i5*strides_[5] + i6];
03466 }
03467
03468 template<typename Numeric, unsigned Len, unsigned Dim>
03469 inline const Numeric& ArrayND<Numeric,Len,Dim>::operator()(
03470 const unsigned i0,
03471 const unsigned i1,
03472 const unsigned i2,
03473 const unsigned i3,
03474 const unsigned i4,
03475 const unsigned i5,
03476 const unsigned i6,
03477 const unsigned i7) const
03478 {
03479 if (8U != dim_) throw npstat::NpstatInvalidArgument(
03480 "In npstat::ArrayND::operator(): wrong # of args (not rank 8 array)");
03481 return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] +
03482 i3*strides_[3] + i4*strides_[4] + i5*strides_[5] +
03483 i6*strides_[6] + i7];
03484 }
03485
03486 template<typename Numeric, unsigned Len, unsigned Dim>
03487 inline const Numeric& ArrayND<Numeric,Len,Dim>::operator()(
03488 const unsigned i0,
03489 const unsigned i1,
03490 const unsigned i2,
03491 const unsigned i3,
03492 const unsigned i4,
03493 const unsigned i5,
03494 const unsigned i6,
03495 const unsigned i7,
03496 const unsigned i8) const
03497 {
03498 if (9U != dim_) throw npstat::NpstatInvalidArgument(
03499 "In npstat::ArrayND::operator(): wrong # of args (not rank 9 array)");
03500 return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] +
03501 i3*strides_[3] + i4*strides_[4] + i5*strides_[5] +
03502 i6*strides_[6] + i7*strides_[7] + i8];
03503 }
03504
03505 template<typename Numeric, unsigned Len, unsigned Dim>
03506 inline const Numeric& ArrayND<Numeric,Len,Dim>::operator()(
03507 const unsigned i0,
03508 const unsigned i1,
03509 const unsigned i2,
03510 const unsigned i3,
03511 const unsigned i4,
03512 const unsigned i5,
03513 const unsigned i6,
03514 const unsigned i7,
03515 const unsigned i8,
03516 const unsigned i9) const
03517 {
03518 if (10U != dim_) throw npstat::NpstatInvalidArgument(
03519 "In npstat::ArrayND::operator(): wrong # of args (not rank 10 array)");
03520 return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] +
03521 i3*strides_[3] + i4*strides_[4] + i5*strides_[5] +
03522 i6*strides_[6] + i7*strides_[7] + i8*strides_[8] + i9];
03523 }
03524
03525 template<typename Numeric, unsigned Len, unsigned Dim>
03526 inline Numeric& ArrayND<Numeric,Len,Dim>::operator()(
03527 const unsigned i0,
03528 const unsigned i1,
03529 const unsigned i2)
03530 {
03531 if (3U != dim_) throw npstat::NpstatInvalidArgument(
03532 "In npstat::ArrayND::operator(): wrong # of args (not rank 3 array)");
03533 return data_[i0*strides_[0] + i1*strides_[1] + i2];
03534 }
03535
03536 template<typename Numeric, unsigned Len, unsigned Dim>
03537 inline Numeric& ArrayND<Numeric,Len,Dim>::operator()(
03538 const unsigned i0,
03539 const unsigned i1,
03540 const unsigned i2,
03541 const unsigned i3)
03542 {
03543 if (4U != dim_) throw npstat::NpstatInvalidArgument(
03544 "In npstat::ArrayND::operator(): wrong # of args (not rank 4 array)");
03545 return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] + i3];
03546 }
03547
03548 template<typename Numeric, unsigned Len, unsigned Dim>
03549 inline Numeric& ArrayND<Numeric,Len,Dim>::operator()(
03550 const unsigned i0,
03551 const unsigned i1,
03552 const unsigned i2,
03553 const unsigned i3,
03554 const unsigned i4)
03555 {
03556 if (5U != dim_) throw npstat::NpstatInvalidArgument(
03557 "In npstat::ArrayND::operator(): wrong # of args (not rank 5 array)");
03558 return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] +
03559 i3*strides_[3] + i4];
03560 }
03561
03562 template<typename Numeric, unsigned Len, unsigned Dim>
03563 inline Numeric& ArrayND<Numeric,Len,Dim>::operator()(
03564 const unsigned i0,
03565 const unsigned i1,
03566 const unsigned i2,
03567 const unsigned i3,
03568 const unsigned i4,
03569 const unsigned i5)
03570 {
03571 if (6U != dim_) throw npstat::NpstatInvalidArgument(
03572 "In npstat::ArrayND::operator(): wrong # of args (not rank 6 array)");
03573 return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] +
03574 i3*strides_[3] + i4*strides_[4] + i5];
03575 }
03576
03577 template<typename Numeric, unsigned Len, unsigned Dim>
03578 inline Numeric& ArrayND<Numeric,Len,Dim>::operator()(
03579 const unsigned i0,
03580 const unsigned i1,
03581 const unsigned i2,
03582 const unsigned i3,
03583 const unsigned i4,
03584 const unsigned i5,
03585 const unsigned i6)
03586 {
03587 if (7U != dim_) throw npstat::NpstatInvalidArgument(
03588 "In npstat::ArrayND::operator(): wrong # of args (not rank 7 array)");
03589 return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] +
03590 i3*strides_[3] + i4*strides_[4] + i5*strides_[5] + i6];
03591 }
03592
03593 template<typename Numeric, unsigned Len, unsigned Dim>
03594 inline Numeric& ArrayND<Numeric,Len,Dim>::operator()(
03595 const unsigned i0,
03596 const unsigned i1,
03597 const unsigned i2,
03598 const unsigned i3,
03599 const unsigned i4,
03600 const unsigned i5,
03601 const unsigned i6,
03602 const unsigned i7)
03603 {
03604 if (8U != dim_) throw npstat::NpstatInvalidArgument(
03605 "In npstat::ArrayND::operator(): wrong # of args (not rank 8 array)");
03606 return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] +
03607 i3*strides_[3] + i4*strides_[4] + i5*strides_[5] +
03608 i6*strides_[6] + i7];
03609 }
03610
03611 template<typename Numeric, unsigned Len, unsigned Dim>
03612 inline Numeric& ArrayND<Numeric,Len,Dim>::operator()(
03613 const unsigned i0,
03614 const unsigned i1,
03615 const unsigned i2,
03616 const unsigned i3,
03617 const unsigned i4,
03618 const unsigned i5,
03619 const unsigned i6,
03620 const unsigned i7,
03621 const unsigned i8)
03622 {
03623 if (9U != dim_) throw npstat::NpstatInvalidArgument(
03624 "In npstat::ArrayND::operator(): wrong # of args (not rank 9 array)");
03625 return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] +
03626 i3*strides_[3] + i4*strides_[4] + i5*strides_[5] +
03627 i6*strides_[6] + i7*strides_[7] + i8];
03628 }
03629
03630 template<typename Numeric, unsigned Len, unsigned Dim>
03631 inline Numeric& ArrayND<Numeric,Len,Dim>::operator()(
03632 const unsigned i0,
03633 const unsigned i1,
03634 const unsigned i2,
03635 const unsigned i3,
03636 const unsigned i4,
03637 const unsigned i5,
03638 const unsigned i6,
03639 const unsigned i7,
03640 const unsigned i8,
03641 const unsigned i9)
03642 {
03643 if (10U != dim_) throw npstat::NpstatInvalidArgument(
03644 "In npstat::ArrayND::operator(): wrong # of args (not rank 10 array)");
03645 return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] +
03646 i3*strides_[3] + i4*strides_[4] + i5*strides_[5] +
03647 i6*strides_[6] + i7*strides_[7] + i8*strides_[8] + i9];
03648 }
03649
03650 template<typename Numeric, unsigned Len, unsigned Dim>
03651 const Numeric& ArrayND<Numeric,Len,Dim>::at(
03652 const unsigned i0,
03653 const unsigned i1,
03654 const unsigned i2) const
03655 {
03656 if (3U != dim_) throw npstat::NpstatInvalidArgument(
03657 "In npstat::ArrayND::at: wrong # of args (not rank 3 array)");
03658 if (i0 >= shape_[0]) throw npstat::NpstatOutOfRange(
03659 "In npstat::ArrayND::at: index 0 out of range (rank 3)");
03660 if (i1 >= shape_[1]) throw npstat::NpstatOutOfRange(
03661 "In npstat::ArrayND::at: index 1 out of range (rank 3)");
03662 if (i2 >= shape_[2]) throw npstat::NpstatOutOfRange(
03663 "In npstat::ArrayND::at: index 2 out of range (rank 3)");
03664 return data_[i0*strides_[0] + i1*strides_[1] + i2];
03665 }
03666
03667 template<typename Numeric, unsigned Len, unsigned Dim>
03668 const Numeric& ArrayND<Numeric,Len,Dim>::at(
03669 const unsigned i0,
03670 const unsigned i1,
03671 const unsigned i2,
03672 const unsigned i3) const
03673 {
03674 if (4U != dim_) throw npstat::NpstatInvalidArgument(
03675 "In npstat::ArrayND::at: wrong # of args (not rank 4 array)");
03676 if (i0 >= shape_[0]) throw npstat::NpstatOutOfRange(
03677 "In npstat::ArrayND::at: index 0 out of range (rank 4)");
03678 if (i1 >= shape_[1]) throw npstat::NpstatOutOfRange(
03679 "In npstat::ArrayND::at: index 1 out of range (rank 4)");
03680 if (i2 >= shape_[2]) throw npstat::NpstatOutOfRange(
03681 "In npstat::ArrayND::at: index 2 out of range (rank 4)");
03682 if (i3 >= shape_[3]) throw npstat::NpstatOutOfRange(
03683 "In npstat::ArrayND::at: index 3 out of range (rank 4)");
03684 return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] + i3];
03685 }
03686
03687 template<typename Numeric, unsigned Len, unsigned Dim>
03688 const Numeric& ArrayND<Numeric,Len,Dim>::at(
03689 const unsigned i0,
03690 const unsigned i1,
03691 const unsigned i2,
03692 const unsigned i3,
03693 const unsigned i4) const
03694 {
03695 if (5U != dim_) throw npstat::NpstatInvalidArgument(
03696 "In npstat::ArrayND::at: wrong # of args (not rank 5 array)");
03697 if (i0 >= shape_[0]) throw npstat::NpstatOutOfRange(
03698 "In npstat::ArrayND::at: index 0 out of range (rank 5)");
03699 if (i1 >= shape_[1]) throw npstat::NpstatOutOfRange(
03700 "In npstat::ArrayND::at: index 1 out of range (rank 5)");
03701 if (i2 >= shape_[2]) throw npstat::NpstatOutOfRange(
03702 "In npstat::ArrayND::at: index 2 out of range (rank 5)");
03703 if (i3 >= shape_[3]) throw npstat::NpstatOutOfRange(
03704 "In npstat::ArrayND::at: index 3 out of range (rank 5)");
03705 if (i4 >= shape_[4]) throw npstat::NpstatOutOfRange(
03706 "In npstat::ArrayND::at: index 4 out of range (rank 5)");
03707 return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] +
03708 i3*strides_[3] + i4];
03709 }
03710
03711 template<typename Numeric, unsigned Len, unsigned Dim>
03712 const Numeric& ArrayND<Numeric,Len,Dim>::at(
03713 const unsigned i0,
03714 const unsigned i1,
03715 const unsigned i2,
03716 const unsigned i3,
03717 const unsigned i4,
03718 const unsigned i5) const
03719 {
03720 if (6U != dim_) throw npstat::NpstatInvalidArgument(
03721 "In npstat::ArrayND::at: wrong # of args (not rank 6 array)");
03722 if (i0 >= shape_[0]) throw npstat::NpstatOutOfRange(
03723 "In npstat::ArrayND::at: index 0 out of range (rank 6)");
03724 if (i1 >= shape_[1]) throw npstat::NpstatOutOfRange(
03725 "In npstat::ArrayND::at: index 1 out of range (rank 6)");
03726 if (i2 >= shape_[2]) throw npstat::NpstatOutOfRange(
03727 "In npstat::ArrayND::at: index 2 out of range (rank 6)");
03728 if (i3 >= shape_[3]) throw npstat::NpstatOutOfRange(
03729 "In npstat::ArrayND::at: index 3 out of range (rank 6)");
03730 if (i4 >= shape_[4]) throw npstat::NpstatOutOfRange(
03731 "In npstat::ArrayND::at: index 4 out of range (rank 6)");
03732 if (i5 >= shape_[5]) throw npstat::NpstatOutOfRange(
03733 "In npstat::ArrayND::at: index 5 out of range (rank 6)");
03734 return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] +
03735 i3*strides_[3] + i4*strides_[4] + i5];
03736 }
03737
03738 template<typename Numeric, unsigned Len, unsigned Dim>
03739 const Numeric& ArrayND<Numeric,Len,Dim>::at(
03740 const unsigned i0,
03741 const unsigned i1,
03742 const unsigned i2,
03743 const unsigned i3,
03744 const unsigned i4,
03745 const unsigned i5,
03746 const unsigned i6) const
03747 {
03748 if (7U != dim_) throw npstat::NpstatInvalidArgument(
03749 "In npstat::ArrayND::at: wrong # of args (not rank 7 array)");
03750 if (i0 >= shape_[0]) throw npstat::NpstatOutOfRange(
03751 "In npstat::ArrayND::at: index 0 out of range (rank 7)");
03752 if (i1 >= shape_[1]) throw npstat::NpstatOutOfRange(
03753 "In npstat::ArrayND::at: index 1 out of range (rank 7)");
03754 if (i2 >= shape_[2]) throw npstat::NpstatOutOfRange(
03755 "In npstat::ArrayND::at: index 2 out of range (rank 7)");
03756 if (i3 >= shape_[3]) throw npstat::NpstatOutOfRange(
03757 "In npstat::ArrayND::at: index 3 out of range (rank 7)");
03758 if (i4 >= shape_[4]) throw npstat::NpstatOutOfRange(
03759 "In npstat::ArrayND::at: index 4 out of range (rank 7)");
03760 if (i5 >= shape_[5]) throw npstat::NpstatOutOfRange(
03761 "In npstat::ArrayND::at: index 5 out of range (rank 7)");
03762 if (i6 >= shape_[6]) throw npstat::NpstatOutOfRange(
03763 "In npstat::ArrayND::at: index 6 out of range (rank 7)");
03764 return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] +
03765 i3*strides_[3] + i4*strides_[4] + i5*strides_[5] + i6];
03766 }
03767
03768 template<typename Numeric, unsigned Len, unsigned Dim>
03769 const Numeric& ArrayND<Numeric,Len,Dim>::at(
03770 const unsigned i0,
03771 const unsigned i1,
03772 const unsigned i2,
03773 const unsigned i3,
03774 const unsigned i4,
03775 const unsigned i5,
03776 const unsigned i6,
03777 const unsigned i7) const
03778 {
03779 if (8U != dim_) throw npstat::NpstatInvalidArgument(
03780 "In npstat::ArrayND::at: wrong # of args (not rank 8 array)");
03781 if (i0 >= shape_[0]) throw npstat::NpstatOutOfRange(
03782 "In npstat::ArrayND::at: index 0 out of range (rank 8)");
03783 if (i1 >= shape_[1]) throw npstat::NpstatOutOfRange(
03784 "In npstat::ArrayND::at: index 1 out of range (rank 8)");
03785 if (i2 >= shape_[2]) throw npstat::NpstatOutOfRange(
03786 "In npstat::ArrayND::at: index 2 out of range (rank 8)");
03787 if (i3 >= shape_[3]) throw npstat::NpstatOutOfRange(
03788 "In npstat::ArrayND::at: index 3 out of range (rank 8)");
03789 if (i4 >= shape_[4]) throw npstat::NpstatOutOfRange(
03790 "In npstat::ArrayND::at: index 4 out of range (rank 8)");
03791 if (i5 >= shape_[5]) throw npstat::NpstatOutOfRange(
03792 "In npstat::ArrayND::at: index 5 out of range (rank 8)");
03793 if (i6 >= shape_[6]) throw npstat::NpstatOutOfRange(
03794 "In npstat::ArrayND::at: index 6 out of range (rank 8)");
03795 if (i7 >= shape_[7]) throw npstat::NpstatOutOfRange(
03796 "In npstat::ArrayND::at: index 7 out of range (rank 8)");
03797 return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] +
03798 i3*strides_[3] + i4*strides_[4] + i5*strides_[5] +
03799 i6*strides_[6] + i7];
03800 }
03801
03802 template<typename Numeric, unsigned Len, unsigned Dim>
03803 const Numeric& ArrayND<Numeric,Len,Dim>::at(
03804 const unsigned i0,
03805 const unsigned i1,
03806 const unsigned i2,
03807 const unsigned i3,
03808 const unsigned i4,
03809 const unsigned i5,
03810 const unsigned i6,
03811 const unsigned i7,
03812 const unsigned i8) const
03813 {
03814 if (9U != dim_) throw npstat::NpstatInvalidArgument(
03815 "In npstat::ArrayND::at: wrong # of args (not rank 9 array)");
03816 if (i0 >= shape_[0]) throw npstat::NpstatOutOfRange(
03817 "In npstat::ArrayND::at: index 0 out of range (rank 9)");
03818 if (i1 >= shape_[1]) throw npstat::NpstatOutOfRange(
03819 "In npstat::ArrayND::at: index 1 out of range (rank 9)");
03820 if (i2 >= shape_[2]) throw npstat::NpstatOutOfRange(
03821 "In npstat::ArrayND::at: index 2 out of range (rank 9)");
03822 if (i3 >= shape_[3]) throw npstat::NpstatOutOfRange(
03823 "In npstat::ArrayND::at: index 3 out of range (rank 9)");
03824 if (i4 >= shape_[4]) throw npstat::NpstatOutOfRange(
03825 "In npstat::ArrayND::at: index 4 out of range (rank 9)");
03826 if (i5 >= shape_[5]) throw npstat::NpstatOutOfRange(
03827 "In npstat::ArrayND::at: index 5 out of range (rank 9)");
03828 if (i6 >= shape_[6]) throw npstat::NpstatOutOfRange(
03829 "In npstat::ArrayND::at: index 6 out of range (rank 9)");
03830 if (i7 >= shape_[7]) throw npstat::NpstatOutOfRange(
03831 "In npstat::ArrayND::at: index 7 out of range (rank 9)");
03832 if (i8 >= shape_[8]) throw npstat::NpstatOutOfRange(
03833 "In npstat::ArrayND::at: index 8 out of range (rank 9)");
03834 return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] +
03835 i3*strides_[3] + i4*strides_[4] + i5*strides_[5] +
03836 i6*strides_[6] + i7*strides_[7] + i8];
03837 }
03838
03839 template<typename Numeric, unsigned Len, unsigned Dim>
03840 const Numeric& ArrayND<Numeric,Len,Dim>::at(
03841 const unsigned i0,
03842 const unsigned i1,
03843 const unsigned i2,
03844 const unsigned i3,
03845 const unsigned i4,
03846 const unsigned i5,
03847 const unsigned i6,
03848 const unsigned i7,
03849 const unsigned i8,
03850 const unsigned i9) const
03851 {
03852 if (10U != dim_) throw npstat::NpstatInvalidArgument(
03853 "In npstat::ArrayND::at: wrong # of args (not rank 10 array)");
03854 if (i0 >= shape_[0]) throw npstat::NpstatOutOfRange(
03855 "In npstat::ArrayND::at: index 0 out of range (rank 10)");
03856 if (i1 >= shape_[1]) throw npstat::NpstatOutOfRange(
03857 "In npstat::ArrayND::at: index 1 out of range (rank 10)");
03858 if (i2 >= shape_[2]) throw npstat::NpstatOutOfRange(
03859 "In npstat::ArrayND::at: index 2 out of range (rank 10)");
03860 if (i3 >= shape_[3]) throw npstat::NpstatOutOfRange(
03861 "In npstat::ArrayND::at: index 3 out of range (rank 10)");
03862 if (i4 >= shape_[4]) throw npstat::NpstatOutOfRange(
03863 "In npstat::ArrayND::at: index 4 out of range (rank 10)");
03864 if (i5 >= shape_[5]) throw npstat::NpstatOutOfRange(
03865 "In npstat::ArrayND::at: index 5 out of range (rank 10)");
03866 if (i6 >= shape_[6]) throw npstat::NpstatOutOfRange(
03867 "In npstat::ArrayND::at: index 6 out of range (rank 10)");
03868 if (i7 >= shape_[7]) throw npstat::NpstatOutOfRange(
03869 "In npstat::ArrayND::at: index 7 out of range (rank 10)");
03870 if (i8 >= shape_[8]) throw npstat::NpstatOutOfRange(
03871 "In npstat::ArrayND::at: index 8 out of range (rank 10)");
03872 if (i9 >= shape_[9]) throw npstat::NpstatOutOfRange(
03873 "In npstat::ArrayND::at: index 9 out of range (rank 10)");
03874 return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] +
03875 i3*strides_[3] + i4*strides_[4] + i5*strides_[5] +
03876 i6*strides_[6] + i7*strides_[7] + i8*strides_[8] + i9];
03877 }
03878
03879 template<typename Numeric, unsigned Len, unsigned Dim>
03880 Numeric& ArrayND<Numeric,Len,Dim>::at(
03881 const unsigned i0,
03882 const unsigned i1,
03883 const unsigned i2)
03884 {
03885 if (3U != dim_) throw npstat::NpstatInvalidArgument(
03886 "In npstat::ArrayND::at: wrong # of args (not rank 3 array)");
03887 if (i0 >= shape_[0]) throw npstat::NpstatOutOfRange(
03888 "In npstat::ArrayND::at: index 0 out of range (rank 3)");
03889 if (i1 >= shape_[1]) throw npstat::NpstatOutOfRange(
03890 "In npstat::ArrayND::at: index 1 out of range (rank 3)");
03891 if (i2 >= shape_[2]) throw npstat::NpstatOutOfRange(
03892 "In npstat::ArrayND::at: index 2 out of range (rank 3)");
03893 return data_[i0*strides_[0] + i1*strides_[1] + i2];
03894 }
03895
03896 template<typename Numeric, unsigned Len, unsigned Dim>
03897 Numeric& ArrayND<Numeric,Len,Dim>::at(
03898 const unsigned i0,
03899 const unsigned i1,
03900 const unsigned i2,
03901 const unsigned i3)
03902 {
03903 if (4U != dim_) throw npstat::NpstatInvalidArgument(
03904 "In npstat::ArrayND::at: wrong # of args (not rank 4 array)");
03905 if (i0 >= shape_[0]) throw npstat::NpstatOutOfRange(
03906 "In npstat::ArrayND::at: index 0 out of range (rank 4)");
03907 if (i1 >= shape_[1]) throw npstat::NpstatOutOfRange(
03908 "In npstat::ArrayND::at: index 1 out of range (rank 4)");
03909 if (i2 >= shape_[2]) throw npstat::NpstatOutOfRange(
03910 "In npstat::ArrayND::at: index 2 out of range (rank 4)");
03911 if (i3 >= shape_[3]) throw npstat::NpstatOutOfRange(
03912 "In npstat::ArrayND::at: index 3 out of range (rank 4)");
03913 return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] + i3];
03914 }
03915
03916 template<typename Numeric, unsigned Len, unsigned Dim>
03917 Numeric& ArrayND<Numeric,Len,Dim>::at(
03918 const unsigned i0,
03919 const unsigned i1,
03920 const unsigned i2,
03921 const unsigned i3,
03922 const unsigned i4)
03923 {
03924 if (5U != dim_) throw npstat::NpstatInvalidArgument(
03925 "In npstat::ArrayND::at: wrong # of args (not rank 5 array)");
03926 if (i0 >= shape_[0]) throw npstat::NpstatOutOfRange(
03927 "In npstat::ArrayND::at: index 0 out of range (rank 5)");
03928 if (i1 >= shape_[1]) throw npstat::NpstatOutOfRange(
03929 "In npstat::ArrayND::at: index 1 out of range (rank 5)");
03930 if (i2 >= shape_[2]) throw npstat::NpstatOutOfRange(
03931 "In npstat::ArrayND::at: index 2 out of range (rank 5)");
03932 if (i3 >= shape_[3]) throw npstat::NpstatOutOfRange(
03933 "In npstat::ArrayND::at: index 3 out of range (rank 5)");
03934 if (i4 >= shape_[4]) throw npstat::NpstatOutOfRange(
03935 "In npstat::ArrayND::at: index 4 out of range (rank 5)");
03936 return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] +
03937 i3*strides_[3] + i4];
03938 }
03939
03940 template<typename Numeric, unsigned Len, unsigned Dim>
03941 Numeric& ArrayND<Numeric,Len,Dim>::at(
03942 const unsigned i0,
03943 const unsigned i1,
03944 const unsigned i2,
03945 const unsigned i3,
03946 const unsigned i4,
03947 const unsigned i5)
03948 {
03949 if (6U != dim_) throw npstat::NpstatInvalidArgument(
03950 "In npstat::ArrayND::at: wrong # of args (not rank 6 array)");
03951 if (i0 >= shape_[0]) throw npstat::NpstatOutOfRange(
03952 "In npstat::ArrayND::at: index 0 out of range (rank 6)");
03953 if (i1 >= shape_[1]) throw npstat::NpstatOutOfRange(
03954 "In npstat::ArrayND::at: index 1 out of range (rank 6)");
03955 if (i2 >= shape_[2]) throw npstat::NpstatOutOfRange(
03956 "In npstat::ArrayND::at: index 2 out of range (rank 6)");
03957 if (i3 >= shape_[3]) throw npstat::NpstatOutOfRange(
03958 "In npstat::ArrayND::at: index 3 out of range (rank 6)");
03959 if (i4 >= shape_[4]) throw npstat::NpstatOutOfRange(
03960 "In npstat::ArrayND::at: index 4 out of range (rank 6)");
03961 if (i5 >= shape_[5]) throw npstat::NpstatOutOfRange(
03962 "In npstat::ArrayND::at: index 5 out of range (rank 6)");
03963 return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] +
03964 i3*strides_[3] + i4*strides_[4] + i5];
03965 }
03966
03967 template<typename Numeric, unsigned Len, unsigned Dim>
03968 Numeric& ArrayND<Numeric,Len,Dim>::at(
03969 const unsigned i0,
03970 const unsigned i1,
03971 const unsigned i2,
03972 const unsigned i3,
03973 const unsigned i4,
03974 const unsigned i5,
03975 const unsigned i6)
03976 {
03977 if (7U != dim_) throw npstat::NpstatInvalidArgument(
03978 "In npstat::ArrayND::at: wrong # of args (not rank 7 array)");
03979 if (i0 >= shape_[0]) throw npstat::NpstatOutOfRange(
03980 "In npstat::ArrayND::at: index 0 out of range (rank 7)");
03981 if (i1 >= shape_[1]) throw npstat::NpstatOutOfRange(
03982 "In npstat::ArrayND::at: index 1 out of range (rank 7)");
03983 if (i2 >= shape_[2]) throw npstat::NpstatOutOfRange(
03984 "In npstat::ArrayND::at: index 2 out of range (rank 7)");
03985 if (i3 >= shape_[3]) throw npstat::NpstatOutOfRange(
03986 "In npstat::ArrayND::at: index 3 out of range (rank 7)");
03987 if (i4 >= shape_[4]) throw npstat::NpstatOutOfRange(
03988 "In npstat::ArrayND::at: index 4 out of range (rank 7)");
03989 if (i5 >= shape_[5]) throw npstat::NpstatOutOfRange(
03990 "In npstat::ArrayND::at: index 5 out of range (rank 7)");
03991 if (i6 >= shape_[6]) throw npstat::NpstatOutOfRange(
03992 "In npstat::ArrayND::at: index 6 out of range (rank 7)");
03993 return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] +
03994 i3*strides_[3] + i4*strides_[4] + i5*strides_[5] + i6];
03995 }
03996
03997 template<typename Numeric, unsigned Len, unsigned Dim>
03998 Numeric& ArrayND<Numeric,Len,Dim>::at(
03999 const unsigned i0,
04000 const unsigned i1,
04001 const unsigned i2,
04002 const unsigned i3,
04003 const unsigned i4,
04004 const unsigned i5,
04005 const unsigned i6,
04006 const unsigned i7)
04007 {
04008 if (8U != dim_) throw npstat::NpstatInvalidArgument(
04009 "In npstat::ArrayND::at: wrong # of args (not rank 8 array)");
04010 if (i0 >= shape_[0]) throw npstat::NpstatOutOfRange(
04011 "In npstat::ArrayND::at: index 0 out of range (rank 8)");
04012 if (i1 >= shape_[1]) throw npstat::NpstatOutOfRange(
04013 "In npstat::ArrayND::at: index 1 out of range (rank 8)");
04014 if (i2 >= shape_[2]) throw npstat::NpstatOutOfRange(
04015 "In npstat::ArrayND::at: index 2 out of range (rank 8)");
04016 if (i3 >= shape_[3]) throw npstat::NpstatOutOfRange(
04017 "In npstat::ArrayND::at: index 3 out of range (rank 8)");
04018 if (i4 >= shape_[4]) throw npstat::NpstatOutOfRange(
04019 "In npstat::ArrayND::at: index 4 out of range (rank 8)");
04020 if (i5 >= shape_[5]) throw npstat::NpstatOutOfRange(
04021 "In npstat::ArrayND::at: index 5 out of range (rank 8)");
04022 if (i6 >= shape_[6]) throw npstat::NpstatOutOfRange(
04023 "In npstat::ArrayND::at: index 6 out of range (rank 8)");
04024 if (i7 >= shape_[7]) throw npstat::NpstatOutOfRange(
04025 "In npstat::ArrayND::at: index 7 out of range (rank 8)");
04026 return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] +
04027 i3*strides_[3] + i4*strides_[4] + i5*strides_[5] +
04028 i6*strides_[6] + i7];
04029 }
04030
04031 template<typename Numeric, unsigned Len, unsigned Dim>
04032 Numeric& ArrayND<Numeric,Len,Dim>::at(
04033 const unsigned i0,
04034 const unsigned i1,
04035 const unsigned i2,
04036 const unsigned i3,
04037 const unsigned i4,
04038 const unsigned i5,
04039 const unsigned i6,
04040 const unsigned i7,
04041 const unsigned i8)
04042 {
04043 if (9U != dim_) throw npstat::NpstatInvalidArgument(
04044 "In npstat::ArrayND::at: wrong # of args (not rank 9 array)");
04045 if (i0 >= shape_[0]) throw npstat::NpstatOutOfRange(
04046 "In npstat::ArrayND::at: index 0 out of range (rank 9)");
04047 if (i1 >= shape_[1]) throw npstat::NpstatOutOfRange(
04048 "In npstat::ArrayND::at: index 1 out of range (rank 9)");
04049 if (i2 >= shape_[2]) throw npstat::NpstatOutOfRange(
04050 "In npstat::ArrayND::at: index 2 out of range (rank 9)");
04051 if (i3 >= shape_[3]) throw npstat::NpstatOutOfRange(
04052 "In npstat::ArrayND::at: index 3 out of range (rank 9)");
04053 if (i4 >= shape_[4]) throw npstat::NpstatOutOfRange(
04054 "In npstat::ArrayND::at: index 4 out of range (rank 9)");
04055 if (i5 >= shape_[5]) throw npstat::NpstatOutOfRange(
04056 "In npstat::ArrayND::at: index 5 out of range (rank 9)");
04057 if (i6 >= shape_[6]) throw npstat::NpstatOutOfRange(
04058 "In npstat::ArrayND::at: index 6 out of range (rank 9)");
04059 if (i7 >= shape_[7]) throw npstat::NpstatOutOfRange(
04060 "In npstat::ArrayND::at: index 7 out of range (rank 9)");
04061 if (i8 >= shape_[8]) throw npstat::NpstatOutOfRange(
04062 "In npstat::ArrayND::at: index 8 out of range (rank 9)");
04063 return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] +
04064 i3*strides_[3] + i4*strides_[4] + i5*strides_[5] +
04065 i6*strides_[6] + i7*strides_[7] + i8];
04066 }
04067
04068 template<typename Numeric, unsigned Len, unsigned Dim>
04069 Numeric& ArrayND<Numeric,Len,Dim>::at(
04070 const unsigned i0,
04071 const unsigned i1,
04072 const unsigned i2,
04073 const unsigned i3,
04074 const unsigned i4,
04075 const unsigned i5,
04076 const unsigned i6,
04077 const unsigned i7,
04078 const unsigned i8,
04079 const unsigned i9)
04080 {
04081 if (10U != dim_) throw npstat::NpstatInvalidArgument(
04082 "In npstat::ArrayND::at: wrong # of args (not rank 10 array)");
04083 if (i0 >= shape_[0]) throw npstat::NpstatOutOfRange(
04084 "In npstat::ArrayND::at: index 0 out of range (rank 10)");
04085 if (i1 >= shape_[1]) throw npstat::NpstatOutOfRange(
04086 "In npstat::ArrayND::at: index 1 out of range (rank 10)");
04087 if (i2 >= shape_[2]) throw npstat::NpstatOutOfRange(
04088 "In npstat::ArrayND::at: index 2 out of range (rank 10)");
04089 if (i3 >= shape_[3]) throw npstat::NpstatOutOfRange(
04090 "In npstat::ArrayND::at: index 3 out of range (rank 10)");
04091 if (i4 >= shape_[4]) throw npstat::NpstatOutOfRange(
04092 "In npstat::ArrayND::at: index 4 out of range (rank 10)");
04093 if (i5 >= shape_[5]) throw npstat::NpstatOutOfRange(
04094 "In npstat::ArrayND::at: index 5 out of range (rank 10)");
04095 if (i6 >= shape_[6]) throw npstat::NpstatOutOfRange(
04096 "In npstat::ArrayND::at: index 6 out of range (rank 10)");
04097 if (i7 >= shape_[7]) throw npstat::NpstatOutOfRange(
04098 "In npstat::ArrayND::at: index 7 out of range (rank 10)");
04099 if (i8 >= shape_[8]) throw npstat::NpstatOutOfRange(
04100 "In npstat::ArrayND::at: index 8 out of range (rank 10)");
04101 if (i9 >= shape_[9]) throw npstat::NpstatOutOfRange(
04102 "In npstat::ArrayND::at: index 9 out of range (rank 10)");
04103 return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] +
04104 i3*strides_[3] + i4*strides_[4] + i5*strides_[5] +
04105 i6*strides_[6] + i7*strides_[7] + i8*strides_[8] + i9];
04106 }
04107
04108 template<typename Numeric, unsigned Len, unsigned Dim>
04109 template<unsigned Len2, unsigned Dim2>
04110 double ArrayND<Numeric,Len,Dim>::maxAbsDifference(
04111 const ArrayND<Numeric,Len2,Dim2>& r) const
04112 {
04113 if (!isShapeCompatible(r)) throw npstat::NpstatInvalidArgument(
04114 "In npstat::ArrayND::maxAbsDifference: "
04115 "incompatible argument array shape");
04116 if (dim_)
04117 {
04118 double maxd = 0.0;
04119 for (unsigned long i=0; i<len_; ++i)
04120 {
04121 const Numeric rval = r.data_[i];
04122 const double d = absDifference(data_[i], rval);
04123 if (d > maxd)
04124 maxd = d;
04125 }
04126 return maxd;
04127 }
04128 else
04129 {
04130 const Numeric rval = r.localData_[0];
04131 return absDifference(localData_[0], rval);
04132 }
04133 }
04134
04135 template<typename Numeric, unsigned Len, unsigned Dim>
04136 template<unsigned Len2, unsigned Dim2>
04137 bool ArrayND<Numeric,Len,Dim>::operator==(
04138 const ArrayND<Numeric,Len2,Dim2>& r) const
04139 {
04140 if (shapeIsKnown_ != r.shapeIsKnown_)
04141 return false;
04142 if (r.dim_ != dim_)
04143 return false;
04144 if (r.len_ != len_)
04145 return false;
04146 for (unsigned i=0; i<dim_; ++i)
04147 if (shape_[i] != r.shape_[i])
04148 return false;
04149 for (unsigned i=0; i<dim_; ++i)
04150 assert(strides_[i] == r.strides_[i]);
04151 for (unsigned long j=0; j<len_; ++j)
04152 if (data_[j] != r.data_[j])
04153 return false;
04154 return true;
04155 }
04156
04157 template<typename Numeric, unsigned Len, unsigned Dim>
04158 template<unsigned Len2, unsigned Dim2>
04159 inline bool ArrayND<Numeric,Len,Dim>::operator!=(
04160 const ArrayND<Numeric,Len2,Dim2>& r) const
04161 {
04162 return !(*this == r);
04163 }
04164
04165 template<typename Numeric, unsigned Len, unsigned Dim>
04166 template<typename Num2>
04167 ArrayND<Numeric,Len,Dim>
04168 ArrayND<Numeric,Len,Dim>::operator*(const Num2& r) const
04169 {
04170 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
04171 "Initialize npstat::ArrayND before calling method \"operator*\"");
04172 ArrayND<Numeric,Len,Dim> result(shape_, dim_);
04173 for (unsigned long i=0; i<len_; ++i)
04174 result.data_[i] = data_[i]*r;
04175 return result;
04176 }
04177
04178 template<typename Numeric, unsigned Len, unsigned Dim>
04179 template<typename Num2>
04180 ArrayND<Numeric,Len,Dim>
04181 ArrayND<Numeric,Len,Dim>::operator/(const Num2& r) const
04182 {
04183 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
04184 "Initialize npstat::ArrayND before calling method \"operator/\"");
04185 if (r == Num2()) throw npstat::NpstatRuntimeError(
04186 "In npstat::ArrayND::operator/: division by zero");
04187 ArrayND<Numeric,Len,Dim> result(shape_, dim_);
04188 for (unsigned long i=0; i<len_; ++i)
04189 result.data_[i] = data_[i]/r;
04190 return result;
04191 }
04192
04193 template<typename Numeric, unsigned Len, unsigned Dim>
04194 template <unsigned Len2, unsigned Dim2>
04195 ArrayND<Numeric,Len,Dim>
04196 ArrayND<Numeric,Len,Dim>::operator+(
04197 const ArrayND<Numeric,Len2,Dim2>& r) const
04198 {
04199 if (!isShapeCompatible(r)) throw npstat::NpstatInvalidArgument(
04200 "In npstat::ArrayND::operator+: "
04201 "incompatible argument array shape");
04202 ArrayND<Numeric,Len,Dim> result(shape_, dim_);
04203 for (unsigned long i=0; i<len_; ++i)
04204 result.data_[i] = data_[i] + r.data_[i];
04205 return result;
04206 }
04207
04208 template<typename Numeric, unsigned Len, unsigned Dim>
04209 template <unsigned Len2, unsigned Dim2>
04210 ArrayND<Numeric,Len,Dim>
04211 ArrayND<Numeric,Len,Dim>::operator-(
04212 const ArrayND<Numeric,Len2,Dim2>& r) const
04213 {
04214 if (!isShapeCompatible(r)) throw npstat::NpstatInvalidArgument(
04215 "In npstat::ArrayND::operator-: "
04216 "incompatible argument array shape");
04217 ArrayND<Numeric,Len,Dim> result(shape_, dim_);
04218 for (unsigned long i=0; i<len_; ++i)
04219 result.data_[i] = data_[i] - r.data_[i];
04220 return result;
04221 }
04222
04223 template<typename Numeric, unsigned Len, unsigned Dim>
04224 inline ArrayND<Numeric,Len,Dim> ArrayND<Numeric,Len,Dim>::operator+() const
04225 {
04226 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
04227 "Initialize npstat::ArrayND before calling method \"operator+\"");
04228 return *this;
04229 }
04230
04231 template<typename Numeric, unsigned Len, unsigned Dim>
04232 ArrayND<Numeric,Len,Dim> ArrayND<Numeric,Len,Dim>::operator-() const
04233 {
04234 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
04235 "Initialize npstat::ArrayND before calling method \"operator-\"");
04236 ArrayND<Numeric,Len,Dim> result(shape_, dim_);
04237 for (unsigned long i=0; i<len_; ++i)
04238 result.data_[i] = -data_[i];
04239 return result;
04240 }
04241
04242 template<typename Numeric, unsigned Len, unsigned Dim>
04243 template <typename Num2>
04244 ArrayND<Numeric,Len,Dim>&
04245 ArrayND<Numeric,Len,Dim>::operator*=(const Num2& r)
04246 {
04247 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
04248 "Initialize npstat::ArrayND before calling method \"operator*=\"");
04249 for (unsigned long i=0; i<len_; ++i)
04250 data_[i] *= r;
04251 return *this;
04252 }
04253
04254 template<typename Numeric, unsigned Len, unsigned Dim>
04255 ArrayND<Numeric,Len,Dim>&
04256 ArrayND<Numeric,Len,Dim>::makeNonNegative()
04257 {
04258 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
04259 "Initialize npstat::ArrayND before calling method \"makeNonNegative\"");
04260 const Numeric zero = Numeric();
04261 if (dim_)
04262 {
04263 for (unsigned long i=0; i<len_; ++i)
04264 if (!(ComplexComparesAbs<Numeric>::more(data_[i], zero)))
04265 data_[i] = zero;
04266 }
04267 else
04268 if (!(ComplexComparesAbs<Numeric>::more(localData_[0], zero)))
04269 localData_[0] = zero;
04270 return *this;
04271 }
04272
04273 template<typename Numeric, unsigned Len, unsigned Dim>
04274 unsigned ArrayND<Numeric,Len,Dim>::makeCopulaSteps(
04275 const double tolerance, const unsigned nCycles)
04276 {
04277 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
04278 "Initialize npstat::ArrayND before calling method \"makeCopulaSteps\"");
04279 if (nCycles == 0U)
04280 return 0U;
04281 if (!dim_) throw npstat::NpstatInvalidArgument(
04282 "npstat::ArrayND::makeCopulaSteps method "
04283 "can not be used with array of 0 rank");
04284
04285 const Numeric zero = Numeric();
04286 for (unsigned long i=0; i<len_; ++i)
04287 if (!(ComplexComparesAbs<Numeric>::more(data_[i], zero)))
04288 data_[i] = zero;
04289
04290 std::vector<Numeric*> axesPtrBuf(dim_);
04291 Numeric** axes = &axesPtrBuf[0];
04292 const Numeric one = static_cast<Numeric>(1);
04293
04294
04295 unsigned idxSum = 0;
04296 for (unsigned i=0; i<dim_; ++i)
04297 idxSum += shape_[i];
04298 std::vector<Numeric> axesBuf(idxSum);
04299 axes[0] = &axesBuf[0];
04300 for (unsigned i=1; i<dim_; ++i)
04301 axes[i] = axes[i-1] + shape_[i-1];
04302
04303
04304 unsigned icycle = 0;
04305 for (; icycle<nCycles; ++icycle)
04306 {
04307 for (unsigned i=0; i<idxSum; ++i)
04308 axesBuf[i] = zero;
04309
04310
04311 for (unsigned long idat=0; idat<len_; ++idat)
04312 {
04313 unsigned long l = idat;
04314 for (unsigned i=0; i<dim_; ++i)
04315 {
04316 const unsigned idx = l / strides_[i];
04317 l -= (idx * strides_[i]);
04318 axes[i][idx] += data_[idat];
04319 }
04320 }
04321
04322
04323 bool withinTolerance = true;
04324 Numeric totalSum = zero;
04325 for (unsigned i=0; i<dim_; ++i)
04326 {
04327 Numeric axisSum = zero;
04328 const unsigned amax = shape_[i];
04329 for (unsigned a=0; a<amax; ++a)
04330 {
04331 if (axes[i][a] == zero)
04332 throw npstat::NpstatRuntimeError(
04333 "In npstat::ArrayND::makeCopulaSteps: "
04334 "marginal density is zero");
04335 axisSum += axes[i][a];
04336 }
04337 totalSum += axisSum;
04338 const Numeric axisAverage = axisSum/static_cast<Numeric>(amax);
04339 for (unsigned a=0; a<amax; ++a)
04340 axes[i][a] /= axisAverage;
04341 for (unsigned a=0; a<amax && withinTolerance; ++a)
04342 {
04343 const double adelta = absDifference(axes[i][a], one);
04344 if (adelta > tolerance)
04345 withinTolerance = false;
04346 }
04347 }
04348
04349 if (withinTolerance)
04350 break;
04351
04352 const Numeric totalAverage = totalSum/
04353 static_cast<Numeric>(len_)/static_cast<Numeric>(dim_);
04354
04355
04356
04357 for (unsigned long idat=0; idat<len_; ++idat)
04358 {
04359 unsigned long l = idat;
04360 for (unsigned i=0; i<dim_; ++i)
04361 {
04362 const unsigned idx = l / strides_[i];
04363 l -= (idx * strides_[i]);
04364 data_[idat] /= axes[i][idx];
04365 }
04366 data_[idat] /= totalAverage;
04367 }
04368 }
04369
04370 return icycle;
04371 }
04372
04373 template<typename Numeric, unsigned Len, unsigned Dim>
04374 template <typename Num2>
04375 ArrayND<Numeric,Len,Dim>&
04376 ArrayND<Numeric,Len,Dim>::operator/=(const Num2& r)
04377 {
04378 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
04379 "Initialize npstat::ArrayND before calling method \"operator/=\"");
04380 if (r == Num2()) throw npstat::NpstatRuntimeError(
04381 "In npstat::ArrayND::operator/=: division by zero");
04382 for (unsigned long i=0; i<len_; ++i)
04383 data_[i] /= r;
04384 return *this;
04385 }
04386
04387 template<typename Numeric, unsigned Len, unsigned Dim>
04388 template <typename Num2, unsigned Len2, unsigned Dim2>
04389 ArrayND<Numeric,Len,Dim>&
04390 ArrayND<Numeric,Len,Dim>::operator+=(const ArrayND<Num2,Len2,Dim2>& r)
04391 {
04392 if (!isShapeCompatible(r)) throw npstat::NpstatInvalidArgument(
04393 "In npstat::ArrayND::operator+=: "
04394 "incompatible argument array shape");
04395 for (unsigned long i=0; i<len_; ++i)
04396 data_[i] += r.data_[i];
04397 return *this;
04398 }
04399
04400 template<typename Numeric, unsigned Len, unsigned Dim>
04401 template <typename Num3, typename Num2, unsigned Len2, unsigned Dim2>
04402 ArrayND<Numeric,Len,Dim>&
04403 ArrayND<Numeric,Len,Dim>::addmul(const ArrayND<Num2,Len2,Dim2>& r,
04404 const Num3& c)
04405 {
04406 if (!isShapeCompatible(r)) throw npstat::NpstatInvalidArgument(
04407 "In npstat::ArrayND::addmul: "
04408 "incompatible argument array shape");
04409 for (unsigned long i=0; i<len_; ++i)
04410 data_[i] += r.data_[i]*c;
04411 return *this;
04412 }
04413
04414 template<typename Numeric, unsigned Len, unsigned Dim>
04415 template <typename Num2, unsigned Len2, unsigned Dim2>
04416 ArrayND<Numeric,Len,Dim>&
04417 ArrayND<Numeric,Len,Dim>::operator-=(const ArrayND<Num2,Len2,Dim2>& r)
04418 {
04419 if (!isShapeCompatible(r)) throw npstat::NpstatInvalidArgument(
04420 "In npstat::ArrayND::operator-=: "
04421 "incompatible argument array shape");
04422 for (unsigned long i=0; i<len_; ++i)
04423 data_[i] -= r.data_[i];
04424 return *this;
04425 }
04426
04427 template<typename Numeric, unsigned Len, unsigned Dim>
04428 Numeric ArrayND<Numeric,Len,Dim>::interpolate1(
04429 const double *coords, const unsigned dim) const
04430 {
04431 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
04432 "Initialize npstat::ArrayND before calling method \"interpolate1\"");
04433 if (dim != dim_) throw npstat::NpstatInvalidArgument(
04434 "In npstat::ArrayND::interpolate1: incompatible coordinate length");
04435 if (dim)
04436 {
04437 const unsigned maxdim = CHAR_BIT*sizeof(unsigned long);
04438 if (dim_ >= maxdim)
04439 throw npstat::NpstatInvalidArgument(
04440 "In npstat::ArrayND::interpolate1: array rank is too large");
04441
04442 double dx[maxdim];
04443 unsigned ix[maxdim];
04444 for (unsigned i=0; i<dim; ++i)
04445 {
04446 const double x = coords[i];
04447 if (x <= 0.0)
04448 {
04449 ix[i] = 0;
04450 dx[i] = 0.0;
04451 }
04452 else if (x >= static_cast<double>(shape_[i] - 1))
04453 {
04454 ix[i] = shape_[i] - 1;
04455 dx[i] = 0.0;
04456 }
04457 else
04458 {
04459 ix[i] = static_cast<unsigned>(std::floor(x));
04460 dx[i] = x - ix[i];
04461 }
04462 }
04463
04464 Numeric sum = Numeric();
04465 const unsigned long maxcycle = 1UL << dim;
04466 for (unsigned long icycle=0UL; icycle<maxcycle; ++icycle)
04467 {
04468 double w = 1.0;
04469 unsigned long icell = 0UL;
04470 for (unsigned i=0; i<dim; ++i)
04471 {
04472 if (icycle & (1UL << i))
04473 {
04474 w *= dx[i];
04475 icell += strides_[i]*(ix[i] + 1U);
04476 }
04477 else
04478 {
04479 w *= (1.0 - dx[i]);
04480 icell += strides_[i]*ix[i];
04481 }
04482 }
04483 if (w > 0.0)
04484 sum += data_[icell]*static_cast<proper_double>(w);
04485 }
04486 return sum;
04487 }
04488 else
04489 return localData_[0];
04490 }
04491
04492 template<typename Numeric, unsigned Len, unsigned Dim>
04493 Numeric ArrayND<Numeric,Len,Dim>::interpolateLoop(
04494 const unsigned level, const double* coords, const Numeric* base) const
04495 {
04496 const unsigned npoints = shape_[level];
04497 const double x = coords[level];
04498
04499 unsigned ix, npt = 1;
04500 double dx = 0.0;
04501 if (x < 0.0)
04502 ix = 0;
04503 else if (x > static_cast<double>(npoints - 1))
04504 ix = npoints - 1;
04505 else
04506 {
04507 ix = static_cast<unsigned>(std::floor(x));
04508 if (ix) --ix;
04509 unsigned imax = ix + 3;
04510 while (imax >= npoints)
04511 {
04512 if (ix) --ix;
04513 --imax;
04514 }
04515 dx = x - ix;
04516 npt = imax + 1 - ix;
04517 }
04518 assert(npt >= 1 && npt <= 4);
04519
04520 Numeric fit[4];
04521 if (level < dim_ - 1)
04522 for (unsigned ipt=0; ipt<npt; ++ipt)
04523 fit[ipt] = interpolateLoop(level + 1, coords,
04524 base + (ix + ipt)*strides_[level]);
04525
04526 const Numeric* const v = (level == dim_ - 1 ? base + ix : fit);
04527 switch (npt)
04528 {
04529 case 1:
04530 return v[0];
04531 case 2:
04532 return interpolate_linear(dx, v[0], v[1]);
04533 case 3:
04534 return interpolate_quadratic(dx, v[0], v[1], v[2]);
04535 case 4:
04536 return interpolate_cubic(dx, v[0], v[1], v[2], v[3]);
04537 default:
04538 assert(0);
04539 return Numeric();
04540 }
04541 }
04542
04543 template<typename Numeric, unsigned Len, unsigned Dim>
04544 inline Numeric ArrayND<Numeric,Len,Dim>::interpolate3(
04545 const double* coords, const unsigned dim) const
04546 {
04547 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
04548 "Initialize npstat::ArrayND before calling method \"interpolate3\"");
04549 if (dim != dim_) throw npstat::NpstatInvalidArgument(
04550 "In npstat::ArrayND::interpolate3: incompatible coordinate length");
04551 if (dim)
04552 {
04553 assert(coords);
04554 return interpolateLoop(0, coords, data_);
04555 }
04556 else
04557 return localData_[0];
04558 }
04559
04560 template<typename Numeric, unsigned Len, unsigned Dim>
04561 template<class Functor>
04562 ArrayND<Numeric,Len,Dim>& ArrayND<Numeric,Len,Dim>::apply(Functor f)
04563 {
04564 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
04565 "Initialize npstat::ArrayND before calling method \"apply\"");
04566 for (unsigned long i=0; i<len_; ++i)
04567 data_[i] = static_cast<Numeric>(f(data_[i]));
04568 return *this;
04569 }
04570
04571 template<typename Numeric, unsigned Len, unsigned Dim>
04572 template<class Functor>
04573 ArrayND<Numeric,Len,Dim>& ArrayND<Numeric,Len,Dim>::scanInPlace(
04574 Functor f)
04575 {
04576 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
04577 "Initialize npstat::ArrayND before calling method \"scanInPlace\"");
04578 for (unsigned long i=0; i<len_; ++i)
04579 f(data_[i]);
04580 return *this;
04581 }
04582
04583 template<typename Numeric, unsigned Len, unsigned Dim>
04584 ArrayND<Numeric,Len,Dim>& ArrayND<Numeric,Len,Dim>::constFill(
04585 const Numeric c)
04586 {
04587 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
04588 "Initialize npstat::ArrayND before calling method \"constFill\"");
04589 for (unsigned long i=0; i<len_; ++i)
04590 data_[i] = c;
04591 return *this;
04592 }
04593
04594 template<typename Numeric, unsigned Len, unsigned Dim>
04595 inline ArrayND<Numeric,Len,Dim>& ArrayND<Numeric,Len,Dim>::clear()
04596 {
04597 return constFill(Numeric());
04598 }
04599
04600 template<typename Numeric, unsigned Len, unsigned Dim>
04601 ArrayND<Numeric,Len,Dim>& ArrayND<Numeric,Len,Dim>::uninitialize()
04602 {
04603 destroyBuffer(data_, localData_);
04604 destroyBuffer(strides_, localStrides_);
04605 destroyBuffer(shape_, localShape_);
04606 localData_[0] = Numeric();
04607 data_ = localData_;
04608 strides_ = 0;
04609 shape_ = 0;
04610 len_ = 0;
04611 dim_ = 0;
04612 shapeIsKnown_ = false;
04613 return *this;
04614 }
04615
04616 template<typename Numeric, unsigned Len, unsigned Dim>
04617 ArrayND<Numeric,Len,Dim>& ArrayND<Numeric,Len,Dim>::makeUnit()
04618 {
04619 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
04620 "Initialize npstat::ArrayND before calling method \"makeUnit\"");
04621 if (dim_ < 2) throw npstat::NpstatInvalidArgument(
04622 "npstat::ArrayND::makeUnit method "
04623 "can not be used with arrays of rank less than 2");
04624 constFill(Numeric());
04625 unsigned long stride = 0UL;
04626 const unsigned dimlen = shape_[0];
04627 for (unsigned i=0; i<dim_; ++i)
04628 {
04629 if (shape_[i] != dimlen) throw npstat::NpstatInvalidArgument(
04630 "npstat::ArrayND::makeUnit method needs "
04631 "the array span to be the same in ech dimension");
04632 stride += strides_[i];
04633 }
04634 const Numeric one(static_cast<Numeric>(1));
04635 for (unsigned i=0; i<dimlen; ++i)
04636 data_[i*stride] = one;
04637 return *this;
04638 }
04639
04640 template<typename Numeric, unsigned Len, unsigned Dim>
04641 void ArrayND<Numeric,Len,Dim>::linearFillLoop(
04642 const unsigned level, const double s0, const unsigned long idx,
04643 const double shift, const double* coeffs)
04644 {
04645 const unsigned imax = shape_[level];
04646 const double c = coeffs[level];
04647 if (level == dim_ - 1)
04648 {
04649 Numeric* d = &data_[idx];
04650 for (unsigned i=0; i<imax; ++i)
04651 {
04652
04653
04654
04655 const double sum = s0 + c*i + shift;
04656 d[i] = static_cast<Numeric>(sum);
04657 }
04658 }
04659 else
04660 {
04661 const unsigned long stride = strides_[level];
04662 for (unsigned i=0; i<imax; ++i)
04663 linearFillLoop(level+1, s0 + c*i, idx + i*stride,
04664 shift, coeffs);
04665 }
04666 }
04667
04668 template<typename Numeric, unsigned Len, unsigned Dim>
04669 ArrayND<Numeric,Len,Dim>& ArrayND<Numeric,Len,Dim>::linearFill(
04670 const double* coeffs, const unsigned dimCoeffs, const double shift)
04671 {
04672
04673 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
04674 "Initialize npstat::ArrayND before calling method \"linearFill\"");
04675 if (dim_ != dimCoeffs) throw npstat::NpstatInvalidArgument(
04676 "In npstat::ArrayND::linearFill: incompatible number of coefficients");
04677 if (dim_)
04678 {
04679 assert(coeffs);
04680 linearFillLoop(0U, 0.0, 0UL, shift, coeffs);
04681 }
04682 else
04683 localData_[0] = static_cast<Numeric>(shift);
04684 return *this;
04685 }
04686
04687 template<typename Numeric, unsigned Len, unsigned Dim>
04688 template<class Functor>
04689 void ArrayND<Numeric,Len,Dim>::functorFillLoop(
04690 const unsigned level, const unsigned long idx,
04691 Functor f, unsigned* farg)
04692 {
04693 const unsigned imax = shape_[level];
04694 if (level == dim_ - 1)
04695 {
04696 Numeric* d = &data_[idx];
04697 const unsigned* myarg = farg;
04698 for (unsigned i = 0; i<imax; ++i)
04699 {
04700 farg[level] = i;
04701 d[i] = static_cast<Numeric>(f(myarg, dim_));
04702 }
04703 }
04704 else
04705 {
04706 const unsigned long stride = strides_[level];
04707 for (unsigned i = 0; i<imax; ++i)
04708 {
04709 farg[level] = i;
04710 functorFillLoop(level+1, idx + i*stride, f, farg);
04711 }
04712 }
04713 }
04714
04715 template<typename Numeric, unsigned Len, unsigned Dim>
04716 template <typename Accumulator>
04717 void ArrayND<Numeric,Len,Dim>::convertToLastDimCdfLoop(
04718 ArrayND* sumSlice, const unsigned level, unsigned long idx0,
04719 unsigned long idxSlice, const bool useTrapezoids)
04720 {
04721 static const proper_double half = 0.5;
04722 const unsigned imax = shape_[level];
04723 if (level == dim_ - 1)
04724 {
04725 Accumulator acc = Accumulator();
04726 Numeric* data = data_ + idx0;
04727 if (useTrapezoids)
04728 {
04729 Numeric oldval = Numeric();
04730 for (unsigned i = 0; i<imax; ++i)
04731 {
04732 acc += (data[i] + oldval)*half;
04733 oldval = data[i];
04734 data[i] = static_cast<Numeric>(acc);
04735 }
04736 acc += oldval*half;
04737 }
04738 else
04739 for (unsigned i = 0; i<imax; ++i)
04740 {
04741 acc += data[i];
04742 data[i] = static_cast<Numeric>(acc);
04743 }
04744 if (sumSlice->dim_)
04745 sumSlice->data_[idxSlice] = static_cast<Numeric>(acc);
04746 else
04747 sumSlice->localData_[0] = static_cast<Numeric>(acc);
04748 }
04749 else
04750 {
04751 const unsigned long stride = strides_[level];
04752 unsigned long sumStride = 0UL;
04753 if (sumSlice->dim_)
04754 sumStride = sumSlice->strides_[level];
04755 for (unsigned i = 0; i<imax; ++i)
04756 {
04757 convertToLastDimCdfLoop<Accumulator>(
04758 sumSlice, level+1, idx0, idxSlice, useTrapezoids);
04759 idx0 += stride;
04760 idxSlice += sumStride;
04761 }
04762 }
04763 }
04764
04765 template<typename Numeric, unsigned Len, unsigned Dim>
04766 template <typename Accumulator>
04767 inline void ArrayND<Numeric,Len,Dim>::convertToLastDimCdf(
04768 ArrayND* sumSlice, const bool useTrapezoids)
04769 {
04770 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
04771 "Initialize npstat::ArrayND before calling "
04772 "method \"convertToLastDimCdf\"");
04773 if (!dim_) throw npstat::NpstatInvalidArgument(
04774 "npstat::ArrayND::convertToLastDimCdf method "
04775 "can not be used with array of 0 rank");
04776 assert(sumSlice);
04777 if (!sumSlice->shapeIsKnown_) throw npstat::NpstatInvalidArgument(
04778 "In npstat::ArrayND::convertToLastDimCdf: "
04779 "uninitialized argument array");
04780 convertToLastDimCdfLoop<Accumulator>(sumSlice, 0U, 0UL, 0UL,
04781 useTrapezoids);
04782 }
04783
04784 template<typename Numeric, unsigned Len, unsigned Dim>
04785 template<class Functor>
04786 ArrayND<Numeric,Len,Dim>& ArrayND<Numeric,Len,Dim>::functorFill(Functor f)
04787 {
04788 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
04789 "Initialize npstat::ArrayND before calling method \"functorFill\"");
04790 if (dim_)
04791 {
04792 unsigned localIndex[Dim];
04793 unsigned* index = makeBuffer(dim_, localIndex, Dim);
04794 functorFillLoop(0U, 0UL, f, index);
04795 destroyBuffer(index, localIndex);
04796 }
04797 else
04798 localData_[0] = static_cast<Numeric>(
04799 f(static_cast<unsigned*>(0), 0U));
04800 return *this;
04801 }
04802
04803 template<typename Numeric, unsigned Len, unsigned Dim>
04804 template <typename Num2, unsigned Len2, unsigned Dim2>
04805 bool ArrayND<Numeric,Len,Dim>::isClose(
04806 const ArrayND<Num2,Len2,Dim2>& r, const double eps) const
04807 {
04808 if (eps < 0.0) throw npstat::NpstatDomainError(
04809 "In npstat::ArrayND::isClose: tolerance must not be negative");
04810 if (!isShapeCompatible(r)) throw npstat::NpstatInvalidArgument(
04811 "In npstat::ArrayND::isClose: incompatible argument array shape");
04812 if (dim_)
04813 {
04814 for (unsigned long i=0; i<len_; ++i)
04815 {
04816 const Numeric rval = r.data_[i];
04817 if (static_cast<double>(absDifference(data_[i], rval)) > eps)
04818 return false;
04819 }
04820 }
04821 else
04822 {
04823 const Numeric rval = r.localData_[0];
04824 if (static_cast<double>(absDifference(localData_[0], rval)) > eps)
04825 return false;
04826 }
04827 return true;
04828 }
04829
04830 template<typename Numeric, unsigned Len, unsigned Dim>
04831 template<typename Num2, unsigned Len2, unsigned Dim2>
04832 ArrayND<Numeric,Len,Dim> ArrayND<Numeric,Len,Dim>::outer(
04833 const ArrayND<Num2,Len2,Dim2>& r) const
04834 {
04835 return ArrayND<Numeric,Len,Dim>(*this, r);
04836 }
04837
04838 template<typename Numeric, unsigned Len, unsigned Dim>
04839 void ArrayND<Numeric,Len,Dim>::contractLoop(
04840 unsigned thisLevel, const unsigned resLevel,
04841 const unsigned pos1, const unsigned pos2,
04842 unsigned long idxThis, unsigned long idxRes, ArrayND& result) const
04843 {
04844 while (thisLevel == pos1 || thisLevel == pos2)
04845 ++thisLevel;
04846 assert(thisLevel < dim_);
04847
04848 if (resLevel == result.dim_ - 1)
04849 {
04850 const unsigned ncontract = shape_[pos1];
04851 const unsigned imax = result.shape_[resLevel];
04852 const unsigned long stride = strides_[pos1] + strides_[pos2];
04853 for (unsigned i=0; i<imax; ++i)
04854 {
04855 const Numeric* tmp = data_ + (idxThis + i*strides_[thisLevel]);
04856 Numeric sum = Numeric();
04857 for (unsigned j=0; j<ncontract; ++j)
04858 sum += tmp[j*stride];
04859 result.data_[idxRes + i] = sum;
04860 }
04861 }
04862 else
04863 {
04864 const unsigned imax = result.shape_[resLevel];
04865 assert(imax == shape_[thisLevel]);
04866 for (unsigned i=0; i<imax; ++i)
04867 {
04868 contractLoop(thisLevel+1, resLevel+1, pos1, pos2,
04869 idxThis, idxRes, result);
04870 idxThis += strides_[thisLevel];
04871 idxRes += result.strides_[resLevel];
04872 }
04873 }
04874 }
04875
04876 template<typename Numeric, unsigned Len, unsigned Dim>
04877 ArrayND<Numeric,Len,Dim> ArrayND<Numeric,Len,Dim>::contract(
04878 const unsigned pos1, const unsigned pos2) const
04879 {
04880 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
04881 "Initialize npstat::ArrayND before calling method \"contract\"");
04882 if (!(pos1 < dim_ && pos2 < dim_ && pos1 != pos2))
04883 throw npstat::NpstatInvalidArgument("In npstat::ArrayND::contract: "
04884 "incompatible contraction indices");
04885 if (shape_[pos1] != shape_[pos2])
04886 throw npstat::NpstatInvalidArgument(
04887 "In npstat::ArrayND::contract: incompatible "
04888 "length of contracted dimensions");
04889
04890
04891 unsigned newshapeBuf[Dim];
04892 unsigned* newshape = makeBuffer(dim_ - 2, newshapeBuf, Dim);
04893 unsigned ishap = 0;
04894 for (unsigned i=0; i<dim_; ++i)
04895 if (i != pos1 && i != pos2)
04896 newshape[ishap++] = shape_[i];
04897
04898
04899 ArrayND<Numeric,Len,Dim> result(newshape, ishap);
04900 if (ishap)
04901 contractLoop(0, 0, pos1, pos2, 0UL, 0UL, result);
04902 else
04903 {
04904
04905 Numeric sum = Numeric();
04906 const unsigned imax = shape_[0];
04907 const unsigned long stride = strides_[0] + strides_[1];
04908 for (unsigned i=0; i<imax; ++i)
04909 sum += data_[i*stride];
04910 result() = sum;
04911 }
04912
04913 destroyBuffer(newshape, newshapeBuf);
04914 return result;
04915 }
04916
04917 template<typename Numeric, unsigned Len, unsigned Dim>
04918 void ArrayND<Numeric,Len,Dim>::transposeLoop(
04919 const unsigned level, const unsigned pos1, const unsigned pos2,
04920 unsigned long idxThis, unsigned long idxRes, ArrayND& result) const
04921 {
04922 const unsigned imax = shape_[level];
04923 const unsigned long mystride = strides_[level];
04924 const unsigned relevel = level == pos1 ? pos2 :
04925 (level == pos2 ? pos1 : level);
04926 const unsigned long restride = result.strides_[relevel];
04927 const bool ready = (level == dim_ - 1);
04928 for (unsigned i=0; i<imax; ++i)
04929 {
04930 if (ready)
04931 result.data_[idxRes] = data_[idxThis];
04932 else
04933 transposeLoop(level+1, pos1, pos2, idxThis, idxRes, result);
04934 idxThis += mystride;
04935 idxRes += restride;
04936 }
04937 }
04938
04939 template<typename Numeric, unsigned Len, unsigned Dim>
04940 template<typename Num2>
04941 Num2 ArrayND<Numeric,Len,Dim>::sum() const
04942 {
04943 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
04944 "Initialize npstat::ArrayND before calling method \"sum\"");
04945 Num2 sum = Num2();
04946 for (unsigned long i=0; i<len_; ++i)
04947 sum += data_[i];
04948 return sum;
04949 }
04950
04951 template<typename Numeric, unsigned Len, unsigned Dim>
04952 template<typename Num2>
04953 Num2 ArrayND<Numeric,Len,Dim>::sumsq() const
04954 {
04955 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
04956 "Initialize npstat::ArrayND before calling method \"sumsq\"");
04957 Num2 sum = Num2();
04958 for (unsigned long i=0; i<len_; ++i)
04959 {
04960 const Num2 absval = absValue(data_[i]);
04961 sum += absval*absval;
04962 }
04963 return sum;
04964 }
04965
04966 template<typename Numeric, unsigned Len, unsigned Dim>
04967 Numeric ArrayND<Numeric,Len,Dim>::min() const
04968 {
04969 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
04970 "Initialize npstat::ArrayND before calling method \"min\"");
04971 if (dim_)
04972 {
04973 Numeric minval(data_[0]);
04974 for (unsigned long i=1UL; i<len_; ++i)
04975 if (ComplexComparesAbs<Numeric>::less(data_[i], minval))
04976 minval = data_[i];
04977 return minval;
04978 }
04979 else
04980 return localData_[0];
04981 }
04982
04983 template<typename Numeric, unsigned Len, unsigned Dim>
04984 Numeric ArrayND<Numeric,Len,Dim>::min(
04985 unsigned *index, const unsigned indexLen) const
04986 {
04987 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
04988 "Initialize npstat::ArrayND before calling method \"min\"");
04989 if (indexLen != dim_) throw npstat::NpstatInvalidArgument(
04990 "In npstat::ArrayND::min: incompatible index length");
04991 if (dim_)
04992 {
04993 unsigned long minind = 0UL;
04994 Numeric minval(data_[0]);
04995 for (unsigned long i=1UL; i<len_; ++i)
04996 if (ComplexComparesAbs<Numeric>::less(data_[i], minval))
04997 {
04998 minval = data_[i];
04999 minind = i;
05000 }
05001 convertLinearIndex(minind, index, indexLen);
05002 return minval;
05003 }
05004 else
05005 return localData_[0];
05006 }
05007
05008 template<typename Numeric, unsigned Len, unsigned Dim>
05009 Numeric ArrayND<Numeric,Len,Dim>::max() const
05010 {
05011 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
05012 "Initialize npstat::ArrayND before calling method \"max\"");
05013 if (dim_)
05014 {
05015 Numeric maxval(data_[0]);
05016 for (unsigned long i=1UL; i<len_; ++i)
05017 if (ComplexComparesAbs<Numeric>::less(maxval, data_[i]))
05018 maxval = data_[i];
05019 return maxval;
05020 }
05021 else
05022 return localData_[0];
05023 }
05024
05025 template<typename Numeric, unsigned Len, unsigned Dim>
05026 Numeric ArrayND<Numeric,Len,Dim>::max(
05027 unsigned *index, const unsigned indexLen) const
05028 {
05029 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
05030 "Initialize npstat::ArrayND before calling method \"max\"");
05031 if (indexLen != dim_) throw npstat::NpstatInvalidArgument(
05032 "In npstat::ArrayND::max: incompatible index length");
05033 if (dim_)
05034 {
05035 unsigned long maxind = 0UL;
05036 Numeric maxval(data_[0]);
05037 for (unsigned long i=1UL; i<len_; ++i)
05038 if (ComplexComparesAbs<Numeric>::less(maxval, data_[i]))
05039 {
05040 maxval = data_[i];
05041 maxind = i;
05042 }
05043 convertLinearIndex(maxind, index, indexLen);
05044 return maxval;
05045 }
05046 else
05047 return localData_[0];
05048 }
05049
05050
05051 template<typename Numeric, unsigned Len, unsigned Dim>
05052 ArrayND<Numeric,Len,Dim> ArrayND<Numeric,Len,Dim>::transpose() const
05053 {
05054 if (dim_ != 2) throw npstat::NpstatInvalidArgument(
05055 "npstat::ArrayND::transpose method "
05056 "can not be used with arrays of rank other than 2");
05057 unsigned newshape[2];
05058 newshape[0] = shape_[1];
05059 newshape[1] = shape_[0];
05060 ArrayND<Numeric,Len,Dim> result(newshape, dim_);
05061 const unsigned imax = shape_[0];
05062 const unsigned jmax = shape_[1];
05063 for (unsigned i=0; i<imax; ++i)
05064 for (unsigned j=0; j<jmax; ++j)
05065 result.data_[j*imax + i] = data_[i*jmax + j];
05066 return result;
05067 }
05068
05069 template<typename Numeric, unsigned Len, unsigned Dim>
05070 template <typename Accumulator>
05071 ArrayND<Numeric,Len,Dim> ArrayND<Numeric,Len,Dim>::derivative(
05072 const double inscale) const
05073 {
05074 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
05075 "Initialize npstat::ArrayND before calling method \"derivative\"");
05076 if (!dim_) throw npstat::NpstatInvalidArgument(
05077 "npstat::ArrayND::derivative method "
05078 "can not be used with array of 0 rank");
05079
05080 const typename ProperDblFromCmpl<Accumulator>::type scale = inscale;
05081 const unsigned maxdim = CHAR_BIT*sizeof(unsigned long);
05082 if (dim_ >= maxdim) throw npstat::NpstatInvalidArgument(
05083 "In npstat::ArrayND::derivative: array rank is too large");
05084 const unsigned long maxcycle = 1UL << dim_;
05085
05086 ArrayShape sh;
05087 sh.reserve(dim_);
05088 for (unsigned i=0; i<dim_; ++i)
05089 {
05090 if (shape_[i] <= 1U)
05091 throw npstat::NpstatInvalidArgument(
05092 "In npstat::ArrayND::derivative: in some dimendions "
05093 "array size is too small");
05094 sh.push_back(shape_[i] - 1U);
05095 }
05096
05097 ArrayND result(sh);
05098 const unsigned long rLen = result.length();
05099 for (unsigned long ilin=0; ilin<rLen; ++ilin)
05100 {
05101 result.convertLinearIndex(ilin, &sh[0], dim_);
05102
05103 Accumulator deriv = Accumulator();
05104 for (unsigned long icycle=0UL; icycle<maxcycle; ++icycle)
05105 {
05106 unsigned long icell = 0UL;
05107 unsigned n1 = 0U;
05108 for (unsigned i=0; i<dim_; ++i)
05109 {
05110 if (icycle & (1UL << i))
05111 {
05112 ++n1;
05113 icell += strides_[i]*(sh[i] + 1);
05114 }
05115 else
05116 icell += strides_[i]*sh[i];
05117 }
05118 if ((dim_ - n1) % 2U)
05119 deriv -= data_[icell];
05120 else
05121 deriv += data_[icell];
05122 }
05123 result.data_[ilin] = static_cast<Numeric>(deriv*scale);
05124 }
05125
05126 return result;
05127 }
05128
05129 template<typename Numeric, unsigned Len, unsigned Dim>
05130 template <typename Accumulator>
05131 Accumulator ArrayND<Numeric,Len,Dim>::sumBelowLoop(
05132 const unsigned level, unsigned long idx0,
05133 const unsigned* limit) const
05134 {
05135 Accumulator cdf = Accumulator();
05136 const unsigned imax = limit[level] + 1U;
05137 if (level == dim_ - 1)
05138 {
05139 Numeric* base = data_ + idx0;
05140 for (unsigned i=0; i<imax; ++i)
05141 cdf += base[i];
05142 }
05143 else
05144 {
05145 const unsigned long stride = strides_[level];
05146 for (unsigned i=0; i<imax; ++i, idx0+=stride)
05147 cdf += sumBelowLoop<Accumulator>(level+1, idx0, limit);
05148 }
05149 return cdf;
05150 }
05151
05152 template<typename Numeric, unsigned Len, unsigned Dim>
05153 template <typename Accumulator>
05154 Accumulator ArrayND<Numeric,Len,Dim>::cdfValue(
05155 const unsigned *index, const unsigned indexLen) const
05156 {
05157 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
05158 "Initialize npstat::ArrayND before calling method \"cdfValue\"");
05159 if (!dim_) throw npstat::NpstatInvalidArgument(
05160 "npstat::ArrayND::cdfValue method "
05161 "can not be used with array of 0 rank");
05162 if (indexLen != dim_) throw npstat::NpstatInvalidArgument(
05163 "In npstat::ArrayND::cdfValue: incompatible index length");
05164 for (unsigned i=0; i<indexLen; ++i)
05165 if (index[i] >= shape_[i])
05166 throw npstat::NpstatOutOfRange(
05167 "In npstat::ArrayND::cdfValue: index out of range");
05168 return sumBelowLoop<Accumulator>(0, 0U, index);
05169 }
05170
05171 template<typename Numeric, unsigned Len, unsigned Dim>
05172 template <typename Accumulator>
05173 ArrayND<Numeric,Len,Dim> ArrayND<Numeric,Len,Dim>::cdfArray(
05174 const double inscale) const
05175 {
05176 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
05177 "Initialize npstat::ArrayND before calling method \"cdfArray\"");
05178 if (!dim_) throw npstat::NpstatInvalidArgument(
05179 "npstat::ArrayND::cdfArray method "
05180 "can not be used with array of 0 rank");
05181
05182 const proper_double scale = inscale;
05183 const unsigned maxdim = CHAR_BIT*sizeof(unsigned long);
05184 if (dim_ >= maxdim)
05185 throw npstat::NpstatInvalidArgument(
05186 "In npstat::ArrayND::cdfArray: array rank is too large");
05187 const unsigned long maxcycle = 1UL << dim_;
05188
05189 ArrayShape sh;
05190 sh.reserve(dim_);
05191 for (unsigned i=0; i<dim_; ++i)
05192 sh.push_back(shape_[i] + 1U);
05193
05194 ArrayND<Accumulator> result(sh);
05195
05196 unsigned* psh = &sh[0];
05197 const unsigned long len = result.length();
05198 for (unsigned long ipre=0; ipre<len; ++ipre)
05199 {
05200 result.convertLinearIndex(ipre, psh, dim_);
05201 Accumulator deriv = Accumulator();
05202 bool has0 = false;
05203 for (unsigned i=0; i<dim_; ++i)
05204 if (psh[i]-- == 0U)
05205 {
05206 has0 = true;
05207 break;
05208 }
05209 if (!has0)
05210 {
05211 for (unsigned long icycle=0UL; icycle<maxcycle; ++icycle)
05212 {
05213 unsigned long icell = 0UL;
05214 unsigned n1 = 0U;
05215 for (unsigned i=0; i<dim_; ++i)
05216 {
05217 if (icycle & (1UL << i))
05218 {
05219 ++n1;
05220 icell += result.strides_[i]*(psh[i] + 1);
05221 }
05222 else
05223 icell += result.strides_[i]*psh[i];
05224 }
05225 if (n1 < dim_)
05226 {
05227 if ((dim_ - n1) % 2U)
05228 deriv += result.data_[icell];
05229 else
05230 deriv -= result.data_[icell];
05231 }
05232 }
05233 deriv += static_cast<Accumulator>(value(psh, dim_)*scale);
05234 }
05235 result.data_[ipre] = deriv;
05236 }
05237
05238
05239 return result;
05240 }
05241
05242 template<typename Numeric, unsigned Len, unsigned Dim>
05243 ArrayND<Numeric,Len,Dim> ArrayND<Numeric,Len,Dim>::transpose(
05244 const unsigned pos1, const unsigned pos2) const
05245 {
05246 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
05247 "Initialize npstat::ArrayND before calling method \"transpose\"");
05248 if (!(pos1 < dim_ && pos2 < dim_ && pos1 != pos2))
05249 throw npstat::NpstatInvalidArgument("In npstat::ArrayND::transpose: "
05250 "incompatible transposition indices");
05251 if (dim_ == 2)
05252 return transpose();
05253 else
05254 {
05255
05256 unsigned newshapeBuf[Dim];
05257 unsigned *newshape = makeBuffer(dim_, newshapeBuf, Dim);
05258 copyBuffer(newshape, shape_, dim_);
05259 std::swap(newshape[pos1], newshape[pos2]);
05260
05261
05262 ArrayND<Numeric,Len,Dim> result(newshape, dim_);
05263
05264
05265 transposeLoop(0, pos1, pos2, 0UL, 0UL, result);
05266
05267 destroyBuffer(newshape, newshapeBuf);
05268 return result;
05269 }
05270 }
05271
05272 template<typename Numeric, unsigned Len, unsigned Dim>
05273 template<typename Num2, unsigned Len2, unsigned Dim2>
05274 void ArrayND<Numeric,Len,Dim>::multiMirror(
05275 ArrayND<Num2, Len2, Dim2>* out) const
05276 {
05277 assert(out);
05278 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
05279 "Initialize npstat::ArrayND before calling method \"multiMirror\"");
05280 if (!out->shapeIsKnown_)
05281 *out = ArrayND<Num2, Len2, Dim2>(doubleShape(shape()));
05282 if (dim_ != out->dim_) throw npstat::NpstatInvalidArgument(
05283 "In npstat::ArrayND::multiMirror: incompatible argument array rank");
05284
05285 if (dim_)
05286 {
05287 const unsigned *dshape = out->shape_;
05288 for (unsigned i=0; i<dim_; ++i)
05289 if (dshape[i] != shape_[i]*2U) throw npstat::NpstatInvalidArgument(
05290 "In npstat::ArrayND::multiMirror: "
05291 "incompatible argument array shape");
05292
05293 if (dim_ >= CHAR_BIT*sizeof(unsigned long))
05294 throw npstat::NpstatInvalidArgument(
05295 "In npstat::ArrayND::multiMirror: "
05296 "array rank is too large");
05297 const unsigned long maxcycle = 1UL << dim_;
05298 std::vector<unsigned> indexbuf(dim_*2U);
05299 unsigned* idx = &indexbuf[0];
05300 unsigned* mirror = idx + dim_;
05301
05302 for (unsigned long ipt=0; ipt<len_; ++ipt)
05303 {
05304 unsigned long l = ipt;
05305 for (unsigned i=0; i<dim_; ++i)
05306 {
05307 idx[i] = l / strides_[i];
05308 l -= (idx[i] * strides_[i]);
05309 }
05310 for (unsigned long icycle=0UL; icycle<maxcycle; ++icycle)
05311 {
05312 for (unsigned i=0; i<dim_; ++i)
05313 {
05314 if (icycle & (1UL << i))
05315 mirror[i] = dshape[i] - idx[i] - 1U;
05316 else
05317 mirror[i] = idx[i];
05318 }
05319 out->value(mirror, dim_) = data_[ipt];
05320 }
05321 }
05322 }
05323 else
05324 out->localData_[0] = static_cast<Num2>(localData_[0]);
05325 }
05326
05327 template<typename Numeric, unsigned Len, unsigned Dim>
05328 template<typename Num2, unsigned Len2, unsigned Dim2>
05329 void ArrayND<Numeric,Len,Dim>::rotate(
05330 const unsigned* shifts, const unsigned lenShifts,
05331 ArrayND<Num2, Len2, Dim2>* rotated) const
05332 {
05333 assert(rotated);
05334 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
05335 "Initialize npstat::ArrayND before calling method \"rotate\"");
05336
05337 if ((void*)rotated == (void*)this) throw npstat::NpstatInvalidArgument(
05338 "In npstat::ArrayND::rotate: can not rotate array into itself");
05339 if (!rotated->shapeIsKnown_)
05340 *rotated = *this;
05341 if (dim_ != rotated->dim_) throw npstat::NpstatInvalidArgument(
05342 "In npstat::ArrayND::rotate: incompatible argument array rank");
05343 if (lenShifts != dim_) throw npstat::NpstatInvalidArgument(
05344 "In npstat::ArrayND::rotate: incompatible dimensionality of shifts");
05345
05346 if (dim_)
05347 {
05348 assert(shifts);
05349 if (dim_ > CHAR_BIT*sizeof(unsigned long))
05350 throw npstat::NpstatInvalidArgument(
05351 "In npstat::ArrayND::rotate: array rank is too large");
05352 unsigned buf[CHAR_BIT*sizeof(unsigned long)];
05353 clearBuffer(buf, dim_);
05354 (const_cast<ArrayND*>(this))->flatCircularLoop(
05355 0U, 0UL, 0UL, buf, shape_, shifts,
05356 *rotated, scast_assign_right<Numeric,Num2>());
05357 }
05358 else
05359 rotated->localData_[0] = static_cast<Num2>(localData_[0]);
05360 }
05361
05362 template<typename Numeric, unsigned Len, unsigned Dim>
05363 template<typename Num2, unsigned Len2, unsigned Dim2>
05364 void ArrayND<Numeric,Len,Dim>::dotProductLoop(
05365 const unsigned level, unsigned long idx0,
05366 unsigned long idx1, unsigned long idx2,
05367 const ArrayND<Num2, Len2, Dim2>& r,
05368 ArrayND& result) const
05369 {
05370
05371
05372
05373 if (level == result.dim_)
05374 {
05375 Numeric sum = Numeric();
05376 const unsigned imax = r.shape_[0];
05377 const unsigned rstride = r.strides_[0];
05378 const Numeric* l = data_ + idx0;
05379 const Num2* ri = r.data_ + idx1;
05380 for (unsigned i=0; i<imax; ++i)
05381 sum += l[i]*ri[i*rstride];
05382 result.data_[idx2] = sum;
05383 }
05384 else
05385 {
05386 const unsigned imax = result.shape_[level];
05387 for (unsigned i=0; i<imax; ++i)
05388 {
05389 dotProductLoop(level+1, idx0, idx1, idx2, r, result);
05390 idx2 += result.strides_[level];
05391 if (level < dim_ - 1)
05392 idx0 += strides_[level];
05393 else
05394 idx1 += r.strides_[level + 2 - dim_];
05395 }
05396 }
05397 }
05398
05399 template<typename Numeric, unsigned Len, unsigned Dim>
05400 template<typename Num2, unsigned Len2, unsigned Dim2>
05401 ArrayND<Numeric,Len,Dim> ArrayND<Numeric,Len,Dim>::dot(
05402 const ArrayND<Num2,Len2,Dim2>& r) const
05403 {
05404 if (!dim_) throw npstat::NpstatInvalidArgument(
05405 "npstat::ArrayND::dot method "
05406 "can not be used with array of 0 rank");
05407 if (!r.dim_) throw npstat::NpstatInvalidArgument(
05408 "npstat::ArrayND::dot method "
05409 "can not be used with argument array of 0 rank");
05410 if (shape_[dim_ - 1] != r.shape_[0]) throw npstat::NpstatInvalidArgument(
05411 "In npstat::ArrayND::dot: incompatible argument array shape");
05412
05413 if (dim_ == 1 && r.dim_ == 1)
05414 {
05415
05416 ArrayND<Numeric,Len,Dim> result(static_cast<unsigned*>(0), 0U);
05417 Numeric sum = Numeric();
05418 const unsigned imax = shape_[0];
05419 for (unsigned i=0; i<imax; ++i)
05420 sum += data_[i]*r.data_[i];
05421 result() = sum;
05422 return result;
05423 }
05424 else
05425 {
05426 unsigned newshapeBuf[2*Dim];
05427 unsigned *newshape = makeBuffer(dim_+r.dim_-2, newshapeBuf, 2*Dim);
05428 copyBuffer(newshape, shape_, dim_-1);
05429 copyBuffer(newshape+(dim_-1), r.shape_+1, r.dim_-1);
05430 ArrayND<Numeric,Len,Dim> result(newshape, dim_+r.dim_-2);
05431
05432 dotProductLoop(0U, 0UL, 0UL, 0UL, r, result);
05433
05434 destroyBuffer(newshape, newshapeBuf);
05435 return result;
05436 }
05437 }
05438
05439 template<typename Numeric, unsigned Len, unsigned Dim>
05440 inline unsigned ArrayND<Numeric,Len,Dim>::span(const unsigned dim) const
05441 {
05442 if (dim >= dim_)
05443 throw npstat::NpstatOutOfRange(
05444 "In npstat::ArrayND::span: dimension number is out of range");
05445 return shape_[dim];
05446 }
05447
05448 template<typename Numeric, unsigned Len, unsigned Dim>
05449 unsigned ArrayND<Numeric,Len,Dim>::maximumSpan() const
05450 {
05451 unsigned maxspan = 0;
05452 for (unsigned i=0; i<dim_; ++i)
05453 if (shape_[i] > maxspan)
05454 maxspan = shape_[i];
05455 return maxspan;
05456 }
05457
05458 template<typename Numeric, unsigned Len, unsigned Dim>
05459 unsigned ArrayND<Numeric,Len,Dim>::minimumSpan() const
05460 {
05461 if (dim_ == 0)
05462 return 0U;
05463 unsigned minspan = shape_[0];
05464 for (unsigned i=1; i<dim_; ++i)
05465 if (shape_[i] < minspan)
05466 minspan = shape_[i];
05467 return minspan;
05468 }
05469
05470 template<typename Numeric, unsigned Len, unsigned Dim>
05471 inline const Numeric& ArrayND<Numeric,Len,Dim>::cl() const
05472 {
05473 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
05474 "Initialize npstat::ArrayND before calling method \"cl\"");
05475 if (dim_) throw npstat::NpstatInvalidArgument(
05476 "In npstat::ArrayND::cl: wrong # of args (not rank 0 array)");
05477 return localData_[0];
05478 }
05479
05480 template<typename Numeric, unsigned Len, unsigned Dim>
05481 inline const Numeric& ArrayND<Numeric,Len,Dim>::cl(
05482 const double i0) const
05483 {
05484 if (1U != dim_) throw npstat::NpstatInvalidArgument(
05485 "In npstat::ArrayND::cl: wrong # of args (not rank 1 array)");
05486 return data_[coordToIndex(i0, 0)];
05487 }
05488
05489 template<typename Numeric, unsigned Len, unsigned Dim>
05490 inline const Numeric& ArrayND<Numeric,Len,Dim>::cl(
05491 const double i0,
05492 const double i1) const
05493 {
05494 if (2U != dim_) throw npstat::NpstatInvalidArgument(
05495 "In npstat::ArrayND::cl: wrong # of args (not rank 2 array)");
05496 return data_[coordToIndex(i0, 0)*strides_[0] +
05497 coordToIndex(i1, 1)];
05498 }
05499
05500 template<typename Numeric, unsigned Len, unsigned Dim>
05501 inline const Numeric& ArrayND<Numeric,Len,Dim>::cl(
05502 const double i0,
05503 const double i1,
05504 const double i2) const
05505 {
05506 if (3U != dim_) throw npstat::NpstatInvalidArgument(
05507 "In npstat::ArrayND::cl: wrong # of args (not rank 3 array)");
05508 return data_[coordToIndex(i0, 0)*strides_[0] +
05509 coordToIndex(i1, 1)*strides_[1] +
05510 coordToIndex(i2, 2)];
05511 }
05512
05513 template<typename Numeric, unsigned Len, unsigned Dim>
05514 inline const Numeric& ArrayND<Numeric,Len,Dim>::cl(
05515 const double i0,
05516 const double i1,
05517 const double i2,
05518 const double i3) const
05519 {
05520 if (4U != dim_) throw npstat::NpstatInvalidArgument(
05521 "In npstat::ArrayND::cl: wrong # of args (not rank 4 array)");
05522 return data_[coordToIndex(i0, 0)*strides_[0] +
05523 coordToIndex(i1, 1)*strides_[1] +
05524 coordToIndex(i2, 2)*strides_[2] +
05525 coordToIndex(i3, 3)];
05526 }
05527
05528 template<typename Numeric, unsigned Len, unsigned Dim>
05529 inline const Numeric& ArrayND<Numeric,Len,Dim>::cl(
05530 const double i0,
05531 const double i1,
05532 const double i2,
05533 const double i3,
05534 const double i4) const
05535 {
05536 if (5U != dim_) throw npstat::NpstatInvalidArgument(
05537 "In npstat::ArrayND::cl: wrong # of args (not rank 5 array)");
05538 return data_[coordToIndex(i0, 0)*strides_[0] +
05539 coordToIndex(i1, 1)*strides_[1] +
05540 coordToIndex(i2, 2)*strides_[2] +
05541 coordToIndex(i3, 3)*strides_[3] +
05542 coordToIndex(i4, 4)];
05543 }
05544
05545 template<typename Numeric, unsigned Len, unsigned Dim>
05546 inline const Numeric& ArrayND<Numeric,Len,Dim>::cl(
05547 const double i0,
05548 const double i1,
05549 const double i2,
05550 const double i3,
05551 const double i4,
05552 const double i5) const
05553 {
05554 if (6U != dim_) throw npstat::NpstatInvalidArgument(
05555 "In npstat::ArrayND::cl: wrong # of args (not rank 6 array)");
05556 return data_[coordToIndex(i0, 0)*strides_[0] +
05557 coordToIndex(i1, 1)*strides_[1] +
05558 coordToIndex(i2, 2)*strides_[2] +
05559 coordToIndex(i3, 3)*strides_[3] +
05560 coordToIndex(i4, 4)*strides_[4] +
05561 coordToIndex(i5, 5)];
05562 }
05563
05564 template<typename Numeric, unsigned Len, unsigned Dim>
05565 inline const Numeric& ArrayND<Numeric,Len,Dim>::cl(
05566 const double i0,
05567 const double i1,
05568 const double i2,
05569 const double i3,
05570 const double i4,
05571 const double i5,
05572 const double i6) const
05573 {
05574 if (7U != dim_) throw npstat::NpstatInvalidArgument(
05575 "In npstat::ArrayND::cl: wrong # of args (not rank 7 array)");
05576 return data_[coordToIndex(i0, 0)*strides_[0] +
05577 coordToIndex(i1, 1)*strides_[1] +
05578 coordToIndex(i2, 2)*strides_[2] +
05579 coordToIndex(i3, 3)*strides_[3] +
05580 coordToIndex(i4, 4)*strides_[4] +
05581 coordToIndex(i5, 5)*strides_[5] +
05582 coordToIndex(i6, 6)];
05583 }
05584
05585 template<typename Numeric, unsigned Len, unsigned Dim>
05586 inline const Numeric& ArrayND<Numeric,Len,Dim>::cl(
05587 const double i0,
05588 const double i1,
05589 const double i2,
05590 const double i3,
05591 const double i4,
05592 const double i5,
05593 const double i6,
05594 const double i7) const
05595 {
05596 if (8U != dim_) throw npstat::NpstatInvalidArgument(
05597 "In npstat::ArrayND::cl: wrong # of args (not rank 8 array)");
05598 return data_[coordToIndex(i0, 0)*strides_[0] +
05599 coordToIndex(i1, 1)*strides_[1] +
05600 coordToIndex(i2, 2)*strides_[2] +
05601 coordToIndex(i3, 3)*strides_[3] +
05602 coordToIndex(i4, 4)*strides_[4] +
05603 coordToIndex(i5, 5)*strides_[5] +
05604 coordToIndex(i6, 6)*strides_[6] +
05605 coordToIndex(i7, 7)];
05606 }
05607
05608 template<typename Numeric, unsigned Len, unsigned Dim>
05609 inline const Numeric& ArrayND<Numeric,Len,Dim>::cl(
05610 const double i0,
05611 const double i1,
05612 const double i2,
05613 const double i3,
05614 const double i4,
05615 const double i5,
05616 const double i6,
05617 const double i7,
05618 const double i8) const
05619 {
05620 if (9U != dim_) throw npstat::NpstatInvalidArgument(
05621 "In npstat::ArrayND::cl: wrong # of args (not rank 9 array)");
05622 return data_[coordToIndex(i0, 0)*strides_[0] +
05623 coordToIndex(i1, 1)*strides_[1] +
05624 coordToIndex(i2, 2)*strides_[2] +
05625 coordToIndex(i3, 3)*strides_[3] +
05626 coordToIndex(i4, 4)*strides_[4] +
05627 coordToIndex(i5, 5)*strides_[5] +
05628 coordToIndex(i6, 6)*strides_[6] +
05629 coordToIndex(i7, 7)*strides_[7] +
05630 coordToIndex(i8, 8)];
05631 }
05632
05633 template<typename Numeric, unsigned Len, unsigned Dim>
05634 inline const Numeric& ArrayND<Numeric,Len,Dim>::cl(
05635 const double i0,
05636 const double i1,
05637 const double i2,
05638 const double i3,
05639 const double i4,
05640 const double i5,
05641 const double i6,
05642 const double i7,
05643 const double i8,
05644 const double i9) const
05645 {
05646 if (10U != dim_) throw npstat::NpstatInvalidArgument(
05647 "In npstat::ArrayND::cl: wrong # of args (not rank 10 array)");
05648 return data_[coordToIndex(i0, 0)*strides_[0] +
05649 coordToIndex(i1, 1)*strides_[1] +
05650 coordToIndex(i2, 2)*strides_[2] +
05651 coordToIndex(i3, 3)*strides_[3] +
05652 coordToIndex(i4, 4)*strides_[4] +
05653 coordToIndex(i5, 5)*strides_[5] +
05654 coordToIndex(i6, 6)*strides_[6] +
05655 coordToIndex(i7, 7)*strides_[7] +
05656 coordToIndex(i8, 8)*strides_[8] +
05657 coordToIndex(i9, 9)];
05658 }
05659
05660 template<typename Numeric, unsigned Len, unsigned Dim>
05661 inline Numeric& ArrayND<Numeric,Len,Dim>::cl()
05662 {
05663 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
05664 "Initialize npstat::ArrayND before calling method \"cl\"");
05665 if (dim_) throw npstat::NpstatInvalidArgument(
05666 "In npstat::ArrayND::cl: wrong # of args (not rank 0 array)");
05667 return localData_[0];
05668 }
05669
05670 template<typename Numeric, unsigned Len, unsigned Dim>
05671 inline Numeric& ArrayND<Numeric,Len,Dim>::cl(
05672 const double i0)
05673 {
05674 if (1U != dim_) throw npstat::NpstatInvalidArgument(
05675 "In npstat::ArrayND::cl: wrong # of args (not rank 1 array)");
05676 return data_[coordToIndex(i0, 0)];
05677 }
05678
05679 template<typename Numeric, unsigned Len, unsigned Dim>
05680 inline Numeric& ArrayND<Numeric,Len,Dim>::cl(
05681 const double i0,
05682 const double i1)
05683 {
05684 if (2U != dim_) throw npstat::NpstatInvalidArgument(
05685 "In npstat::ArrayND::cl: wrong # of args (not rank 2 array)");
05686 return data_[coordToIndex(i0, 0)*strides_[0] +
05687 coordToIndex(i1, 1)];
05688 }
05689
05690 template<typename Numeric, unsigned Len, unsigned Dim>
05691 inline Numeric& ArrayND<Numeric,Len,Dim>::cl(
05692 const double i0,
05693 const double i1,
05694 const double i2)
05695 {
05696 if (3U != dim_) throw npstat::NpstatInvalidArgument(
05697 "In npstat::ArrayND::cl: wrong # of args (not rank 3 array)");
05698 return data_[coordToIndex(i0, 0)*strides_[0] +
05699 coordToIndex(i1, 1)*strides_[1] +
05700 coordToIndex(i2, 2)];
05701 }
05702
05703 template<typename Numeric, unsigned Len, unsigned Dim>
05704 inline Numeric& ArrayND<Numeric,Len,Dim>::cl(
05705 const double i0,
05706 const double i1,
05707 const double i2,
05708 const double i3)
05709 {
05710 if (4U != dim_) throw npstat::NpstatInvalidArgument(
05711 "In npstat::ArrayND::cl: wrong # of args (not rank 4 array)");
05712 return data_[coordToIndex(i0, 0)*strides_[0] +
05713 coordToIndex(i1, 1)*strides_[1] +
05714 coordToIndex(i2, 2)*strides_[2] +
05715 coordToIndex(i3, 3)];
05716 }
05717
05718 template<typename Numeric, unsigned Len, unsigned Dim>
05719 inline Numeric& ArrayND<Numeric,Len,Dim>::cl(
05720 const double i0,
05721 const double i1,
05722 const double i2,
05723 const double i3,
05724 const double i4)
05725 {
05726 if (5U != dim_) throw npstat::NpstatInvalidArgument(
05727 "In npstat::ArrayND::cl: wrong # of args (not rank 5 array)");
05728 return data_[coordToIndex(i0, 0)*strides_[0] +
05729 coordToIndex(i1, 1)*strides_[1] +
05730 coordToIndex(i2, 2)*strides_[2] +
05731 coordToIndex(i3, 3)*strides_[3] +
05732 coordToIndex(i4, 4)];
05733 }
05734
05735 template<typename Numeric, unsigned Len, unsigned Dim>
05736 inline Numeric& ArrayND<Numeric,Len,Dim>::cl(
05737 const double i0,
05738 const double i1,
05739 const double i2,
05740 const double i3,
05741 const double i4,
05742 const double i5)
05743 {
05744 if (6U != dim_) throw npstat::NpstatInvalidArgument(
05745 "In npstat::ArrayND::cl: wrong # of args (not rank 6 array)");
05746 return data_[coordToIndex(i0, 0)*strides_[0] +
05747 coordToIndex(i1, 1)*strides_[1] +
05748 coordToIndex(i2, 2)*strides_[2] +
05749 coordToIndex(i3, 3)*strides_[3] +
05750 coordToIndex(i4, 4)*strides_[4] +
05751 coordToIndex(i5, 5)];
05752 }
05753
05754 template<typename Numeric, unsigned Len, unsigned Dim>
05755 inline Numeric& ArrayND<Numeric,Len,Dim>::cl(
05756 const double i0,
05757 const double i1,
05758 const double i2,
05759 const double i3,
05760 const double i4,
05761 const double i5,
05762 const double i6)
05763 {
05764 if (7U != dim_) throw npstat::NpstatInvalidArgument(
05765 "In npstat::ArrayND::cl: wrong # of args (not rank 7 array)");
05766 return data_[coordToIndex(i0, 0)*strides_[0] +
05767 coordToIndex(i1, 1)*strides_[1] +
05768 coordToIndex(i2, 2)*strides_[2] +
05769 coordToIndex(i3, 3)*strides_[3] +
05770 coordToIndex(i4, 4)*strides_[4] +
05771 coordToIndex(i5, 5)*strides_[5] +
05772 coordToIndex(i6, 6)];
05773 }
05774
05775 template<typename Numeric, unsigned Len, unsigned Dim>
05776 inline Numeric& ArrayND<Numeric,Len,Dim>::cl(
05777 const double i0,
05778 const double i1,
05779 const double i2,
05780 const double i3,
05781 const double i4,
05782 const double i5,
05783 const double i6,
05784 const double i7)
05785 {
05786 if (8U != dim_) throw npstat::NpstatInvalidArgument(
05787 "In npstat::ArrayND::cl: wrong # of args (not rank 8 array)");
05788 return data_[coordToIndex(i0, 0)*strides_[0] +
05789 coordToIndex(i1, 1)*strides_[1] +
05790 coordToIndex(i2, 2)*strides_[2] +
05791 coordToIndex(i3, 3)*strides_[3] +
05792 coordToIndex(i4, 4)*strides_[4] +
05793 coordToIndex(i5, 5)*strides_[5] +
05794 coordToIndex(i6, 6)*strides_[6] +
05795 coordToIndex(i7, 7)];
05796 }
05797
05798 template<typename Numeric, unsigned Len, unsigned Dim>
05799 inline Numeric& ArrayND<Numeric,Len,Dim>::cl(
05800 const double i0,
05801 const double i1,
05802 const double i2,
05803 const double i3,
05804 const double i4,
05805 const double i5,
05806 const double i6,
05807 const double i7,
05808 const double i8)
05809 {
05810 if (9U != dim_) throw npstat::NpstatInvalidArgument(
05811 "In npstat::ArrayND::cl: wrong # of args (not rank 9 array)");
05812 return data_[coordToIndex(i0, 0)*strides_[0] +
05813 coordToIndex(i1, 1)*strides_[1] +
05814 coordToIndex(i2, 2)*strides_[2] +
05815 coordToIndex(i3, 3)*strides_[3] +
05816 coordToIndex(i4, 4)*strides_[4] +
05817 coordToIndex(i5, 5)*strides_[5] +
05818 coordToIndex(i6, 6)*strides_[6] +
05819 coordToIndex(i7, 7)*strides_[7] +
05820 coordToIndex(i8, 8)];
05821 }
05822
05823 template<typename Numeric, unsigned Len, unsigned Dim>
05824 inline Numeric& ArrayND<Numeric,Len,Dim>::cl(
05825 const double i0,
05826 const double i1,
05827 const double i2,
05828 const double i3,
05829 const double i4,
05830 const double i5,
05831 const double i6,
05832 const double i7,
05833 const double i8,
05834 const double i9)
05835 {
05836 if (10U != dim_) throw npstat::NpstatInvalidArgument(
05837 "In npstat::ArrayND::cl: wrong # of args (not rank 10 array)");
05838 return data_[coordToIndex(i0, 0)*strides_[0] +
05839 coordToIndex(i1, 1)*strides_[1] +
05840 coordToIndex(i2, 2)*strides_[2] +
05841 coordToIndex(i3, 3)*strides_[3] +
05842 coordToIndex(i4, 4)*strides_[4] +
05843 coordToIndex(i5, 5)*strides_[5] +
05844 coordToIndex(i6, 6)*strides_[6] +
05845 coordToIndex(i7, 7)*strides_[7] +
05846 coordToIndex(i8, 8)*strides_[8] +
05847 coordToIndex(i9, 9)];
05848 }
05849
05850 template<typename Numeric, unsigned StackLen, unsigned StackDim>
05851 const char* ArrayND<Numeric,StackLen,StackDim>::classname()
05852 {
05853 static const std::string name(
05854 gs::template_class_name<Numeric>("npstat::ArrayND"));
05855 return name.c_str();
05856 }
05857
05858 template<typename Numeric, unsigned StackLen, unsigned StackDim>
05859 bool ArrayND<Numeric,StackLen,StackDim>::write(std::ostream& os) const
05860 {
05861 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
05862 "Initialize npstat::ArrayND before calling method \"write\"");
05863 gs::write_pod_vector(os, shape());
05864 return !os.fail() &&
05865 (dim_ ? gs::write_array(os, data_, len_) :
05866 gs::write_item(os, localData_[0], false));
05867 }
05868
05869 template<typename Numeric, unsigned Len, unsigned Dim>
05870 void ArrayND<Numeric,Len,Dim>::restore(
05871 const gs::ClassId& id, std::istream& in, ArrayND* array)
05872 {
05873 static const gs::ClassId current(gs::ClassId::makeId<ArrayND<Numeric,Len,Dim> >());
05874 current.ensureSameId(id);
05875
05876 ArrayShape rshape;
05877 gs::read_pod_vector(in, &rshape);
05878 if (in.fail()) throw gs::IOReadFailure(
05879 "In npstat::ArrayND::restore: input stream failure (checkpoint 0)");
05880
05881 assert(array);
05882 array->uninitialize();
05883 array->dim_ = rshape.size();
05884 array->shapeIsKnown_ = true;
05885 array->len_ = 1UL;
05886 if (array->dim_)
05887 {
05888 array->shape_ = makeBuffer(array->dim_, array->localShape_, Dim);
05889 for (unsigned i=0; i<array->dim_; ++i)
05890 {
05891 array->shape_[i] = rshape[i];
05892 assert(array->shape_[i]);
05893 array->len_ *= array->shape_[i];
05894 }
05895 array->buildStrides();
05896 array->data_ = makeBuffer(array->len_, array->localData_, Len);
05897 gs::read_array(in, array->data_, array->len_);
05898 }
05899 else
05900 gs::restore_item(in, array->localData_, false);
05901 if (in.fail()) throw gs::IOReadFailure(
05902 "In npstat::ArrayND::restore: input stream failure (checkpoint 1)");
05903 }
05904
05905 template<typename Numeric, unsigned Len, unsigned StackDim>
05906 template <typename Num2, unsigned Len2, unsigned Dim2>
05907 void ArrayND<Numeric,Len,StackDim>::exportSubrange(
05908 const unsigned* corner, const unsigned lenCorner,
05909 ArrayND<Num2, Len2, Dim2>* out) const
05910 {
05911 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
05912 "Initialize npstat::ArrayND before calling method \"exportSubrange\"");
05913 if (dim_ != lenCorner) throw npstat::NpstatInvalidArgument(
05914 "In npstat::ArrayND::exportSubrange: incompatible corner index length");
05915 assert(out);
05916 if (!out->shapeIsKnown_) throw npstat::NpstatInvalidArgument(
05917 "In npstat::ArrayND::exportSubrange: uninitialized argument array");
05918 if (out->dim_ != dim_) throw npstat::NpstatInvalidArgument(
05919 "In npstat::ArrayND::exportSubrange: incompatible argument array rank");
05920
05921 if (dim_)
05922 {
05923 assert(corner);
05924 if (dim_ > CHAR_BIT*sizeof(unsigned long))
05925 throw npstat::NpstatInvalidArgument(
05926 "In npstat::ArrayND::exportSubrange: "
05927 "array rank is too large");
05928 unsigned toBuf[CHAR_BIT*sizeof(unsigned long)];
05929 clearBuffer(toBuf, dim_);
05930 (const_cast<ArrayND*>(this))->commonSubrangeLoop(
05931 0U, 0UL, 0UL, corner, out->shape_, toBuf, *out,
05932 scast_assign_right<Numeric,Num2>());
05933 }
05934 else
05935 out->localData_[0] = static_cast<Num2>(localData_[0]);
05936 }
05937
05938 template<typename Numeric, unsigned Len, unsigned StackDim>
05939 template <typename Num2, unsigned Len2, unsigned Dim2>
05940 void ArrayND<Numeric,Len,StackDim>::importSubrange(
05941 const unsigned* corner, const unsigned lenCorner,
05942 const ArrayND<Num2, Len2, Dim2>& from)
05943 {
05944 if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument(
05945 "Initialize npstat::ArrayND before calling method \"importSubrange\"");
05946 if (dim_ != lenCorner) throw npstat::NpstatInvalidArgument(
05947 "In npstat::ArrayND::importSubrange: incompatible corner index length");
05948 if (!from.shapeIsKnown_) throw npstat::NpstatInvalidArgument(
05949 "In npstat::ArrayND::importSubrange: uninitialized argument array");
05950 if (from.dim_ != dim_) throw npstat::NpstatInvalidArgument(
05951 "In npstat::ArrayND::importSubrange: incompatible argument array rank");
05952
05953 if (dim_)
05954 {
05955 assert(corner);
05956 if (dim_ > CHAR_BIT*sizeof(unsigned long))
05957 throw npstat::NpstatInvalidArgument(
05958 "In npstat::ArrayND::importSubrange: "
05959 "array rank is too large");
05960 unsigned toBuf[CHAR_BIT*sizeof(unsigned long)];
05961 clearBuffer(toBuf, dim_);
05962 commonSubrangeLoop(0U, 0UL, 0UL, corner, from.shape_, toBuf,
05963 const_cast<ArrayND<Num2, Len2, Dim2>&>(from),
05964 scast_assign_left<Numeric,Num2>());
05965 }
05966 else
05967 localData_[0] = static_cast<Numeric>(from.localData_[0]);
05968 }
05969 }
05970
05971
05972 #endif // NPSTAT_ARRAYND_HH_
05973
1.7.3