CMS 3D CMS Logo

ArrayND.h
Go to the documentation of this file.
1 #ifndef NPSTAT_ARRAYND_HH_
2 #define NPSTAT_ARRAYND_HH_
3 
14 #include <cassert>
15 
16 #include "Alignment/Geners/interface/ClassId.hh"
17 
24 
25 namespace npstat {
47  template <typename Numeric, unsigned StackLen=1U, unsigned StackDim=10U>
48  class ArrayND
49  {
50  template <typename Num2, unsigned Len2, unsigned Dim2>
51  friend class ArrayND;
52 
53  public:
54  typedef Numeric value_type;
56 
71  ArrayND();
72 
74 
82  explicit ArrayND(const ArrayShape& shape);
83  ArrayND(const unsigned* shape, unsigned dim);
85 
87  ArrayND(const ArrayND&);
88 
97  template <typename Num2, unsigned Len2, unsigned Dim2>
99 
104  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
105  ArrayND(const ArrayND<Num2, Len2, Dim2>&, Functor f);
106 
108  template <typename Num2, unsigned Len2, unsigned Dim2>
110  const ArrayRange& fromRange);
111 
113  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
115  const ArrayRange& fromRange, Functor f);
116 
128  template <typename Num2, unsigned Len2, unsigned Dim2>
129  ArrayND(const ArrayND<Num2, Len2, Dim2>& slicedArray,
130  const unsigned *indices, unsigned nIndices);
131 
133  template <typename Num1, unsigned Len1, unsigned Dim1,
134  typename Num2, unsigned Len2, unsigned Dim2>
136  const ArrayND<Num2, Len2, Dim2>& a2);
137 
139 
143  explicit ArrayND(unsigned n0);
144  ArrayND(unsigned n0, unsigned n1);
145  ArrayND(unsigned n0, unsigned n1, unsigned n2);
146  ArrayND(unsigned n0, unsigned n1, unsigned n2, unsigned n3);
147  ArrayND(unsigned n0, unsigned n1, unsigned n2, unsigned n3,
148  unsigned n4);
149  ArrayND(unsigned n0, unsigned n1, unsigned n2, unsigned n3,
150  unsigned n4, unsigned n5);
151  ArrayND(unsigned n0, unsigned n1, unsigned n2, unsigned n3,
152  unsigned n4, unsigned n5, unsigned n6);
153  ArrayND(unsigned n0, unsigned n1, unsigned n2, unsigned n3,
154  unsigned n4, unsigned n5, unsigned n6, unsigned n7);
155  ArrayND(unsigned n0, unsigned n1, unsigned n2, unsigned n3,
156  unsigned n4, unsigned n5, unsigned n6, unsigned n7,
157  unsigned n8);
158  ArrayND(unsigned n0, unsigned n1, unsigned n2, unsigned n3,
159  unsigned n4, unsigned n5, unsigned n6, unsigned n7,
160  unsigned n8, unsigned n9);
162 
164  ~ArrayND();
165 
174  ArrayND& operator=(const ArrayND&);
175 
177  template <typename Num2, unsigned Len2, unsigned Dim2>
179 
181  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
182  ArrayND& assign(const ArrayND<Num2, Len2, Dim2>&, Functor f);
183 
190 
192 
197  Numeric& value(const unsigned *index, unsigned indexLen);
198  const Numeric& value(const unsigned *index, unsigned indexLen) const;
200 
202 
206  Numeric& valueAt(const unsigned *index, unsigned indexLen);
207  const Numeric& valueAt(const unsigned *index, unsigned indexLen) const;
209 
211 
212  Numeric& linearValue(unsigned long index);
213  const Numeric& linearValue(unsigned long index) const;
215 
217 
218  Numeric& linearValueAt(unsigned long index);
219  const Numeric& linearValueAt(unsigned long index) const;
221 
223  void convertLinearIndex(unsigned long l, unsigned* index,
224  unsigned indexLen) const;
225 
227  unsigned long linearIndex(const unsigned* idx, unsigned idxLen) const;
228 
229  // Some inspectors
231  inline unsigned long length() const {return len_;}
232 
234  inline const Numeric* data() const {return data_;}
235 
237  inline bool isShapeKnown() const {return shapeIsKnown_;}
238 
240  inline unsigned rank() const {return dim_;}
241 
243  ArrayShape shape() const;
244 
246  inline const unsigned *shapeData() const {return shape_;}
247 
249  ArrayRange fullRange() const;
250 
252  unsigned span(unsigned dim) const;
253 
255  unsigned maximumSpan() const;
256 
258  unsigned minimumSpan() const;
259 
261  inline const unsigned long* strides() const {return strides_;}
262 
264  bool isZero() const;
265 
271  bool isDensity() const;
272 
274  template <typename Num2>
275  ArrayND& setData(const Num2* data, unsigned long dataLength);
276 
278  template <unsigned Len2, unsigned Dim2>
279  bool operator==(const ArrayND<Numeric,Len2,Dim2>&) const;
280 
282  template <unsigned Len2, unsigned Dim2>
283  bool operator!=(const ArrayND<Numeric,Len2,Dim2>&) const;
284 
286  template <unsigned Len2, unsigned Dim2>
287  double maxAbsDifference(const ArrayND<Numeric,Len2,Dim2>&) const;
288 
290  ArrayND operator+() const;
291 
293  ArrayND operator-() const;
294 
296  template <unsigned Len2, unsigned Dim2>
298 
300  template <unsigned Len2, unsigned Dim2>
302 
304  template <typename Num2>
305  ArrayND operator*(const Num2& r) const;
306 
308  template <typename Num2>
309  ArrayND operator/(const Num2& r) const;
310 
312 
316  template <typename Num2>
317  ArrayND& operator*=(const Num2& r);
318 
319  template <typename Num2>
320  ArrayND& operator/=(const Num2& r);
321 
322  template <typename Num2, unsigned Len2, unsigned Dim2>
324 
325  template <typename Num2, unsigned Len2, unsigned Dim2>
328 
330  template <typename Num3, typename Num2, unsigned Len2, unsigned Dim2>
331  ArrayND& addmul(const ArrayND<Num2,Len2,Dim2>& r, const Num3& c);
332 
334  template <typename Num2, unsigned Len2, unsigned Dim2>
335  ArrayND outer(const ArrayND<Num2,Len2,Dim2>& r) const;
336 
341  ArrayND contract(unsigned pos1, unsigned pos2) const;
342 
348  template <typename Num2, unsigned Len2, unsigned Dim2>
349  ArrayND dot(const ArrayND<Num2,Len2,Dim2>& r) const;
350 
365  template <typename Num2, unsigned Len2, unsigned Dim2>
367  const unsigned* indexMap, unsigned mapLen) const;
368 
370  ArrayND transpose(unsigned pos1, unsigned pos2) const;
371 
373  ArrayND transpose() const;
374 
381  template <typename Num2>
382  Num2 sum() const;
383 
388  template <typename Num2>
389  Num2 sumsq() const;
390 
399  template <typename Num2>
400  ArrayND derivative(double scale=1.0) const;
401 
406  template <typename Num2>
407  ArrayND cdfArray(double scale=1.0) const;
408 
413  template <typename Num2>
414  Num2 cdfValue(const unsigned *index, unsigned indexLen) const;
415 
427  template <typename Num2>
428  void convertToLastDimCdf(ArrayND* sumSlice, bool useTrapezoids);
429 
431  Numeric min() const;
432 
434  Numeric min(unsigned *index, unsigned indexLen) const;
435 
437  Numeric max() const;
438 
440  Numeric max(unsigned *index, unsigned indexLen) const;
441 
443 
451  Numeric& closest(const double *x, unsigned xDim);
452  const Numeric& closest(const double *x, unsigned xDim) const;
454 
463  Numeric interpolate1(const double *x, unsigned xDim) const;
464 
471  Numeric interpolate3(const double *x, unsigned xDim) const;
472 
482  template <class Functor>
483  ArrayND& apply(Functor f);
484 
492  template <class Functor>
493  ArrayND& scanInPlace(Functor f);
494 
496  ArrayND& constFill(Numeric c);
497 
499  ArrayND& clear();
500 
507  ArrayND& linearFill(const double* coeff, unsigned coeffLen, double c);
508 
515  template <class Functor>
516  ArrayND& functorFill(Functor f);
517 
524  ArrayND& makeUnit();
525 
528 
538  unsigned makeCopulaSteps(double tolerance, unsigned maxIterations);
539 
544  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
545  void jointScan(ArrayND<Num2, Len2, Dim2>& other, Functor binaryFunct);
546 
549  template <typename Num2, unsigned Len2, unsigned Dim2>
551  {
552  jointScan(const_cast<ArrayND<Num2,Len2,Dim2>&>(r),
554  return *this;
555  }
556 
574  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
576  const unsigned* thisCorner,
577  const unsigned* range,
578  const unsigned* otherCorner,
579  unsigned arrLen,
580  Functor binaryFunct);
581 
587  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
589  const unsigned* thisCorner,
590  const unsigned* range,
591  const unsigned* otherCorner,
592  unsigned arrLen,
593  Functor binaryFunct);
594 
599  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
601  const unsigned* thisCorner,
602  const unsigned* range,
603  const unsigned* otherCorner,
604  unsigned arrLen,
605  Functor binaryFunct);
606 
611  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
613  const unsigned* thisCorner,
614  const unsigned* range,
615  const unsigned* otherCorner,
616  unsigned arrLen,
617  Functor binaryFunct);
618 
625  template <typename Num2, typename Integer>
627  const BoxND<Integer>& subrange) const;
628 
636  template <typename Num2, unsigned Len2, unsigned Dim2>
637  void exportSubrange(const unsigned* fromCorner, unsigned lenCorner,
639 
641  template <typename Num2, unsigned Len2, unsigned Dim2>
642  void importSubrange(const unsigned* fromCorner, unsigned lenCorner,
643  const ArrayND<Num2, Len2, Dim2>& from);
644 
650  template <typename Num2, unsigned Len2, unsigned Dim2>
651  bool isClose(const ArrayND<Num2,Len2,Dim2>& r, double eps) const;
652 
654  bool isCompatible(const ArrayShape& shape) const;
655 
660  template <typename Num2, unsigned Len2, unsigned Dim2>
661  bool isShapeCompatible(const ArrayND<Num2,Len2,Dim2>& r) const;
662 
672  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
674  const unsigned *fixedIndices,
675  const unsigned *fixedIndexValues,
676  unsigned nFixedIndices,
677  Functor binaryFunct);
678 
687  template <typename Num2, class Functor>
688  void jointMemSliceScan(Num2* buffer, unsigned long bufLen,
689  const unsigned *fixedIndices,
690  const unsigned *fixedIndexValues,
691  unsigned nFixedIndices,
692  Functor binaryFunct);
693 
695  ArrayShape sliceShape(const unsigned *fixedIndices,
696  unsigned nFixedIndices) const;
697 
699  template <typename Num2, unsigned Len2, unsigned Dim2>
701  const unsigned *fixedIndices,
702  const unsigned *fixedIndexValues,
703  unsigned nFixedIndices) const
704  {
705  assert(slice);
706  (const_cast<ArrayND*>(this))->jointSliceScan(
707  *slice, fixedIndices, fixedIndexValues, nFixedIndices,
709  }
710 
715  template <typename Num2>
716  inline void exportMemSlice(Num2* buffer, unsigned long bufLen,
717  const unsigned *fixedIndices,
718  const unsigned *fixedIndexValues,
719  unsigned nFixedIndices) const
720  {
721  (const_cast<ArrayND*>(this))->jointMemSliceScan(
722  buffer, bufLen, fixedIndices, fixedIndexValues,
723  nFixedIndices, scast_assign_right<Numeric,Num2>());
724  }
725 
727  template <typename Num2, unsigned Len2, unsigned Dim2>
728  inline void importSlice(const ArrayND<Num2,Len2,Dim2>& slice,
729  const unsigned *fixedIndices,
730  const unsigned *fixedIndexValues,
731  unsigned nFixedIndices)
732  {
733  jointSliceScan(const_cast<ArrayND<Num2,Len2,Dim2>&>(slice),
734  fixedIndices, fixedIndexValues, nFixedIndices,
736  }
737 
742  template <typename Num2>
743  inline void importMemSlice(const Num2* buffer, unsigned long bufLen,
744  const unsigned *fixedIndices,
745  const unsigned *fixedIndexValues,
746  unsigned nFixedIndices)
747  {
748  jointMemSliceScan(const_cast<Num2*>(buffer), bufLen,
749  fixedIndices, fixedIndexValues, nFixedIndices,
751  }
752 
759  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
761  const unsigned *fixedIndices, unsigned nFixedIndices,
762  Functor binaryFunct);
763 
768  template <typename Num2, unsigned Len2, unsigned Dim2>
770  const unsigned *fixedIndices,
771  unsigned nFixedIndices)
772  {
773  applySlice(const_cast<ArrayND<Num2,Len2,Dim2>&>(slice),
774  fixedIndices, nFixedIndices,
776  return *this;
777  }
778 
780 
787  template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
788  void project(ArrayND<Num2,Len2,Dim2>* projection,
790  const unsigned *projectedIndices,
791  unsigned nProjectedIndices) const;
792 
793  template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
794  void project(ArrayND<Num2,Len2,Dim2>* projection,
795  AbsVisitor<Numeric,Num3>& projector,
796  const unsigned *projectedIndices,
797  unsigned nProjectedIndices) const;
799 
801 
806  template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
807  void addToProjection(ArrayND<Num2,Len2,Dim2>* projection,
809  const unsigned *projectedIndices,
810  unsigned nProjectedIndices) const;
811 
812  template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
815  const unsigned *projectedIndices,
816  unsigned nProjectedIndices) const;
817 
818  template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
819  void addToProjection(ArrayND<Num2,Len2,Dim2>* projection,
820  AbsVisitor<Numeric,Num3>& projector,
821  const unsigned *projectedIndices,
822  unsigned nProjectedIndices) const;
823 
824  template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
826  AbsVisitor<Numeric,Num3>& projector,
827  const unsigned *projectedIndices,
828  unsigned nProjectedIndices) const;
830 
836  template <typename Num2, unsigned Len2, unsigned Dim2>
837  void rotate(const unsigned* shifts, unsigned lenShifts,
838  ArrayND<Num2, Len2, Dim2>* rotated) const;
839 
845  template <typename Num2, unsigned Len2, unsigned Dim2>
847 
849 
853  Numeric& operator()();
854  const Numeric& operator()() const;
855 
856  Numeric& operator()(unsigned i0);
857  const Numeric& operator()(unsigned i0) const;
858 
859  Numeric& operator()(unsigned i0, unsigned i1);
860  const Numeric& operator()(unsigned i0, unsigned i1) const;
861 
862  Numeric& operator()(unsigned i0, unsigned i1, unsigned i2);
863  const Numeric& operator()(unsigned i0, unsigned i1, unsigned i2) const;
864 
865  Numeric& operator()(unsigned i0, unsigned i1,
866  unsigned i2, unsigned i3);
867  const Numeric& operator()(unsigned i0, unsigned i1,
868  unsigned i2, unsigned i3) const;
869 
870  Numeric& operator()(unsigned i0, unsigned i1,
871  unsigned i2, unsigned i3, unsigned i4);
872  const Numeric& operator()(unsigned i0, unsigned i1,
873  unsigned i2, unsigned i3, unsigned i4) const;
874 
875  Numeric& operator()(unsigned i0, unsigned i1, unsigned i2,
876  unsigned i3, unsigned i4, unsigned i5);
877  const Numeric& operator()(unsigned i0, unsigned i1, unsigned i2,
878  unsigned i3, unsigned i4, unsigned i5) const;
879 
880  Numeric& operator()(unsigned i0, unsigned i1, unsigned i2,
881  unsigned i3, unsigned i4, unsigned i5,
882  unsigned i6);
883  const Numeric& operator()(unsigned i0, unsigned i1, unsigned i2,
884  unsigned i3, unsigned i4, unsigned i5,
885  unsigned i6) const;
886 
887  Numeric& operator()(unsigned i0, unsigned i1, unsigned i2,
888  unsigned i3, unsigned i4, unsigned i5,
889  unsigned i6, unsigned i7);
890  const Numeric& operator()(unsigned i0, unsigned i1, unsigned i2,
891  unsigned i3, unsigned i4, unsigned i5,
892  unsigned i6, unsigned i7) const;
893 
894  Numeric& operator()(unsigned i0, unsigned i1, unsigned i2,
895  unsigned i3, unsigned i4, unsigned i5,
896  unsigned i6, unsigned i7, unsigned i8);
897  const Numeric& operator()(unsigned i0, unsigned i1, unsigned i2,
898  unsigned i3, unsigned i4, unsigned i5,
899  unsigned i6, unsigned i7, unsigned i8) const;
900 
901  Numeric& operator()(unsigned i0, unsigned i1, unsigned i2,
902  unsigned i3, unsigned i4, unsigned i5,
903  unsigned i6, unsigned i7, unsigned i8,
904  unsigned i9);
905  const Numeric& operator()(unsigned i0, unsigned i1, unsigned i2,
906  unsigned i3, unsigned i4, unsigned i5,
907  unsigned i6, unsigned i7, unsigned i8,
908  unsigned i9) const;
910 
912 
916  Numeric& at();
917  const Numeric& at() const;
918 
919  Numeric& at(unsigned i0);
920  const Numeric& at(unsigned i0) const;
921 
922  Numeric& at(unsigned i0, unsigned i1);
923  const Numeric& at(unsigned i0, unsigned i1) const;
924 
925  Numeric& at(unsigned i0, unsigned i1, unsigned i2);
926  const Numeric& at(unsigned i0, unsigned i1, unsigned i2) const;
927 
928  Numeric& at(unsigned i0, unsigned i1,
929  unsigned i2, unsigned i3);
930  const Numeric& at(unsigned i0, unsigned i1,
931  unsigned i2, unsigned i3) const;
932 
933  Numeric& at(unsigned i0, unsigned i1,
934  unsigned i2, unsigned i3, unsigned i4);
935  const Numeric& at(unsigned i0, unsigned i1,
936  unsigned i2, unsigned i3, unsigned i4) const;
937 
938  Numeric& at(unsigned i0, unsigned i1, unsigned i2,
939  unsigned i3, unsigned i4, unsigned i5);
940  const Numeric& at(unsigned i0, unsigned i1, unsigned i2,
941  unsigned i3, unsigned i4, unsigned i5) const;
942 
943  Numeric& at(unsigned i0, unsigned i1, unsigned i2,
944  unsigned i3, unsigned i4, unsigned i5,
945  unsigned i6);
946  const Numeric& at(unsigned i0, unsigned i1, unsigned i2,
947  unsigned i3, unsigned i4, unsigned i5,
948  unsigned i6) const;
949 
950  Numeric& at(unsigned i0, unsigned i1, unsigned i2,
951  unsigned i3, unsigned i4, unsigned i5,
952  unsigned i6, unsigned i7);
953  const Numeric& at(unsigned i0, unsigned i1, unsigned i2,
954  unsigned i3, unsigned i4, unsigned i5,
955  unsigned i6, unsigned i7) const;
956 
957  Numeric& at(unsigned i0, unsigned i1, unsigned i2,
958  unsigned i3, unsigned i4, unsigned i5,
959  unsigned i6, unsigned i7, unsigned i8);
960  const Numeric& at(unsigned i0, unsigned i1, unsigned i2,
961  unsigned i3, unsigned i4, unsigned i5,
962  unsigned i6, unsigned i7, unsigned i8) const;
963 
964  Numeric& at(unsigned i0, unsigned i1, unsigned i2,
965  unsigned i3, unsigned i4, unsigned i5,
966  unsigned i6, unsigned i7, unsigned i8,
967  unsigned i9);
968  const Numeric& at(unsigned i0, unsigned i1, unsigned i2,
969  unsigned i3, unsigned i4, unsigned i5,
970  unsigned i6, unsigned i7, unsigned i8,
971  unsigned i9) const;
973 
975 
979  Numeric& cl();
980  const Numeric& cl() const;
981 
982  Numeric& cl(double x0);
983  const Numeric& cl(double x0) const;
984 
985  Numeric& cl(double x0, double x1);
986  const Numeric& cl(double x0, double x1) const;
987 
988  Numeric& cl(double x0, double x1, double x2);
989  const Numeric& cl(double x0, double x1, double x2) const;
990 
991  Numeric& cl(double x0, double x1,
992  double x2, double x3);
993  const Numeric& cl(double x0, double x1,
994  double x2, double x3) const;
995 
996  Numeric& cl(double x0, double x1,
997  double x2, double x3, double x4);
998  const Numeric& cl(double x0, double x1,
999  double x2, double x3, double x4) const;
1000 
1001  Numeric& cl(double x0, double x1, double x2,
1002  double x3, double x4, double x5);
1003  const Numeric& cl(double x0, double x1, double x2,
1004  double x3, double x4, double x5) const;
1005 
1006  Numeric& cl(double x0, double x1, double x2,
1007  double x3, double x4, double x5,
1008  double x6);
1009  const Numeric& cl(double x0, double x1, double x2,
1010  double x3, double x4, double x5,
1011  double x6) const;
1012 
1013  Numeric& cl(double x0, double x1, double x2,
1014  double x3, double x4, double x5,
1015  double x6, double x7);
1016  const Numeric& cl(double x0, double x1, double x2,
1017  double x3, double x4, double x5,
1018  double x6, double x7) const;
1019 
1020  Numeric& cl(double x0, double x1, double x2,
1021  double x3, double x4, double x5,
1022  double x6, double x7, double x8);
1023  const Numeric& cl(double x0, double x1, double x2,
1024  double x3, double x4, double x5,
1025  double x6, double x7, double x8) const;
1026 
1027  Numeric& cl(double x0, double x1, double x2,
1028  double x3, double x4, double x5,
1029  double x6, double x7, double x8,
1030  double x9);
1031  const Numeric& cl(double x0, double x1, double x2,
1032  double x3, double x4, double x5,
1033  double x6, double x7, double x8,
1034  double x9) const;
1036 
1038 
1039  inline gs::ClassId classId() const {return gs::ClassId(*this);}
1040  bool write(std::ostream& of) const;
1042 
1043  static const char* classname();
1044  static inline unsigned version() {return 1;}
1045  static void restore(const gs::ClassId& id, std::istream& in,
1046  ArrayND* array);
1047  private:
1048  Numeric localData_[StackLen];
1049  Numeric* data_;
1050 
1051  unsigned long localStrides_[StackDim];
1052  unsigned long *strides_;
1053 
1054  unsigned localShape_[StackDim];
1055  unsigned *shape_;
1056 
1057  unsigned long len_;
1058  unsigned dim_;
1059 
1061 
1062  // Basic initialization from unsigned* shape and dimensionality
1063  void buildFromShapePtr(const unsigned*, unsigned);
1064 
1065  // Build strides_ array out of the shape_ array
1066  void buildStrides();
1067 
1068  // Recursive implementation of nested loops for "linearFill"
1069  void linearFillLoop(unsigned level, double s0,
1070  unsigned long idx, double shift,
1071  const double* coeffs);
1072 
1073  // Recursive implementation of nested loops for "functorFill"
1074  template <class Functor>
1075  void functorFillLoop(unsigned level, unsigned long idx,
1076  Functor f, unsigned* farg);
1077 
1078  // Recursive implementation of nested loops for "interpolate3"
1079  Numeric interpolateLoop(unsigned level, const double *x,
1080  const Numeric* base) const;
1081 
1082  // Recursive implementation of nested loops for the outer product
1083  template <typename Num1, unsigned Len1, unsigned Dim1,
1084  typename Num2, unsigned Len2, unsigned Dim2>
1085  void outerProductLoop(unsigned level, unsigned long idx0,
1086  unsigned long idx1, unsigned long idx2,
1087  const ArrayND<Num1, Len1, Dim1>& a1,
1088  const ArrayND<Num2, Len2, Dim2>& a2);
1089 
1090  // Recursive implementation of nested loops for contraction
1091  void contractLoop(unsigned thisLevel, unsigned resLevel,
1092  unsigned pos1, unsigned pos2,
1093  unsigned long idxThis, unsigned long idxRes,
1094  ArrayND& result) const;
1095 
1096  // Recursive implementation of nested loops for transposition
1097  void transposeLoop(unsigned level, unsigned pos1, unsigned pos2,
1098  unsigned long idxThis, unsigned long idxRes,
1099  ArrayND& result) const;
1100 
1101  // Recursive implementation of nested loops for the dot product
1102  template <typename Num2, unsigned Len2, unsigned Dim2>
1103  void dotProductLoop(unsigned level, unsigned long idx0,
1104  unsigned long idx1, unsigned long idx2,
1105  const ArrayND<Num2, Len2, Dim2>& r,
1106  ArrayND& result) const;
1107 
1108  // Recursive implementation of nested loops for marginalization
1109  template <typename Num2, unsigned Len2, unsigned Dim2>
1110  Numeric marginalizeInnerLoop(unsigned long idx,
1111  unsigned levelPr, unsigned long idxPr,
1112  const ArrayND<Num2,Len2,Dim2>& prior,
1113  const unsigned* indexMap) const;
1114  template <typename Num2, unsigned Len2, unsigned Dim2>
1115  void marginalizeLoop(unsigned level, unsigned long idx,
1116  unsigned levelRes, unsigned long idxRes,
1117  const ArrayND<Num2,Len2,Dim2>& prior,
1118  const unsigned* indexMap, ArrayND& res) const;
1119 
1120  // Recursive implementation of nested loops for range copy
1121  // with functor modification of elements
1122  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
1123  void copyRangeLoopFunct(unsigned level, unsigned long idx0,
1124  unsigned long idx1,
1125  const ArrayND<Num2, Len2, Dim2>& r,
1126  const ArrayRange& range, Functor f);
1127 
1128  // Loop over subrange in such a way that the functor is called
1129  // only if indices on both sides are valid. The topology of both
1130  // arrays is that of the hyperplane (flat).
1131  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
1132  void commonSubrangeLoop(unsigned level, unsigned long idx0,
1133  unsigned long idx1,
1134  const unsigned* thisCorner,
1135  const unsigned* range,
1136  const unsigned* otherCorner,
1138  Functor binaryFunct);
1139 
1140  // Similar loop with the topology of the hypertorus for both
1141  // arrays (all indices of both arrays are wrapped around)
1142  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
1143  void dualCircularLoop(unsigned level, unsigned long idx0,
1144  unsigned long idx1,
1145  const unsigned* thisCorner,
1146  const unsigned* range,
1147  const unsigned* otherCorner,
1149  Functor binaryFunct);
1150 
1151  // Similar loop in which the topology of this array is assumed
1152  // to be flat and the topology of the destination array is that
1153  // of the hypertorus. Due to the asymmetry of the topologies,
1154  // "const" function is not provided (use const_cast as appropriate).
1155  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
1156  void flatCircularLoop(unsigned level, unsigned long idx0,
1157  unsigned long idx1,
1158  const unsigned* thisCorner,
1159  const unsigned* range,
1160  const unsigned* otherCorner,
1162  Functor binaryFunct);
1163 
1164  // Similar loop in which the topology of this array is assumed
1165  // to be hypertoroidal and the topology of the destination array
1166  // is flat.
1167  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
1168  void circularFlatLoop(unsigned level, unsigned long idx0,
1169  unsigned long idx1,
1170  const unsigned* thisCorner,
1171  const unsigned* range,
1172  const unsigned* otherCorner,
1174  Functor binaryFunct);
1175 
1176  // Slice compatibility verification
1177  template <typename Num2, unsigned Len2, unsigned Dim2>
1178  unsigned long verifySliceCompatibility(
1179  const ArrayND<Num2,Len2,Dim2>& slice,
1180  const unsigned *fixedIndices,
1181  const unsigned *fixedIndexValues,
1182  unsigned nFixedIndices) const;
1183 
1184  // Buffer compatibility verification with a slice
1185  unsigned long verifyBufferSliceCompatibility(
1186  unsigned long bufLen,
1187  const unsigned *fixedIndices,
1188  const unsigned *fixedIndexValues,
1189  unsigned nFixedIndices,
1190  unsigned long* sliceStrides) const;
1191 
1192  // Recursive implementation of nested loops for slice operations
1193  template <typename Num2, class Functor>
1194  void jointSliceLoop(unsigned level, unsigned long idx0,
1195  unsigned level1, unsigned long idx1,
1196  Num2* sliceData, const unsigned long* sliceStrides,
1197  const unsigned *fixedIndices,
1198  const unsigned *fixedIndexValues,
1199  unsigned nFixedIndices, Functor binaryFunctor);
1200 
1201  // Recursive implementation of nested loops for "applySlice"
1202  template <typename Num2, class Functor>
1203  void scaleBySliceInnerLoop(unsigned level, unsigned long idx0,
1204  Num2& scale,
1205  const unsigned* projectedIndices,
1206  unsigned nProjectedIndices,
1207  Functor binaryFunct);
1208 
1209  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
1210  void scaleBySliceLoop(unsigned level, unsigned long idx0,
1211  unsigned level1, unsigned long idx1,
1212  ArrayND<Num2,Len2,Dim2>& slice,
1213  const unsigned *fixedIndices,
1214  unsigned nFixedIndices,
1215  Functor binaryFunct);
1216 
1217  // Recursive implementation of nested loops for projections
1218  template <typename Num2>
1219  void projectInnerLoop(unsigned level, unsigned long idx0,
1220  unsigned* currentIndex,
1222  const unsigned* projectedIndices,
1223  unsigned nProjectedIndices) const;
1224 
1225  template <typename Num2, unsigned Len2, unsigned Dim2,
1226  typename Num3, class Op>
1227  void projectLoop(unsigned level, unsigned long idx0,
1228  unsigned level1, unsigned long idx1,
1229  unsigned* currentIndex,
1230  ArrayND<Num2,Len2,Dim2>* projection,
1232  const unsigned* projectedIndices,
1233  unsigned nProjectedIndices, Op fcn) const;
1234 
1235  // Note that "projectLoop2" is almost identical to "projectLoop"
1236  // while "projectInnerLoop2" is almost identical to "projectInnerLoop".
1237  // It would make a lot of sense to combine these functions into
1238  // the same code and then partially specialize the little piece
1239  // where the "AbsVisitor" or "AbsArrayProjector" is actually called.
1240  // Unfortunately, "AbsVisitor" and "AbsArrayProjector" are
1241  // templates themselves, and it is not possible in C++ to partially
1242  // specialize a function template (that is, even if we can specialize
1243  // on "AbsVisitor" vs. "AbsArrayProjector", there is no way to
1244  // specialize on their parameter types).
1245  template <typename Num2, unsigned Len2, unsigned Dim2,
1246  typename Num3, class Op>
1247  void projectLoop2(unsigned level, unsigned long idx0,
1248  unsigned level1, unsigned long idx1,
1249  ArrayND<Num2,Len2,Dim2>* projection,
1250  AbsVisitor<Numeric,Num3>& projector,
1251  const unsigned* projectedIndices,
1252  unsigned nProjectedIndices, Op fcn) const;
1253 
1254  template <typename Num2>
1255  void projectInnerLoop2(unsigned level, unsigned long idx0,
1256  AbsVisitor<Numeric,Num2>& projector,
1257  const unsigned* projectedIndices,
1258  unsigned nProjectedIndices) const;
1259 
1260  template <typename Num2, typename Integer>
1261  void processSubrangeLoop(unsigned level, unsigned long idx0,
1262  unsigned* currentIndex,
1264  const BoxND<Integer>& subrange) const;
1265 
1266  // Sum of all points with the given index and below
1267  template <typename Accumulator>
1268  Accumulator sumBelowLoop(unsigned level, unsigned long idx0,
1269  const unsigned* limit) const;
1270 
1271  // Loop for "convertToLastDimCdf"
1272  template <typename Accumulator>
1273  void convertToLastDimCdfLoop(ArrayND* sumSlice, unsigned level,
1274  unsigned long idx0,
1275  unsigned long idxSlice,
1276  bool useTrapezoids);
1277 
1278  // Convert a coordinate into index.
1279  // No checking whether "idim" is within limits.
1280  unsigned coordToIndex(double coord, unsigned idim) const;
1281 
1282  // Verify that projection array is compatible with this one
1283  template <typename Num2, unsigned Len2, unsigned Dim2>
1285  const ArrayND<Num2,Len2,Dim2>& projection,
1286  const unsigned *projectedIndices,
1287  unsigned nProjectedIndices) const;
1288 
1289  };
1290 }
1291 
1292 #include <cmath>
1293 #include <climits>
1294 #include <algorithm>
1295 #include <sstream>
1297 
1298 #include "Alignment/Geners/interface/GenericIO.hh"
1299 #include "Alignment/Geners/interface/IOIsUnsigned.hh"
1300 
1302 
1307 
1308 #define me_macro_check_loop_prerequisites(METHOD, INNERLOOP) \
1309  template<typename Numeric, unsigned Len, unsigned Dim> \
1310  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor> \
1311  void ArrayND<Numeric,Len,Dim>:: METHOD ( \
1312  ArrayND<Num2, Len2, Dim2>& other, \
1313  const unsigned* thisCorner, \
1314  const unsigned* range, \
1315  const unsigned* otherCorner, \
1316  const unsigned arrLen, \
1317  Functor binaryFunct) \
1318  { \
1319  if (!shapeIsKnown_) throw npstat::NpstatInvalidArgument( \
1320  "Initialize npstat::ArrayND before calling method \"" \
1321  #METHOD "\""); \
1322  if (!other.shapeIsKnown_) throw npstat::NpstatInvalidArgument( \
1323  "In npstat::ArrayND::" #METHOD ": uninitialized argument array");\
1324  if (dim_ != other.dim_) throw npstat::NpstatInvalidArgument( \
1325  "In npstat::ArrayND::" #METHOD ": incompatible argument array rank");\
1326  if (arrLen != dim_) throw npstat::NpstatInvalidArgument( \
1327  "In npstat::ArrayND::" #METHOD ": incompatible index length"); \
1328  if (dim_) \
1329  { \
1330  assert(thisCorner); \
1331  assert(range); \
1332  assert(otherCorner); \
1333  INNERLOOP (0U, 0UL, 0UL, thisCorner, range, \
1334  otherCorner, other, binaryFunct); \
1335  } \
1336  else \
1337  binaryFunct(localData_[0], other.localData_[0]); \
1338  }
1339 
1340 namespace npstat {
1341  template<typename Numeric, unsigned Len, unsigned Dim>
1342  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
1344  unsigned level, unsigned long idx0,
1345  unsigned long idx1,
1346  const unsigned* thisCorner,
1347  const unsigned* range,
1348  const unsigned* otherCorner,
1350  Functor binaryFunct)
1351  {
1352  const unsigned imax = range[level];
1353 
1354  if (level == dim_ - 1)
1355  {
1356  Numeric* left = data_ + (idx0 + thisCorner[level]);
1357  Numeric* const lMax = data_ + (idx0 + shape_[level]);
1358  Num2* right = r.data_ + (idx1 + otherCorner[level]);
1359  Num2* const rMax = r.data_ + (idx1 + r.shape_[level]);
1360 
1361  for (unsigned i=0; i<imax && left<lMax && right<rMax; ++i)
1362  binaryFunct(*left++, *right++);
1363  }
1364  else
1365  {
1366  const unsigned long leftStride = strides_[level];
1367  const unsigned long leftMax = idx0 + shape_[level]*leftStride;
1368  idx0 += thisCorner[level]*leftStride;
1369  const unsigned long rightStride = r.strides_[level];
1370  const unsigned long rightMax = idx1 + r.shape_[level]*rightStride;
1371  idx1 += otherCorner[level]*rightStride;
1372 
1373  for (unsigned i=0; i<imax && idx0 < leftMax && idx1 < rightMax;
1374  ++i, idx0 += leftStride, idx1 += rightStride)
1375  commonSubrangeLoop(level+1, idx0, idx1, thisCorner, range,
1376  otherCorner, r, binaryFunct);
1377  }
1378  }
1379 
1380  template<typename Numeric, unsigned Len, unsigned Dim>
1381  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
1383  unsigned level, unsigned long idx0,
1384  unsigned long idx1,
1385  const unsigned* thisCorner,
1386  const unsigned* range,
1387  const unsigned* otherCorner,
1389  Functor binaryFunct)
1390  {
1391  const unsigned imax = range[level];
1392  const unsigned leftShift = thisCorner[level];
1393  const unsigned leftPeriod = shape_[level];
1394  const unsigned rightShift = otherCorner[level];
1395  const unsigned rightPeriod = r.shape_[level];
1396 
1397  if (level == dim_ - 1)
1398  {
1399  Numeric* left = data_ + idx0;
1400  Num2* right = r.data_ + idx1;
1401  for (unsigned i=0; i<imax; ++i)
1402  binaryFunct(left[(i+leftShift)%leftPeriod],
1403  right[(i+rightShift)%rightPeriod]);
1404  }
1405  else
1406  {
1407  const unsigned long leftStride = strides_[level];
1408  const unsigned long rightStride = r.strides_[level];
1409  for (unsigned i=0; i<imax; ++i)
1411  level+1, idx0+((i+leftShift)%leftPeriod)*leftStride,
1412  idx1+((i+rightShift)%rightPeriod)*rightStride,
1413  thisCorner, range, otherCorner, r, binaryFunct);
1414  }
1415  }
1416 
1417  template<typename Numeric, unsigned Len, unsigned Dim>
1418  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
1420  unsigned level, unsigned long idx0,
1421  unsigned long idx1,
1422  const unsigned* thisCorner,
1423  const unsigned* range,
1424  const unsigned* otherCorner,
1426  Functor binaryFunct)
1427  {
1428  const unsigned imax = range[level];
1429  const unsigned rightShift = otherCorner[level];
1430  const unsigned rightPeriod = r.shape_[level];
1431 
1432  if (level == dim_ - 1)
1433  {
1434  Numeric* left = data_ + (idx0 + thisCorner[level]);
1435  Numeric* const lMax = data_ + (idx0 + shape_[level]);
1436  Num2* right = r.data_ + idx1;
1437 
1438  for (unsigned i=0; i<imax && left<lMax; ++i)
1439  binaryFunct(*left++, right[(i+rightShift)%rightPeriod]);
1440  }
1441  else
1442  {
1443  const unsigned long leftStride = strides_[level];
1444  const unsigned long leftMax = idx0 + shape_[level]*leftStride;
1445  idx0 += thisCorner[level]*leftStride;
1446  const unsigned long rightStride = r.strides_[level];
1447 
1448  for (unsigned i=0; i<imax && idx0 < leftMax; ++i, idx0+=leftStride)
1450  level+1, idx0,
1451  idx1+((i+rightShift)%rightPeriod)*rightStride,
1452  thisCorner, range, otherCorner, r, binaryFunct);
1453  }
1454  }
1455 
1456  template<typename Numeric, unsigned Len, unsigned Dim>
1457  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
1459  unsigned level, unsigned long idx0,
1460  unsigned long idx1,
1461  const unsigned* thisCorner,
1462  const unsigned* range,
1463  const unsigned* otherCorner,
1465  Functor binaryFunct)
1466  {
1467  const unsigned imax = range[level];
1468  const unsigned leftShift = thisCorner[level];
1469  const unsigned leftPeriod = shape_[level];
1470 
1471  if (level == dim_ - 1)
1472  {
1473  Numeric* left = data_ + idx0;
1474  Num2* right = r.data_ + (idx1 + otherCorner[level]);
1475  Num2* const rMax = r.data_ + (idx1 + r.shape_[level]);
1476 
1477  for (unsigned i=0; i<imax && right<rMax; ++i)
1478  binaryFunct(left[(i+leftShift)%leftPeriod], *right++);
1479  }
1480  else
1481  {
1482  const unsigned long leftStride = strides_[level];
1483  const unsigned long rightStride = r.strides_[level];
1484  const unsigned long rightMax = idx1 + r.shape_[level]*rightStride;
1485  idx1 += otherCorner[level]*rightStride;
1486 
1487  for (unsigned i=0; i<imax && idx1<rightMax; ++i, idx1+=rightStride)
1489  level+1, idx0+((i+leftShift)%leftPeriod)*leftStride,
1490  idx1, thisCorner, range, otherCorner, r, binaryFunct);
1491  }
1492  }
1493 
1498 
1499  template <typename Numeric, unsigned StackLen, unsigned StackDim>
1500  template <typename Num2, unsigned Len2, unsigned Dim2>
1501  Numeric ArrayND<Numeric,StackLen,StackDim>::marginalizeInnerLoop(
1502  unsigned long idx, const unsigned levelPr, unsigned long idxPr,
1503  const ArrayND<Num2,Len2,Dim2>& prior,
1504  const unsigned* indexMap) const
1505  {
1506  Numeric sum = Numeric();
1507  const unsigned long myStride = strides_[indexMap[levelPr]];
1508  const unsigned imax = prior.shape_[levelPr];
1509  assert(imax == shape_[indexMap[levelPr]]);
1510  if (levelPr == prior.dim_ - 1)
1511  {
1512  for (unsigned i=0; i<imax; ++i)
1513  sum += data_[idx+i*myStride]*prior.data_[idxPr++];
1514  }
1515  else
1516  {
1517  const unsigned long priorStride = prior.strides_[levelPr];
1518  for (unsigned i=0; i<imax; ++i)
1519  {
1520  sum += marginalizeInnerLoop(idx, levelPr+1U, idxPr,
1521  prior, indexMap);
1522  idx += myStride;
1523  idxPr += priorStride;
1524  }
1525  }
1526  return sum;
1527  }
1528 
1529  template <typename Numeric, unsigned StackLen, unsigned StackDim>
1530  template <typename Num2, unsigned Len2, unsigned Dim2>
1532  const unsigned level, unsigned long idx,
1533  const unsigned levelRes, unsigned long idxRes,
1535  const unsigned* indexMap, ArrayND& result) const
1536  {
1537  if (level == dim_)
1538  {
1539  const Numeric res = marginalizeInnerLoop(
1540  idx, 0U, 0UL, prior, indexMap);
1541  if (result.dim_)
1542  result.data_[idxRes] = res;
1543  else
1544  result.localData_[0] = res;
1545  }
1546  else
1547  {
1548  // Check if this level is mapped or not
1549  bool mapped = false;
1550  for (unsigned i=0; i<prior.dim_; ++i)
1551  if (level == indexMap[i])
1552  {
1553  mapped = true;
1554  break;
1555  }
1556  if (mapped)
1557  marginalizeLoop(level+1U, idx, levelRes, idxRes,
1558  prior, indexMap, result);
1559  else
1560  {
1561  const unsigned imax = shape_[level];
1562  const unsigned long myStride = strides_[level];
1563  const unsigned long resStride = result.strides_[levelRes];
1564  for (unsigned i=0; i<imax; ++i)
1565  {
1566  marginalizeLoop(level+1U, idx, levelRes+1U, idxRes,
1567  prior, indexMap, result);
1568  idx += myStride;
1569  idxRes += resStride;
1570  }
1571  }
1572  }
1573  }
1574 
1575  template<typename Numeric, unsigned StackLen, unsigned StackDim>
1576  template <typename Num2, unsigned Len2, unsigned Dim2>
1580  const unsigned* indexMap, const unsigned mapLen) const
1581  {
1583  "Initialize npstat::ArrayND before calling method \"marginalize\"");
1584  if (!(prior.dim_ && prior.dim_ <= dim_)) throw npstat::NpstatInvalidArgument(
1585  "In npstat::ArrayND::marginalize: incompatible argument array rank");
1586  const unsigned resultDim = dim_ - prior.dim_;
1587 
1588  // Check that the index map is reasonable
1589  if (mapLen != prior.dim_) throw npstat::NpstatInvalidArgument(
1590  "In npstat::ArrayND::marginalize: incompatible index map length");
1591  assert(indexMap);
1592  for (unsigned i=0; i<mapLen; ++i)
1593  {
1594  const unsigned thisInd = indexMap[i];
1595  if (shape_[thisInd] != prior.shape_[i]) throw npstat::NpstatInvalidArgument(
1596  "In npstat::ArrayND::marginalize: "
1597  "incompatible argument array dimensions");
1598  if (thisInd >= dim_) throw npstat::NpstatOutOfRange(
1599  "In npstat::ArrayND::marginalize: index map entry out of range");
1600  for (unsigned j=0; j<i; ++j)
1601  if (indexMap[j] == thisInd) throw npstat::NpstatInvalidArgument(
1602  "In npstat::ArrayND::marginalize: "
1603  "duplicate entry in the index map");
1604  }
1605 
1606  // Build the shape for the array of results
1607  ArrayShape newShape;
1608  newShape.reserve(resultDim);
1609  for (unsigned i=0; i<dim_; ++i)
1610  {
1611  bool mapped = false;
1612  for (unsigned j=0; j<mapLen; ++j)
1613  if (indexMap[j] == i)
1614  {
1615  mapped = true;
1616  break;
1617  }
1618  if (!mapped)
1619  newShape.push_back(shape_[i]);
1620  }
1621 
1622  ArrayND result(newShape);
1623  assert(result.dim_ == resultDim);
1624  bool calculated = false;
1625  if (resultDim == 0)
1626  {
1627  calculated = true;
1628  for (unsigned i=0; i<dim_; ++i)
1629  if (indexMap[i] != i)
1630  {
1631  calculated = false;
1632  break;
1633  }
1634  if (calculated)
1635  {
1636  Numeric sum = Numeric();
1637  for (unsigned long i=0; i<len_; ++i)
1638  sum += data_[i]*prior.data_[i];
1639  result.localData_[0] = sum;
1640  }
1641  }
1642 
1643  if (!calculated)
1644  marginalizeLoop(0U, 0UL, 0U, 0UL, prior, indexMap, result);
1645 
1646  return result;
1647  }
1648 
1649  template<typename Numeric, unsigned Len, unsigned Dim>
1651  const unsigned* sizes, const unsigned dim)
1652  {
1653  dim_ = dim;
1654  if (dim_)
1655  {
1656  assert(sizes);
1657  for (unsigned i=0; i<dim_; ++i)
1658  if (sizes[i] == 0)
1660  "In npstat::ArrayND::buildFromShapePtr: "
1661  "detected span of zero");
1662 
1663  // Copy the array shape and figure out the array length
1664  shape_ = makeBuffer(dim_, localShape_, Dim);
1665  for (unsigned i=0; i<dim_; ++i)
1666  {
1667  shape_[i] = sizes[i];
1668  len_ *= shape_[i];
1669  }
1670 
1671  // Figure out the array strides
1672  buildStrides();
1673 
1674  // Allocate the data array
1675  data_ = makeBuffer(len_, localData_, Len);
1676  }
1677  else
1678  {
1679  localData_[0] = Numeric();
1680  data_ = localData_;
1681  }
1682  }
1683 
1684  template<typename Numeric, unsigned StackLen, unsigned StackDim>
1685  template <typename Num2, unsigned Len2, unsigned Dim2>
1687  const ArrayND<Num2, Len2, Dim2>& slicedArray,
1688  const unsigned *fixedIndices, const unsigned nFixedIndices)
1689  : data_(0), strides_(nullptr), shape_(nullptr),
1690  len_(1UL), dim_(slicedArray.dim_ - nFixedIndices),
1692  {
1693  if (nFixedIndices)
1694  {
1695  assert(fixedIndices);
1696  if (nFixedIndices > slicedArray.dim_) throw npstat::NpstatInvalidArgument(
1697  "In npstat::ArrayND slicing constructor: too many fixed indices");
1698  if (!slicedArray.shapeIsKnown_) throw npstat::NpstatInvalidArgument(
1699  "In npstat::ArrayND slicing constructor: "
1700  "uninitialized argument array");
1701 
1702  // Check that the fixed indices are within range
1703  for (unsigned j=0; j<nFixedIndices; ++j)
1704  if (fixedIndices[j] >= slicedArray.dim_)
1705  throw npstat::NpstatOutOfRange("In npstat::ArrayND slicing "
1706  "constructor: fixed index out of range");
1707 
1708  // Build the shape for the slice
1709  shape_ = makeBuffer(dim_, localShape_, StackDim);
1710  unsigned idim = 0;
1711  for (unsigned i=0; i<slicedArray.dim_; ++i)
1712  {
1713  bool fixed = false;
1714  for (unsigned j=0; j<nFixedIndices; ++j)
1715  if (fixedIndices[j] == i)
1716  {
1717  fixed = true;
1718  break;
1719  }
1720  if (!fixed)
1721  {
1722  assert(idim < dim_);
1723  shape_[idim++] = slicedArray.shape_[i];
1724  }
1725  }
1726  assert(idim == dim_);
1727 
1728  if (dim_)
1729  {
1730  // Copy the array shape and figure out the array length
1731  for (unsigned i=0; i<dim_; ++i)
1732  len_ *= shape_[i];
1733 
1734  // Figure out the array strides
1735  buildStrides();
1736 
1737  // Allocate the data array
1738  data_ = makeBuffer(len_, localData_, StackLen);
1739  }
1740  else
1741  {
1742  localData_[0] = Numeric();
1743  data_ = localData_;
1744  }
1745  }
1746  else
1747  {
1748  new (this) ArrayND(slicedArray);
1749  }
1750  }
1751 
1752  template<typename Numeric, unsigned StackLen, unsigned StackDim>
1754  const unsigned *fixedIndices, const unsigned nFixedIndices) const
1755  {
1757  "Initialize npstat::ArrayND before calling method \"sliceShape\"");
1758  if (nFixedIndices)
1759  {
1760  assert(fixedIndices);
1761  if (nFixedIndices > dim_)
1763  "In npstat::ArrayND::sliceShape: too many fixed indices");
1764  for (unsigned j=0; j<nFixedIndices; ++j)
1765  if (fixedIndices[j] >= dim_)
1766  throw npstat::NpstatOutOfRange("In npstat::ArrayND::sliceShape: "
1767  "fixed index out of range");
1768  ArrayShape sh;
1769  sh.reserve(dim_ - nFixedIndices);
1770  for (unsigned i=0; i<dim_; ++i)
1771  {
1772  bool fixed = false;
1773  for (unsigned j=0; j<nFixedIndices; ++j)
1774  if (fixedIndices[j] == i)
1775  {
1776  fixed = true;
1777  break;
1778  }
1779  if (!fixed)
1780  sh.push_back(shape_[i]);
1781  }
1782  return sh;
1783  }
1784  else
1785  return ArrayShape(shape_, shape_+dim_);
1786  }
1787 
1788  template<typename Numeric, unsigned StackLen, unsigned StackDim>
1789  template <typename Num2, unsigned Len2, unsigned Dim2>
1791  const ArrayND<Num2,Len2,Dim2>& slice,
1792  const unsigned *fixedIndices,
1793  const unsigned *fixedIndexValues,
1794  const unsigned nFixedIndices) const
1795  {
1796  if (!(nFixedIndices && nFixedIndices <= dim_))
1798  "In npstat::ArrayND::verifySliceCompatibility: "
1799  "invalid number of fixed indices");
1801  "Initialize npstat::ArrayND before calling "
1802  "method \"verifySliceCompatibility\"");
1803  if (!slice.shapeIsKnown_) throw npstat::NpstatInvalidArgument(
1804  "In npstat::ArrayND::verifySliceCompatibility: "
1805  "uninitialized argument array");
1806  if (slice.dim_ != dim_ - nFixedIndices) throw npstat::NpstatInvalidArgument(
1807  "In npstat::ArrayND::verifySliceCompatibility: "
1808  "incompatible argument array rank");
1809  assert(fixedIndices);
1810  assert(fixedIndexValues);
1811 
1812  for (unsigned j=0; j<nFixedIndices; ++j)
1813  if (fixedIndices[j] >= dim_) throw npstat::NpstatOutOfRange(
1814  "In npstat::ArrayND::verifySliceCompatibility: "
1815  "fixed index out of range");
1816 
1817  // Check slice shape compatibility
1818  unsigned long idx = 0UL;
1819  unsigned sliceDim = 0U;
1820  for (unsigned i=0; i<dim_; ++i)
1821  {
1822  bool fixed = false;
1823  for (unsigned j=0; j<nFixedIndices; ++j)
1824  if (fixedIndices[j] == i)
1825  {
1826  fixed = true;
1827  if (fixedIndexValues[j] >= shape_[i])
1829  "In npstat::ArrayND::verifySliceCompatibility: "
1830  "fixed index value out of range");
1831  idx += fixedIndexValues[j]*strides_[i];
1832  break;
1833  }
1834  if (!fixed)
1835  {
1836  if (shape_[i] != slice.shape_[sliceDim])
1838  "In npstat::ArrayND::verifySliceCompatibility: "
1839  "incompatible argument array dimensions");
1840  ++sliceDim;
1841  }
1842  }
1843  assert(sliceDim == slice.dim_);
1844  return idx;
1845  }
1846 
1847  template<typename Numeric, unsigned StackLen, unsigned StackDim>
1848  unsigned long
1850  const unsigned long bufLen,
1851  const unsigned *fixedIndices,
1852  const unsigned *fixedIndexValues,
1853  const unsigned nFixedIndices,
1854  unsigned long* sliceStrides) const
1855  {
1856  if (!(nFixedIndices && nFixedIndices <= dim_))
1858  "In npstat::ArrayND::verifyBufferSliceCompatibility: "
1859  "invalid number of fixed indices");
1861  "Initialize npstat::ArrayND before calling "
1862  "method \"verifyBufferSliceCompatibility\"");
1863  assert(fixedIndices);
1864  assert(fixedIndexValues);
1865 
1866  for (unsigned j=0; j<nFixedIndices; ++j)
1867  if (fixedIndices[j] >= dim_) throw npstat::NpstatOutOfRange(
1868  "In npstat::ArrayND::verifyBufferSliceCompatibility: "
1869  "fixed index out of range");
1870 
1871  // Figure out slice shape. Store it, temporarily, in sliceStrides.
1872  unsigned long idx = 0UL;
1873  unsigned sliceDim = 0U;
1874  for (unsigned i=0; i<dim_; ++i)
1875  {
1876  bool fixed = false;
1877  for (unsigned j=0; j<nFixedIndices; ++j)
1878  if (fixedIndices[j] == i)
1879  {
1880  fixed = true;
1881  if (fixedIndexValues[j] >= shape_[i])
1883  "In npstat::ArrayND::verifyBufferSliceCompatibility:"
1884  " fixed index value out of range");
1885  idx += fixedIndexValues[j]*strides_[i];
1886  break;
1887  }
1888  if (!fixed)
1889  sliceStrides[sliceDim++] = shape_[i];
1890  }
1891  assert(sliceDim + nFixedIndices == dim_);
1892 
1893  // Convert the slice shape into slice strides
1894  unsigned long expectedBufLen = 1UL;
1895  if (sliceDim)
1896  {
1897  unsigned long shapeJ = sliceStrides[sliceDim - 1];
1898  sliceStrides[sliceDim - 1] = 1UL;
1899  for (unsigned j=sliceDim - 1; j>0; --j)
1900  {
1901  const unsigned long nextStride = sliceStrides[j]*shapeJ;
1902  shapeJ = sliceStrides[j - 1];
1903  sliceStrides[j - 1] = nextStride;
1904  }
1905  expectedBufLen = sliceStrides[0]*shapeJ;
1906  }
1907  if (expectedBufLen != bufLen)
1909  "In npstat::ArrayND::verifyBufferSliceCompatibility: "
1910  "invalid memory buffer length");
1911 
1912  return idx;
1913  }
1914 
1915  template<typename Numeric, unsigned StackLen, unsigned StackDim>
1916  template <typename Num2, class Op>
1918  const unsigned level, const unsigned long idx0,
1919  const unsigned level1, const unsigned long idx1,
1920  Num2* sliceData, const unsigned long* sliceStrides,
1921  const unsigned *fixedIndices,
1922  const unsigned *fixedIndexValues,
1923  const unsigned nFixedIndices,
1924  Op fcn)
1925  {
1926  bool fixed = false;
1927  for (unsigned j=0; j<nFixedIndices; ++j)
1928  if (fixedIndices[j] == level)
1929  {
1930  fixed = true;
1931  break;
1932  }
1933  if (fixed)
1934  jointSliceLoop(level+1, idx0, level1, idx1,
1935  sliceData, sliceStrides, fixedIndices,
1936  fixedIndexValues, nFixedIndices, fcn);
1937  else
1938  {
1939  const unsigned imax = shape_[level];
1940  const unsigned long stride = strides_[level];
1941 
1942  if (level1 == dim_ - nFixedIndices - 1)
1943  {
1944  sliceData += idx1;
1945  Numeric* localData = data_ + idx0;
1946  for (unsigned i = 0; i<imax; ++i)
1947  fcn(localData[i*stride], sliceData[i]);
1948  }
1949  else
1950  {
1951  const unsigned long stride2 = sliceStrides[level1];
1952  for (unsigned i = 0; i<imax; ++i)
1953  jointSliceLoop(level+1, idx0+i*stride,
1954  level1+1, idx1+i*stride2,
1955  sliceData, sliceStrides, fixedIndices,
1956  fixedIndexValues, nFixedIndices, fcn);
1957  }
1958  }
1959  }
1960 
1961  template<typename Numeric, unsigned StackLen, unsigned StackDim>
1962  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
1964  ArrayND<Num2,Len2,Dim2>& slice,
1965  const unsigned *fixedIndices,
1966  const unsigned *fixedIndexValues,
1967  const unsigned nFixedIndices,
1968  Functor binaryFunct)
1969  {
1970  const unsigned long idx = verifySliceCompatibility(
1971  slice, fixedIndices, fixedIndexValues, nFixedIndices);
1972  if (slice.dim_)
1973  jointSliceLoop(0U, idx, 0U, 0UL, slice.data_, slice.strides_,
1974  fixedIndices, fixedIndexValues,
1975  nFixedIndices, binaryFunct);
1976  else
1977  binaryFunct(data_[idx], slice.localData_[0]);
1978  }
1979 
1980  template<typename Numeric, unsigned StackLen, unsigned StackDim>
1981  template <typename Num2, class Functor>
1983  Num2* slice, const unsigned long len,
1984  const unsigned *fixedIndices, const unsigned *fixedIndexValues,
1985  unsigned nFixedIndices, Functor binaryFunct)
1986  {
1987  assert(slice);
1988  if (dim_ > CHAR_BIT*sizeof(unsigned long))
1990  "In npstat::ArrayND::jointMemSliceScan: "
1991  "rank of this array is too large");
1992  unsigned long sliceStrides[CHAR_BIT*sizeof(unsigned long)];
1993  const unsigned long idx = verifyBufferSliceCompatibility(
1994  len, fixedIndices, fixedIndexValues, nFixedIndices, sliceStrides);
1995  if (dim_ > nFixedIndices)
1996  jointSliceLoop(0U, idx, 0U, 0UL, slice, sliceStrides,
1997  fixedIndices, fixedIndexValues,
1998  nFixedIndices, binaryFunct);
1999  else
2000  binaryFunct(data_[idx], *slice);
2001  }
2002 
2003  template<typename Numeric, unsigned StackLen, unsigned StackDim>
2004  template <typename Num2>
2006  const unsigned level, unsigned long idx0,
2007  unsigned* currentIndex,
2009  const unsigned *projectedIndices,
2010  const unsigned nProjectedIndices) const
2011  {
2012  // level : dimension number among indices which are being iterated
2013  const unsigned idx = projectedIndices[level];
2014  const unsigned imax = shape_[idx];
2015  const unsigned long stride = strides_[idx];
2016  const bool last = (level == nProjectedIndices - 1);
2017 
2018  for (unsigned i = 0; i<imax; ++i)
2019  {
2020  currentIndex[idx] = i;
2021  if (last)
2022  projector.process(currentIndex, dim_, idx0, data_[idx0]);
2023  else
2024  projectInnerLoop(level+1, idx0, currentIndex, projector,
2025  projectedIndices, nProjectedIndices);
2026  idx0 += stride;
2027  }
2028  }
2029 
2030  template<typename Numeric, unsigned StackLen, unsigned StackDim>
2031  template <typename Num2, unsigned Len2, unsigned Dim2,
2032  typename Num3, class Op>
2034  const unsigned level, const unsigned long idx0,
2035  const unsigned level1, const unsigned long idx1,
2036  unsigned* currentIndex,
2037  ArrayND<Num2,Len2,Dim2>* projection,
2039  const unsigned *projectedIndices,
2040  const unsigned nProjectedIndices, Op fcn) const
2041  {
2042  // level : dimension number in this array
2043  // level1 : dimension number in the projection
2044  // idx0 : linear index in this array
2045  // idx1 : linear index in the projection
2046  // currentIndex : cycled over in this loop with the exception of the
2047  // dimensions which are iterated over to build the
2048  // projection
2049  if (level == dim_)
2050  {
2051  assert(level1 == projection->dim_);
2052  projector.clear();
2053  projectInnerLoop(0U, idx0, currentIndex, projector,
2054  projectedIndices, nProjectedIndices);
2055  if (projection->dim_)
2056  fcn(projection->data_[idx1], projector.result());
2057  else
2058  fcn(projection->localData_[0], projector.result());
2059  }
2060  else
2061  {
2062  bool iterated = false;
2063  for (unsigned j=0; j<nProjectedIndices; ++j)
2064  if (projectedIndices[j] == level)
2065  {
2066  iterated = true;
2067  break;
2068  }
2069  if (iterated)
2070  {
2071  // This index will be iterated over inside "projectInnerLoop"
2072  projectLoop(level+1, idx0, level1, idx1,
2073  currentIndex, projection, projector,
2074  projectedIndices, nProjectedIndices, fcn);
2075  }
2076  else
2077  {
2078  const unsigned imax = shape_[level];
2079  const unsigned long stride = strides_[level];
2080  // We will not be able to get here if projection->dim_ is 0.
2081  // Therefore, it is safe to access projection->strides_.
2082  const unsigned long stride2 = projection->strides_[level1];
2083  for (unsigned i = 0; i<imax; ++i)
2084  {
2085  currentIndex[level] = i;
2086  projectLoop(level+1, idx0+i*stride,
2087  level1+1, idx1+i*stride2,
2088  currentIndex, projection, projector,
2089  projectedIndices, nProjectedIndices, fcn);
2090  }
2091  }
2092  }
2093  }
2094 
2095  template<typename Numeric, unsigned StackLen, unsigned StackDim>
2096  template <typename Num2, unsigned Len2, unsigned Dim2>
2098  const ArrayND<Num2,Len2,Dim2>& projection,
2099  const unsigned *projectedIndices,
2100  const unsigned nProjectedIndices) const
2101  {
2102  if (!(nProjectedIndices && nProjectedIndices <= dim_))
2104  "In npstat::ArrayND::verifyProjectionCompatibility: "
2105  "invalid number of projected indices");
2107  "Initialize npstat::ArrayND before calling "
2108  "method \"verifyProjectionCompatibility\"");
2109  if (!projection.shapeIsKnown_) throw npstat::NpstatInvalidArgument(
2110  "In npstat::ArrayND::verifyProjectionCompatibility: "
2111  "uninitialized argument array");
2112  if (projection.dim_ != dim_ - nProjectedIndices)
2114  "In npstat::ArrayND::verifyProjectionCompatibility: "
2115  "incompatible argument array rank");
2116  assert(projectedIndices);
2117 
2118  for (unsigned j=0; j<nProjectedIndices; ++j)
2119  if (projectedIndices[j] >= dim_) throw npstat::NpstatOutOfRange(
2120  "In npstat::ArrayND::verifyProjectionCompatibility: "
2121  "projected index out of range");
2122 
2123  // Check projection shape compatibility
2124  unsigned sliceDim = 0U;
2125  for (unsigned i=0; i<dim_; ++i)
2126  {
2127  bool fixed = false;
2128  for (unsigned j=0; j<nProjectedIndices; ++j)
2129  if (projectedIndices[j] == i)
2130  {
2131  fixed = true;
2132  break;
2133  }
2134  if (!fixed)
2135  {
2136  if (shape_[i] != projection.shape_[sliceDim])
2138  "In npstat::ArrayND::verifyProjectionCompatibility: "
2139  "incompatible argument array dimensions");
2140  ++sliceDim;
2141  }
2142  }
2143  assert(sliceDim == projection.dim_);
2144  }
2145 
2146  template<typename Numeric, unsigned StackLen, unsigned StackDim>
2147  template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
2149  ArrayND<Num2,Len2,Dim2>* projection,
2151  const unsigned *projectedIndices,
2152  const unsigned nProjectedIndices) const
2153  {
2154  assert(projection);
2155  verifyProjectionCompatibility(*projection, projectedIndices,
2156  nProjectedIndices);
2157  unsigned ibuf[StackDim];
2158  unsigned* buf = makeBuffer(dim_, ibuf, StackDim);
2159  for (unsigned i=0; i<dim_; ++i)
2160  buf[i] = 0U;
2161  projectLoop(0U, 0UL, 0U, 0UL, buf, projection,
2162  projector, projectedIndices, nProjectedIndices,
2164  destroyBuffer(buf, ibuf);
2165  }
2166 
2167  template<typename Numeric, unsigned StackLen, unsigned StackDim>
2168  template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
2170  ArrayND<Num2,Len2,Dim2>* projection,
2172  const unsigned *projectedIndices,
2173  const unsigned nProjectedIndices) const
2174  {
2175  assert(projection);
2176  verifyProjectionCompatibility(*projection, projectedIndices,
2177  nProjectedIndices);
2178  unsigned ibuf[StackDim];
2179  unsigned* buf = makeBuffer(dim_, ibuf, StackDim);
2180  for (unsigned i=0; i<dim_; ++i)
2181  buf[i] = 0U;
2182  projectLoop(0U, 0UL, 0U, 0UL, buf, projection,
2183  projector, projectedIndices, nProjectedIndices,
2185  destroyBuffer(buf, ibuf);
2186  }
2187 
2188  template<typename Numeric, unsigned StackLen, unsigned StackDim>
2189  template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
2191  ArrayND<Num2,Len2,Dim2>* projection,
2193  const unsigned *projectedIndices,
2194  const unsigned nProjectedIndices) const
2195  {
2196  assert(projection);
2197  verifyProjectionCompatibility(*projection, projectedIndices,
2198  nProjectedIndices);
2199  unsigned ibuf[StackDim];
2200  unsigned* buf = makeBuffer(dim_, ibuf, StackDim);
2201  for (unsigned i=0; i<dim_; ++i)
2202  buf[i] = 0U;
2203  projectLoop(0U, 0UL, 0U, 0UL, buf, projection,
2204  projector, projectedIndices, nProjectedIndices,
2206  destroyBuffer(buf, ibuf);
2207  }
2208 
2209  template<typename Numeric, unsigned StackLen, unsigned StackDim>
2210  template <typename Num2>
2212  const unsigned level, const unsigned long idx0,
2213  AbsVisitor<Numeric,Num2>& projector,
2214  const unsigned *projectedIndices,
2215  const unsigned nProjectedIndices) const
2216  {
2217  const unsigned idx = projectedIndices[level];
2218  const unsigned imax = shape_[idx];
2219  const unsigned long stride = strides_[idx];
2220  const bool last = (level == nProjectedIndices - 1);
2221 
2222  for (unsigned i = 0; i<imax; ++i)
2223  {
2224  if (last)
2225  projector.process(data_[idx0+i*stride]);
2226  else
2227  projectInnerLoop2(level+1, idx0+i*stride, projector,
2228  projectedIndices, nProjectedIndices);
2229  }
2230  }
2231 
2232  template<typename Numeric, unsigned StackLen, unsigned StackDim>
2233  template <typename Num2, unsigned Len2, unsigned Dim2,
2234  typename Num3, class Op>
2236  const unsigned level, const unsigned long idx0,
2237  const unsigned level1, const unsigned long idx1,
2238  ArrayND<Num2,Len2,Dim2>* projection,
2239  AbsVisitor<Numeric,Num3>& projector,
2240  const unsigned *projectedIndices,
2241  const unsigned nProjectedIndices, Op fcn) const
2242  {
2243  if (level == dim_)
2244  {
2245  assert(level1 == projection->dim_);
2246  projector.clear();
2247  projectInnerLoop2(0U, idx0, projector,
2248  projectedIndices, nProjectedIndices);
2249  if (projection->dim_)
2250  fcn(projection->data_[idx1], projector.result());
2251  else
2252  fcn(projection->localData_[0], projector.result());
2253  }
2254  else
2255  {
2256  bool fixed = false;
2257  for (unsigned j=0; j<nProjectedIndices; ++j)
2258  if (projectedIndices[j] == level)
2259  {
2260  fixed = true;
2261  break;
2262  }
2263  if (fixed)
2264  {
2265  projectLoop2(level+1, idx0, level1, idx1,
2266  projection, projector,
2267  projectedIndices, nProjectedIndices, fcn);
2268  }
2269  else
2270  {
2271  const unsigned imax = shape_[level];
2272  const unsigned long stride = strides_[level];
2273  const unsigned long stride2 = projection->strides_[level1];
2274  for (unsigned i = 0; i<imax; ++i)
2275  projectLoop2(level+1, idx0+i*stride,
2276  level1+1, idx1+i*stride2,
2277  projection, projector,
2278  projectedIndices, nProjectedIndices, fcn);
2279  }
2280  }
2281  }
2282 
2283  template<typename Numeric, unsigned StackLen, unsigned StackDim>
2284  template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
2286  ArrayND<Num2,Len2,Dim2>* projection,
2287  AbsVisitor<Numeric,Num3>& projector,
2288  const unsigned *projectedIndices,
2289  const unsigned nProjectedIndices) const
2290  {
2291  assert(projection);
2292  verifyProjectionCompatibility(*projection, projectedIndices,
2293  nProjectedIndices);
2294  projectLoop2(0U, 0UL, 0U, 0UL, projection,
2295  projector, projectedIndices, nProjectedIndices,
2297  }
2298 
2299  template<typename Numeric, unsigned StackLen, unsigned StackDim>
2300  template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
2302  ArrayND<Num2,Len2,Dim2>* projection,
2303  AbsVisitor<Numeric,Num3>& projector,
2304  const unsigned *projectedIndices,
2305  const unsigned nProjectedIndices) const
2306  {
2307  assert(projection);
2308  verifyProjectionCompatibility(*projection, projectedIndices,
2309  nProjectedIndices);
2310  projectLoop2(0U, 0UL, 0U, 0UL, projection,
2311  projector, projectedIndices, nProjectedIndices,
2313  }
2314 
2315  template<typename Numeric, unsigned StackLen, unsigned StackDim>
2316  template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
2318  ArrayND<Num2,Len2,Dim2>* projection,
2319  AbsVisitor<Numeric,Num3>& projector,
2320  const unsigned *projectedIndices,
2321  const unsigned nProjectedIndices) const
2322  {
2323  assert(projection);
2324  verifyProjectionCompatibility(*projection, projectedIndices,
2325  nProjectedIndices);
2326  projectLoop2(0U, 0UL, 0U, 0UL, projection,
2327  projector, projectedIndices, nProjectedIndices,
2329  }
2330 
2331  template<typename Numeric, unsigned StackLen, unsigned StackDim>
2332  template <typename Num2, class Functor>
2334  const unsigned level, const unsigned long idx0,
2335  Num2& scale, const unsigned *projectedIndices,
2336  const unsigned nProjectedIndices, Functor binaryFunct)
2337  {
2338  const unsigned idx = projectedIndices[level];
2339  const unsigned imax = shape_[idx];
2340  const unsigned long stride = strides_[idx];
2341 
2342  if (level == nProjectedIndices - 1)
2343  {
2344  Numeric* data = data_ + idx0;
2345  for (unsigned i = 0; i<imax; ++i)
2346  binaryFunct(data[i*stride], scale);
2347  }
2348  else
2349  for (unsigned i = 0; i<imax; ++i)
2350  scaleBySliceInnerLoop(level+1, idx0+i*stride, scale,
2351  projectedIndices, nProjectedIndices,
2352  binaryFunct);
2353  }
2354 
2355  template<typename Numeric, unsigned StackLen, unsigned StackDim>
2356  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
2358  unsigned level, unsigned long idx0,
2359  unsigned level1, unsigned long idx1,
2360  ArrayND<Num2,Len2,Dim2>& slice,
2361  const unsigned *projectedIndices,
2362  const unsigned nProjectedIndices,
2363  Functor binaryFunct)
2364  {
2365  if (level == dim_)
2366  {
2367  assert(level1 == slice.dim_);
2368  Num2& scaleFactor = slice.dim_ ? slice.data_[idx1] :
2369  slice.localData_[0];
2370  scaleBySliceInnerLoop(0U, idx0, scaleFactor, projectedIndices,
2371  nProjectedIndices, binaryFunct);
2372  }
2373  else
2374  {
2375  bool fixed = false;
2376  for (unsigned j=0; j<nProjectedIndices; ++j)
2377  if (projectedIndices[j] == level)
2378  {
2379  fixed = true;
2380  break;
2381  }
2382  if (fixed)
2383  {
2384  scaleBySliceLoop(level+1, idx0, level1, idx1, slice,
2385  projectedIndices, nProjectedIndices,
2386  binaryFunct);
2387  }
2388  else
2389  {
2390  const unsigned imax = shape_[level];
2391  const unsigned long stride = strides_[level];
2392  const unsigned long stride2 = slice.strides_[level1];
2393  for (unsigned i = 0; i<imax; ++i)
2394  scaleBySliceLoop(level+1, idx0+i*stride, level1+1,
2395  idx1+i*stride2, slice, projectedIndices,
2396  nProjectedIndices, binaryFunct);
2397  }
2398  }
2399  }
2400 
2401  template<typename Numeric, unsigned StackLen, unsigned StackDim>
2402  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
2404  ArrayND<Num2, Len2, Dim2>& r, Functor binaryFunct)
2405  {
2407  "In npstat::ArrayND::jointScan: incompatible argument array shape");
2408  if (dim_)
2409  for (unsigned long i=0; i<len_; ++i)
2410  binaryFunct(data_[i], r.data_[i]);
2411  else
2412  binaryFunct(localData_[0], r.localData_[0]);
2413  }
2414 
2415  template<typename Numeric, unsigned StackLen, unsigned StackDim>
2416  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
2418  ArrayND<Num2,Len2,Dim2>& slice,
2419  const unsigned *fixedIndices, const unsigned nFixedIndices,
2420  Functor binaryFunct)
2421  {
2422  if (nFixedIndices)
2423  {
2424  verifyProjectionCompatibility(slice, fixedIndices, nFixedIndices);
2425  if (slice.dim_ == 0U)
2426  for (unsigned long i=0; i<len_; ++i)
2427  binaryFunct(data_[i], slice.localData_[0]);
2428  else
2429  scaleBySliceLoop(0U, 0UL, 0U, 0UL, slice,
2430  fixedIndices, nFixedIndices, binaryFunct);
2431  }
2432  else
2433  jointScan(slice, binaryFunct);
2434  }
2435 
2436  template<typename Numeric, unsigned StackLen, unsigned StackDim>
2438  const ArrayShape& shape) const
2439  {
2440  if (!shapeIsKnown_)
2441  return false;
2442  if (dim_ != shape.size())
2443  return false;
2444  if (dim_)
2445  {
2446  for (unsigned i=0; i<dim_; ++i)
2447  if (shape_[i] != shape[i])
2448  return false;
2449  }
2450  else
2451  assert(len_ == 1UL);
2452  return true;
2453  }
2454 
2455  template<typename Numeric, unsigned StackLen, unsigned StackDim>
2456  template <typename Num2, unsigned Len2, unsigned Dim2>
2458  const ArrayND<Num2,Len2,Dim2>& r) const
2459  {
2460  if (!shapeIsKnown_)
2461  return false;
2462  if (!r.shapeIsKnown_)
2463  return false;
2464  if (dim_ != r.dim_)
2465  return false;
2466  if (len_ != r.len_)
2467  return false;
2468  if (dim_)
2469  {
2470  assert(shape_);
2471  assert(r.shape_);
2472  for (unsigned i=0; i<dim_; ++i)
2473  if (shape_[i] != r.shape_[i])
2474  return false;
2475  }
2476  else
2477  assert(len_ == 1UL);
2478  return true;
2479  }
2480 
2481  template<typename Numeric, unsigned Len, unsigned Dim>
2482  template<typename Num2, typename Integer>
2484  const unsigned level, unsigned long idx0,
2485  unsigned* currentIndex,
2487  const BoxND<Integer>& subrange) const
2488  {
2489  // Deal with possible negative limits first
2490  const Interval<Integer>& levelRange(subrange[level]);
2491  long long int iminl = static_cast<long long int>(levelRange.min());
2492  if (iminl < 0LL) iminl = 0LL;
2493  long long int imaxl = static_cast<long long int>(levelRange.max());
2494  if (imaxl < 0LL) imaxl = 0LL;
2495 
2496  // Now deal with possible out-of-range limits
2497  const unsigned imin = static_cast<unsigned>(iminl);
2498  unsigned imax = static_cast<unsigned>(imaxl);
2499  if (imax > shape_[level])
2500  imax = shape_[level];
2501 
2502  if (level == dim_ - 1)
2503  {
2504  idx0 += imin;
2505  for (unsigned i=imin; i<imax; ++i, ++idx0)
2506  {
2507  currentIndex[level] = i;
2508  f.process(currentIndex, dim_, idx0, data_[idx0]);
2509  }
2510  }
2511  else
2512  {
2513  const unsigned long stride = strides_[level];
2514  idx0 += imin*stride;
2515  for (unsigned i=imin; i<imax; ++i)
2516  {
2517  currentIndex[level] = i;
2518  processSubrangeLoop(level+1U, idx0, currentIndex, f, subrange);
2519  idx0 += stride;
2520  }
2521  }
2522  }
2523 
2524  template<typename Numeric, unsigned Len, unsigned StackDim>
2525  template <typename Num2, typename Integer>
2528  const BoxND<Integer>& subrange) const
2529  {
2531  "Initialize npstat::ArrayND before calling method \"processSubrange\"");
2532  if (!dim_) throw npstat::NpstatInvalidArgument(
2533  "npstat::ArrayND::processSubrange method "
2534  "can not be used with array of 0 rank");
2535  if (dim_ != subrange.dim()) throw npstat::NpstatInvalidArgument(
2536  "In npstat::ArrayND::processSubrange: incompatible subrange rank");
2537  unsigned ibuf[StackDim];
2538  unsigned* buf = makeBuffer(dim_, ibuf, StackDim);
2539  for (unsigned i=0; i<dim_; ++i)
2540  buf[i] = 0U;
2541  processSubrangeLoop(0U, 0UL, buf, f, subrange);
2542  destroyBuffer(buf, ibuf);
2543  }
2544 
2545  template<typename Numeric, unsigned Len, unsigned Dim>
2546  template<typename Num2>
2548  const Num2* data, const unsigned long dataLength)
2549  {
2551  "Initialize npstat::ArrayND before calling method \"setData\"");
2552  if (dataLength != len_) throw npstat::NpstatInvalidArgument(
2553  "In npstat::ArrayND::setData: incompatible input data length");
2554  copyBuffer(data_, data, dataLength);
2555  return *this;
2556  }
2557 
2558  template<typename Numeric, unsigned Len, unsigned Dim>
2560  {
2561  assert(dim_);
2562  if (strides_ == nullptr)
2564  strides_[dim_ - 1] = 1UL;
2565  for (unsigned j=dim_ - 1; j>0; --j)
2566  strides_[j - 1] = strides_[j]*shape_[j];
2567  }
2568 
2569  template<typename Numeric, unsigned StackLen, unsigned StackDim>
2572  len_(0UL), dim_(0U), shapeIsKnown_(false)
2573  {
2574  localData_[0] = Numeric();
2575  data_ = localData_;
2576  }
2577 
2578  template<typename Numeric, unsigned Len, unsigned Dim>
2582  {
2583  if (dim_)
2584  {
2585  // Copy the shape
2586  shape_ = makeBuffer(dim_, localShape_, Dim);
2588 
2589  // Copy the strides
2592 
2593  // Copy the data
2594  data_ = makeBuffer(len_, localData_, Len);
2595  copyBuffer(data_, r.data_, len_);
2596  }
2597  else
2598  {
2599  assert(len_ == 1UL);
2600  localData_[0] = r.localData_[0];
2601  data_ = localData_;
2602  }
2603  }
2604 
2605  template<typename Numeric, unsigned Len, unsigned Dim>
2606  template <typename Num2, unsigned Len2, unsigned Dim2>
2608  : data_(0), strides_(nullptr), shape_(nullptr),
2610  {
2611  if (dim_)
2612  {
2613  // Copy the shape
2614  shape_ = makeBuffer(dim_, localShape_, Dim);
2616 
2617  // Copy the strides
2620 
2621  // Copy the data
2622  data_ = makeBuffer(len_, localData_, Len);
2623  copyBuffer(data_, r.data_, len_);
2624  }
2625  else
2626  {
2627  assert(len_ == 1UL);
2628  localData_[0] = static_cast<Numeric>(r.localData_[0]);
2629  data_ = localData_;
2630  }
2631  }
2632 
2633  template<typename Numeric, unsigned Len, unsigned Dim>
2634  template<typename Num2, unsigned Len2, unsigned Dim2, class Functor>
2636  Functor f)
2637  : data_(0), strides_(nullptr), shape_(nullptr),
2639  {
2640  if (dim_)
2641  {
2642  // Copy the shape
2643  shape_ = makeBuffer(dim_, localShape_, Dim);
2645 
2646  // Copy the strides
2649 
2650  // Copy the data
2651  data_ = makeBuffer(len_, localData_, Len);
2652  for (unsigned long i=0; i<len_; ++i)
2653  data_[i] = static_cast<Numeric>(f(r.data_[i]));
2654  }
2655  else
2656  {
2657  assert(len_ == 1UL);
2658  localData_[0] = static_cast<Numeric>(f(r.localData_[0]));
2659  data_ = localData_;
2660  }
2661  }
2662 
2663  template<typename Numeric, unsigned Len, unsigned Dim>
2664  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
2666  const unsigned level, unsigned long idx0,
2667  unsigned long idx1,
2669  const ArrayRange& range, Functor f)
2670  {
2671  const unsigned imax = shape_[level];
2672  if (level == dim_ - 1)
2673  {
2674  Numeric* to = data_ + idx0;
2675  const Num2* from = r.data_ + (idx1 + range[level].min());
2676  for (unsigned i=0; i<imax; ++i)
2677  *to++ = static_cast<Numeric>(f(*from++));
2678  }
2679  else
2680  {
2681  const unsigned long fromstride = r.strides_[level];
2682  const unsigned long tostride = strides_[level];
2683  idx1 += range[level].min()*fromstride;
2684  for (unsigned i=0; i<imax; ++i)
2685  {
2686  copyRangeLoopFunct(level+1, idx0, idx1, r, range, f);
2687  idx0 += tostride;
2688  idx1 += fromstride;
2689  }
2690  }
2691  }
2692 
2693  template<typename Numeric, unsigned Len, unsigned Dim>
2694  template <typename Num2, unsigned Len2, unsigned Dim2>
2696  const ArrayND<Num2, Len2, Dim2>& r, const ArrayRange& range)
2697  : data_(0), strides_(nullptr), shape_(nullptr),
2699  {
2700  if (!range.isCompatible(r.shape_, r.dim_))
2702  "In npstat::ArrayND subrange constructor: invalid subrange");
2703  if (dim_)
2704  {
2705  len_ = range.rangeSize();
2706  if (!len_)
2708  "In npstat::ArrayND subrange constructor: empty subrange");
2709 
2710  // Figure out the shape
2711  shape_ = makeBuffer(dim_, localShape_, Dim);
2712  range.rangeLength(shape_, dim_);
2713 
2714  // Figure out the strides
2715  buildStrides();
2716 
2717  // Allocate the data array
2718  data_ = makeBuffer(len_, localData_, Len);
2719 
2720  // Copy the data
2721  if (dim_ > CHAR_BIT*sizeof(unsigned long))
2723  "In npstat::ArrayND subrange constructor: "
2724  "input array rank is too large");
2725  unsigned lolim[CHAR_BIT*sizeof(unsigned long)];
2726  range.lowerLimits(lolim, dim_);
2727  unsigned toBuf[CHAR_BIT*sizeof(unsigned long)];
2728  clearBuffer(toBuf, dim_);
2729  (const_cast<ArrayND<Num2, Len2, Dim2>&>(r)).commonSubrangeLoop(
2730  0U, 0UL, 0UL, lolim, shape_, toBuf, *this,
2732  }
2733  else
2734  {
2735  assert(len_ == 1UL);
2736  localData_[0] = static_cast<Numeric>(r.localData_[0]);
2737  data_ = localData_;
2738  }
2739  }
2740 
2741  template<typename Numeric, unsigned Len, unsigned Dim>
2742  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
2744  const ArrayND<Num2, Len2, Dim2>& r, const ArrayRange& range,
2745  Functor f)
2746  : data_(0), strides_(nullptr), shape_(nullptr),
2748  {
2749  if (!range.isCompatible(r.shape_, r.dim_))
2751  "In npstat::ArrayND transforming subrange constructor: "
2752  "incompatible subrange");
2753  if (dim_)
2754  {
2755  len_ = range.rangeSize();
2756  if (!len_)
2758  "In npstat::ArrayND transforming subrange constructor: "
2759  "empty subrange");
2760 
2761  // Figure out the shape
2762  shape_ = makeBuffer(dim_, localShape_, Dim);
2763  for (unsigned i=0; i<dim_; ++i)
2764  shape_[i] = range[i].length();
2765 
2766  // Figure out the strides
2767  buildStrides();
2768 
2769  // Allocate the data array
2770  data_ = makeBuffer(len_, localData_, Len);
2771 
2772  // Transform the data
2773  copyRangeLoopFunct(0U, 0UL, 0UL, r, range, f);
2774  }
2775  else
2776  {
2777  assert(len_ == 1UL);
2778  localData_[0] = static_cast<Numeric>(f(r.localData_[0]));
2779  data_ = localData_;
2780  }
2781  }
2782 
2783  template<typename Numeric, unsigned Len, unsigned Dim>
2786  {
2787  const unsigned sz = sh.size();
2788  buildFromShapePtr(sz ? &sh[0] : nullptr, sz);
2789  }
2790 
2791  template<typename Numeric, unsigned Len, unsigned Dim>
2792  ArrayND<Numeric,Len,Dim>::ArrayND(const unsigned* sizes,
2793  const unsigned dim)
2795  {
2796  buildFromShapePtr(sizes, dim);
2797  }
2798 
2799  template<typename Numeric, unsigned Len, unsigned Dim>
2802  {
2803  const unsigned dim = 1U;
2804  unsigned sizes[dim];
2805  sizes[0] = n0;
2806  buildFromShapePtr(sizes, dim);
2807  }
2808 
2809  template<typename Numeric, unsigned Len, unsigned Dim>
2811  const unsigned n1)
2813  {
2814  const unsigned dim = 2U;
2815  unsigned sizes[dim];
2816  sizes[0] = n0;
2817  sizes[1] = n1;
2818  buildFromShapePtr(sizes, dim);
2819  }
2820 
2821  template<typename Numeric, unsigned Len, unsigned Dim>
2823  const unsigned n1,
2824  const unsigned n2)
2826  {
2827  const unsigned dim = 3U;
2828  unsigned sizes[dim];
2829  sizes[0] = n0;
2830  sizes[1] = n1;
2831  sizes[2] = n2;
2832  buildFromShapePtr(sizes, dim);
2833  }
2834 
2835  template<typename Numeric, unsigned Len, unsigned Dim>
2837  const unsigned n1,
2838  const unsigned n2,
2839  const unsigned n3)
2841  {
2842  const unsigned dim = 4U;
2843  unsigned sizes[dim];
2844  sizes[0] = n0;
2845  sizes[1] = n1;
2846  sizes[2] = n2;
2847  sizes[3] = n3;
2848  buildFromShapePtr(sizes, dim);
2849  }
2850 
2851  template<typename Numeric, unsigned Len, unsigned Dim>
2853  const unsigned n1,
2854  const unsigned n2,
2855  const unsigned n3,
2856  const unsigned n4)
2858  {
2859  const unsigned dim = 5U;
2860  unsigned sizes[dim];
2861  sizes[0] = n0;
2862  sizes[1] = n1;
2863  sizes[2] = n2;
2864  sizes[3] = n3;
2865  sizes[4] = n4;
2866  buildFromShapePtr(sizes, dim);
2867  }
2868 
2869  template<typename Numeric, unsigned Len, unsigned Dim>
2871  const unsigned n1,
2872  const unsigned n2,
2873  const unsigned n3,
2874  const unsigned n4,
2875  const unsigned n5)
2877  {
2878  const unsigned dim = 6U;
2879  unsigned sizes[dim];
2880  sizes[0] = n0;
2881  sizes[1] = n1;
2882  sizes[2] = n2;
2883  sizes[3] = n3;
2884  sizes[4] = n4;
2885  sizes[5] = n5;
2886  buildFromShapePtr(sizes, dim);
2887  }
2888 
2889  template<typename Numeric, unsigned Len, unsigned Dim>
2891  const unsigned n1,
2892  const unsigned n2,
2893  const unsigned n3,
2894  const unsigned n4,
2895  const unsigned n5,
2896  const unsigned n6)
2898  {
2899  const unsigned dim = 7U;
2900  unsigned sizes[dim];
2901  sizes[0] = n0;
2902  sizes[1] = n1;
2903  sizes[2] = n2;
2904  sizes[3] = n3;
2905  sizes[4] = n4;
2906  sizes[5] = n5;
2907  sizes[6] = n6;
2908  buildFromShapePtr(sizes, dim);
2909  }
2910 
2911  template<typename Numeric, unsigned Len, unsigned Dim>
2913  const unsigned n1,
2914  const unsigned n2,
2915  const unsigned n3,
2916  const unsigned n4,
2917  const unsigned n5,
2918  const unsigned n6,
2919  const unsigned n7)
2921  {
2922  const unsigned dim = 8U;
2923  unsigned sizes[dim];
2924  sizes[0] = n0;
2925  sizes[1] = n1;
2926  sizes[2] = n2;
2927  sizes[3] = n3;
2928  sizes[4] = n4;
2929  sizes[5] = n5;
2930  sizes[6] = n6;
2931  sizes[7] = n7;
2932  buildFromShapePtr(sizes, dim);
2933  }
2934 
2935  template<typename Numeric, unsigned Len, unsigned Dim>
2937  const unsigned n1,
2938  const unsigned n2,
2939  const unsigned n3,
2940  const unsigned n4,
2941  const unsigned n5,
2942  const unsigned n6,
2943  const unsigned n7,
2944  const unsigned n8)
2946  {
2947  const unsigned dim = 9U;
2948  unsigned sizes[dim];
2949  sizes[0] = n0;
2950  sizes[1] = n1;
2951  sizes[2] = n2;
2952  sizes[3] = n3;
2953  sizes[4] = n4;
2954  sizes[5] = n5;
2955  sizes[6] = n6;
2956  sizes[7] = n7;
2957  sizes[8] = n8;
2958  buildFromShapePtr(sizes, dim);
2959  }
2960 
2961  template<typename Numeric, unsigned Len, unsigned Dim>
2963  const unsigned n1,
2964  const unsigned n2,
2965  const unsigned n3,
2966  const unsigned n4,
2967  const unsigned n5,
2968  const unsigned n6,
2969  const unsigned n7,
2970  const unsigned n8,
2971  const unsigned n9)
2973  {
2974  const unsigned dim = 10U;
2975  unsigned sizes[dim];
2976  sizes[0] = n0;
2977  sizes[1] = n1;
2978  sizes[2] = n2;
2979  sizes[3] = n3;
2980  sizes[4] = n4;
2981  sizes[5] = n5;
2982  sizes[6] = n6;
2983  sizes[7] = n7;
2984  sizes[8] = n8;
2985  sizes[9] = n9;
2986  buildFromShapePtr(sizes, dim);
2987  }
2988 
2989  template<typename Numeric, unsigned Len, unsigned Dim>
2990  template<typename Num1, unsigned Len1, unsigned Dim1,
2991  typename Num2, unsigned Len2, unsigned Dim2>
2993  const unsigned level, unsigned long idx0,
2994  unsigned long idx1, unsigned long idx2,
2995  const ArrayND<Num1, Len1, Dim1>& a1,
2996  const ArrayND<Num2, Len2, Dim2>& a2)
2997  {
2998  const unsigned imax = shape_[level];
2999  if (level == dim_ - 1)
3000  {
3001  for (unsigned i=0; i<imax; ++i)
3002  data_[idx0 + i] = a1.data_[idx1]*a2.data_[idx2 + i];
3003  }
3004  else
3005  {
3006  for (unsigned i=0; i<imax; ++i)
3007  {
3008  outerProductLoop(level+1, idx0, idx1, idx2, a1, a2);
3009  idx0 += strides_[level];
3010  if (level < a1.dim_)
3011  idx1 += a1.strides_[level];
3012  else
3013  idx2 += a2.strides_[level - a1.dim_];
3014  }
3015  }
3016  }
3017 
3018  template<typename Numeric, unsigned Len, unsigned Dim>
3019  template<typename Num1, unsigned Len1, unsigned Dim1,
3020  typename Num2, unsigned Len2, unsigned Dim2>
3022  const ArrayND<Num2, Len2, Dim2>& a2)
3023  : data_(0), strides_(nullptr), shape_(nullptr),
3024  len_(1UL), dim_(a1.dim_ + a2.dim_), shapeIsKnown_(true)
3025  {
3026  if (!(a1.shapeIsKnown_ && a2.shapeIsKnown_))
3028  "In npstat::ArrayND outer product constructor: "
3029  "uninitialized argument array");
3030  if (dim_)
3031  {
3032  shape_ = makeBuffer(dim_, localShape_, Dim);
3033  copyBuffer(shape_, a1.shape_, a1.dim_);
3034  copyBuffer(shape_+a1.dim_, a2.shape_, a2.dim_);
3035 
3036  for (unsigned i=0; i<dim_; ++i)
3037  {
3038  assert(shape_[i]);
3039  len_ *= shape_[i];
3040  }
3041 
3042  // Figure out the array strides
3043  buildStrides();
3044 
3045  // Allocate the data array
3046  data_ = makeBuffer(len_, localData_, Len);
3047 
3048  // Fill the data array
3049  if (a1.dim_ == 0)
3050  {
3051  for (unsigned long i=0; i<len_; ++i)
3052  data_[i] = a1.localData_[0] * a2.data_[i];
3053  }
3054  else if (a2.dim_ == 0)
3055  {
3056  for (unsigned long i=0; i<len_; ++i)
3057  data_[i] = a1.data_[i] * a2.localData_[0];
3058  }
3059  else
3060  outerProductLoop(0U, 0UL, 0UL, 0UL, a1, a2);
3061  }
3062  else
3063  {
3064  localData_[0] = a1.localData_[0] * a2.localData_[0];
3065  data_ = localData_;
3066  }
3067  }
3068 
3069  template<typename Numeric, unsigned Len, unsigned Dim>
3071  {
3075  }
3076 
3077  template<typename Numeric, unsigned Len, unsigned Dim>
3080  {
3081  if (this == &r)
3082  return *this;
3083  if (shapeIsKnown_)
3084  {
3086  "In npstat::ArrayND assignment operator: "
3087  "uninitialized argument array");
3089  "In npstat::ArrayND assignment operator: "
3090  "incompatible argument array shape");
3091  if (dim_)
3092  copyBuffer(data_, r.data_, len_);
3093  else
3094  localData_[0] = r.localData_[0];
3095  }
3096  else
3097  {
3098  // This object is uninitialized. If the object on the
3099  // right is itself initialized, make an in-place copy.
3100  if (r.shapeIsKnown_)
3101  new (this) ArrayND(r);
3102  }
3103  return *this;
3104  }
3105 
3106  template<typename Numeric, unsigned Len, unsigned Dim>
3107  template <typename Num2, unsigned Len2, unsigned Dim2>
3110  {
3111  if ((void*)this == (void*)(&r))
3112  return *this;
3113  if (shapeIsKnown_)
3114  {
3116  "In npstat::ArrayND assignment operator: "
3117  "uninitialized argument array");
3119  "In npstat::ArrayND assignment operator: "
3120  "incompatible argument array shape");
3121  if (dim_)
3122  copyBuffer(data_, r.data_, len_);
3123  else
3124  localData_[0] = static_cast<Numeric>(r.localData_[0]);
3125  }
3126  else
3127  {
3128  // This object is uninitialized. If the object on the
3129  // right is itself initialized, make an in-place copy.
3130  if (r.shapeIsKnown_)
3131  new (this) ArrayND(r);
3132  }
3133  return *this;
3134  }
3135 
3136  template<typename Numeric, unsigned Len, unsigned Dim>
3137  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
3140  Functor f)
3141  {
3142  if (shapeIsKnown_)
3143  {
3145  "In npstat::ArrayND::assign: uninitialized argument array");
3147  "In npstat::ArrayND::assign: incompatible argument array shape");
3148  if (dim_)
3149  for (unsigned long i=0; i<len_; ++i)
3150  data_[i] = static_cast<Numeric>(f(r.data_[i]));
3151  else
3152  localData_[0] = static_cast<Numeric>(f(r.localData_[0]));
3153  }
3154  else
3155  {
3156  // This object is uninitialized. If the object on the
3157  // right is itself initialized, build new array in place.
3158  if (r.shapeIsKnown_)
3159  new (this) ArrayND(r, f);
3160  }
3161  return *this;
3162  }
3163 
3164  template<typename Numeric, unsigned Len, unsigned Dim>
3166  {
3168  "Initialize npstat::ArrayND before calling method \"shape\"");
3169  return ArrayShape(shape_, shape_+dim_);
3170  }
3171 
3172  template<typename Numeric, unsigned Len, unsigned Dim>
3174  {
3176  "Initialize npstat::ArrayND before calling method \"fullRange\"");
3177  ArrayRange range;
3178  if (dim_)
3179  {
3180  range.reserve(dim_);
3181  for (unsigned i=0; i<dim_; ++i)
3182  range.push_back(Interval<unsigned>(0U, shape_[i]));
3183  }
3184  return range;
3185  }
3186 
3187  template<typename Numeric, unsigned Len, unsigned Dim>
3189  {
3191  "Initialize npstat::ArrayND before calling method \"isDensity\"");
3192  const Numeric zero = Numeric();
3193  bool hasPositive = false;
3194  if (dim_)
3195  for (unsigned long i=0; i<len_; ++i)
3196  {
3197  // Don't make comparisons whose result can be
3198  // determined in advance by assuming that Numeric
3199  // is an unsigned type. Some compilers will
3200  // complain about it when this template is
3201  // instantiated with such a type.
3202  if (data_[i] == zero)
3203  continue;
3205  hasPositive = true;
3206  else
3207  return false;
3208  }
3209  else
3211  zero, localData_[0]);
3212  return hasPositive;
3213  }
3214 
3215  template<typename Numeric, unsigned Len, unsigned Dim>
3217  {
3219  "Initialize npstat::ArrayND before calling method \"isZero\"");
3220  const Numeric zero = Numeric();
3221  if (dim_)
3222  {
3223  for (unsigned long i=0; i<len_; ++i)
3224  if (data_[i] != zero)
3225  return false;
3226  }
3227  else
3228  if (localData_[0] != zero)
3229  return false;
3230  return true;
3231  }
3232 
3233  template<typename Numeric, unsigned Len, unsigned Dim>
3235  unsigned long l, unsigned* idx, const unsigned idxLen) const
3236  {
3238  "Initialize npstat::ArrayND before calling "
3239  "method \"convertLinearIndex\"");
3240  if (!dim_) throw npstat::NpstatInvalidArgument(
3241  "npstat::ArrayND::convertLinearIndex method "
3242  "can not be used with array of 0 rank");
3243  if (idxLen != dim_) throw npstat::NpstatInvalidArgument(
3244  "In npstat::ArrayND::convertLinearIndex: incompatible index length");
3245  if (l >= len_) throw npstat::NpstatOutOfRange(
3246  "In npstat::ArrayND::convertLinearIndex: linear index out of range");
3247  assert(idx);
3248 
3249  for (unsigned i=0; i<dim_; ++i)
3250  {
3251  idx[i] = l / strides_[i];
3252  l -= (idx[i] * strides_[i]);
3253  }
3254  }
3255 
3256  template<typename Numeric, unsigned Len, unsigned Dim>
3258  const unsigned* index, unsigned idxLen) const
3259  {
3261  "Initialize npstat::ArrayND before calling method \"linearIndex\"");
3262  if (!dim_) throw npstat::NpstatInvalidArgument(
3263  "npstat::ArrayND::linearIndex method "
3264  "can not be used with array of 0 rank");
3265  if (idxLen != dim_) throw npstat::NpstatInvalidArgument(
3266  "In npstat::ArrayND::linearIndex: incompatible index length");
3267  assert(index);
3268 
3269  unsigned long idx = 0UL;
3270  for (unsigned i=0; i<dim_; ++i)
3271  {
3272  if (index[i] >= shape_[i])
3274  "In npstat::ArrayND::linearIndex: index out of range");
3275  idx += index[i]*strides_[i];
3276  }
3277  return idx;
3278  }
3279 
3280  template<typename Numeric, unsigned Len, unsigned Dim>
3282  const unsigned *index, const unsigned dim)
3283  {
3285  "Initialize npstat::ArrayND before calling method \"value\"");
3286  if (dim != dim_) throw npstat::NpstatInvalidArgument(
3287  "In npstat::ArrayND::value: incompatible index length");
3288  if (dim)
3289  {
3290  assert(index);
3291  unsigned long idx = 0UL;
3292  for (unsigned i=0; i<dim_; ++i)
3293  idx += index[i]*strides_[i];
3294  return data_[idx];
3295  }
3296  else
3297  return localData_[0];
3298  }
3299 
3300  template<typename Numeric, unsigned Len, unsigned Dim>
3301  inline const Numeric& ArrayND<Numeric,Len,Dim>::value(
3302  const unsigned *index, const unsigned dim) const
3303  {
3305  "Initialize npstat::ArrayND before calling method \"value\"");
3306  if (dim != dim_) throw npstat::NpstatInvalidArgument(
3307  "In npstat::ArrayND::value: incompatible index length");
3308  if (dim)
3309  {
3310  assert(index);
3311  unsigned long idx = 0UL;
3312  for (unsigned i=0; i<dim_; ++i)
3313  idx += index[i]*strides_[i];
3314  return data_[idx];
3315  }
3316  else
3317  return localData_[0];
3318  }
3319 
3320  template<typename Numeric, unsigned Len, unsigned Dim>
3322  const unsigned long index)
3323  {
3324  return data_[index];
3325  }
3326 
3327  template<typename Numeric, unsigned Len, unsigned Dim>
3329  const unsigned long index) const
3330  {
3331  return data_[index];
3332  }
3333 
3334  template<typename Numeric, unsigned Len, unsigned Dim>
3336  const unsigned long index)
3337  {
3338  if (index >= len_)
3340  "In npstat::ArrayND::linearValueAt: linear index out of range");
3341  return data_[index];
3342  }
3343 
3344  template<typename Numeric, unsigned Len, unsigned Dim>
3346  const unsigned long index) const
3347  {
3348  if (index >= len_)
3350  "In npstat::ArrayND::linearValueAt: linear index out of range");
3351  return data_[index];
3352  }
3353 
3354  template<typename Numeric, unsigned Len, unsigned Dim>
3356  const double x, const unsigned idim) const
3357  {
3358  if (x <= 0.0)
3359  return 0;
3360  else if (x >= static_cast<double>(shape_[idim] - 1))
3361  return shape_[idim] - 1;
3362  else
3363  return static_cast<unsigned>(std::floor(x + 0.5));
3364  }
3365 
3366  template<typename Numeric, unsigned Len, unsigned Dim>
3367  inline const Numeric& ArrayND<Numeric,Len,Dim>::closest(
3368  const double *x, const unsigned dim) const
3369  {
3371  "Initialize npstat::ArrayND before calling method \"closest\"");
3372  if (dim != dim_) throw npstat::NpstatInvalidArgument(
3373  "In npstat::ArrayND::closest: incompatible data length");
3374  if (dim)
3375  {
3376  assert(x);
3377  unsigned long idx = 0UL;
3378  for (unsigned i=0; i<dim_; ++i)
3379  idx += coordToIndex(x[i], i)*strides_[i];
3380  return data_[idx];
3381  }
3382  else
3383  return localData_[0];
3384  }
3385 
3386  template<typename Numeric, unsigned Len, unsigned Dim>
3388  const double *x, const unsigned dim)
3389  {
3391  "Initialize npstat::ArrayND before calling method \"closest\"");
3392  if (dim != dim_) throw npstat::NpstatInvalidArgument(
3393  "In npstat::ArrayND::closest: incompatible data length");
3394  if (dim)
3395  {
3396  assert(x);
3397  unsigned long idx = 0UL;
3398  for (unsigned i=0; i<dim_; ++i)
3399  idx += coordToIndex(x[i], i)*strides_[i];
3400  return data_[idx];
3401  }
3402  else
3403  return localData_[0];
3404  }
3405 
3406  template<typename Numeric, unsigned Len, unsigned Dim>
3407  inline const Numeric& ArrayND<Numeric,Len,Dim>::valueAt(
3408  const unsigned *index, const unsigned dim) const
3409  {
3411  "Initialize npstat::ArrayND before calling method \"valueAt\"");
3412  if (dim != dim_) throw npstat::NpstatInvalidArgument(
3413  "In npstat::ArrayND::valueAt: incompatible index length");
3414  if (dim)
3415  {
3416  assert(index);
3417  unsigned long idx = 0UL;
3418  for (unsigned i=0; i<dim_; ++i)
3419  {
3420  if (index[i] >= shape_[i]) throw npstat::NpstatOutOfRange(
3421  "In npstat::ArrayND::valueAt: index out of range");
3422  idx += index[i]*strides_[i];
3423  }
3424  return data_[idx];
3425  }
3426  else
3427  return localData_[0];
3428  }
3429 
3430  template<typename Numeric, unsigned Len, unsigned Dim>
3432  const unsigned *index, const unsigned dim)
3433  {
3435  "Initialize npstat::ArrayND before calling method \"valueAt\"");
3436  if (dim != dim_) throw npstat::NpstatInvalidArgument(
3437  "In npstat::ArrayND::valueAt: incompatible index length");
3438  if (dim)
3439  {
3440  assert(index);
3441  unsigned long idx = 0UL;
3442  for (unsigned i=0; i<dim_; ++i)
3443  {
3444  if (index[i] >= shape_[i]) throw npstat::NpstatOutOfRange(
3445  "In npstat::ArrayND::valueAt: index out of range");
3446  idx += index[i]*strides_[i];
3447  }
3448  return data_[idx];
3449  }
3450  else
3451  return localData_[0];
3452  }
3453 
3454  template<typename Numeric, unsigned Len, unsigned Dim>
3456  {
3458  "Initialize npstat::ArrayND before calling method \"operator()\"");
3460  "In npstat::ArrayND::operator(): wrong # of args (not rank 0 array)");
3461  return localData_[0];
3462  }
3463 
3464  template<typename Numeric, unsigned Len, unsigned Dim>
3465  inline const Numeric& ArrayND<Numeric,Len,Dim>::operator()() const
3466  {
3468  "Initialize npstat::ArrayND before calling method \"operator()\"");
3470  "In npstat::ArrayND::operator(): wrong # of args (not rank 0 array)");
3471  return localData_[0];
3472  }
3473 
3474  template<typename Numeric, unsigned Len, unsigned Dim>
3476  const unsigned i)
3477  {
3478  if (1U != dim_) throw npstat::NpstatInvalidArgument(
3479  "In npstat::ArrayND::operator(): wrong # of args (not rank 1 array)");
3480  return data_[i];
3481  }
3482 
3483  template<typename Numeric, unsigned Len, unsigned Dim>
3485  const unsigned i) const
3486  {
3487  if (1U != dim_) throw npstat::NpstatInvalidArgument(
3488  "In npstat::ArrayND::operator(): wrong # of args (not rank 1 array)");
3489  return data_[i];
3490  }
3491 
3492  template<typename Numeric, unsigned Len, unsigned Dim>
3493  const Numeric& ArrayND<Numeric,Len,Dim>::at() const
3494  {
3496  "Initialize npstat::ArrayND before calling method \"at\"");
3498  "In npstat::ArrayND::at: wrong # of args (not rank 0 array)");
3499  return localData_[0];
3500  }
3501 
3502  template<typename Numeric, unsigned Len, unsigned Dim>
3504  {
3506  "Initialize npstat::ArrayND before calling method \"at\"");
3508  "In npstat::ArrayND::at: wrong # of args (not rank 0 array)");
3509  return localData_[0];
3510  }
3511 
3512  template<typename Numeric, unsigned Len, unsigned Dim>
3514  const unsigned i0) const
3515  {
3516  if (1U != dim_) throw npstat::NpstatInvalidArgument(
3517  "In npstat::ArrayND::at: wrong # of args (not rank 1 array)");
3518  if (i0 >= shape_[0]) throw npstat::NpstatOutOfRange(
3519  "In npstat::ArrayND::at: index 0 out of range (rank 1)");
3520  return data_[i0];
3521  }
3522 
3523  template<typename Numeric, unsigned Len, unsigned Dim>
3525  const unsigned i0)
3526  {
3527  if (1U != dim_) throw npstat::NpstatInvalidArgument(
3528  "In npstat::ArrayND::at: wrong # of args (not rank 1 array)");
3529  if (i0 >= shape_[0]) throw npstat::NpstatOutOfRange(
3530  "In npstat::ArrayND::at: index 0 out of range (rank 1)");
3531  return data_[i0];
3532  }
3533 
3534  template<typename Numeric, unsigned Len, unsigned Dim>
3536  const unsigned i0,
3537  const unsigned i1)
3538  {
3539  if (2U != dim_) throw npstat::NpstatInvalidArgument(
3540  "In npstat::ArrayND::operator(): wrong # of args (not rank 2 array)");
3541  return data_[i0*strides_[0] + i1];
3542  }
3543 
3544  template<typename Numeric, unsigned Len, unsigned Dim>
3546  const unsigned i0,
3547  const unsigned i1) const
3548  {
3549  if (2U != dim_) throw npstat::NpstatInvalidArgument(
3550  "In npstat::ArrayND::operator(): wrong # of args (not rank 2 array)");
3551  return data_[i0*strides_[0] + i1];
3552  }
3553 
3554  template<typename Numeric, unsigned Len, unsigned Dim>
3556  const unsigned i0,
3557  const unsigned i1) const
3558  {
3559  if (2U != dim_) throw npstat::NpstatInvalidArgument(
3560  "In npstat::ArrayND::at: wrong # of args (not rank 2 array)");
3561  if (i0 >= shape_[0]) throw npstat::NpstatOutOfRange(
3562  "In npstat::ArrayND::at: index 0 out of range (rank 2)");
3563  if (i1 >= shape_[1]) throw npstat::NpstatOutOfRange(
3564  "In npstat::ArrayND::at: index 1 out of range (rank 2)");
3565  return data_[i0*strides_[0] + i1];
3566  }
3567 
3568  template<typename Numeric, unsigned Len, unsigned Dim>
3570  const unsigned i0,
3571  const unsigned i1)
3572  {
3573  if (2U != dim_) throw npstat::NpstatInvalidArgument(
3574  "In npstat::ArrayND::at: wrong # of args (not rank 2 array)");
3575  if (i0 >= shape_[0]) throw npstat::NpstatOutOfRange(
3576  "In npstat::ArrayND::at: index 0 out of range (rank 2)");
3577  if (i1 >= shape_[1]) throw npstat::NpstatOutOfRange(
3578  "In npstat::ArrayND::at: index 1 out of range (rank 2)");
3579  return data_[i0*strides_[0] + i1];
3580  }
3581 
3582  template<typename Numeric, unsigned Len, unsigned Dim>
3584  const unsigned i0,
3585  const unsigned i1,
3586  const unsigned i2) const
3587  {
3588  if (3U != dim_) throw npstat::NpstatInvalidArgument(
3589  "In npstat::ArrayND::operator(): wrong # of args (not rank 3 array)");
3590  return data_[i0*strides_[0] + i1*strides_[1] + i2];
3591  }
3592 
3593  template<typename Numeric, unsigned Len, unsigned Dim>
3595  const unsigned i0,
3596  const unsigned i1,
3597  const unsigned i2,
3598  const unsigned i3) const
3599  {
3600  if (4U != dim_) throw npstat::NpstatInvalidArgument(
3601  "In npstat::ArrayND::operator(): wrong # of args (not rank 4 array)");
3602  return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] + i3];
3603  }
3604 
3605  template<typename Numeric, unsigned Len, unsigned Dim>
3607  const unsigned i0,
3608  const unsigned i1,
3609  const unsigned i2,
3610  const unsigned i3,
3611  const unsigned i4) const
3612  {
3613  if (5U != dim_) throw npstat::NpstatInvalidArgument(
3614  "In npstat::ArrayND::operator(): wrong # of args (not rank 5 array)");
3615  return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] +
3616  i3*strides_[3] + i4];
3617  }
3618 
3619  template<typename Numeric, unsigned Len, unsigned Dim>
3621  const unsigned i0,
3622  const unsigned i1,
3623  const unsigned i2,
3624  const unsigned i3,
3625  const unsigned i4,
3626  const unsigned i5) const
3627  {
3628  if (6U != dim_) throw npstat::NpstatInvalidArgument(
3629  "In npstat::ArrayND::operator(): wrong # of args (not rank 6 array)");
3630  return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] +
3631  i3*strides_[3] + i4*strides_[4] + i5];
3632  }
3633 
3634  template<typename Numeric, unsigned Len, unsigned Dim>
3636  const unsigned i0,
3637  const unsigned i1,
3638  const unsigned i2,
3639  const unsigned i3,
3640  const unsigned i4,
3641  const unsigned i5,
3642  const unsigned i6) const
3643  {
3644  if (7U != dim_) throw npstat::NpstatInvalidArgument(
3645  "In npstat::ArrayND::operator(): wrong # of args (not rank 7 array)");
3646  return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] +
3647  i3*strides_[3] + i4*strides_[4] + i5*strides_[5] + i6];
3648  }
3649 
3650  template<typename Numeric, unsigned Len, unsigned Dim>
3652  const unsigned i0,
3653  const unsigned i1,
3654  const unsigned i2,
3655  const unsigned i3,
3656  const unsigned i4,
3657  const unsigned i5,
3658  const unsigned i6,
3659  const unsigned i7) const
3660  {
3661  if (8U != dim_) throw npstat::NpstatInvalidArgument(
3662  "In npstat::ArrayND::operator(): wrong # of args (not rank 8 array)");
3663  return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] +
3664  i3*strides_[3] + i4*strides_[4] + i5*strides_[5] +
3665  i6*strides_[6] + i7];
3666  }
3667 
3668  template<typename Numeric, unsigned Len, unsigned Dim>
3670  const unsigned i0,
3671  const unsigned i1,
3672  const unsigned i2,
3673  const unsigned i3,
3674  const unsigned i4,
3675  const unsigned i5,
3676  const unsigned i6,
3677  const unsigned i7,
3678  const unsigned i8) const
3679  {
3680  if (9U != dim_) throw npstat::NpstatInvalidArgument(
3681  "In npstat::ArrayND::operator(): wrong # of args (not rank 9 array)");
3682  return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] +
3683  i3*strides_[3] + i4*strides_[4] + i5*strides_[5] +
3684  i6*strides_[6] + i7*strides_[7] + i8];
3685  }
3686 
3687  template<typename Numeric, unsigned Len, unsigned Dim>
3689  const unsigned i0,
3690  const unsigned i1,
3691  const unsigned i2,
3692  const unsigned i3,
3693  const unsigned i4,
3694  const unsigned i5,
3695  const unsigned i6,
3696  const unsigned i7,
3697  const unsigned i8,
3698  const unsigned i9) const
3699  {
3700  if (10U != dim_) throw npstat::NpstatInvalidArgument(
3701  "In npstat::ArrayND::operator(): wrong # of args (not rank 10 array)");
3702  return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] +
3703  i3*strides_[3] + i4*strides_[4] + i5*strides_[5] +
3704  i6*strides_[6] + i7*strides_[7] + i8*strides_[8] + i9];
3705  }
3706 
3707  template<typename Numeric, unsigned Len, unsigned Dim>
3709  const unsigned i0,
3710  const unsigned i1,
3711  const unsigned i2)
3712  {
3713  if (3U != dim_) throw npstat::NpstatInvalidArgument(
3714  "In npstat::ArrayND::operator(): wrong # of args (not rank 3 array)");
3715  return data_[i0*strides_[0] + i1*strides_[1] + i2];
3716  }
3717 
3718  template<typename Numeric, unsigned Len, unsigned Dim>
3720  const unsigned i0,
3721  const unsigned i1,
3722  const unsigned i2,
3723  const unsigned i3)
3724  {
3725  if (4U != dim_) throw npstat::NpstatInvalidArgument(
3726  "In npstat::ArrayND::operator(): wrong # of args (not rank 4 array)");
3727  return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] + i3];
3728  }
3729 
3730  template<typename Numeric, unsigned Len, unsigned Dim>
3732  const unsigned i0,
3733  const unsigned i1,
3734  const unsigned i2,
3735  const unsigned i3,
3736  const unsigned i4)
3737  {
3738  if (5U != dim_) throw npstat::NpstatInvalidArgument(
3739  "In npstat::ArrayND::operator(): wrong # of args (not rank 5 array)");
3740  return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] +
3741  i3*strides_[3] + i4];
3742  }
3743 
3744  template<typename Numeric, unsigned Len, unsigned Dim>
3746  const unsigned i0,
3747  const unsigned i1,
3748  const unsigned i2,
3749  const unsigned i3,
3750  const unsigned i4,
3751  const unsigned i5)
3752  {
3753  if (6U != dim_) throw npstat::NpstatInvalidArgument(
3754  "In npstat::ArrayND::operator(): wrong # of args (not rank 6 array)");
3755  return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] +
3756  i3*strides_[3] + i4*strides_[4] + i5];
3757  }
3758 
3759  template<typename Numeric, unsigned Len, unsigned Dim>
3761  const unsigned i0,
3762  const unsigned i1,
3763  const unsigned i2,
3764  const unsigned i3,
3765  const unsigned i4,
3766  const unsigned i5,
3767  const unsigned i6)
3768  {
3769  if (7U != dim_) throw npstat::NpstatInvalidArgument(
3770  "In npstat::ArrayND::operator(): wrong # of args (not rank 7 array)");
3771  return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] +
3772  i3*strides_[3] + i4*strides_[4] + i5*strides_[5] + i6];
3773  }
3774 
3775  template<typename Numeric, unsigned Len, unsigned Dim>
3777  const unsigned i0,
3778  const unsigned i1,
3779  const unsigned i2,
3780  const unsigned i3,
3781  const unsigned i4,
3782  const unsigned i5,
3783  const unsigned i6,
3784  const unsigned i7)
3785  {
3786  if (8U != dim_) throw npstat::NpstatInvalidArgument(
3787  "In npstat::ArrayND::operator(): wrong # of args (not rank 8 array)");
3788  return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] +
3789  i3*strides_[3] + i4*strides_[4] + i5*strides_[5] +
3790  i6*strides_[6] + i7];
3791  }
3792 
3793  template<typename Numeric, unsigned Len, unsigned Dim>
3795  const unsigned i0,
3796  const unsigned i1,
3797  const unsigned i2,
3798  const unsigned i3,
3799  const unsigned i4,
3800  const unsigned i5,
3801  const unsigned i6,
3802  const unsigned i7,
3803  const unsigned i8)
3804  {
3805  if (9U != dim_) throw npstat::NpstatInvalidArgument(
3806  "In npstat::ArrayND::operator(): wrong # of args (not rank 9 array)");
3807  return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] +
3808  i3*strides_[3] + i4*strides_[4] + i5*strides_[5] +
3809  i6*strides_[6] + i7*strides_[7] + i8];
3810  }
3811 
3812  template<typename Numeric, unsigned Len, unsigned Dim>
3814  const unsigned i0,
3815  const unsigned i1,
3816  const unsigned i2,
3817  const unsigned i3,
3818  const unsigned i4,
3819  const unsigned i5,
3820  const unsigned i6,
3821  const unsigned i7,
3822  const unsigned i8,
3823  const unsigned i9)
3824  {
3825  if (10U != dim_) throw npstat::NpstatInvalidArgument(
3826  "In npstat::ArrayND::operator(): wrong # of args (not rank 10 array)");
3827  return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] +
3828  i3*strides_[3] + i4*strides_[4] + i5*strides_[5] +
3829  i6*strides_[6] + i7*strides_[7] + i8*strides_[8] + i9];
3830  }
3831 
3832  template<typename Numeric, unsigned Len, unsigned Dim>
3834  const unsigned i0,
3835  const unsigned i1,
3836  const unsigned i2) const
3837  {
3838  if (3U != dim_) throw npstat::NpstatInvalidArgument(
3839  "In npstat::ArrayND::at: wrong # of args (not rank 3 array)");
3840  if (i0 >= shape_[0]) throw npstat::NpstatOutOfRange(
3841  "In npstat::ArrayND::at: index 0 out of range (rank 3)");
3842  if (i1 >= shape_[1]) throw npstat::NpstatOutOfRange(
3843  "In npstat::ArrayND::at: index 1 out of range (rank 3)");
3844  if (i2 >= shape_[2]) throw npstat::NpstatOutOfRange(
3845  "In npstat::ArrayND::at: index 2 out of range (rank 3)");
3846  return data_[i0*strides_[0] + i1*strides_[1] + i2];
3847  }
3848 
3849  template<typename Numeric, unsigned Len, unsigned Dim>
3851  const unsigned i0,
3852  const unsigned i1,
3853  const unsigned i2,
3854  const unsigned i3) const
3855  {
3856  if (4U != dim_) throw npstat::NpstatInvalidArgument(
3857  "In npstat::ArrayND::at: wrong # of args (not rank 4 array)");
3858  if (i0 >= shape_[0]) throw npstat::NpstatOutOfRange(
3859  "In npstat::ArrayND::at: index 0 out of range (rank 4)");
3860  if (i1 >= shape_[1]) throw npstat::NpstatOutOfRange(
3861  "In npstat::ArrayND::at: index 1 out of range (rank 4)");
3862  if (i2 >= shape_[2]) throw npstat::NpstatOutOfRange(
3863  "In npstat::ArrayND::at: index 2 out of range (rank 4)");
3864  if (i3 >= shape_[3]) throw npstat::NpstatOutOfRange(
3865  "In npstat::ArrayND::at: index 3 out of range (rank 4)");
3866  return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] + i3];
3867  }
3868 
3869  template<typename Numeric, unsigned Len, unsigned Dim>
3871  const unsigned i0,
3872  const unsigned i1,
3873  const unsigned i2,
3874  const unsigned i3,
3875  const unsigned i4) const
3876  {
3877  if (5U != dim_) throw npstat::NpstatInvalidArgument(
3878  "In npstat::ArrayND::at: wrong # of args (not rank 5 array)");
3879  if (i0 >= shape_[0]) throw npstat::NpstatOutOfRange(
3880  "In npstat::ArrayND::at: index 0 out of range (rank 5)");
3881  if (i1 >= shape_[1]) throw npstat::NpstatOutOfRange(
3882  "In npstat::ArrayND::at: index 1 out of range (rank 5)");
3883  if (i2 >= shape_[2]) throw npstat::NpstatOutOfRange(
3884  "In npstat::ArrayND::at: index 2 out of range (rank 5)");
3885  if (i3 >= shape_[3]) throw npstat::NpstatOutOfRange(
3886  "In npstat::ArrayND::at: index 3 out of range (rank 5)");
3887  if (i4 >= shape_[4]) throw npstat::NpstatOutOfRange(
3888  "In npstat::ArrayND::at: index 4 out of range (rank 5)");
3889  return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] +
3890  i3*strides_[3] + i4];
3891  }
3892 
3893  template<typename Numeric, unsigned Len, unsigned Dim>
3895  const unsigned i0,
3896  const unsigned i1,
3897  const unsigned i2,
3898  const unsigned i3,
3899  const unsigned i4,
3900  const unsigned i5) const
3901  {
3902  if (6U != dim_) throw npstat::NpstatInvalidArgument(
3903  "In npstat::ArrayND::at: wrong # of args (not rank 6 array)");
3904  if (i0 >= shape_[0]) throw npstat::NpstatOutOfRange(
3905  "In npstat::ArrayND::at: index 0 out of range (rank 6)");
3906  if (i1 >= shape_[1]) throw npstat::NpstatOutOfRange(
3907  "In npstat::ArrayND::at: index 1 out of range (rank 6)");
3908  if (i2 >= shape_[2]) throw npstat::NpstatOutOfRange(
3909  "In npstat::ArrayND::at: index 2 out of range (rank 6)");
3910  if (i3 >= shape_[3]) throw npstat::NpstatOutOfRange(
3911  "In npstat::ArrayND::at: index 3 out of range (rank 6)");
3912  if (i4 >= shape_[4]) throw npstat::NpstatOutOfRange(
3913  "In npstat::ArrayND::at: index 4 out of range (rank 6)");
3914  if (i5 >= shape_[5]) throw npstat::NpstatOutOfRange(
3915  "In npstat::ArrayND::at: index 5 out of range (rank 6)");
3916  return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] +
3917  i3*strides_[3] + i4*strides_[4] + i5];
3918  }
3919 
3920  template<typename Numeric, unsigned Len, unsigned Dim>
3922  const unsigned i0,
3923  const unsigned i1,
3924  const unsigned i2,
3925  const unsigned i3,
3926  const unsigned i4,
3927  const unsigned i5,
3928  const unsigned i6) const
3929  {
3930  if (7U != dim_) throw npstat::NpstatInvalidArgument(
3931  "In npstat::ArrayND::at: wrong # of args (not rank 7 array)");
3932  if (i0 >= shape_[0]) throw npstat::NpstatOutOfRange(
3933  "In npstat::ArrayND::at: index 0 out of range (rank 7)");
3934  if (i1 >= shape_[1]) throw npstat::NpstatOutOfRange(
3935  "In npstat::ArrayND::at: index 1 out of range (rank 7)");
3936  if (i2 >= shape_[2]) throw npstat::NpstatOutOfRange(
3937  "In npstat::ArrayND::at: index 2 out of range (rank 7)");
3938  if (i3 >= shape_[3]) throw npstat::NpstatOutOfRange(
3939  "In npstat::ArrayND::at: index 3 out of range (rank 7)");
3940  if (i4 >= shape_[4]) throw npstat::NpstatOutOfRange(
3941  "In npstat::ArrayND::at: index 4 out of range (rank 7)");
3942  if (i5 >= shape_[5]) throw npstat::NpstatOutOfRange(
3943  "In npstat::ArrayND::at: index 5 out of range (rank 7)");
3944  if (i6 >= shape_[6]) throw npstat::NpstatOutOfRange(
3945  "In npstat::ArrayND::at: index 6 out of range (rank 7)");
3946  return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] +
3947  i3*strides_[3] + i4*strides_[4] + i5*strides_[5] + i6];
3948  }
3949 
3950  template<typename Numeric, unsigned Len, unsigned Dim>
3952  const unsigned i0,
3953  const unsigned i1,
3954  const unsigned i2,
3955  const unsigned i3,
3956  const unsigned i4,
3957  const unsigned i5,
3958  const unsigned i6,
3959  const unsigned i7) const
3960  {
3961  if (8U != dim_) throw npstat::NpstatInvalidArgument(
3962  "In npstat::ArrayND::at: wrong # of args (not rank 8 array)");
3963  if (i0 >= shape_[0]) throw npstat::NpstatOutOfRange(
3964  "In npstat::ArrayND::at: index 0 out of range (rank 8)");
3965  if (i1 >= shape_[1]) throw npstat::NpstatOutOfRange(
3966  "In npstat::ArrayND::at: index 1 out of range (rank 8)");
3967  if (i2 >= shape_[2]) throw npstat::NpstatOutOfRange(
3968  "In npstat::ArrayND::at: index 2 out of range (rank 8)");
3969  if (i3 >= shape_[3]) throw npstat::NpstatOutOfRange(
3970  "In npstat::ArrayND::at: index 3 out of range (rank 8)");
3971  if (i4 >= shape_[4]) throw npstat::NpstatOutOfRange(
3972  "In npstat::ArrayND::at: index 4 out of range (rank 8)");
3973  if (i5 >= shape_[5]) throw npstat::NpstatOutOfRange(
3974  "In npstat::ArrayND::at: index 5 out of range (rank 8)");
3975  if (i6 >= shape_[6]) throw npstat::NpstatOutOfRange(
3976  "In npstat::ArrayND::at: index 6 out of range (rank 8)");
3977  if (i7 >= shape_[7]) throw npstat::NpstatOutOfRange(
3978  "In npstat::ArrayND::at: index 7 out of range (rank 8)");
3979  return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] +
3980  i3*strides_[3] + i4*strides_[4] + i5*strides_[5] +
3981  i6*strides_[6] + i7];
3982  }
3983 
3984  template<typename Numeric, unsigned Len, unsigned Dim>
3986  const unsigned i0,
3987  const unsigned i1,
3988  const unsigned i2,
3989  const unsigned i3,
3990  const unsigned i4,
3991  const unsigned i5,
3992  const unsigned i6,
3993  const unsigned i7,
3994  const unsigned i8) const
3995  {
3996  if (9U != dim_) throw npstat::NpstatInvalidArgument(
3997  "In npstat::ArrayND::at: wrong # of args (not rank 9 array)");
3998  if (i0 >= shape_[0]) throw npstat::NpstatOutOfRange(
3999  "In npstat::ArrayND::at: index 0 out of range (rank 9)");
4000  if (i1 >= shape_[1]) throw npstat::NpstatOutOfRange(
4001  "In npstat::ArrayND::at: index 1 out of range (rank 9)");
4002  if (i2 >= shape_[2]) throw npstat::NpstatOutOfRange(
4003  "In npstat::ArrayND::at: index 2 out of range (rank 9)");
4004  if (i3 >= shape_[3]) throw npstat::NpstatOutOfRange(
4005  "In npstat::ArrayND::at: index 3 out of range (rank 9)");
4006  if (i4 >= shape_[4]) throw npstat::NpstatOutOfRange(
4007  "In npstat::ArrayND::at: index 4 out of range (rank 9)");
4008  if (i5 >= shape_[5]) throw npstat::NpstatOutOfRange(
4009  "In npstat::ArrayND::at: index 5 out of range (rank 9)");
4010  if (i6 >= shape_[6]) throw npstat::NpstatOutOfRange(
4011  "In npstat::ArrayND::at: index 6 out of range (rank 9)");
4012  if (i7 >= shape_[7]) throw npstat::NpstatOutOfRange(
4013  "In npstat::ArrayND::at: index 7 out of range (rank 9)");
4014  if (i8 >= shape_[8]) throw npstat::NpstatOutOfRange(
4015  "In npstat::ArrayND::at: index 8 out of range (rank 9)");
4016  return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] +
4017  i3*strides_[3] + i4*strides_[4] + i5*strides_[5] +
4018  i6*strides_[6] + i7*strides_[7] + i8];
4019  }
4020 
4021  template<typename Numeric, unsigned Len, unsigned Dim>
4023  const unsigned i0,
4024  const unsigned i1,
4025  const unsigned i2,
4026  const unsigned i3,
4027  const unsigned i4,
4028  const unsigned i5,
4029  const unsigned i6,
4030  const unsigned i7,
4031  const unsigned i8,
4032  const unsigned i9) const
4033  {
4034  if (10U != dim_) throw npstat::NpstatInvalidArgument(
4035  "In npstat::ArrayND::at: wrong # of args (not rank 10 array)");
4036  if (i0 >= shape_[0]) throw npstat::NpstatOutOfRange(
4037  "In npstat::ArrayND::at: index 0 out of range (rank 10)");
4038  if (i1 >= shape_[1]) throw npstat::NpstatOutOfRange(
4039  "In npstat::ArrayND::at: index 1 out of range (rank 10)");
4040  if (i2 >= shape_[2]) throw npstat::NpstatOutOfRange(
4041  "In npstat::ArrayND::at: index 2 out of range (rank 10)");
4042  if (i3 >= shape_[3]) throw npstat::NpstatOutOfRange(
4043  "In npstat::ArrayND::at: index 3 out of range (rank 10)");
4044  if (i4 >= shape_[4]) throw npstat::NpstatOutOfRange(
4045  "In npstat::ArrayND::at: index 4 out of range (rank 10)");
4046  if (i5 >= shape_[5]) throw npstat::NpstatOutOfRange(
4047  "In npstat::ArrayND::at: index 5 out of range (rank 10)");
4048  if (i6 >= shape_[6]) throw npstat::NpstatOutOfRange(
4049  "In npstat::ArrayND::at: index 6 out of range (rank 10)");
4050  if (i7 >= shape_[7]) throw npstat::NpstatOutOfRange(
4051  "In npstat::ArrayND::at: index 7 out of range (rank 10)");
4052  if (i8 >= shape_[8]) throw npstat::NpstatOutOfRange(
4053  "In npstat::ArrayND::at: index 8 out of range (rank 10)");
4054  if (i9 >= shape_[9]) throw npstat::NpstatOutOfRange(
4055  "In npstat::ArrayND::at: index 9 out of range (rank 10)");
4056  return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] +
4057  i3*strides_[3] + i4*strides_[4] + i5*strides_[5] +
4058  i6*strides_[6] + i7*strides_[7] + i8*strides_[8] + i9];
4059  }
4060 
4061  template<typename Numeric, unsigned Len, unsigned Dim>
4063  const unsigned i0,
4064  const unsigned i1,
4065  const unsigned i2)
4066  {
4067  if (3U != dim_) throw npstat::NpstatInvalidArgument(
4068  "In npstat::ArrayND::at: wrong # of args (not rank 3 array)");
4069  if (i0 >= shape_[0]) throw npstat::NpstatOutOfRange(
4070  "In npstat::ArrayND::at: index 0 out of range (rank 3)");
4071  if (i1 >= shape_[1]) throw npstat::NpstatOutOfRange(
4072  "In npstat::ArrayND::at: index 1 out of range (rank 3)");
4073  if (i2 >= shape_[2]) throw npstat::NpstatOutOfRange(
4074  "In npstat::ArrayND::at: index 2 out of range (rank 3)");
4075  return data_[i0*strides_[0] + i1*strides_[1] + i2];
4076  }
4077 
4078  template<typename Numeric, unsigned Len, unsigned Dim>
4080  const unsigned i0,
4081  const unsigned i1,
4082  const unsigned i2,
4083  const unsigned i3)
4084  {
4085  if (4U != dim_) throw npstat::NpstatInvalidArgument(
4086  "In npstat::ArrayND::at: wrong # of args (not rank 4 array)");
4087  if (i0 >= shape_[0]) throw npstat::NpstatOutOfRange(
4088  "In npstat::ArrayND::at: index 0 out of range (rank 4)");
4089  if (i1 >= shape_[1]) throw npstat::NpstatOutOfRange(
4090  "In npstat::ArrayND::at: index 1 out of range (rank 4)");
4091  if (i2 >= shape_[2]) throw npstat::NpstatOutOfRange(
4092  "In npstat::ArrayND::at: index 2 out of range (rank 4)");
4093  if (i3 >= shape_[3]) throw npstat::NpstatOutOfRange(
4094  "In npstat::ArrayND::at: index 3 out of range (rank 4)");
4095  return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] + i3];
4096  }
4097 
4098  template<typename Numeric, unsigned Len, unsigned Dim>
4100  const unsigned i0,
4101  const unsigned i1,
4102  const unsigned i2,
4103  const unsigned i3,
4104  const unsigned i4)
4105  {
4106  if (5U != dim_) throw npstat::NpstatInvalidArgument(
4107  "In npstat::ArrayND::at: wrong # of args (not rank 5 array)");
4108  if (i0 >= shape_[0]) throw npstat::NpstatOutOfRange(
4109  "In npstat::ArrayND::at: index 0 out of range (rank 5)");
4110  if (i1 >= shape_[1]) throw npstat::NpstatOutOfRange(
4111  "In npstat::ArrayND::at: index 1 out of range (rank 5)");
4112  if (i2 >= shape_[2]) throw npstat::NpstatOutOfRange(
4113  "In npstat::ArrayND::at: index 2 out of range (rank 5)");
4114  if (i3 >= shape_[3]) throw npstat::NpstatOutOfRange(
4115  "In npstat::ArrayND::at: index 3 out of range (rank 5)");
4116  if (i4 >= shape_[4]) throw npstat::NpstatOutOfRange(
4117  "In npstat::ArrayND::at: index 4 out of range (rank 5)");
4118  return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] +
4119  i3*strides_[3] + i4];
4120  }
4121 
4122  template<typename Numeric, unsigned Len, unsigned Dim>
4124  const unsigned i0,
4125  const unsigned i1,
4126  const unsigned i2,
4127  const unsigned i3,
4128  const unsigned i4,
4129  const unsigned i5)
4130  {
4131  if (6U != dim_) throw npstat::NpstatInvalidArgument(
4132  "In npstat::ArrayND::at: wrong # of args (not rank 6 array)");
4133  if (i0 >= shape_[0]) throw npstat::NpstatOutOfRange(
4134  "In npstat::ArrayND::at: index 0 out of range (rank 6)");
4135  if (i1 >= shape_[1]) throw npstat::NpstatOutOfRange(
4136  "In npstat::ArrayND::at: index 1 out of range (rank 6)");
4137  if (i2 >= shape_[2]) throw npstat::NpstatOutOfRange(
4138  "In npstat::ArrayND::at: index 2 out of range (rank 6)");
4139  if (i3 >= shape_[3]) throw npstat::NpstatOutOfRange(
4140  "In npstat::ArrayND::at: index 3 out of range (rank 6)");
4141  if (i4 >= shape_[4]) throw npstat::NpstatOutOfRange(
4142  "In npstat::ArrayND::at: index 4 out of range (rank 6)");
4143  if (i5 >= shape_[5]) throw npstat::NpstatOutOfRange(
4144  "In npstat::ArrayND::at: index 5 out of range (rank 6)");
4145  return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] +
4146  i3*strides_[3] + i4*strides_[4] + i5];
4147  }
4148 
4149  template<typename Numeric, unsigned Len, unsigned Dim>
4151  const unsigned i0,
4152  const unsigned i1,
4153  const unsigned i2,
4154  const unsigned i3,
4155  const unsigned i4,
4156  const unsigned i5,
4157  const unsigned i6)
4158  {
4159  if (7U != dim_) throw npstat::NpstatInvalidArgument(
4160  "In npstat::ArrayND::at: wrong # of args (not rank 7 array)");
4161  if (i0 >= shape_[0]) throw npstat::NpstatOutOfRange(
4162  "In npstat::ArrayND::at: index 0 out of range (rank 7)");
4163  if (i1 >= shape_[1]) throw npstat::NpstatOutOfRange(
4164  "In npstat::ArrayND::at: index 1 out of range (rank 7)");
4165  if (i2 >= shape_[2]) throw npstat::NpstatOutOfRange(
4166  "In npstat::ArrayND::at: index 2 out of range (rank 7)");
4167  if (i3 >= shape_[3]) throw npstat::NpstatOutOfRange(
4168  "In npstat::ArrayND::at: index 3 out of range (rank 7)");
4169  if (i4 >= shape_[4]) throw npstat::NpstatOutOfRange(
4170  "In npstat::ArrayND::at: index 4 out of range (rank 7)");
4171  if (i5 >= shape_[5]) throw npstat::NpstatOutOfRange(
4172  "In npstat::ArrayND::at: index 5 out of range (rank 7)");
4173  if (i6 >= shape_[6]) throw npstat::NpstatOutOfRange(
4174  "In npstat::ArrayND::at: index 6 out of range (rank 7)");
4175  return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] +
4176  i3*strides_[3] + i4*strides_[4] + i5*strides_[5] + i6];
4177  }
4178 
4179  template<typename Numeric, unsigned Len, unsigned Dim>
4181  const unsigned i0,
4182  const unsigned i1,
4183  const unsigned i2,
4184  const unsigned i3,
4185  const unsigned i4,
4186  const unsigned i5,
4187  const unsigned i6,
4188  const unsigned i7)
4189  {
4190  if (8U != dim_) throw npstat::NpstatInvalidArgument(
4191  "In npstat::ArrayND::at: wrong # of args (not rank 8 array)");
4192  if (i0 >= shape_[0]) throw npstat::NpstatOutOfRange(
4193  "In npstat::ArrayND::at: index 0 out of range (rank 8)");
4194  if (i1 >= shape_[1]) throw npstat::NpstatOutOfRange(
4195  "In npstat::ArrayND::at: index 1 out of range (rank 8)");
4196  if (i2 >= shape_[2]) throw npstat::NpstatOutOfRange(
4197  "In npstat::ArrayND::at: index 2 out of range (rank 8)");
4198  if (i3 >= shape_[3]) throw npstat::NpstatOutOfRange(
4199  "In npstat::ArrayND::at: index 3 out of range (rank 8)");
4200  if (i4 >= shape_[4]) throw npstat::NpstatOutOfRange(
4201  "In npstat::ArrayND::at: index 4 out of range (rank 8)");
4202  if (i5 >= shape_[5]) throw npstat::NpstatOutOfRange(
4203  "In npstat::ArrayND::at: index 5 out of range (rank 8)");
4204  if (i6 >= shape_[6]) throw npstat::NpstatOutOfRange(
4205  "In npstat::ArrayND::at: index 6 out of range (rank 8)");
4206  if (i7 >= shape_[7]) throw npstat::NpstatOutOfRange(
4207  "In npstat::ArrayND::at: index 7 out of range (rank 8)");
4208  return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] +
4209  i3*strides_[3] + i4*strides_[4] + i5*strides_[5] +
4210  i6*strides_[6] + i7];
4211  }
4212 
4213  template<typename Numeric, unsigned Len, unsigned Dim>
4215  const unsigned i0,
4216  const unsigned i1,
4217  const unsigned i2,
4218  const unsigned i3,
4219  const unsigned i4,
4220  const unsigned i5,
4221  const unsigned i6,
4222  const unsigned i7,
4223  const unsigned i8)
4224  {
4225  if (9U != dim_) throw npstat::NpstatInvalidArgument(
4226  "In npstat::ArrayND::at: wrong # of args (not rank 9 array)");
4227  if (i0 >= shape_[0]) throw npstat::NpstatOutOfRange(
4228  "In npstat::ArrayND::at: index 0 out of range (rank 9)");
4229  if (i1 >= shape_[1]) throw npstat::NpstatOutOfRange(
4230  "In npstat::ArrayND::at: index 1 out of range (rank 9)");
4231  if (i2 >= shape_[2]) throw npstat::NpstatOutOfRange(
4232  "In npstat::ArrayND::at: index 2 out of range (rank 9)");
4233  if (i3 >= shape_[3]) throw npstat::NpstatOutOfRange(
4234  "In npstat::ArrayND::at: index 3 out of range (rank 9)");
4235  if (i4 >= shape_[4]) throw npstat::NpstatOutOfRange(
4236  "In npstat::ArrayND::at: index 4 out of range (rank 9)");
4237  if (i5 >= shape_[5]) throw npstat::NpstatOutOfRange(
4238  "In npstat::ArrayND::at: index 5 out of range (rank 9)");
4239  if (i6 >= shape_[6]) throw npstat::NpstatOutOfRange(
4240  "In npstat::ArrayND::at: index 6 out of range (rank 9)");
4241  if (i7 >= shape_[7]) throw npstat::NpstatOutOfRange(
4242  "In npstat::ArrayND::at: index 7 out of range (rank 9)");
4243  if (i8 >= shape_[8]) throw npstat::NpstatOutOfRange(
4244  "In npstat::ArrayND::at: index 8 out of range (rank 9)");
4245  return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] +
4246  i3*strides_[3] + i4*strides_[4] + i5*strides_[5] +
4247  i6*strides_[6] + i7*strides_[7] + i8];
4248  }
4249 
4250  template<typename Numeric, unsigned Len, unsigned Dim>
4252  const unsigned i0,
4253  const unsigned i1,
4254  const unsigned i2,
4255  const unsigned i3,
4256  const unsigned i4,
4257  const unsigned i5,
4258  const unsigned i6,
4259  const unsigned i7,
4260  const unsigned i8,
4261  const unsigned i9)
4262  {
4263  if (10U != dim_) throw npstat::NpstatInvalidArgument(
4264  "In npstat::ArrayND::at: wrong # of args (not rank 10 array)");
4265  if (i0 >= shape_[0]) throw npstat::NpstatOutOfRange(
4266  "In npstat::ArrayND::at: index 0 out of range (rank 10)");
4267  if (i1 >= shape_[1]) throw npstat::NpstatOutOfRange(
4268  "In npstat::ArrayND::at: index 1 out of range (rank 10)");
4269  if (i2 >= shape_[2]) throw npstat::NpstatOutOfRange(
4270  "In npstat::ArrayND::at: index 2 out of range (rank 10)");
4271  if (i3 >= shape_[3]) throw npstat::NpstatOutOfRange(
4272  "In npstat::ArrayND::at: index 3 out of range (rank 10)");
4273  if (i4 >= shape_[4]) throw npstat::NpstatOutOfRange(
4274  "In npstat::ArrayND::at: index 4 out of range (rank 10)");
4275  if (i5 >= shape_[5]) throw npstat::NpstatOutOfRange(
4276  "In npstat::ArrayND::at: index 5 out of range (rank 10)");
4277  if (i6 >= shape_[6]) throw npstat::NpstatOutOfRange(
4278  "In npstat::ArrayND::at: index 6 out of range (rank 10)");
4279  if (i7 >= shape_[7]) throw npstat::NpstatOutOfRange(
4280  "In npstat::ArrayND::at: index 7 out of range (rank 10)");
4281  if (i8 >= shape_[8]) throw npstat::NpstatOutOfRange(
4282  "In npstat::ArrayND::at: index 8 out of range (rank 10)");
4283  if (i9 >= shape_[9]) throw npstat::NpstatOutOfRange(
4284  "In npstat::ArrayND::at: index 9 out of range (rank 10)");
4285  return data_[i0*strides_[0] + i1*strides_[1] + i2*strides_[2] +
4286  i3*strides_[3] + i4*strides_[4] + i5*strides_[5] +
4287  i6*strides_[6] + i7*strides_[7] + i8*strides_[8] + i9];
4288  }
4289 
4290  template<typename Numeric, unsigned Len, unsigned Dim>
4291  template<unsigned Len2, unsigned Dim2>
4293  const ArrayND<Numeric,Len2,Dim2>& r) const
4294  {
4296  "In npstat::ArrayND::maxAbsDifference: "
4297  "incompatible argument array shape");
4298  if (dim_)
4299  {
4300  double maxd = 0.0;
4301  for (unsigned long i=0; i<len_; ++i)
4302  {
4303  const Numeric rval = r.data_[i];
4304  const double d = absDifference(data_[i], rval);
4305  if (d > maxd)
4306  maxd = d;
4307  }
4308  return maxd;
4309  }
4310  else
4311  {
4312  const Numeric rval = r.localData_[0];
4313  return absDifference(localData_[0], rval);
4314  }
4315  }
4316 
4317  template<typename Numeric, unsigned Len, unsigned Dim>
4318  template<unsigned Len2, unsigned Dim2>
4320  const ArrayND<Numeric,Len2,Dim2>& r) const
4321  {
4322  if (shapeIsKnown_ != r.shapeIsKnown_)
4323  return false;
4324  if (r.dim_ != dim_)
4325  return false;
4326  if (r.len_ != len_)
4327  return false;
4328  for (unsigned i=0; i<dim_; ++i)
4329  if (shape_[i] != r.shape_[i])
4330  return false;
4331  for (unsigned i=0; i<dim_; ++i)
4332  assert(strides_[i] == r.strides_[i]);
4333  for (unsigned long j=0; j<len_; ++j)
4334  if (data_[j] != r.data_[j])
4335  return false;
4336  return true;
4337  }
4338 
4339  template<typename Numeric, unsigned Len, unsigned Dim>
4340  template<unsigned Len2, unsigned Dim2>
4342  const ArrayND<Numeric,Len2,Dim2>& r) const
4343  {
4344  return !(*this == r);
4345  }
4346 
4347  template<typename Numeric, unsigned Len, unsigned Dim>
4348  template<typename Num2>
4351  {
4353  "Initialize npstat::ArrayND before calling method \"operator*\"");
4355  for (unsigned long i=0; i<len_; ++i)
4356  result.data_[i] = data_[i]*r;
4357  return result;
4358  }
4359 
4360  template<typename Numeric, unsigned Len, unsigned Dim>
4361  template<typename Num2>
4364  {
4366  "Initialize npstat::ArrayND before calling method \"operator/\"");
4367  if (r == Num2()) throw npstat::NpstatRuntimeError(
4368  "In npstat::ArrayND::operator/: division by zero");
4370  for (unsigned long i=0; i<len_; ++i)
4371  result.data_[i] = data_[i]/r;
4372  return result;
4373  }
4374 
4375  template<typename Numeric, unsigned Len, unsigned Dim>
4376  template <unsigned Len2, unsigned Dim2>
4379  const ArrayND<Numeric,Len2,Dim2>& r) const
4380  {
4382  "In npstat::ArrayND::operator+: "
4383  "incompatible argument array shape");
4385  for (unsigned long i=0; i<len_; ++i)
4386  result.data_[i] = data_[i] + r.data_[i];
4387  return result;
4388  }
4389 
4390  template<typename Numeric, unsigned Len, unsigned Dim>
4391  template <unsigned Len2, unsigned Dim2>
4394  const ArrayND<Numeric,Len2,Dim2>& r) const
4395  {
4397  "In npstat::ArrayND::operator-: "
4398  "incompatible argument array shape");
4400  for (unsigned long i=0; i<len_; ++i)
4401  result.data_[i] = data_[i] - r.data_[i];
4402  return result;
4403  }
4404 
4405  template<typename Numeric, unsigned Len, unsigned Dim>
4407  {
4409  "Initialize npstat::ArrayND before calling method \"operator+\"");
4410  return *this;
4411  }
4412 
4413  template<typename Numeric, unsigned Len, unsigned Dim>
4415  {
4417  "Initialize npstat::ArrayND before calling method \"operator-\"");
4419  for (unsigned long i=0; i<len_; ++i)
4420  result.data_[i] = -data_[i];
4421  return result;
4422  }
4423 
4424  template<typename Numeric, unsigned Len, unsigned Dim>
4425  template <typename Num2>
4428  {
4430  "Initialize npstat::ArrayND before calling method \"operator*=\"");
4431  for (unsigned long i=0; i<len_; ++i)
4432  data_[i] *= r;
4433  return *this;
4434  }
4435 
4436  template<typename Numeric, unsigned Len, unsigned Dim>
4439  {
4441  "Initialize npstat::ArrayND before calling method \"makeNonNegative\"");
4442  const Numeric zero = Numeric();
4443  if (dim_)
4444  {
4445  for (unsigned long i=0; i<len_; ++i)
4446  if (!(ComplexComparesAbs<Numeric>::more(data_[i], zero)))
4447  data_[i] = zero;
4448  }
4449  else
4451  localData_[0] = zero;
4452  return *this;
4453  }
4454 
4455  template<typename Numeric, unsigned Len, unsigned Dim>
4457  const double tolerance, const unsigned nCycles)
4458  {
4460  "Initialize npstat::ArrayND before calling method \"makeCopulaSteps\"");
4461  if (nCycles == 0U)
4462  return 0U;
4463  if (!dim_) throw npstat::NpstatInvalidArgument(
4464  "npstat::ArrayND::makeCopulaSteps method "
4465  "can not be used with array of 0 rank");
4466 
4467  const Numeric zero = Numeric();
4468  for (unsigned long i=0; i<len_; ++i)
4469  if (!(ComplexComparesAbs<Numeric>::more(data_[i], zero)))
4470  data_[i] = zero;
4471 
4472  std::vector<Numeric*> axesPtrBuf(dim_);
4473  Numeric** axes = &axesPtrBuf[0];
4474  const Numeric one = static_cast<Numeric>(1);
4475 
4476  // Memory for the axis accumulators
4477  unsigned idxSum = 0;
4478  for (unsigned i=0; i<dim_; ++i)
4479  idxSum += shape_[i];
4480  std::vector<Numeric> axesBuf(idxSum);
4481  axes[0] = &axesBuf[0];
4482  for (unsigned i=1; i<dim_; ++i)
4483  axes[i] = axes[i-1] + shape_[i-1];
4484 
4485  // Accumulate axis projections
4486  unsigned icycle = 0;
4487  for (; icycle<nCycles; ++icycle)
4488  {
4489  for (unsigned i=0; i<idxSum; ++i)
4490  axesBuf[i] = zero;
4491 
4492  // Accumulate sums for each axis
4493  for (unsigned long idat=0; idat<len_; ++idat)
4494  {
4495  unsigned long l = idat;
4496  for (unsigned i=0; i<dim_; ++i)
4497  {
4498  const unsigned idx = l / strides_[i];
4499  l -= (idx * strides_[i]);
4500  axes[i][idx] += data_[idat];
4501  }
4502  }
4503 
4504  // Make averages out of sums
4505  bool withinTolerance = true;
4506  Numeric totalSum = zero;
4507  for (unsigned i=0; i<dim_; ++i)
4508  {
4509  Numeric axisSum = zero;
4510  const unsigned amax = shape_[i];
4511  for (unsigned a=0; a<amax; ++a)
4512  {
4513  if (axes[i][a] == zero)
4515  "In npstat::ArrayND::makeCopulaSteps: "
4516  "marginal density is zero");
4517  axisSum += axes[i][a];
4518  }
4519  totalSum += axisSum;
4520  const Numeric axisAverage = axisSum/static_cast<Numeric>(amax);
4521  for (unsigned a=0; a<amax; ++a)
4522  axes[i][a] /= axisAverage;
4523  for (unsigned a=0; a<amax && withinTolerance; ++a)
4524  {
4525  const double adelta = absDifference(axes[i][a], one);
4526  if (adelta > tolerance)
4527  withinTolerance = false;
4528  }
4529  }
4530 
4531  if (withinTolerance)
4532  break;
4533 
4534  const Numeric totalAverage = totalSum/
4535  static_cast<Numeric>(len_)/static_cast<Numeric>(dim_);
4536 
4537  // Run over all points again and divide by
4538  // the product of marginals
4539  for (unsigned long idat=0; idat<len_; ++idat)
4540  {
4541  unsigned long l = idat;
4542  for (unsigned i=0; i<dim_; ++i)
4543  {
4544  const unsigned idx = l / strides_[i];
4545  l -= (idx * strides_[i]);
4546  data_[idat] /= axes[i][idx];
4547  }
4548  data_[idat] /= totalAverage;
4549  }
4550  }
4551 
4552  return icycle;
4553  }
4554 
4555  template<typename Numeric, unsigned Len, unsigned Dim>
4556  template <typename Num2>
4559  {
4561  "Initialize npstat::ArrayND before calling method \"operator/=\"");
4562  if (r == Num2()) throw npstat::NpstatRuntimeError(
4563  "In npstat::ArrayND::operator/=: division by zero");
4564  for (unsigned long i=0; i<len_; ++i)
4565  data_[i] /= r;
4566  return *this;
4567  }
4568 
4569  template<typename Numeric, unsigned Len, unsigned Dim>
4570  template <typename Num2, unsigned Len2, unsigned Dim2>
4573  {
4575  "In npstat::ArrayND::operator+=: "
4576  "incompatible argument array shape");
4577  for (unsigned long i=0; i<len_; ++i)
4578  data_[i] += r.data_[i];
4579  return *this;
4580  }
4581 
4582  template<typename Numeric, unsigned Len, unsigned Dim>
4583  template <typename Num3, typename Num2, unsigned Len2, unsigned Dim2>
4586  const Num3& c)
4587  {
4589  "In npstat::ArrayND::addmul: "
4590  "incompatible argument array shape");
4591  for (unsigned long i=0; i<len_; ++i)
4592  data_[i] += r.data_[i]*c;
4593  return *this;
4594  }
4595 
4596  template<typename Numeric, unsigned Len, unsigned Dim>
4597  template <typename Num2, unsigned Len2, unsigned Dim2>
4600  {
4602  "In npstat::ArrayND::operator-=: "
4603  "incompatible argument array shape");
4604  for (unsigned long i=0; i<len_; ++i)
4605  data_[i] -= r.data_[i];
4606  return *this;
4607  }
4608 
4609  template<typename Numeric, unsigned Len, unsigned Dim>
4611  const double *coords, const unsigned dim) const
4612  {
4614  "Initialize npstat::ArrayND before calling method \"interpolate1\"");
4615  if (dim != dim_) throw npstat::NpstatInvalidArgument(
4616  "In npstat::ArrayND::interpolate1: incompatible coordinate length");
4617  if (dim)
4618  {
4619  const unsigned maxdim = CHAR_BIT*sizeof(unsigned long);
4620  if (dim_ >= maxdim)
4622  "In npstat::ArrayND::interpolate1: array rank is too large");
4623 
4624  double dx[maxdim];
4625  unsigned ix[maxdim];
4626  for (unsigned i=0; i<dim; ++i)
4627  {
4628  const double x = coords[i];
4629  if (x <= 0.0)
4630  {
4631  ix[i] = 0;
4632  dx[i] = 0.0;
4633  }
4634  else if (x >= static_cast<double>(shape_[i] - 1))
4635  {
4636  ix[i] = shape_[i] - 1;
4637  dx[i] = 0.0;
4638  }
4639  else
4640  {
4641  ix[i] = static_cast<unsigned>(std::floor(x));
4642  dx[i] = x - ix[i];
4643  }
4644  }
4645 
4646  Numeric sum = Numeric();
4647  const unsigned long maxcycle = 1UL << dim;
4648  for (unsigned long icycle=0UL; icycle<maxcycle; ++icycle)
4649  {
4650  double w = 1.0;
4651  unsigned long icell = 0UL;
4652  for (unsigned i=0; i<dim; ++i)
4653  {
4654  if (icycle & (1UL << i))
4655  {
4656  w *= dx[i];
4657  icell += strides_[i]*(ix[i] + 1U);
4658  }
4659  else
4660  {
4661  w *= (1.0 - dx[i]);
4662  icell += strides_[i]*ix[i];
4663  }
4664  }
4665  if (w > 0.0)
4666  sum += data_[icell]*static_cast<proper_double>(w);
4667  }
4668  return sum;
4669  }
4670  else
4671  return localData_[0];
4672  }
4673 
4674  template<typename Numeric, unsigned Len, unsigned Dim>
4676  const unsigned level, const double* coords, const Numeric* base) const
4677  {
4678  const unsigned npoints = shape_[level];
4679  const double x = coords[level];
4680 
4681  unsigned ix, npt = 1;
4682  double dx = 0.0;
4683  if (x < 0.0)
4684  ix = 0;
4685  else if (x > static_cast<double>(npoints - 1))
4686  ix = npoints - 1;
4687  else
4688  {
4689  ix = static_cast<unsigned>(std::floor(x));
4690  if (ix) --ix;
4691  unsigned imax = ix + 3;
4692  while (imax >= npoints)
4693  {
4694  if (ix) --ix;
4695  --imax;
4696  }
4697  dx = x - ix;
4698  npt = imax + 1 - ix;
4699  }
4700  assert(npt >= 1 && npt <= 4);
4701 
4702  Numeric fit[4];
4703  if (level < dim_ - 1)
4704  for (unsigned ipt=0; ipt<npt; ++ipt)
4705  fit[ipt] = interpolateLoop(level + 1, coords,
4706  base + (ix + ipt)*strides_[level]);
4707 
4708  const Numeric* const v = (level == dim_ - 1 ? base + ix : fit);
4709  switch (npt)
4710  {
4711  case 1:
4712  return v[0];
4713  case 2:
4714  return interpolate_linear(dx, v[0], v[1]);
4715  case 3:
4716  return interpolate_quadratic(dx, v[0], v[1], v[2]);
4717  case 4:
4718  return interpolate_cubic(dx, v[0], v[1], v[2], v[3]);
4719  default:
4720  assert(0);
4721  return Numeric();
4722  }
4723  }
4724 
4725  template<typename Numeric, unsigned Len, unsigned Dim>
4727  const double* coords, const unsigned dim) const
4728  {
4730  "Initialize npstat::ArrayND before calling method \"interpolate3\"");
4731  if (dim != dim_) throw npstat::NpstatInvalidArgument(
4732  "In npstat::ArrayND::interpolate3: incompatible coordinate length");
4733  if (dim)
4734  {
4735  assert(coords);
4736  return interpolateLoop(0, coords, data_);
4737  }
4738  else
4739  return localData_[0];
4740  }
4741 
4742  template<typename Numeric, unsigned Len, unsigned Dim>
4743  template<class Functor>
4745  {
4747  "Initialize npstat::ArrayND before calling method \"apply\"");
4748  for (unsigned long i=0; i<len_; ++i)
4749  data_[i] = static_cast<Numeric>(f(data_[i]));
4750  return *this;
4751  }
4752 
4753  template<typename Numeric, unsigned Len, unsigned Dim>
4754  template<class Functor>
4756  Functor f)
4757  {
4759  "Initialize npstat::ArrayND before calling method \"scanInPlace\"");
4760  for (unsigned long i=0; i<len_; ++i)
4761  f(data_[i]);
4762  return *this;