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 
18 #include "JetMETCorrections/InterpolationTables/interface/SimpleFunctors.h"
19 #include "JetMETCorrections/InterpolationTables/interface/ArrayRange.h"
20 #include "JetMETCorrections/InterpolationTables/interface/AbsArrayProjector.h"
21 #include "JetMETCorrections/InterpolationTables/interface/AbsVisitor.h"
22 #include "JetMETCorrections/InterpolationTables/interface/PreciseType.h"
23 #include "JetMETCorrections/InterpolationTables/interface/ProperDblFromCmpl.h"
24 
25 namespace npstat {
47  template <typename Numeric, unsigned StackLen = 1U, unsigned StackDim = 10U>
48  class ArrayND {
49  template <typename Num2, unsigned Len2, unsigned Dim2>
50  friend class ArrayND;
51 
52  public:
53  typedef Numeric value_type;
54  typedef typename ProperDblFromCmpl<Numeric>::type proper_double;
55 
70  ArrayND();
71 
73 
81  explicit ArrayND(const ArrayShape& shape);
82  ArrayND(const unsigned* shape, unsigned dim);
84 
86  ArrayND(const ArrayND&);
87 
96  template <typename Num2, unsigned Len2, unsigned Dim2>
97  ArrayND(const ArrayND<Num2, Len2, Dim2>&);
98 
103  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
104  ArrayND(const ArrayND<Num2, Len2, Dim2>&, Functor f);
105 
107  template <typename Num2, unsigned Len2, unsigned Dim2>
108  ArrayND(const ArrayND<Num2, Len2, Dim2>& from, const ArrayRange& fromRange);
109 
111  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
112  ArrayND(const ArrayND<Num2, Len2, Dim2>& from, const ArrayRange& fromRange, Functor f);
113 
125  template <typename Num2, unsigned Len2, unsigned Dim2>
126  ArrayND(const ArrayND<Num2, Len2, Dim2>& slicedArray, const unsigned* indices, unsigned nIndices);
127 
129  template <typename Num1, unsigned Len1, unsigned Dim1, typename Num2, unsigned Len2, unsigned Dim2>
130  ArrayND(const ArrayND<Num1, Len1, Dim1>& a1, const ArrayND<Num2, Len2, Dim2>& a2);
131 
133 
137  explicit ArrayND(unsigned n0);
138  ArrayND(unsigned n0, unsigned n1);
139  ArrayND(unsigned n0, unsigned n1, unsigned n2);
140  ArrayND(unsigned n0, unsigned n1, unsigned n2, unsigned n3);
141  ArrayND(unsigned n0, unsigned n1, unsigned n2, unsigned n3, unsigned n4);
142  ArrayND(unsigned n0, unsigned n1, unsigned n2, unsigned n3, unsigned n4, unsigned n5);
143  ArrayND(unsigned n0, unsigned n1, unsigned n2, unsigned n3, unsigned n4, unsigned n5, unsigned n6);
144  ArrayND(unsigned n0, unsigned n1, unsigned n2, unsigned n3, unsigned n4, unsigned n5, unsigned n6, unsigned n7);
145  ArrayND(unsigned n0,
146  unsigned n1,
147  unsigned n2,
148  unsigned n3,
149  unsigned n4,
150  unsigned n5,
151  unsigned n6,
152  unsigned n7,
153  unsigned n8);
154  ArrayND(unsigned n0,
155  unsigned n1,
156  unsigned n2,
157  unsigned n3,
158  unsigned n4,
159  unsigned n5,
160  unsigned n6,
161  unsigned n7,
162  unsigned n8,
163  unsigned n9);
165 
167  ~ArrayND();
168 
177  ArrayND& operator=(const ArrayND&);
178 
180  template <typename Num2, unsigned Len2, unsigned Dim2>
181  ArrayND& operator=(const ArrayND<Num2, Len2, Dim2>&);
182 
184  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
185  ArrayND& assign(const ArrayND<Num2, Len2, Dim2>&, Functor f);
186 
192  ArrayND& uninitialize();
193 
195 
200  Numeric& value(const unsigned* index, unsigned indexLen);
201  const Numeric& value(const unsigned* index, unsigned indexLen) const;
203 
205 
209  Numeric& valueAt(const unsigned* index, unsigned indexLen);
210  const Numeric& valueAt(const unsigned* index, unsigned indexLen) const;
212 
214 
215  Numeric& linearValue(unsigned long index);
216  const Numeric& linearValue(unsigned long index) const;
218 
220 
221  Numeric& linearValueAt(unsigned long index);
222  const Numeric& linearValueAt(unsigned long index) const;
224 
226  void convertLinearIndex(unsigned long l, unsigned* index, unsigned indexLen) const;
227 
229  unsigned long linearIndex(const unsigned* idx, unsigned idxLen) const;
230 
231  // Some inspectors
233  inline unsigned long length() const { return len_; }
234 
236  inline const Numeric* data() const { return data_; }
237 
239  inline bool isShapeKnown() const { return shapeIsKnown_; }
240 
242  inline unsigned rank() const { return dim_; }
243 
245  ArrayShape shape() const;
246 
248  inline const unsigned* shapeData() const { return shape_; }
249 
251  ArrayRange fullRange() const;
252 
254  unsigned span(unsigned dim) const;
255 
257  unsigned maximumSpan() const;
258 
260  unsigned minimumSpan() const;
261 
263  inline const unsigned long* strides() const { return strides_; }
264 
266  bool isZero() const;
267 
273  bool isDensity() const;
274 
276  template <typename Num2>
277  ArrayND& setData(const Num2* data, unsigned long dataLength);
278 
280  template <unsigned Len2, unsigned Dim2>
281  bool operator==(const ArrayND<Numeric, Len2, Dim2>&) const;
282 
284  template <unsigned Len2, unsigned Dim2>
285  bool operator!=(const ArrayND<Numeric, Len2, Dim2>&) const;
286 
288  template <unsigned Len2, unsigned Dim2>
289  double maxAbsDifference(const ArrayND<Numeric, Len2, Dim2>&) const;
290 
292  ArrayND operator+() const;
293 
295  ArrayND operator-() const;
296 
298  template <unsigned Len2, unsigned Dim2>
299  ArrayND operator+(const ArrayND<Numeric, Len2, Dim2>& r) const;
300 
302  template <unsigned Len2, unsigned Dim2>
303  ArrayND operator-(const ArrayND<Numeric, Len2, Dim2>& r) const;
304 
306  template <typename Num2>
307  ArrayND operator*(const Num2& r) const;
308 
310  template <typename Num2>
311  ArrayND operator/(const Num2& r) const;
312 
314 
318  template <typename Num2>
319  ArrayND& operator*=(const Num2& r);
320 
321  template <typename Num2>
322  ArrayND& operator/=(const Num2& r);
323 
324  template <typename Num2, unsigned Len2, unsigned Dim2>
325  ArrayND& operator+=(const ArrayND<Num2, Len2, Dim2>& r);
326 
327  template <typename Num2, unsigned Len2, unsigned Dim2>
328  ArrayND& operator-=(const ArrayND<Num2, Len2, Dim2>& r);
330 
332  template <typename Num3, typename Num2, unsigned Len2, unsigned Dim2>
333  ArrayND& addmul(const ArrayND<Num2, Len2, Dim2>& r, const Num3& c);
334 
336  template <typename Num2, unsigned Len2, unsigned Dim2>
337  ArrayND outer(const ArrayND<Num2, Len2, Dim2>& r) const;
338 
343  ArrayND contract(unsigned pos1, unsigned pos2) const;
344 
350  template <typename Num2, unsigned Len2, unsigned Dim2>
351  ArrayND dot(const ArrayND<Num2, Len2, Dim2>& r) const;
352 
367  template <typename Num2, unsigned Len2, unsigned Dim2>
368  ArrayND marginalize(const ArrayND<Num2, Len2, Dim2>& prior, const unsigned* indexMap, unsigned mapLen) const;
369 
371  ArrayND transpose(unsigned pos1, unsigned pos2) const;
372 
374  ArrayND transpose() const;
375 
382  template <typename Num2>
383  Num2 sum() const;
384 
389  template <typename Num2>
390  Num2 sumsq() const;
391 
400  template <typename Num2>
401  ArrayND derivative(double scale = 1.0) const;
402 
407  template <typename Num2>
408  ArrayND cdfArray(double scale = 1.0) const;
409 
414  template <typename Num2>
415  Num2 cdfValue(const unsigned* index, unsigned indexLen) const;
416 
428  template <typename Num2>
429  void convertToLastDimCdf(ArrayND* sumSlice, bool useTrapezoids);
430 
432  Numeric min() const;
433 
435  Numeric min(unsigned* index, unsigned indexLen) const;
436 
438  Numeric max() const;
439 
441  Numeric max(unsigned* index, unsigned indexLen) const;
442 
444 
452  Numeric& closest(const double* x, unsigned xDim);
453  const Numeric& closest(const double* x, unsigned xDim) const;
455 
464  Numeric interpolate1(const double* x, unsigned xDim) const;
465 
472  Numeric interpolate3(const double* x, unsigned xDim) const;
473 
483  template <class Functor>
484  ArrayND& apply(Functor f);
485 
493  template <class Functor>
494  ArrayND& scanInPlace(Functor f);
495 
497  ArrayND& constFill(Numeric c);
498 
500  ArrayND& clear();
501 
508  ArrayND& linearFill(const double* coeff, unsigned coeffLen, double c);
509 
516  template <class Functor>
517  ArrayND& functorFill(Functor f);
518 
525  ArrayND& makeUnit();
526 
528  ArrayND& makeNonNegative();
529 
539  unsigned makeCopulaSteps(double tolerance, unsigned maxIterations);
540 
545  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
546  void jointScan(ArrayND<Num2, Len2, Dim2>& other, Functor binaryFunct);
547 
550  template <typename Num2, unsigned Len2, unsigned Dim2>
551  inline ArrayND& inPlaceMul(const ArrayND<Num2, Len2, Dim2>& r) {
552  jointScan(const_cast<ArrayND<Num2, Len2, Dim2>&>(r), multeq_left<Numeric, Num2>());
553  return *this;
554  }
555 
573  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
574  void jointSubrangeScan(ArrayND<Num2, Len2, Dim2>& other,
575  const unsigned* thisCorner,
576  const unsigned* range,
577  const unsigned* otherCorner,
578  unsigned arrLen,
579  Functor binaryFunct);
580 
586  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
587  void dualCircularScan(ArrayND<Num2, Len2, Dim2>& other,
588  const unsigned* thisCorner,
589  const unsigned* range,
590  const unsigned* otherCorner,
591  unsigned arrLen,
592  Functor binaryFunct);
593 
598  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
599  void flatCircularScan(ArrayND<Num2, Len2, Dim2>& other,
600  const unsigned* thisCorner,
601  const unsigned* range,
602  const unsigned* otherCorner,
603  unsigned arrLen,
604  Functor binaryFunct);
605 
610  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
611  void circularFlatScan(ArrayND<Num2, Len2, Dim2>& other,
612  const unsigned* thisCorner,
613  const unsigned* range,
614  const unsigned* otherCorner,
615  unsigned arrLen,
616  Functor binaryFunct);
617 
624  template <typename Num2, typename Integer>
625  void processSubrange(AbsArrayProjector<Numeric, Num2>& f, const BoxND<Integer>& subrange) const;
626 
634  template <typename Num2, unsigned Len2, unsigned Dim2>
635  void exportSubrange(const unsigned* fromCorner, unsigned lenCorner, ArrayND<Num2, Len2, Dim2>* dest) const;
636 
638  template <typename Num2, unsigned Len2, unsigned Dim2>
639  void importSubrange(const unsigned* fromCorner, unsigned lenCorner, const ArrayND<Num2, Len2, Dim2>& from);
640 
646  template <typename Num2, unsigned Len2, unsigned Dim2>
647  bool isClose(const ArrayND<Num2, Len2, Dim2>& r, double eps) const;
648 
650  bool isCompatible(const ArrayShape& shape) const;
651 
656  template <typename Num2, unsigned Len2, unsigned Dim2>
657  bool isShapeCompatible(const ArrayND<Num2, Len2, Dim2>& r) const;
658 
668  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
669  void jointSliceScan(ArrayND<Num2, Len2, Dim2>& slice,
670  const unsigned* fixedIndices,
671  const unsigned* fixedIndexValues,
672  unsigned nFixedIndices,
673  Functor binaryFunct);
674 
683  template <typename Num2, class Functor>
684  void jointMemSliceScan(Num2* buffer,
685  unsigned long bufLen,
686  const unsigned* fixedIndices,
687  const unsigned* fixedIndexValues,
688  unsigned nFixedIndices,
689  Functor binaryFunct);
690 
692  ArrayShape sliceShape(const unsigned* fixedIndices, unsigned nFixedIndices) const;
693 
695  template <typename Num2, unsigned Len2, unsigned Dim2>
696  inline void exportSlice(ArrayND<Num2, Len2, Dim2>* slice,
697  const unsigned* fixedIndices,
698  const unsigned* fixedIndexValues,
699  unsigned nFixedIndices) const {
700  assert(slice);
701  (const_cast<ArrayND*>(this))
702  ->jointSliceScan(*slice, fixedIndices, fixedIndexValues, nFixedIndices, scast_assign_right<Numeric, Num2>());
703  }
704 
709  template <typename Num2>
710  inline void exportMemSlice(Num2* buffer,
711  unsigned long bufLen,
712  const unsigned* fixedIndices,
713  const unsigned* fixedIndexValues,
714  unsigned nFixedIndices) const {
715  (const_cast<ArrayND*>(this))
716  ->jointMemSliceScan(
717  buffer, bufLen, fixedIndices, fixedIndexValues, nFixedIndices, scast_assign_right<Numeric, Num2>());
718  }
719 
721  template <typename Num2, unsigned Len2, unsigned Dim2>
722  inline void importSlice(const ArrayND<Num2, Len2, Dim2>& slice,
723  const unsigned* fixedIndices,
724  const unsigned* fixedIndexValues,
725  unsigned nFixedIndices) {
726  jointSliceScan(const_cast<ArrayND<Num2, Len2, Dim2>&>(slice),
727  fixedIndices,
728  fixedIndexValues,
729  nFixedIndices,
730  scast_assign_left<Numeric, Num2>());
731  }
732 
737  template <typename Num2>
738  inline void importMemSlice(const Num2* buffer,
739  unsigned long bufLen,
740  const unsigned* fixedIndices,
741  const unsigned* fixedIndexValues,
742  unsigned nFixedIndices) {
743  jointMemSliceScan(const_cast<Num2*>(buffer),
744  bufLen,
745  fixedIndices,
746  fixedIndexValues,
747  nFixedIndices,
748  scast_assign_left<Numeric, Num2>());
749  }
750 
757  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
758  void applySlice(ArrayND<Num2, Len2, Dim2>& slice,
759  const unsigned* fixedIndices,
760  unsigned nFixedIndices,
761  Functor binaryFunct);
762 
767  template <typename Num2, unsigned Len2, unsigned Dim2>
768  inline ArrayND& multiplyBySlice(const ArrayND<Num2, Len2, Dim2>& slice,
769  const unsigned* fixedIndices,
770  unsigned nFixedIndices) {
771  applySlice(
772  const_cast<ArrayND<Num2, Len2, Dim2>&>(slice), fixedIndices, nFixedIndices, multeq_left<Numeric, Num2>());
773  return *this;
774  }
775 
777 
784  template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
785  void project(ArrayND<Num2, Len2, Dim2>* projection,
786  AbsArrayProjector<Numeric, Num3>& projector,
787  const unsigned* projectedIndices,
788  unsigned nProjectedIndices) const;
789 
790  template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
791  void project(ArrayND<Num2, Len2, Dim2>* projection,
792  AbsVisitor<Numeric, Num3>& projector,
793  const unsigned* projectedIndices,
794  unsigned nProjectedIndices) const;
796 
798 
803  template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
804  void addToProjection(ArrayND<Num2, Len2, Dim2>* projection,
805  AbsArrayProjector<Numeric, Num3>& projector,
806  const unsigned* projectedIndices,
807  unsigned nProjectedIndices) const;
808 
809  template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
810  void subtractFromProjection(ArrayND<Num2, Len2, Dim2>* projection,
811  AbsArrayProjector<Numeric, Num3>& projector,
812  const unsigned* projectedIndices,
813  unsigned nProjectedIndices) const;
814 
815  template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
816  void addToProjection(ArrayND<Num2, Len2, Dim2>* projection,
817  AbsVisitor<Numeric, Num3>& projector,
818  const unsigned* projectedIndices,
819  unsigned nProjectedIndices) const;
820 
821  template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
822  void subtractFromProjection(ArrayND<Num2, Len2, Dim2>* projection,
823  AbsVisitor<Numeric, Num3>& projector,
824  const unsigned* projectedIndices,
825  unsigned nProjectedIndices) const;
827 
833  template <typename Num2, unsigned Len2, unsigned Dim2>
834  void rotate(const unsigned* shifts, unsigned lenShifts, ArrayND<Num2, Len2, Dim2>* rotated) const;
835 
841  template <typename Num2, unsigned Len2, unsigned Dim2>
842  void multiMirror(ArrayND<Num2, Len2, Dim2>* out) const;
843 
845 
849  Numeric& operator()();
850  const Numeric& operator()() const;
851 
852  Numeric& operator()(unsigned i0);
853  const Numeric& operator()(unsigned i0) const;
854 
855  Numeric& operator()(unsigned i0, unsigned i1);
856  const Numeric& operator()(unsigned i0, unsigned i1) const;
857 
858  Numeric& operator()(unsigned i0, unsigned i1, unsigned i2);
859  const Numeric& operator()(unsigned i0, unsigned i1, unsigned i2) const;
860 
861  Numeric& operator()(unsigned i0, unsigned i1, unsigned i2, unsigned i3);
862  const Numeric& operator()(unsigned i0, unsigned i1, unsigned i2, unsigned i3) const;
863 
864  Numeric& operator()(unsigned i0, unsigned i1, unsigned i2, unsigned i3, unsigned i4);
865  const Numeric& operator()(unsigned i0, unsigned i1, unsigned i2, unsigned i3, unsigned i4) const;
866 
867  Numeric& operator()(unsigned i0, unsigned i1, unsigned i2, unsigned i3, unsigned i4, unsigned i5);
868  const Numeric& operator()(unsigned i0, unsigned i1, unsigned i2, unsigned i3, unsigned i4, unsigned i5) const;
869 
870  Numeric& operator()(unsigned i0, unsigned i1, unsigned i2, unsigned i3, unsigned i4, unsigned i5, unsigned i6);
871  const Numeric& operator()(
872  unsigned i0, unsigned i1, unsigned i2, unsigned i3, unsigned i4, unsigned i5, unsigned i6) const;
873 
874  Numeric& operator()(
875  unsigned i0, unsigned i1, unsigned i2, unsigned i3, unsigned i4, unsigned i5, unsigned i6, unsigned i7);
876  const Numeric& operator()(
877  unsigned i0, unsigned i1, unsigned i2, unsigned i3, unsigned i4, unsigned i5, unsigned i6, unsigned i7) const;
878 
879  Numeric& operator()(unsigned i0,
880  unsigned i1,
881  unsigned i2,
882  unsigned i3,
883  unsigned i4,
884  unsigned i5,
885  unsigned i6,
886  unsigned i7,
887  unsigned i8);
888  const Numeric& operator()(unsigned i0,
889  unsigned i1,
890  unsigned i2,
891  unsigned i3,
892  unsigned i4,
893  unsigned i5,
894  unsigned i6,
895  unsigned i7,
896  unsigned i8) const;
897 
898  Numeric& operator()(unsigned i0,
899  unsigned i1,
900  unsigned i2,
901  unsigned i3,
902  unsigned i4,
903  unsigned i5,
904  unsigned i6,
905  unsigned i7,
906  unsigned i8,
907  unsigned i9);
908  const Numeric& operator()(unsigned i0,
909  unsigned i1,
910  unsigned i2,
911  unsigned i3,
912  unsigned i4,
913  unsigned i5,
914  unsigned i6,
915  unsigned i7,
916  unsigned i8,
917  unsigned i9) const;
919 
921 
925  Numeric& at();
926  const Numeric& at() const;
927 
928  Numeric& at(unsigned i0);
929  const Numeric& at(unsigned i0) const;
930 
931  Numeric& at(unsigned i0, unsigned i1);
932  const Numeric& at(unsigned i0, unsigned i1) const;
933 
934  Numeric& at(unsigned i0, unsigned i1, unsigned i2);
935  const Numeric& at(unsigned i0, unsigned i1, unsigned i2) const;
936 
937  Numeric& at(unsigned i0, unsigned i1, unsigned i2, unsigned i3);
938  const Numeric& at(unsigned i0, unsigned i1, unsigned i2, unsigned i3) const;
939 
940  Numeric& at(unsigned i0, unsigned i1, unsigned i2, unsigned i3, unsigned i4);
941  const Numeric& at(unsigned i0, unsigned i1, unsigned i2, unsigned i3, unsigned i4) const;
942 
943  Numeric& at(unsigned i0, unsigned i1, unsigned i2, unsigned i3, unsigned i4, unsigned i5);
944  const Numeric& at(unsigned i0, unsigned i1, unsigned i2, unsigned i3, unsigned i4, unsigned i5) const;
945 
946  Numeric& at(unsigned i0, unsigned i1, unsigned i2, unsigned i3, unsigned i4, unsigned i5, unsigned i6);
947  const Numeric& at(unsigned i0, unsigned i1, unsigned i2, unsigned i3, unsigned i4, unsigned i5, unsigned i6) const;
948 
949  Numeric& at(unsigned i0, unsigned i1, unsigned i2, unsigned i3, unsigned i4, unsigned i5, unsigned i6, unsigned i7);
950  const Numeric& at(
951  unsigned i0, unsigned i1, unsigned i2, unsigned i3, unsigned i4, unsigned i5, unsigned i6, unsigned i7) const;
952 
953  Numeric& at(unsigned i0,
954  unsigned i1,
955  unsigned i2,
956  unsigned i3,
957  unsigned i4,
958  unsigned i5,
959  unsigned i6,
960  unsigned i7,
961  unsigned i8);
962  const Numeric& at(unsigned i0,
963  unsigned i1,
964  unsigned i2,
965  unsigned i3,
966  unsigned i4,
967  unsigned i5,
968  unsigned i6,
969  unsigned i7,
970  unsigned i8) const;
971 
972  Numeric& at(unsigned i0,
973  unsigned i1,
974  unsigned i2,
975  unsigned i3,
976  unsigned i4,
977  unsigned i5,
978  unsigned i6,
979  unsigned i7,
980  unsigned i8,
981  unsigned i9);
982  const Numeric& at(unsigned i0,
983  unsigned i1,
984  unsigned i2,
985  unsigned i3,
986  unsigned i4,
987  unsigned i5,
988  unsigned i6,
989  unsigned i7,
990  unsigned i8,
991  unsigned i9) const;
993 
995 
999  Numeric& cl();
1000  const Numeric& cl() const;
1001 
1002  Numeric& cl(double x0);
1003  const Numeric& cl(double x0) const;
1004 
1005  Numeric& cl(double x0, double x1);
1006  const Numeric& cl(double x0, double x1) const;
1007 
1008  Numeric& cl(double x0, double x1, double x2);
1009  const Numeric& cl(double x0, double x1, double x2) const;
1010 
1011  Numeric& cl(double x0, double x1, double x2, double x3);
1012  const Numeric& cl(double x0, double x1, double x2, double x3) const;
1013 
1014  Numeric& cl(double x0, double x1, double x2, double x3, double x4);
1015  const Numeric& cl(double x0, double x1, double x2, double x3, double x4) const;
1016 
1017  Numeric& cl(double x0, double x1, double x2, double x3, double x4, double x5);
1018  const Numeric& cl(double x0, double x1, double x2, double x3, double x4, double x5) const;
1019 
1020  Numeric& cl(double x0, double x1, double x2, double x3, double x4, double x5, double x6);
1021  const Numeric& cl(double x0, double x1, double x2, double x3, double x4, double x5, double x6) const;
1022 
1023  Numeric& cl(double x0, double x1, double x2, double x3, double x4, double x5, double x6, double x7);
1024  const Numeric& cl(double x0, double x1, double x2, double x3, double x4, double x5, double x6, double x7) const;
1025 
1026  Numeric& cl(double x0, double x1, double x2, double x3, double x4, double x5, double x6, double x7, double x8);
1027  const Numeric& cl(
1028  double x0, double x1, double x2, double x3, double x4, double x5, double x6, double x7, double x8) const;
1029 
1030  Numeric& cl(
1031  double x0, double x1, double x2, double x3, double x4, double x5, double x6, double x7, double x8, double x9);
1032  const Numeric& cl(
1033  double x0, double x1, double x2, double x3, double x4, double x5, double x6, double x7, double x8, double x9)
1034  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, ArrayND* array);
1046 
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 
1060  bool shapeIsKnown_;
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, unsigned long idx, double shift, const double* coeffs);
1070 
1071  // Recursive implementation of nested loops for "functorFill"
1072  template <class Functor>
1073  void functorFillLoop(unsigned level, unsigned long idx, Functor f, unsigned* farg);
1074 
1075  // Recursive implementation of nested loops for "interpolate3"
1076  Numeric interpolateLoop(unsigned level, const double* x, const Numeric* base) const;
1077 
1078  // Recursive implementation of nested loops for the outer product
1079  template <typename Num1, unsigned Len1, unsigned Dim1, typename Num2, unsigned Len2, unsigned Dim2>
1080  void outerProductLoop(unsigned level,
1081  unsigned long idx0,
1082  unsigned long idx1,
1083  unsigned long idx2,
1084  const ArrayND<Num1, Len1, Dim1>& a1,
1085  const ArrayND<Num2, Len2, Dim2>& a2);
1086 
1087  // Recursive implementation of nested loops for contraction
1088  void contractLoop(unsigned thisLevel,
1089  unsigned resLevel,
1090  unsigned pos1,
1091  unsigned pos2,
1092  unsigned long idxThis,
1093  unsigned long idxRes,
1094  ArrayND& result) const;
1095 
1096  // Recursive implementation of nested loops for transposition
1097  void transposeLoop(unsigned level,
1098  unsigned pos1,
1099  unsigned pos2,
1100  unsigned long idxThis,
1101  unsigned long idxRes,
1102  ArrayND& result) const;
1103 
1104  // Recursive implementation of nested loops for the dot product
1105  template <typename Num2, unsigned Len2, unsigned Dim2>
1106  void dotProductLoop(unsigned level,
1107  unsigned long idx0,
1108  unsigned long idx1,
1109  unsigned long idx2,
1110  const ArrayND<Num2, Len2, Dim2>& r,
1111  ArrayND& result) const;
1112 
1113  // Recursive implementation of nested loops for marginalization
1114  template <typename Num2, unsigned Len2, unsigned Dim2>
1115  Numeric marginalizeInnerLoop(unsigned long idx,
1116  unsigned levelPr,
1117  unsigned long idxPr,
1118  const ArrayND<Num2, Len2, Dim2>& prior,
1119  const unsigned* indexMap) const;
1120  template <typename Num2, unsigned Len2, unsigned Dim2>
1121  void marginalizeLoop(unsigned level,
1122  unsigned long idx,
1123  unsigned levelRes,
1124  unsigned long idxRes,
1125  const ArrayND<Num2, Len2, Dim2>& prior,
1126  const unsigned* indexMap,
1127  ArrayND& res) const;
1128 
1129  // Recursive implementation of nested loops for range copy
1130  // with functor modification of elements
1131  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
1132  void copyRangeLoopFunct(unsigned level,
1133  unsigned long idx0,
1134  unsigned long idx1,
1135  const ArrayND<Num2, Len2, Dim2>& r,
1136  const ArrayRange& range,
1137  Functor f);
1138 
1139  // Loop over subrange in such a way that the functor is called
1140  // only if indices on both sides are valid. The topology of both
1141  // arrays is that of the hyperplane (flat).
1142  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
1143  void commonSubrangeLoop(unsigned level,
1144  unsigned long idx0,
1145  unsigned long idx1,
1146  const unsigned* thisCorner,
1147  const unsigned* range,
1148  const unsigned* otherCorner,
1149  ArrayND<Num2, Len2, Dim2>& other,
1150  Functor binaryFunct);
1151 
1152  // Similar loop with the topology of the hypertorus for both
1153  // arrays (all indices of both arrays are wrapped around)
1154  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
1155  void dualCircularLoop(unsigned level,
1156  unsigned long idx0,
1157  unsigned long idx1,
1158  const unsigned* thisCorner,
1159  const unsigned* range,
1160  const unsigned* otherCorner,
1161  ArrayND<Num2, Len2, Dim2>& other,
1162  Functor binaryFunct);
1163 
1164  // Similar loop in which the topology of this array is assumed
1165  // to be flat and the topology of the destination array is that
1166  // of the hypertorus. Due to the asymmetry of the topologies,
1167  // "const" function is not provided (use const_cast as appropriate).
1168  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
1169  void flatCircularLoop(unsigned level,
1170  unsigned long idx0,
1171  unsigned long idx1,
1172  const unsigned* thisCorner,
1173  const unsigned* range,
1174  const unsigned* otherCorner,
1175  ArrayND<Num2, Len2, Dim2>& other,
1176  Functor binaryFunct);
1177 
1178  // Similar loop in which the topology of this array is assumed
1179  // to be hypertoroidal and the topology of the destination array
1180  // is flat.
1181  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
1182  void circularFlatLoop(unsigned level,
1183  unsigned long idx0,
1184  unsigned long idx1,
1185  const unsigned* thisCorner,
1186  const unsigned* range,
1187  const unsigned* otherCorner,
1188  ArrayND<Num2, Len2, Dim2>& other,
1189  Functor binaryFunct);
1190 
1191  // Slice compatibility verification
1192  template <typename Num2, unsigned Len2, unsigned Dim2>
1193  unsigned long verifySliceCompatibility(const ArrayND<Num2, Len2, Dim2>& slice,
1194  const unsigned* fixedIndices,
1195  const unsigned* fixedIndexValues,
1196  unsigned nFixedIndices) const;
1197 
1198  // Buffer compatibility verification with a slice
1199  unsigned long verifyBufferSliceCompatibility(unsigned long bufLen,
1200  const unsigned* fixedIndices,
1201  const unsigned* fixedIndexValues,
1202  unsigned nFixedIndices,
1203  unsigned long* sliceStrides) const;
1204 
1205  // Recursive implementation of nested loops for slice operations
1206  template <typename Num2, class Functor>
1207  void jointSliceLoop(unsigned level,
1208  unsigned long idx0,
1209  unsigned level1,
1210  unsigned long idx1,
1211  Num2* sliceData,
1212  const unsigned long* sliceStrides,
1213  const unsigned* fixedIndices,
1214  const unsigned* fixedIndexValues,
1215  unsigned nFixedIndices,
1216  Functor binaryFunctor);
1217 
1218  // Recursive implementation of nested loops for "applySlice"
1219  template <typename Num2, class Functor>
1220  void scaleBySliceInnerLoop(unsigned level,
1221  unsigned long idx0,
1222  Num2& scale,
1223  const unsigned* projectedIndices,
1224  unsigned nProjectedIndices,
1225  Functor binaryFunct);
1226 
1227  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
1228  void scaleBySliceLoop(unsigned level,
1229  unsigned long idx0,
1230  unsigned level1,
1231  unsigned long idx1,
1232  ArrayND<Num2, Len2, Dim2>& slice,
1233  const unsigned* fixedIndices,
1234  unsigned nFixedIndices,
1235  Functor binaryFunct);
1236 
1237  // Recursive implementation of nested loops for projections
1238  template <typename Num2>
1239  void projectInnerLoop(unsigned level,
1240  unsigned long idx0,
1241  unsigned* currentIndex,
1242  AbsArrayProjector<Numeric, Num2>& projector,
1243  const unsigned* projectedIndices,
1244  unsigned nProjectedIndices) const;
1245 
1246  template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3, class Op>
1247  void projectLoop(unsigned level,
1248  unsigned long idx0,
1249  unsigned level1,
1250  unsigned long idx1,
1251  unsigned* currentIndex,
1252  ArrayND<Num2, Len2, Dim2>* projection,
1253  AbsArrayProjector<Numeric, Num3>& projector,
1254  const unsigned* projectedIndices,
1255  unsigned nProjectedIndices,
1256  Op fcn) const;
1257 
1258  // Note that "projectLoop2" is almost identical to "projectLoop"
1259  // while "projectInnerLoop2" is almost identical to "projectInnerLoop".
1260  // It would make a lot of sense to combine these functions into
1261  // the same code and then partially specialize the little piece
1262  // where the "AbsVisitor" or "AbsArrayProjector" is actually called.
1263  // Unfortunately, "AbsVisitor" and "AbsArrayProjector" are
1264  // templates themselves, and it is not possible in C++ to partially
1265  // specialize a function template (that is, even if we can specialize
1266  // on "AbsVisitor" vs. "AbsArrayProjector", there is no way to
1267  // specialize on their parameter types).
1268  template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3, class Op>
1269  void projectLoop2(unsigned level,
1270  unsigned long idx0,
1271  unsigned level1,
1272  unsigned long idx1,
1273  ArrayND<Num2, Len2, Dim2>* projection,
1274  AbsVisitor<Numeric, Num3>& projector,
1275  const unsigned* projectedIndices,
1276  unsigned nProjectedIndices,
1277  Op fcn) const;
1278 
1279  template <typename Num2>
1280  void projectInnerLoop2(unsigned level,
1281  unsigned long idx0,
1282  AbsVisitor<Numeric, Num2>& projector,
1283  const unsigned* projectedIndices,
1284  unsigned nProjectedIndices) const;
1285 
1286  template <typename Num2, typename Integer>
1287  void processSubrangeLoop(unsigned level,
1288  unsigned long idx0,
1289  unsigned* currentIndex,
1290  AbsArrayProjector<Numeric, Num2>& f,
1291  const BoxND<Integer>& subrange) const;
1292 
1293  // Sum of all points with the given index and below
1294  template <typename Accumulator>
1295  Accumulator sumBelowLoop(unsigned level, unsigned long idx0, const unsigned* limit) const;
1296 
1297  // Loop for "convertToLastDimCdf"
1298  template <typename Accumulator>
1299  void convertToLastDimCdfLoop(
1300  ArrayND* sumSlice, unsigned level, unsigned long idx0, unsigned long idxSlice, bool useTrapezoids);
1301 
1302  // Convert a coordinate into index.
1303  // No checking whether "idim" is within limits.
1304  unsigned coordToIndex(double coord, unsigned idim) const;
1305 
1306  // Verify that projection array is compatible with this one
1307  template <typename Num2, unsigned Len2, unsigned Dim2>
1308  void verifyProjectionCompatibility(const ArrayND<Num2, Len2, Dim2>& projection,
1309  const unsigned* projectedIndices,
1310  unsigned nProjectedIndices) const;
1311  };
1312 } // namespace npstat
1313 
1314 #include <cmath>
1315 #include <climits>
1316 #include <algorithm>
1317 #include <sstream>
1318 #include "JetMETCorrections/InterpolationTables/interface/NpstatException.h"
1319 
1320 #include "Alignment/Geners/interface/GenericIO.hh"
1321 #include "Alignment/Geners/interface/IOIsUnsigned.hh"
1322 
1323 #include "JetMETCorrections/InterpolationTables/interface/allocators.h"
1324 
1325 #include "JetMETCorrections/InterpolationTables/interface/interpolate.h"
1326 #include "JetMETCorrections/InterpolationTables/interface/absDifference.h"
1327 #include "JetMETCorrections/InterpolationTables/interface/ComplexComparesFalse.h"
1328 #include "JetMETCorrections/InterpolationTables/interface/ComplexComparesAbs.h"
1329 
1330 #define me_macro_check_loop_prerequisites(METHOD, INNERLOOP) \
1331  template <typename Numeric, unsigned Len, unsigned Dim> \
1332  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor> \
1333  void ArrayND<Numeric, Len, Dim>::METHOD(ArrayND<Num2, Len2, Dim2>& other, \
1334  const unsigned* thisCorner, \
1335  const unsigned* range, \
1336  const unsigned* otherCorner, \
1337  const unsigned arrLen, \
1338  Functor binaryFunct) { \
1339  if (!shapeIsKnown_) \
1340  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"" #METHOD "\""); \
1341  if (!other.shapeIsKnown_) \
1342  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::" #METHOD ": uninitialized argument array"); \
1343  if (dim_ != other.dim_) \
1344  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::" #METHOD ": incompatible argument array rank"); \
1345  if (arrLen != dim_) \
1346  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::" #METHOD ": incompatible index length"); \
1347  if (dim_) { \
1348  assert(thisCorner); \
1349  assert(range); \
1350  assert(otherCorner); \
1351  INNERLOOP(0U, 0UL, 0UL, thisCorner, range, otherCorner, other, binaryFunct); \
1352  } else \
1353  binaryFunct(localData_[0], other.localData_[0]); \
1354  }
1355 
1356 namespace npstat {
1357  template <typename Numeric, unsigned Len, unsigned Dim>
1358  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
1359  void ArrayND<Numeric, Len, Dim>::commonSubrangeLoop(unsigned level,
1360  unsigned long idx0,
1361  unsigned long idx1,
1362  const unsigned* thisCorner,
1363  const unsigned* range,
1364  const unsigned* otherCorner,
1365  ArrayND<Num2, Len2, Dim2>& r,
1366  Functor binaryFunct) {
1367  const unsigned imax = range[level];
1368 
1369  if (level == dim_ - 1) {
1370  Numeric* left = data_ + (idx0 + thisCorner[level]);
1371  Numeric* const lMax = data_ + (idx0 + shape_[level]);
1372  Num2* right = r.data_ + (idx1 + otherCorner[level]);
1373  Num2* const rMax = r.data_ + (idx1 + r.shape_[level]);
1374 
1375  for (unsigned i = 0; i < imax && left < lMax && right < rMax; ++i)
1376  binaryFunct(*left++, *right++);
1377  } else {
1378  const unsigned long leftStride = strides_[level];
1379  const unsigned long leftMax = idx0 + shape_[level] * leftStride;
1380  idx0 += thisCorner[level] * leftStride;
1381  const unsigned long rightStride = r.strides_[level];
1382  const unsigned long rightMax = idx1 + r.shape_[level] * rightStride;
1383  idx1 += otherCorner[level] * rightStride;
1384 
1385  for (unsigned i = 0; i < imax && idx0 < leftMax && idx1 < rightMax; ++i, idx0 += leftStride, idx1 += rightStride)
1386  commonSubrangeLoop(level + 1, idx0, idx1, thisCorner, range, otherCorner, r, binaryFunct);
1387  }
1388  }
1389 
1390  template <typename Numeric, unsigned Len, unsigned Dim>
1391  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
1392  void ArrayND<Numeric, Len, Dim>::dualCircularLoop(unsigned level,
1393  unsigned long idx0,
1394  unsigned long idx1,
1395  const unsigned* thisCorner,
1396  const unsigned* range,
1397  const unsigned* otherCorner,
1398  ArrayND<Num2, Len2, Dim2>& r,
1399  Functor binaryFunct) {
1400  const unsigned imax = range[level];
1401  const unsigned leftShift = thisCorner[level];
1402  const unsigned leftPeriod = shape_[level];
1403  const unsigned rightShift = otherCorner[level];
1404  const unsigned rightPeriod = r.shape_[level];
1405 
1406  if (level == dim_ - 1) {
1407  Numeric* left = data_ + idx0;
1408  Num2* right = r.data_ + idx1;
1409  for (unsigned i = 0; i < imax; ++i)
1410  binaryFunct(left[(i + leftShift) % leftPeriod], right[(i + rightShift) % rightPeriod]);
1411  } else {
1412  const unsigned long leftStride = strides_[level];
1413  const unsigned long rightStride = r.strides_[level];
1414  for (unsigned i = 0; i < imax; ++i)
1415  dualCircularLoop(level + 1,
1416  idx0 + ((i + leftShift) % leftPeriod) * leftStride,
1417  idx1 + ((i + rightShift) % rightPeriod) * rightStride,
1418  thisCorner,
1419  range,
1420  otherCorner,
1421  r,
1422  binaryFunct);
1423  }
1424  }
1425 
1426  template <typename Numeric, unsigned Len, unsigned Dim>
1427  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
1428  void ArrayND<Numeric, Len, Dim>::flatCircularLoop(unsigned level,
1429  unsigned long idx0,
1430  unsigned long idx1,
1431  const unsigned* thisCorner,
1432  const unsigned* range,
1433  const unsigned* otherCorner,
1434  ArrayND<Num2, Len2, Dim2>& r,
1435  Functor binaryFunct) {
1436  const unsigned imax = range[level];
1437  const unsigned rightShift = otherCorner[level];
1438  const unsigned rightPeriod = r.shape_[level];
1439 
1440  if (level == dim_ - 1) {
1441  Numeric* left = data_ + (idx0 + thisCorner[level]);
1442  Numeric* const lMax = data_ + (idx0 + shape_[level]);
1443  Num2* right = r.data_ + idx1;
1444 
1445  for (unsigned i = 0; i < imax && left < lMax; ++i)
1446  binaryFunct(*left++, right[(i + rightShift) % rightPeriod]);
1447  } else {
1448  const unsigned long leftStride = strides_[level];
1449  const unsigned long leftMax = idx0 + shape_[level] * leftStride;
1450  idx0 += thisCorner[level] * leftStride;
1451  const unsigned long rightStride = r.strides_[level];
1452 
1453  for (unsigned i = 0; i < imax && idx0 < leftMax; ++i, idx0 += leftStride)
1454  flatCircularLoop(level + 1,
1455  idx0,
1456  idx1 + ((i + rightShift) % rightPeriod) * rightStride,
1457  thisCorner,
1458  range,
1459  otherCorner,
1460  r,
1461  binaryFunct);
1462  }
1463  }
1464 
1465  template <typename Numeric, unsigned Len, unsigned Dim>
1466  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
1467  void ArrayND<Numeric, Len, Dim>::circularFlatLoop(unsigned level,
1468  unsigned long idx0,
1469  unsigned long idx1,
1470  const unsigned* thisCorner,
1471  const unsigned* range,
1472  const unsigned* otherCorner,
1473  ArrayND<Num2, Len2, Dim2>& r,
1474  Functor binaryFunct) {
1475  const unsigned imax = range[level];
1476  const unsigned leftShift = thisCorner[level];
1477  const unsigned leftPeriod = shape_[level];
1478 
1479  if (level == dim_ - 1) {
1480  Numeric* left = data_ + idx0;
1481  Num2* right = r.data_ + (idx1 + otherCorner[level]);
1482  Num2* const rMax = r.data_ + (idx1 + r.shape_[level]);
1483 
1484  for (unsigned i = 0; i < imax && right < rMax; ++i)
1485  binaryFunct(left[(i + leftShift) % leftPeriod], *right++);
1486  } else {
1487  const unsigned long leftStride = strides_[level];
1488  const unsigned long rightStride = r.strides_[level];
1489  const unsigned long rightMax = idx1 + r.shape_[level] * rightStride;
1490  idx1 += otherCorner[level] * rightStride;
1491 
1492  for (unsigned i = 0; i < imax && idx1 < rightMax; ++i, idx1 += rightStride)
1493  circularFlatLoop(level + 1,
1494  idx0 + ((i + leftShift) % leftPeriod) * leftStride,
1495  idx1,
1496  thisCorner,
1497  range,
1498  otherCorner,
1499  r,
1500  binaryFunct);
1501  }
1502  }
1503 
1504  me_macro_check_loop_prerequisites(jointSubrangeScan, commonSubrangeLoop)
1505  me_macro_check_loop_prerequisites(dualCircularScan, dualCircularLoop)
1506  me_macro_check_loop_prerequisites(flatCircularScan,
1507  flatCircularLoop) me_macro_check_loop_prerequisites(circularFlatScan,
1508  circularFlatLoop)
1509 
1510  template <typename Numeric, unsigned StackLen, unsigned StackDim>
1511  template <typename Num2, unsigned Len2, unsigned Dim2>
1512  Numeric ArrayND<Numeric, StackLen, StackDim>::marginalizeInnerLoop(unsigned long idx,
1513  const unsigned levelPr,
1514  unsigned long idxPr,
1515  const ArrayND<Num2, Len2, Dim2>& prior,
1516  const unsigned* indexMap) const {
1517  Numeric sum = Numeric();
1518  const unsigned long myStride = strides_[indexMap[levelPr]];
1519  const unsigned imax = prior.shape_[levelPr];
1520  assert(imax == shape_[indexMap[levelPr]]);
1521  if (levelPr == prior.dim_ - 1) {
1522  for (unsigned i = 0; i < imax; ++i)
1523  sum += data_[idx + i * myStride] * prior.data_[idxPr++];
1524  } else {
1525  const unsigned long priorStride = prior.strides_[levelPr];
1526  for (unsigned i = 0; i < imax; ++i) {
1527  sum += marginalizeInnerLoop(idx, levelPr + 1U, idxPr, prior, indexMap);
1528  idx += myStride;
1529  idxPr += priorStride;
1530  }
1531  }
1532  return sum;
1533  }
1534 
1535  template <typename Numeric, unsigned StackLen, unsigned StackDim>
1536  template <typename Num2, unsigned Len2, unsigned Dim2>
1537  void ArrayND<Numeric, StackLen, StackDim>::marginalizeLoop(const unsigned level,
1538  unsigned long idx,
1539  const unsigned levelRes,
1540  unsigned long idxRes,
1541  const ArrayND<Num2, Len2, Dim2>& prior,
1542  const unsigned* indexMap,
1543  ArrayND& result) const {
1544  if (level == dim_) {
1545  const Numeric res = marginalizeInnerLoop(idx, 0U, 0UL, prior, indexMap);
1546  if (result.dim_)
1547  result.data_[idxRes] = res;
1548  else
1549  result.localData_[0] = res;
1550  } else {
1551  // Check if this level is mapped or not
1552  bool mapped = false;
1553  for (unsigned i = 0; i < prior.dim_; ++i)
1554  if (level == indexMap[i]) {
1555  mapped = true;
1556  break;
1557  }
1558  if (mapped)
1559  marginalizeLoop(level + 1U, idx, levelRes, idxRes, prior, indexMap, result);
1560  else {
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  marginalizeLoop(level + 1U, idx, levelRes + 1U, idxRes, prior, indexMap, result);
1566  idx += myStride;
1567  idxRes += resStride;
1568  }
1569  }
1570  }
1571  }
1572 
1573  template <typename Numeric, unsigned StackLen, unsigned StackDim>
1574  template <typename Num2, unsigned Len2, unsigned Dim2>
1575  ArrayND<Numeric, StackLen, StackDim> ArrayND<Numeric, StackLen, StackDim>::marginalize(
1576  const ArrayND<Num2, Len2, Dim2>& prior, const unsigned* indexMap, const unsigned mapLen) const {
1577  if (!shapeIsKnown_)
1578  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"marginalize\"");
1579  if (!(prior.dim_ && prior.dim_ <= dim_))
1580  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::marginalize: incompatible argument array rank");
1581  const unsigned resultDim = dim_ - prior.dim_;
1582 
1583  // Check that the index map is reasonable
1584  if (mapLen != prior.dim_)
1585  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::marginalize: incompatible index map length");
1586  assert(indexMap);
1587  for (unsigned i = 0; i < mapLen; ++i) {
1588  const unsigned thisInd = indexMap[i];
1589  if (shape_[thisInd] != prior.shape_[i])
1590  throw npstat::NpstatInvalidArgument(
1591  "In npstat::ArrayND::marginalize: "
1592  "incompatible argument array dimensions");
1593  if (thisInd >= dim_)
1594  throw npstat::NpstatOutOfRange("In npstat::ArrayND::marginalize: index map entry out of range");
1595  for (unsigned j = 0; j < i; ++j)
1596  if (indexMap[j] == thisInd)
1597  throw npstat::NpstatInvalidArgument(
1598  "In npstat::ArrayND::marginalize: "
1599  "duplicate entry in the index map");
1600  }
1601 
1602  // Build the shape for the array of results
1603  ArrayShape newShape;
1604  newShape.reserve(resultDim);
1605  for (unsigned i = 0; i < dim_; ++i) {
1606  bool mapped = false;
1607  for (unsigned j = 0; j < mapLen; ++j)
1608  if (indexMap[j] == i) {
1609  mapped = true;
1610  break;
1611  }
1612  if (!mapped)
1613  newShape.push_back(shape_[i]);
1614  }
1615 
1616  ArrayND result(newShape);
1617  assert(result.dim_ == resultDim);
1618  bool calculated = false;
1619  if (resultDim == 0) {
1620  calculated = true;
1621  for (unsigned i = 0; i < dim_; ++i)
1622  if (indexMap[i] != i) {
1623  calculated = false;
1624  break;
1625  }
1626  if (calculated) {
1627  Numeric sum = Numeric();
1628  for (unsigned long i = 0; i < len_; ++i)
1629  sum += data_[i] * prior.data_[i];
1630  result.localData_[0] = sum;
1631  }
1632  }
1633 
1634  if (!calculated)
1635  marginalizeLoop(0U, 0UL, 0U, 0UL, prior, indexMap, result);
1636 
1637  return result;
1638  }
1639 
1640  template <typename Numeric, unsigned Len, unsigned Dim>
1641  void ArrayND<Numeric, Len, Dim>::buildFromShapePtr(const unsigned* sizes, const unsigned dim) {
1642  dim_ = dim;
1643  if (dim_) {
1644  assert(sizes);
1645  for (unsigned i = 0; i < dim_; ++i)
1646  if (sizes[i] == 0)
1647  throw npstat::NpstatInvalidArgument(
1648  "In npstat::ArrayND::buildFromShapePtr: "
1649  "detected span of zero");
1650 
1651  // Copy the array shape and figure out the array length
1652  shape_ = makeBuffer(dim_, localShape_, Dim);
1653  for (unsigned i = 0; i < dim_; ++i) {
1654  shape_[i] = sizes[i];
1655  len_ *= shape_[i];
1656  }
1657 
1658  // Figure out the array strides
1659  buildStrides();
1660 
1661  // Allocate the data array
1662  data_ = makeBuffer(len_, localData_, Len);
1663  } else {
1664  localData_[0] = Numeric();
1665  data_ = localData_;
1666  }
1667  }
1668 
1669  template <typename Numeric, unsigned StackLen, unsigned StackDim>
1670  template <typename Num2, unsigned Len2, unsigned Dim2>
1671  ArrayND<Numeric, StackLen, StackDim>::ArrayND(const ArrayND<Num2, Len2, Dim2>& slicedArray,
1672  const unsigned* fixedIndices,
1673  const unsigned nFixedIndices)
1674  : data_(0),
1675  strides_(nullptr),
1676  shape_(nullptr),
1677  len_(1UL),
1678  dim_(slicedArray.dim_ - nFixedIndices),
1679  shapeIsKnown_(true) {
1680  if (nFixedIndices) {
1681  assert(fixedIndices);
1682  if (nFixedIndices > slicedArray.dim_)
1683  throw npstat::NpstatInvalidArgument("In npstat::ArrayND slicing constructor: too many fixed indices");
1684  if (!slicedArray.shapeIsKnown_)
1685  throw npstat::NpstatInvalidArgument(
1686  "In npstat::ArrayND slicing constructor: "
1687  "uninitialized argument array");
1688 
1689  // Check that the fixed indices are within range
1690  for (unsigned j = 0; j < nFixedIndices; ++j)
1691  if (fixedIndices[j] >= slicedArray.dim_)
1692  throw npstat::NpstatOutOfRange(
1693  "In npstat::ArrayND slicing "
1694  "constructor: fixed index out of range");
1695 
1696  // Build the shape for the slice
1697  shape_ = makeBuffer(dim_, localShape_, StackDim);
1698  unsigned idim = 0;
1699  for (unsigned i = 0; i < slicedArray.dim_; ++i) {
1700  bool fixed = false;
1701  for (unsigned j = 0; j < nFixedIndices; ++j)
1702  if (fixedIndices[j] == i) {
1703  fixed = true;
1704  break;
1705  }
1706  if (!fixed) {
1707  assert(idim < dim_);
1708  shape_[idim++] = slicedArray.shape_[i];
1709  }
1710  }
1711  assert(idim == dim_);
1712 
1713  if (dim_) {
1714  // Copy the array shape and figure out the array length
1715  for (unsigned i = 0; i < dim_; ++i)
1716  len_ *= shape_[i];
1717 
1718  // Figure out the array strides
1719  buildStrides();
1720 
1721  // Allocate the data array
1722  data_ = makeBuffer(len_, localData_, StackLen);
1723  } else {
1724  localData_[0] = Numeric();
1725  data_ = localData_;
1726  }
1727  } else {
1728  new (this) ArrayND(slicedArray);
1729  }
1730  }
1731 
1732  template <typename Numeric, unsigned StackLen, unsigned StackDim>
1733  ArrayShape ArrayND<Numeric, StackLen, StackDim>::sliceShape(const unsigned* fixedIndices,
1734  const unsigned nFixedIndices) const {
1735  if (!shapeIsKnown_)
1736  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"sliceShape\"");
1737  if (nFixedIndices) {
1738  assert(fixedIndices);
1739  if (nFixedIndices > dim_)
1740  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::sliceShape: too many fixed indices");
1741  for (unsigned j = 0; j < nFixedIndices; ++j)
1742  if (fixedIndices[j] >= dim_)
1743  throw npstat::NpstatOutOfRange(
1744  "In npstat::ArrayND::sliceShape: "
1745  "fixed index out of range");
1746  ArrayShape sh;
1747  sh.reserve(dim_ - nFixedIndices);
1748  for (unsigned i = 0; i < dim_; ++i) {
1749  bool fixed = false;
1750  for (unsigned j = 0; j < nFixedIndices; ++j)
1751  if (fixedIndices[j] == i) {
1752  fixed = true;
1753  break;
1754  }
1755  if (!fixed)
1756  sh.push_back(shape_[i]);
1757  }
1758  return sh;
1759  } else
1760  return ArrayShape(shape_, shape_ + dim_);
1761  }
1762 
1763  template <typename Numeric, unsigned StackLen, unsigned StackDim>
1764  template <typename Num2, unsigned Len2, unsigned Dim2>
1765  unsigned long ArrayND<Numeric, StackLen, StackDim>::verifySliceCompatibility(const ArrayND<Num2, Len2, Dim2>& slice,
1766  const unsigned* fixedIndices,
1767  const unsigned* fixedIndexValues,
1768  const unsigned nFixedIndices) const {
1769  if (!(nFixedIndices && nFixedIndices <= dim_))
1770  throw npstat::NpstatInvalidArgument(
1771  "In npstat::ArrayND::verifySliceCompatibility: "
1772  "invalid number of fixed indices");
1773  if (!shapeIsKnown_)
1774  throw npstat::NpstatInvalidArgument(
1775  "Initialize npstat::ArrayND before calling "
1776  "method \"verifySliceCompatibility\"");
1777  if (!slice.shapeIsKnown_)
1778  throw npstat::NpstatInvalidArgument(
1779  "In npstat::ArrayND::verifySliceCompatibility: "
1780  "uninitialized argument array");
1781  if (slice.dim_ != dim_ - nFixedIndices)
1782  throw npstat::NpstatInvalidArgument(
1783  "In npstat::ArrayND::verifySliceCompatibility: "
1784  "incompatible argument array rank");
1785  assert(fixedIndices);
1786  assert(fixedIndexValues);
1787 
1788  for (unsigned j = 0; j < nFixedIndices; ++j)
1789  if (fixedIndices[j] >= dim_)
1790  throw npstat::NpstatOutOfRange(
1791  "In npstat::ArrayND::verifySliceCompatibility: "
1792  "fixed index out of range");
1793 
1794  // Check slice shape compatibility
1795  unsigned long idx = 0UL;
1796  unsigned sliceDim = 0U;
1797  for (unsigned i = 0; i < dim_; ++i) {
1798  bool fixed = false;
1799  for (unsigned j = 0; j < nFixedIndices; ++j)
1800  if (fixedIndices[j] == i) {
1801  fixed = true;
1802  if (fixedIndexValues[j] >= shape_[i])
1803  throw npstat::NpstatOutOfRange(
1804  "In npstat::ArrayND::verifySliceCompatibility: "
1805  "fixed index value out of range");
1806  idx += fixedIndexValues[j] * strides_[i];
1807  break;
1808  }
1809  if (!fixed) {
1810  if (shape_[i] != slice.shape_[sliceDim])
1811  throw npstat::NpstatInvalidArgument(
1812  "In npstat::ArrayND::verifySliceCompatibility: "
1813  "incompatible argument array dimensions");
1814  ++sliceDim;
1815  }
1816  }
1817  assert(sliceDim == slice.dim_);
1818  return idx;
1819  }
1820 
1821  template <typename Numeric, unsigned StackLen, unsigned StackDim>
1822  unsigned long ArrayND<Numeric, StackLen, StackDim>::verifyBufferSliceCompatibility(const unsigned long bufLen,
1823  const unsigned* fixedIndices,
1824  const unsigned* fixedIndexValues,
1825  const unsigned nFixedIndices,
1826  unsigned long* sliceStrides) const {
1827  if (!(nFixedIndices && nFixedIndices <= dim_))
1828  throw npstat::NpstatInvalidArgument(
1829  "In npstat::ArrayND::verifyBufferSliceCompatibility: "
1830  "invalid number of fixed indices");
1831  if (!shapeIsKnown_)
1832  throw npstat::NpstatInvalidArgument(
1833  "Initialize npstat::ArrayND before calling "
1834  "method \"verifyBufferSliceCompatibility\"");
1835  assert(fixedIndices);
1836  assert(fixedIndexValues);
1837 
1838  for (unsigned j = 0; j < nFixedIndices; ++j)
1839  if (fixedIndices[j] >= dim_)
1840  throw npstat::NpstatOutOfRange(
1841  "In npstat::ArrayND::verifyBufferSliceCompatibility: "
1842  "fixed index out of range");
1843 
1844  // Figure out slice shape. Store it, temporarily, in sliceStrides.
1845  unsigned long idx = 0UL;
1846  unsigned sliceDim = 0U;
1847  for (unsigned i = 0; i < dim_; ++i) {
1848  bool fixed = false;
1849  for (unsigned j = 0; j < nFixedIndices; ++j)
1850  if (fixedIndices[j] == i) {
1851  fixed = true;
1852  if (fixedIndexValues[j] >= shape_[i])
1853  throw npstat::NpstatOutOfRange(
1854  "In npstat::ArrayND::verifyBufferSliceCompatibility:"
1855  " fixed index value out of range");
1856  idx += fixedIndexValues[j] * strides_[i];
1857  break;
1858  }
1859  if (!fixed)
1860  sliceStrides[sliceDim++] = shape_[i];
1861  }
1862  assert(sliceDim + nFixedIndices == dim_);
1863 
1864  // Convert the slice shape into slice strides
1865  unsigned long expectedBufLen = 1UL;
1866  if (sliceDim) {
1867  unsigned long shapeJ = sliceStrides[sliceDim - 1];
1868  sliceStrides[sliceDim - 1] = 1UL;
1869  for (unsigned j = sliceDim - 1; j > 0; --j) {
1870  const unsigned long nextStride = sliceStrides[j] * shapeJ;
1871  shapeJ = sliceStrides[j - 1];
1872  sliceStrides[j - 1] = nextStride;
1873  }
1874  expectedBufLen = sliceStrides[0] * shapeJ;
1875  }
1876  if (expectedBufLen != bufLen)
1877  throw npstat::NpstatInvalidArgument(
1878  "In npstat::ArrayND::verifyBufferSliceCompatibility: "
1879  "invalid memory buffer length");
1880 
1881  return idx;
1882  }
1883 
1884  template <typename Numeric, unsigned StackLen, unsigned StackDim>
1885  template <typename Num2, class Op>
1886  void ArrayND<Numeric, StackLen, StackDim>::jointSliceLoop(const unsigned level,
1887  const unsigned long idx0,
1888  const unsigned level1,
1889  const unsigned long idx1,
1890  Num2* sliceData,
1891  const unsigned long* sliceStrides,
1892  const unsigned* fixedIndices,
1893  const unsigned* fixedIndexValues,
1894  const unsigned nFixedIndices,
1895  Op fcn) {
1896  bool fixed = false;
1897  for (unsigned j = 0; j < nFixedIndices; ++j)
1898  if (fixedIndices[j] == level) {
1899  fixed = true;
1900  break;
1901  }
1902  if (fixed)
1903  jointSliceLoop(
1904  level + 1, idx0, level1, idx1, sliceData, sliceStrides, fixedIndices, fixedIndexValues, nFixedIndices, fcn);
1905  else {
1906  const unsigned imax = shape_[level];
1907  const unsigned long stride = strides_[level];
1908 
1909  if (level1 == dim_ - nFixedIndices - 1) {
1910  sliceData += idx1;
1911  Numeric* localData = data_ + idx0;
1912  for (unsigned i = 0; i < imax; ++i)
1913  fcn(localData[i * stride], sliceData[i]);
1914  } else {
1915  const unsigned long stride2 = sliceStrides[level1];
1916  for (unsigned i = 0; i < imax; ++i)
1917  jointSliceLoop(level + 1,
1918  idx0 + i * stride,
1919  level1 + 1,
1920  idx1 + i * stride2,
1921  sliceData,
1922  sliceStrides,
1923  fixedIndices,
1924  fixedIndexValues,
1925  nFixedIndices,
1926  fcn);
1927  }
1928  }
1929  }
1930 
1931  template <typename Numeric, unsigned StackLen, unsigned StackDim>
1932  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
1933  void ArrayND<Numeric, StackLen, StackDim>::jointSliceScan(ArrayND<Num2, Len2, Dim2>& slice,
1934  const unsigned* fixedIndices,
1935  const unsigned* fixedIndexValues,
1936  const unsigned nFixedIndices,
1937  Functor binaryFunct) {
1938  const unsigned long idx = verifySliceCompatibility(slice, fixedIndices, fixedIndexValues, nFixedIndices);
1939  if (slice.dim_)
1940  jointSliceLoop(
1941  0U, idx, 0U, 0UL, slice.data_, slice.strides_, fixedIndices, fixedIndexValues, nFixedIndices, binaryFunct);
1942  else
1943  binaryFunct(data_[idx], slice.localData_[0]);
1944  }
1945 
1946  template <typename Numeric, unsigned StackLen, unsigned StackDim>
1947  template <typename Num2, class Functor>
1948  void ArrayND<Numeric, StackLen, StackDim>::jointMemSliceScan(Num2* slice,
1949  const unsigned long len,
1950  const unsigned* fixedIndices,
1951  const unsigned* fixedIndexValues,
1952  unsigned nFixedIndices,
1953  Functor binaryFunct) {
1954  assert(slice);
1955  if (dim_ > CHAR_BIT * sizeof(unsigned long))
1956  throw npstat::NpstatInvalidArgument(
1957  "In npstat::ArrayND::jointMemSliceScan: "
1958  "rank of this array is too large");
1959  unsigned long sliceStrides[CHAR_BIT * sizeof(unsigned long)];
1960  const unsigned long idx =
1961  verifyBufferSliceCompatibility(len, fixedIndices, fixedIndexValues, nFixedIndices, sliceStrides);
1962  if (dim_ > nFixedIndices)
1963  jointSliceLoop(0U, idx, 0U, 0UL, slice, sliceStrides, fixedIndices, fixedIndexValues, nFixedIndices, binaryFunct);
1964  else
1965  binaryFunct(data_[idx], *slice);
1966  }
1967 
1968  template <typename Numeric, unsigned StackLen, unsigned StackDim>
1969  template <typename Num2>
1970  void ArrayND<Numeric, StackLen, StackDim>::projectInnerLoop(const unsigned level,
1971  unsigned long idx0,
1972  unsigned* currentIndex,
1973  AbsArrayProjector<Numeric, Num2>& projector,
1974  const unsigned* projectedIndices,
1975  const unsigned nProjectedIndices) const {
1976  // level : dimension number among indices which are being iterated
1977  const unsigned idx = projectedIndices[level];
1978  const unsigned imax = shape_[idx];
1979  const unsigned long stride = strides_[idx];
1980  const bool last = (level == nProjectedIndices - 1);
1981 
1982  for (unsigned i = 0; i < imax; ++i) {
1983  currentIndex[idx] = i;
1984  if (last)
1985  projector.process(currentIndex, dim_, idx0, data_[idx0]);
1986  else
1987  projectInnerLoop(level + 1, idx0, currentIndex, projector, projectedIndices, nProjectedIndices);
1988  idx0 += stride;
1989  }
1990  }
1991 
1992  template <typename Numeric, unsigned StackLen, unsigned StackDim>
1993  template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3, class Op>
1994  void ArrayND<Numeric, StackLen, StackDim>::projectLoop(const unsigned level,
1995  const unsigned long idx0,
1996  const unsigned level1,
1997  const unsigned long idx1,
1998  unsigned* currentIndex,
1999  ArrayND<Num2, Len2, Dim2>* projection,
2000  AbsArrayProjector<Numeric, Num3>& projector,
2001  const unsigned* projectedIndices,
2002  const unsigned nProjectedIndices,
2003  Op fcn) const {
2004  // level : dimension number in this array
2005  // level1 : dimension number in the projection
2006  // idx0 : linear index in this array
2007  // idx1 : linear index in the projection
2008  // currentIndex : cycled over in this loop with the exception of the
2009  // dimensions which are iterated over to build the
2010  // projection
2011  if (level == dim_) {
2012  assert(level1 == projection->dim_);
2013  projector.clear();
2014  projectInnerLoop(0U, idx0, currentIndex, projector, projectedIndices, nProjectedIndices);
2015  if (projection->dim_)
2016  fcn(projection->data_[idx1], projector.result());
2017  else
2018  fcn(projection->localData_[0], projector.result());
2019  } else {
2020  bool iterated = false;
2021  for (unsigned j = 0; j < nProjectedIndices; ++j)
2022  if (projectedIndices[j] == level) {
2023  iterated = true;
2024  break;
2025  }
2026  if (iterated) {
2027  // This index will be iterated over inside "projectInnerLoop"
2028  projectLoop(
2029  level + 1, idx0, level1, idx1, currentIndex, projection, projector, projectedIndices, nProjectedIndices, fcn);
2030  } else {
2031  const unsigned imax = shape_[level];
2032  const unsigned long stride = strides_[level];
2033  // We will not be able to get here if projection->dim_ is 0.
2034  // Therefore, it is safe to access projection->strides_.
2035  const unsigned long stride2 = projection->strides_[level1];
2036  for (unsigned i = 0; i < imax; ++i) {
2037  currentIndex[level] = i;
2038  projectLoop(level + 1,
2039  idx0 + i * stride,
2040  level1 + 1,
2041  idx1 + i * stride2,
2042  currentIndex,
2043  projection,
2044  projector,
2045  projectedIndices,
2046  nProjectedIndices,
2047  fcn);
2048  }
2049  }
2050  }
2051  }
2052 
2053  template <typename Numeric, unsigned StackLen, unsigned StackDim>
2054  template <typename Num2, unsigned Len2, unsigned Dim2>
2055  void ArrayND<Numeric, StackLen, StackDim>::verifyProjectionCompatibility(const ArrayND<Num2, Len2, Dim2>& projection,
2056  const unsigned* projectedIndices,
2057  const unsigned nProjectedIndices) const {
2058  if (!(nProjectedIndices && nProjectedIndices <= dim_))
2059  throw npstat::NpstatInvalidArgument(
2060  "In npstat::ArrayND::verifyProjectionCompatibility: "
2061  "invalid number of projected indices");
2062  if (!shapeIsKnown_)
2063  throw npstat::NpstatInvalidArgument(
2064  "Initialize npstat::ArrayND before calling "
2065  "method \"verifyProjectionCompatibility\"");
2066  if (!projection.shapeIsKnown_)
2067  throw npstat::NpstatInvalidArgument(
2068  "In npstat::ArrayND::verifyProjectionCompatibility: "
2069  "uninitialized argument array");
2070  if (projection.dim_ != dim_ - nProjectedIndices)
2071  throw npstat::NpstatInvalidArgument(
2072  "In npstat::ArrayND::verifyProjectionCompatibility: "
2073  "incompatible argument array rank");
2074  assert(projectedIndices);
2075 
2076  for (unsigned j = 0; j < nProjectedIndices; ++j)
2077  if (projectedIndices[j] >= dim_)
2078  throw npstat::NpstatOutOfRange(
2079  "In npstat::ArrayND::verifyProjectionCompatibility: "
2080  "projected index out of range");
2081 
2082  // Check projection shape compatibility
2083  unsigned sliceDim = 0U;
2084  for (unsigned i = 0; i < dim_; ++i) {
2085  bool fixed = false;
2086  for (unsigned j = 0; j < nProjectedIndices; ++j)
2087  if (projectedIndices[j] == i) {
2088  fixed = true;
2089  break;
2090  }
2091  if (!fixed) {
2092  if (shape_[i] != projection.shape_[sliceDim])
2093  throw npstat::NpstatInvalidArgument(
2094  "In npstat::ArrayND::verifyProjectionCompatibility: "
2095  "incompatible argument array dimensions");
2096  ++sliceDim;
2097  }
2098  }
2099  assert(sliceDim == projection.dim_);
2100  }
2101 
2102  template <typename Numeric, unsigned StackLen, unsigned StackDim>
2103  template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
2104  void ArrayND<Numeric, StackLen, StackDim>::project(ArrayND<Num2, Len2, Dim2>* projection,
2105  AbsArrayProjector<Numeric, Num3>& projector,
2106  const unsigned* projectedIndices,
2107  const unsigned nProjectedIndices) const {
2108  assert(projection);
2109  verifyProjectionCompatibility(*projection, projectedIndices, nProjectedIndices);
2110  unsigned ibuf[StackDim];
2111  unsigned* buf = makeBuffer(dim_, ibuf, StackDim);
2112  for (unsigned i = 0; i < dim_; ++i)
2113  buf[i] = 0U;
2114  projectLoop(0U,
2115  0UL,
2116  0U,
2117  0UL,
2118  buf,
2119  projection,
2120  projector,
2121  projectedIndices,
2122  nProjectedIndices,
2123  scast_assign_left<Num2, Num3>());
2124  destroyBuffer(buf, ibuf);
2125  }
2126 
2127  template <typename Numeric, unsigned StackLen, unsigned StackDim>
2128  template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
2129  void ArrayND<Numeric, StackLen, StackDim>::addToProjection(ArrayND<Num2, Len2, Dim2>* projection,
2130  AbsArrayProjector<Numeric, Num3>& projector,
2131  const unsigned* projectedIndices,
2132  const unsigned nProjectedIndices) const {
2133  assert(projection);
2134  verifyProjectionCompatibility(*projection, projectedIndices, nProjectedIndices);
2135  unsigned ibuf[StackDim];
2136  unsigned* buf = makeBuffer(dim_, ibuf, StackDim);
2137  for (unsigned i = 0; i < dim_; ++i)
2138  buf[i] = 0U;
2139  projectLoop(0U,
2140  0UL,
2141  0U,
2142  0UL,
2143  buf,
2144  projection,
2145  projector,
2146  projectedIndices,
2147  nProjectedIndices,
2148  scast_pluseq_left<Num2, Num3>());
2149  destroyBuffer(buf, ibuf);
2150  }
2151 
2152  template <typename Numeric, unsigned StackLen, unsigned StackDim>
2153  template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
2154  void ArrayND<Numeric, StackLen, StackDim>::subtractFromProjection(ArrayND<Num2, Len2, Dim2>* projection,
2155  AbsArrayProjector<Numeric, Num3>& projector,
2156  const unsigned* projectedIndices,
2157  const unsigned nProjectedIndices) const {
2158  assert(projection);
2159  verifyProjectionCompatibility(*projection, projectedIndices, nProjectedIndices);
2160  unsigned ibuf[StackDim];
2161  unsigned* buf = makeBuffer(dim_, ibuf, StackDim);
2162  for (unsigned i = 0; i < dim_; ++i)
2163  buf[i] = 0U;
2164  projectLoop(0U,
2165  0UL,
2166  0U,
2167  0UL,
2168  buf,
2169  projection,
2170  projector,
2171  projectedIndices,
2172  nProjectedIndices,
2173  scast_minuseq_left<Num2, Num3>());
2174  destroyBuffer(buf, ibuf);
2175  }
2176 
2177  template <typename Numeric, unsigned StackLen, unsigned StackDim>
2178  template <typename Num2>
2179  void ArrayND<Numeric, StackLen, StackDim>::projectInnerLoop2(const unsigned level,
2180  const unsigned long idx0,
2181  AbsVisitor<Numeric, Num2>& projector,
2182  const unsigned* projectedIndices,
2183  const unsigned nProjectedIndices) const {
2184  const unsigned idx = projectedIndices[level];
2185  const unsigned imax = shape_[idx];
2186  const unsigned long stride = strides_[idx];
2187  const bool last = (level == nProjectedIndices - 1);
2188 
2189  for (unsigned i = 0; i < imax; ++i) {
2190  if (last)
2191  projector.process(data_[idx0 + i * stride]);
2192  else
2193  projectInnerLoop2(level + 1, idx0 + i * stride, projector, projectedIndices, nProjectedIndices);
2194  }
2195  }
2196 
2197  template <typename Numeric, unsigned StackLen, unsigned StackDim>
2198  template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3, class Op>
2199  void ArrayND<Numeric, StackLen, StackDim>::projectLoop2(const unsigned level,
2200  const unsigned long idx0,
2201  const unsigned level1,
2202  const unsigned long idx1,
2203  ArrayND<Num2, Len2, Dim2>* projection,
2204  AbsVisitor<Numeric, Num3>& projector,
2205  const unsigned* projectedIndices,
2206  const unsigned nProjectedIndices,
2207  Op fcn) const {
2208  if (level == dim_) {
2209  assert(level1 == projection->dim_);
2210  projector.clear();
2211  projectInnerLoop2(0U, idx0, projector, projectedIndices, nProjectedIndices);
2212  if (projection->dim_)
2213  fcn(projection->data_[idx1], projector.result());
2214  else
2215  fcn(projection->localData_[0], projector.result());
2216  } else {
2217  bool fixed = false;
2218  for (unsigned j = 0; j < nProjectedIndices; ++j)
2219  if (projectedIndices[j] == level) {
2220  fixed = true;
2221  break;
2222  }
2223  if (fixed) {
2224  projectLoop2(level + 1, idx0, level1, idx1, projection, projector, projectedIndices, nProjectedIndices, fcn);
2225  } else {
2226  const unsigned imax = shape_[level];
2227  const unsigned long stride = strides_[level];
2228  const unsigned long stride2 = projection->strides_[level1];
2229  for (unsigned i = 0; i < imax; ++i)
2230  projectLoop2(level + 1,
2231  idx0 + i * stride,
2232  level1 + 1,
2233  idx1 + i * stride2,
2234  projection,
2235  projector,
2236  projectedIndices,
2237  nProjectedIndices,
2238  fcn);
2239  }
2240  }
2241  }
2242 
2243  template <typename Numeric, unsigned StackLen, unsigned StackDim>
2244  template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
2245  void ArrayND<Numeric, StackLen, StackDim>::project(ArrayND<Num2, Len2, Dim2>* projection,
2246  AbsVisitor<Numeric, Num3>& projector,
2247  const unsigned* projectedIndices,
2248  const unsigned nProjectedIndices) const {
2249  assert(projection);
2250  verifyProjectionCompatibility(*projection, projectedIndices, nProjectedIndices);
2251  projectLoop2(
2252  0U, 0UL, 0U, 0UL, projection, projector, projectedIndices, nProjectedIndices, scast_assign_left<Num2, Num3>());
2253  }
2254 
2255  template <typename Numeric, unsigned StackLen, unsigned StackDim>
2256  template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
2257  void ArrayND<Numeric, StackLen, StackDim>::addToProjection(ArrayND<Num2, Len2, Dim2>* projection,
2258  AbsVisitor<Numeric, Num3>& projector,
2259  const unsigned* projectedIndices,
2260  const unsigned nProjectedIndices) const {
2261  assert(projection);
2262  verifyProjectionCompatibility(*projection, projectedIndices, nProjectedIndices);
2263  projectLoop2(
2264  0U, 0UL, 0U, 0UL, projection, projector, projectedIndices, nProjectedIndices, scast_pluseq_left<Num2, Num3>());
2265  }
2266 
2267  template <typename Numeric, unsigned StackLen, unsigned StackDim>
2268  template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
2269  void ArrayND<Numeric, StackLen, StackDim>::subtractFromProjection(ArrayND<Num2, Len2, Dim2>* projection,
2270  AbsVisitor<Numeric, Num3>& projector,
2271  const unsigned* projectedIndices,
2272  const unsigned nProjectedIndices) const {
2273  assert(projection);
2274  verifyProjectionCompatibility(*projection, projectedIndices, nProjectedIndices);
2275  projectLoop2(
2276  0U, 0UL, 0U, 0UL, projection, projector, projectedIndices, nProjectedIndices, scast_minuseq_left<Num2, Num3>());
2277  }
2278 
2279  template <typename Numeric, unsigned StackLen, unsigned StackDim>
2280  template <typename Num2, class Functor>
2281  void ArrayND<Numeric, StackLen, StackDim>::scaleBySliceInnerLoop(const unsigned level,
2282  const unsigned long idx0,
2283  Num2& scale,
2284  const unsigned* projectedIndices,
2285  const unsigned nProjectedIndices,
2286  Functor binaryFunct) {
2287  const unsigned idx = projectedIndices[level];
2288  const unsigned imax = shape_[idx];
2289  const unsigned long stride = strides_[idx];
2290 
2291  if (level == nProjectedIndices - 1) {
2292  Numeric* data = data_ + idx0;
2293  for (unsigned i = 0; i < imax; ++i)
2294  binaryFunct(data[i * stride], scale);
2295  } else
2296  for (unsigned i = 0; i < imax; ++i)
2297  scaleBySliceInnerLoop(level + 1, idx0 + i * stride, scale, projectedIndices, nProjectedIndices, binaryFunct);
2298  }
2299 
2300  template <typename Numeric, unsigned StackLen, unsigned StackDim>
2301  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
2302  void ArrayND<Numeric, StackLen, StackDim>::scaleBySliceLoop(unsigned level,
2303  unsigned long idx0,
2304  unsigned level1,
2305  unsigned long idx1,
2306  ArrayND<Num2, Len2, Dim2>& slice,
2307  const unsigned* projectedIndices,
2308  const unsigned nProjectedIndices,
2309  Functor binaryFunct) {
2310  if (level == dim_) {
2311  assert(level1 == slice.dim_);
2312  Num2& scaleFactor = slice.dim_ ? slice.data_[idx1] : slice.localData_[0];
2313  scaleBySliceInnerLoop(0U, idx0, scaleFactor, projectedIndices, nProjectedIndices, binaryFunct);
2314  } else {
2315  bool fixed = false;
2316  for (unsigned j = 0; j < nProjectedIndices; ++j)
2317  if (projectedIndices[j] == level) {
2318  fixed = true;
2319  break;
2320  }
2321  if (fixed) {
2322  scaleBySliceLoop(level + 1, idx0, level1, idx1, slice, projectedIndices, nProjectedIndices, binaryFunct);
2323  } else {
2324  const unsigned imax = shape_[level];
2325  const unsigned long stride = strides_[level];
2326  const unsigned long stride2 = slice.strides_[level1];
2327  for (unsigned i = 0; i < imax; ++i)
2328  scaleBySliceLoop(level + 1,
2329  idx0 + i * stride,
2330  level1 + 1,
2331  idx1 + i * stride2,
2332  slice,
2333  projectedIndices,
2334  nProjectedIndices,
2335  binaryFunct);
2336  }
2337  }
2338  }
2339 
2340  template <typename Numeric, unsigned StackLen, unsigned StackDim>
2341  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
2342  void ArrayND<Numeric, StackLen, StackDim>::jointScan(ArrayND<Num2, Len2, Dim2>& r, Functor binaryFunct) {
2343  if (!isShapeCompatible(r))
2344  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::jointScan: incompatible argument array shape");
2345  if (dim_)
2346  for (unsigned long i = 0; i < len_; ++i)
2347  binaryFunct(data_[i], r.data_[i]);
2348  else
2349  binaryFunct(localData_[0], r.localData_[0]);
2350  }
2351 
2352  template <typename Numeric, unsigned StackLen, unsigned StackDim>
2353  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
2354  void ArrayND<Numeric, StackLen, StackDim>::applySlice(ArrayND<Num2, Len2, Dim2>& slice,
2355  const unsigned* fixedIndices,
2356  const unsigned nFixedIndices,
2357  Functor binaryFunct) {
2358  if (nFixedIndices) {
2359  verifyProjectionCompatibility(slice, fixedIndices, nFixedIndices);
2360  if (slice.dim_ == 0U)
2361  for (unsigned long i = 0; i < len_; ++i)
2362  binaryFunct(data_[i], slice.localData_[0]);
2363  else
2364  scaleBySliceLoop(0U, 0UL, 0U, 0UL, slice, fixedIndices, nFixedIndices, binaryFunct);
2365  } else
2366  jointScan(slice, binaryFunct);
2367  }
2368 
2369  template <typename Numeric, unsigned StackLen, unsigned StackDim>
2370  inline bool ArrayND<Numeric, StackLen, StackDim>::isCompatible(const ArrayShape& shape) const {
2371  if (!shapeIsKnown_)
2372  return false;
2373  if (dim_ != shape.size())
2374  return false;
2375  if (dim_) {
2376  for (unsigned i = 0; i < dim_; ++i)
2377  if (shape_[i] != shape[i])
2378  return false;
2379  } else
2380  assert(len_ == 1UL);
2381  return true;
2382  }
2383 
2384  template <typename Numeric, unsigned StackLen, unsigned StackDim>
2385  template <typename Num2, unsigned Len2, unsigned Dim2>
2386  inline bool ArrayND<Numeric, StackLen, StackDim>::isShapeCompatible(const ArrayND<Num2, Len2, Dim2>& r) const {
2387  if (!shapeIsKnown_)
2388  return false;
2389  if (!r.shapeIsKnown_)
2390  return false;
2391  if (dim_ != r.dim_)
2392  return false;
2393  if (len_ != r.len_)
2394  return false;
2395  if (dim_) {
2396  assert(shape_);
2397  assert(r.shape_);
2398  for (unsigned i = 0; i < dim_; ++i)
2399  if (shape_[i] != r.shape_[i])
2400  return false;
2401  } else
2402  assert(len_ == 1UL);
2403  return true;
2404  }
2405 
2406  template <typename Numeric, unsigned Len, unsigned Dim>
2407  template <typename Num2, typename Integer>
2408  void ArrayND<Numeric, Len, Dim>::processSubrangeLoop(const unsigned level,
2409  unsigned long idx0,
2410  unsigned* currentIndex,
2411  AbsArrayProjector<Numeric, Num2>& f,
2412  const BoxND<Integer>& subrange) const {
2413  // Deal with possible negative limits first
2414  const Interval<Integer>& levelRange(subrange[level]);
2415  long long int iminl = static_cast<long long int>(levelRange.min());
2416  if (iminl < 0LL)
2417  iminl = 0LL;
2418  long long int imaxl = static_cast<long long int>(levelRange.max());
2419  if (imaxl < 0LL)
2420  imaxl = 0LL;
2421 
2422  // Now deal with possible out-of-range limits
2423  const unsigned imin = static_cast<unsigned>(iminl);
2424  unsigned imax = static_cast<unsigned>(imaxl);
2425  if (imax > shape_[level])
2426  imax = shape_[level];
2427 
2428  if (level == dim_ - 1) {
2429  idx0 += imin;
2430  for (unsigned i = imin; i < imax; ++i, ++idx0) {
2431  currentIndex[level] = i;
2432  f.process(currentIndex, dim_, idx0, data_[idx0]);
2433  }
2434  } else {
2435  const unsigned long stride = strides_[level];
2436  idx0 += imin * stride;
2437  for (unsigned i = imin; i < imax; ++i) {
2438  currentIndex[level] = i;
2439  processSubrangeLoop(level + 1U, idx0, currentIndex, f, subrange);
2440  idx0 += stride;
2441  }
2442  }
2443  }
2444 
2445  template <typename Numeric, unsigned Len, unsigned StackDim>
2446  template <typename Num2, typename Integer>
2447  void ArrayND<Numeric, Len, StackDim>::processSubrange(AbsArrayProjector<Numeric, Num2>& f,
2448  const BoxND<Integer>& subrange) const {
2449  if (!shapeIsKnown_)
2450  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"processSubrange\"");
2451  if (!dim_)
2452  throw npstat::NpstatInvalidArgument(
2453  "npstat::ArrayND::processSubrange method "
2454  "can not be used with array of 0 rank");
2455  if (dim_ != subrange.dim())
2456  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::processSubrange: incompatible subrange rank");
2457  unsigned ibuf[StackDim];
2458  unsigned* buf = makeBuffer(dim_, ibuf, StackDim);
2459  for (unsigned i = 0; i < dim_; ++i)
2460  buf[i] = 0U;
2461  processSubrangeLoop(0U, 0UL, buf, f, subrange);
2462  destroyBuffer(buf, ibuf);
2463  }
2464 
2465  template <typename Numeric, unsigned Len, unsigned Dim>
2466  template <typename Num2>
2467  inline ArrayND<Numeric, Len, Dim>& ArrayND<Numeric, Len, Dim>::setData(const Num2* data,
2468  const unsigned long dataLength) {
2469  if (!shapeIsKnown_)
2470  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"setData\"");
2471  if (dataLength != len_)
2472  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::setData: incompatible input data length");
2473  copyBuffer(data_, data, dataLength);
2474  return *this;
2475  }
2476 
2477  template <typename Numeric, unsigned Len, unsigned Dim>
2478  void ArrayND<Numeric, Len, Dim>::buildStrides() {
2479  assert(dim_);
2480  if (strides_ == nullptr)
2481  strides_ = makeBuffer(dim_, localStrides_, Dim);
2482  strides_[dim_ - 1] = 1UL;
2483  for (unsigned j = dim_ - 1; j > 0; --j)
2484  strides_[j - 1] = strides_[j] * shape_[j];
2485  }
2486 
2487  template <typename Numeric, unsigned StackLen, unsigned StackDim>
2488  inline ArrayND<Numeric, StackLen, StackDim>::ArrayND()
2489  : data_(nullptr), strides_(nullptr), shape_(nullptr), len_(0UL), dim_(0U), shapeIsKnown_(false) {
2490  localData_[0] = Numeric();
2491  data_ = localData_;
2492  }
2493 
2494  template <typename Numeric, unsigned Len, unsigned Dim>
2495  ArrayND<Numeric, Len, Dim>::ArrayND(const ArrayND& r)
2496  : data_(nullptr), strides_(nullptr), shape_(nullptr), len_(r.len_), dim_(r.dim_), shapeIsKnown_(r.shapeIsKnown_) {
2497  if (dim_) {
2498  // Copy the shape
2499  shape_ = makeBuffer(dim_, localShape_, Dim);
2500  copyBuffer(shape_, r.shape_, dim_);
2501 
2502  // Copy the strides
2503  strides_ = makeBuffer(dim_, localStrides_, Dim);
2504  copyBuffer(strides_, r.strides_, dim_);
2505 
2506  // Copy the data
2507  data_ = makeBuffer(len_, localData_, Len);
2508  copyBuffer(data_, r.data_, len_);
2509  } else {
2510  assert(len_ == 1UL);
2511  localData_[0] = r.localData_[0];
2512  data_ = localData_;
2513  }
2514  }
2515 
2516  template <typename Numeric, unsigned Len, unsigned Dim>
2517  template <typename Num2, unsigned Len2, unsigned Dim2>
2518  ArrayND<Numeric, Len, Dim>::ArrayND(const ArrayND<Num2, Len2, Dim2>& r)
2519  : data_(0), strides_(nullptr), shape_(nullptr), len_(r.len_), dim_(r.dim_), shapeIsKnown_(r.shapeIsKnown_) {
2520  if (dim_) {
2521  // Copy the shape
2522  shape_ = makeBuffer(dim_, localShape_, Dim);
2523  copyBuffer(shape_, r.shape_, dim_);
2524 
2525  // Copy the strides
2526  strides_ = makeBuffer(dim_, localStrides_, Dim);
2527  copyBuffer(strides_, r.strides_, dim_);
2528 
2529  // Copy the data
2530  data_ = makeBuffer(len_, localData_, Len);
2531  copyBuffer(data_, r.data_, len_);
2532  } else {
2533  assert(len_ == 1UL);
2534  localData_[0] = static_cast<Numeric>(r.localData_[0]);
2535  data_ = localData_;
2536  }
2537  }
2538 
2539  template <typename Numeric, unsigned Len, unsigned Dim>
2540  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
2541  ArrayND<Numeric, Len, Dim>::ArrayND(const ArrayND<Num2, Len2, Dim2>& r, Functor f)
2542  : data_(0), strides_(nullptr), shape_(nullptr), len_(r.len_), dim_(r.dim_), shapeIsKnown_(r.shapeIsKnown_) {
2543  if (dim_) {
2544  // Copy the shape
2545  shape_ = makeBuffer(dim_, localShape_, Dim);
2546  copyBuffer(shape_, r.shape_, dim_);
2547 
2548  // Copy the strides
2549  strides_ = makeBuffer(dim_, localStrides_, Dim);
2550  copyBuffer(strides_, r.strides_, dim_);
2551 
2552  // Copy the data
2553  data_ = makeBuffer(len_, localData_, Len);
2554  for (unsigned long i = 0; i < len_; ++i)
2555  data_[i] = static_cast<Numeric>(f(r.data_[i]));
2556  } else {
2557  assert(len_ == 1UL);
2558  localData_[0] = static_cast<Numeric>(f(r.localData_[0]));
2559  data_ = localData_;
2560  }
2561  }
2562 
2563  template <typename Numeric, unsigned Len, unsigned Dim>
2564  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
2565  void ArrayND<Numeric, Len, Dim>::copyRangeLoopFunct(const unsigned level,
2566  unsigned long idx0,
2567  unsigned long idx1,
2568  const ArrayND<Num2, Len2, Dim2>& r,
2569  const ArrayRange& range,
2570  Functor f) {
2571  const unsigned imax = shape_[level];
2572  if (level == dim_ - 1) {
2573  Numeric* to = data_ + idx0;
2574  const Num2* from = r.data_ + (idx1 + range[level].min());
2575  for (unsigned i = 0; i < imax; ++i)
2576  *to++ = static_cast<Numeric>(f(*from++));
2577  } else {
2578  const unsigned long fromstride = r.strides_[level];
2579  const unsigned long tostride = strides_[level];
2580  idx1 += range[level].min() * fromstride;
2581  for (unsigned i = 0; i < imax; ++i) {
2582  copyRangeLoopFunct(level + 1, idx0, idx1, r, range, f);
2583  idx0 += tostride;
2584  idx1 += fromstride;
2585  }
2586  }
2587  }
2588 
2589  template <typename Numeric, unsigned Len, unsigned Dim>
2590  template <typename Num2, unsigned Len2, unsigned Dim2>
2591  ArrayND<Numeric, Len, Dim>::ArrayND(const ArrayND<Num2, Len2, Dim2>& r, const ArrayRange& range)
2592  : data_(0), strides_(nullptr), shape_(nullptr), len_(r.len_), dim_(r.dim_), shapeIsKnown_(r.shapeIsKnown_) {
2593  if (!range.isCompatible(r.shape_, r.dim_))
2594  throw npstat::NpstatInvalidArgument("In npstat::ArrayND subrange constructor: invalid subrange");
2595  if (dim_) {
2596  len_ = range.rangeSize();
2597  if (!len_)
2598  throw npstat::NpstatInvalidArgument("In npstat::ArrayND subrange constructor: empty subrange");
2599 
2600  // Figure out the shape
2601  shape_ = makeBuffer(dim_, localShape_, Dim);
2602  range.rangeLength(shape_, dim_);
2603 
2604  // Figure out the strides
2605  buildStrides();
2606 
2607  // Allocate the data array
2608  data_ = makeBuffer(len_, localData_, Len);
2609 
2610  // Copy the data
2611  if (dim_ > CHAR_BIT * sizeof(unsigned long))
2612  throw npstat::NpstatInvalidArgument(
2613  "In npstat::ArrayND subrange constructor: "
2614  "input array rank is too large");
2615  unsigned lolim[CHAR_BIT * sizeof(unsigned long)];
2616  range.lowerLimits(lolim, dim_);
2617  unsigned toBuf[CHAR_BIT * sizeof(unsigned long)];
2618  clearBuffer(toBuf, dim_);
2619  (const_cast<ArrayND<Num2, Len2, Dim2>&>(r))
2620  .commonSubrangeLoop(0U, 0UL, 0UL, lolim, shape_, toBuf, *this, scast_assign_right<Num2, Numeric>());
2621  } else {
2622  assert(len_ == 1UL);
2623  localData_[0] = static_cast<Numeric>(r.localData_[0]);
2624  data_ = localData_;
2625  }
2626  }
2627 
2628  template <typename Numeric, unsigned Len, unsigned Dim>
2629  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
2630  ArrayND<Numeric, Len, Dim>::ArrayND(const ArrayND<Num2, Len2, Dim2>& r, const ArrayRange& range, Functor f)
2631  : data_(0), strides_(nullptr), shape_(nullptr), len_(r.len_), dim_(r.dim_), shapeIsKnown_(r.shapeIsKnown_) {
2632  if (!range.isCompatible(r.shape_, r.dim_))
2633  throw npstat::NpstatInvalidArgument(
2634  "In npstat::ArrayND transforming subrange constructor: "
2635  "incompatible subrange");
2636  if (dim_) {
2637  len_ = range.rangeSize();
2638  if (!len_)
2639  throw npstat::NpstatInvalidArgument(
2640  "In npstat::ArrayND transforming subrange constructor: "
2641  "empty subrange");
2642 
2643  // Figure out the shape
2644  shape_ = makeBuffer(dim_, localShape_, Dim);
2645  for (unsigned i = 0; i < dim_; ++i)
2646  shape_[i] = range[i].length();
2647 
2648  // Figure out the strides
2649  buildStrides();
2650 
2651  // Allocate the data array
2652  data_ = makeBuffer(len_, localData_, Len);
2653 
2654  // Transform the data
2655  copyRangeLoopFunct(0U, 0UL, 0UL, r, range, f);
2656  } else {
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  ArrayND<Numeric, Len, Dim>::ArrayND(const ArrayShape& sh)
2665  : data_(nullptr), strides_(nullptr), shape_(nullptr), len_(1UL), shapeIsKnown_(true) {
2666  const unsigned sz = sh.size();
2667  buildFromShapePtr(sz ? &sh[0] : nullptr, sz);
2668  }
2669 
2670  template <typename Numeric, unsigned Len, unsigned Dim>
2671  ArrayND<Numeric, Len, Dim>::ArrayND(const unsigned* sizes, const unsigned dim)
2672  : data_(0), strides_(nullptr), shape_(nullptr), len_(1UL), shapeIsKnown_(true) {
2673  buildFromShapePtr(sizes, dim);
2674  }
2675 
2676  template <typename Numeric, unsigned Len, unsigned Dim>
2677  ArrayND<Numeric, Len, Dim>::ArrayND(const unsigned n0)
2678  : data_(0), strides_(nullptr), shape_(nullptr), len_(1UL), shapeIsKnown_(true) {
2679  const unsigned dim = 1U;
2680  unsigned sizes[dim];
2681  sizes[0] = n0;
2682  buildFromShapePtr(sizes, dim);
2683  }
2684 
2685  template <typename Numeric, unsigned Len, unsigned Dim>
2686  ArrayND<Numeric, Len, Dim>::ArrayND(const unsigned n0, const unsigned n1)
2687  : data_(0), strides_(nullptr), shape_(nullptr), len_(1UL), shapeIsKnown_(true) {
2688  const unsigned dim = 2U;
2689  unsigned sizes[dim];
2690  sizes[0] = n0;
2691  sizes[1] = n1;
2692  buildFromShapePtr(sizes, dim);
2693  }
2694 
2695  template <typename Numeric, unsigned Len, unsigned Dim>
2696  ArrayND<Numeric, Len, Dim>::ArrayND(const unsigned n0, const unsigned n1, const unsigned n2)
2697  : data_(0), strides_(nullptr), shape_(nullptr), len_(1UL), shapeIsKnown_(true) {
2698  const unsigned dim = 3U;
2699  unsigned sizes[dim];
2700  sizes[0] = n0;
2701  sizes[1] = n1;
2702  sizes[2] = n2;
2703  buildFromShapePtr(sizes, dim);
2704  }
2705 
2706  template <typename Numeric, unsigned Len, unsigned Dim>
2707  ArrayND<Numeric, Len, Dim>::ArrayND(const unsigned n0, const unsigned n1, const unsigned n2, const unsigned n3)
2708  : data_(0), strides_(nullptr), shape_(nullptr), len_(1UL), shapeIsKnown_(true) {
2709  const unsigned dim = 4U;
2710  unsigned sizes[dim];
2711  sizes[0] = n0;
2712  sizes[1] = n1;
2713  sizes[2] = n2;
2714  sizes[3] = n3;
2715  buildFromShapePtr(sizes, dim);
2716  }
2717 
2718  template <typename Numeric, unsigned Len, unsigned Dim>
2719  ArrayND<Numeric, Len, Dim>::ArrayND(
2720  const unsigned n0, const unsigned n1, const unsigned n2, const unsigned n3, const unsigned n4)
2721  : data_(0), strides_(nullptr), shape_(nullptr), len_(1UL), shapeIsKnown_(true) {
2722  const unsigned dim = 5U;
2723  unsigned sizes[dim];
2724  sizes[0] = n0;
2725  sizes[1] = n1;
2726  sizes[2] = n2;
2727  sizes[3] = n3;
2728  sizes[4] = n4;
2729  buildFromShapePtr(sizes, dim);
2730  }
2731 
2732  template <typename Numeric, unsigned Len, unsigned Dim>
2733  ArrayND<Numeric, Len, Dim>::ArrayND(
2734  const unsigned n0, const unsigned n1, const unsigned n2, const unsigned n3, const unsigned n4, const unsigned n5)
2735  : data_(0), strides_(nullptr), shape_(nullptr), len_(1UL), shapeIsKnown_(true) {
2736  const unsigned dim = 6U;
2737  unsigned sizes[dim];
2738  sizes[0] = n0;
2739  sizes[1] = n1;
2740  sizes[2] = n2;
2741  sizes[3] = n3;
2742  sizes[4] = n4;
2743  sizes[5] = n5;
2744  buildFromShapePtr(sizes, dim);
2745  }
2746 
2747  template <typename Numeric, unsigned Len, unsigned Dim>
2748  ArrayND<Numeric, Len, Dim>::ArrayND(const unsigned n0,
2749  const unsigned n1,
2750  const unsigned n2,
2751  const unsigned n3,
2752  const unsigned n4,
2753  const unsigned n5,
2754  const unsigned n6)
2755  : data_(0), strides_(nullptr), shape_(nullptr), len_(1UL), shapeIsKnown_(true) {
2756  const unsigned dim = 7U;
2757  unsigned sizes[dim];
2758  sizes[0] = n0;
2759  sizes[1] = n1;
2760  sizes[2] = n2;
2761  sizes[3] = n3;
2762  sizes[4] = n4;
2763  sizes[5] = n5;
2764  sizes[6] = n6;
2765  buildFromShapePtr(sizes, dim);
2766  }
2767 
2768  template <typename Numeric, unsigned Len, unsigned Dim>
2769  ArrayND<Numeric, Len, Dim>::ArrayND(const unsigned n0,
2770  const unsigned n1,
2771  const unsigned n2,
2772  const unsigned n3,
2773  const unsigned n4,
2774  const unsigned n5,
2775  const unsigned n6,
2776  const unsigned n7)
2777  : data_(0), strides_(nullptr), shape_(nullptr), len_(1UL), shapeIsKnown_(true) {
2778  const unsigned dim = 8U;
2779  unsigned sizes[dim];
2780  sizes[0] = n0;
2781  sizes[1] = n1;
2782  sizes[2] = n2;
2783  sizes[3] = n3;
2784  sizes[4] = n4;
2785  sizes[5] = n5;
2786  sizes[6] = n6;
2787  sizes[7] = n7;
2788  buildFromShapePtr(sizes, dim);
2789  }
2790 
2791  template <typename Numeric, unsigned Len, unsigned Dim>
2792  ArrayND<Numeric, Len, Dim>::ArrayND(const unsigned n0,
2793  const unsigned n1,
2794  const unsigned n2,
2795  const unsigned n3,
2796  const unsigned n4,
2797  const unsigned n5,
2798  const unsigned n6,
2799  const unsigned n7,
2800  const unsigned n8)
2801  : data_(0), strides_(nullptr), shape_(nullptr), len_(1UL), shapeIsKnown_(true) {
2802  const unsigned dim = 9U;
2803  unsigned sizes[dim];
2804  sizes[0] = n0;
2805  sizes[1] = n1;
2806  sizes[2] = n2;
2807  sizes[3] = n3;
2808  sizes[4] = n4;
2809  sizes[5] = n5;
2810  sizes[6] = n6;
2811  sizes[7] = n7;
2812  sizes[8] = n8;
2813  buildFromShapePtr(sizes, dim);
2814  }
2815 
2816  template <typename Numeric, unsigned Len, unsigned Dim>
2817  ArrayND<Numeric, Len, Dim>::ArrayND(const unsigned n0,
2818  const unsigned n1,
2819  const unsigned n2,
2820  const unsigned n3,
2821  const unsigned n4,
2822  const unsigned n5,
2823  const unsigned n6,
2824  const unsigned n7,
2825  const unsigned n8,
2826  const unsigned n9)
2827  : data_(0), strides_(nullptr), shape_(nullptr), len_(1UL), shapeIsKnown_(true) {
2828  const unsigned dim = 10U;
2829  unsigned sizes[dim];
2830  sizes[0] = n0;
2831  sizes[1] = n1;
2832  sizes[2] = n2;
2833  sizes[3] = n3;
2834  sizes[4] = n4;
2835  sizes[5] = n5;
2836  sizes[6] = n6;
2837  sizes[7] = n7;
2838  sizes[8] = n8;
2839  sizes[9] = n9;
2840  buildFromShapePtr(sizes, dim);
2841  }
2842 
2843  template <typename Numeric, unsigned Len, unsigned Dim>
2844  template <typename Num1, unsigned Len1, unsigned Dim1, typename Num2, unsigned Len2, unsigned Dim2>
2845  void ArrayND<Numeric, Len, Dim>::outerProductLoop(const unsigned level,
2846  unsigned long idx0,
2847  unsigned long idx1,
2848  unsigned long idx2,
2849  const ArrayND<Num1, Len1, Dim1>& a1,
2850  const ArrayND<Num2, Len2, Dim2>& a2) {
2851  const unsigned imax = shape_[level];
2852  if (level == dim_ - 1) {
2853  for (unsigned i = 0; i < imax; ++i)
2854  data_[idx0 + i] = a1.data_[idx1] * a2.data_[idx2 + i];
2855  } else {
2856  for (unsigned i = 0; i < imax; ++i) {
2857  outerProductLoop(level + 1, idx0, idx1, idx2, a1, a2);
2858  idx0 += strides_[level];
2859  if (level < a1.dim_)
2860  idx1 += a1.strides_[level];
2861  else
2862  idx2 += a2.strides_[level - a1.dim_];
2863  }
2864  }
2865  }
2866 
2867  template <typename Numeric, unsigned Len, unsigned Dim>
2868  template <typename Num1, unsigned Len1, unsigned Dim1, typename Num2, unsigned Len2, unsigned Dim2>
2869  ArrayND<Numeric, Len, Dim>::ArrayND(const ArrayND<Num1, Len1, Dim1>& a1, const ArrayND<Num2, Len2, Dim2>& a2)
2870  : data_(0), strides_(nullptr), shape_(nullptr), len_(1UL), dim_(a1.dim_ + a2.dim_), shapeIsKnown_(true) {
2871  if (!(a1.shapeIsKnown_ && a2.shapeIsKnown_))
2872  throw npstat::NpstatInvalidArgument(
2873  "In npstat::ArrayND outer product constructor: "
2874  "uninitialized argument array");
2875  if (dim_) {
2876  shape_ = makeBuffer(dim_, localShape_, Dim);
2877  copyBuffer(shape_, a1.shape_, a1.dim_);
2878  copyBuffer(shape_ + a1.dim_, a2.shape_, a2.dim_);
2879 
2880  for (unsigned i = 0; i < dim_; ++i) {
2881  assert(shape_[i]);
2882  len_ *= shape_[i];
2883  }
2884 
2885  // Figure out the array strides
2886  buildStrides();
2887 
2888  // Allocate the data array
2889  data_ = makeBuffer(len_, localData_, Len);
2890 
2891  // Fill the data array
2892  if (a1.dim_ == 0) {
2893  for (unsigned long i = 0; i < len_; ++i)
2894  data_[i] = a1.localData_[0] * a2.data_[i];
2895  } else if (a2.dim_ == 0) {
2896  for (unsigned long i = 0; i < len_; ++i)
2897  data_[i] = a1.data_[i] * a2.localData_[0];
2898  } else
2899  outerProductLoop(0U, 0UL, 0UL, 0UL, a1, a2);
2900  } else {
2901  localData_[0] = a1.localData_[0] * a2.localData_[0];
2902  data_ = localData_;
2903  }
2904  }
2905 
2906  template <typename Numeric, unsigned Len, unsigned Dim>
2907  inline ArrayND<Numeric, Len, Dim>::~ArrayND() {
2908  destroyBuffer(data_, localData_);
2909  destroyBuffer(strides_, localStrides_);
2910  destroyBuffer(shape_, localShape_);
2911  }
2912 
2913  template <typename Numeric, unsigned Len, unsigned Dim>
2914  ArrayND<Numeric, Len, Dim>& ArrayND<Numeric, Len, Dim>::operator=(const ArrayND& r) {
2915  if (this == &r)
2916  return *this;
2917  if (shapeIsKnown_) {
2918  if (!r.shapeIsKnown_)
2919  throw npstat::NpstatInvalidArgument(
2920  "In npstat::ArrayND assignment operator: "
2921  "uninitialized argument array");
2922  if (!isShapeCompatible(r))
2923  throw npstat::NpstatInvalidArgument(
2924  "In npstat::ArrayND assignment operator: "
2925  "incompatible argument array shape");
2926  if (dim_)
2927  copyBuffer(data_, r.data_, len_);
2928  else
2929  localData_[0] = r.localData_[0];
2930  } else {
2931  // This object is uninitialized. If the object on the
2932  // right is itself initialized, make an in-place copy.
2933  if (r.shapeIsKnown_)
2934  new (this) ArrayND(r);
2935  }
2936  return *this;
2937  }
2938 
2939  template <typename Numeric, unsigned Len, unsigned Dim>
2940  template <typename Num2, unsigned Len2, unsigned Dim2>
2941  ArrayND<Numeric, Len, Dim>& ArrayND<Numeric, Len, Dim>::operator=(const ArrayND<Num2, Len2, Dim2>& r) {
2942  if ((void*)this == (void*)(&r))
2943  return *this;
2944  if (shapeIsKnown_) {
2945  if (!r.shapeIsKnown_)
2946  throw npstat::NpstatInvalidArgument(
2947  "In npstat::ArrayND assignment operator: "
2948  "uninitialized argument array");
2949  if (!isShapeCompatible(r))
2950  throw npstat::NpstatInvalidArgument(
2951  "In npstat::ArrayND assignment operator: "
2952  "incompatible argument array shape");
2953  if (dim_)
2954  copyBuffer(data_, r.data_, len_);
2955  else
2956  localData_[0] = static_cast<Numeric>(r.localData_[0]);
2957  } else {
2958  // This object is uninitialized. If the object on the
2959  // right is itself initialized, make an in-place copy.
2960  if (r.shapeIsKnown_)
2961  new (this) ArrayND(r);
2962  }
2963  return *this;
2964  }
2965 
2966  template <typename Numeric, unsigned Len, unsigned Dim>
2967  template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
2968  ArrayND<Numeric, Len, Dim>& ArrayND<Numeric, Len, Dim>::assign(const ArrayND<Num2, Len2, Dim2>& r, Functor f) {
2969  if (shapeIsKnown_) {
2970  if (!r.shapeIsKnown_)
2971  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::assign: uninitialized argument array");
2972  if (!isShapeCompatible(r))
2973  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::assign: incompatible argument array shape");
2974  if (dim_)
2975  for (unsigned long i = 0; i < len_; ++i)
2976  data_[i] = static_cast<Numeric>(f(r.data_[i]));
2977  else
2978  localData_[0] = static_cast<Numeric>(f(r.localData_[0]));
2979  } else {
2980  // This object is uninitialized. If the object on the
2981  // right is itself initialized, build new array in place.
2982  if (r.shapeIsKnown_)
2983  new (this) ArrayND(r, f);
2984  }
2985  return *this;
2986  }
2987 
2988  template <typename Numeric, unsigned Len, unsigned Dim>
2989  inline ArrayShape ArrayND<Numeric, Len, Dim>::shape() const {
2990  if (!shapeIsKnown_)
2991  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"shape\"");
2992  return ArrayShape(shape_, shape_ + dim_);
2993  }
2994 
2995  template <typename Numeric, unsigned Len, unsigned Dim>
2996  inline ArrayRange ArrayND<Numeric, Len, Dim>::fullRange() const {
2997  if (!shapeIsKnown_)
2998  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"fullRange\"");
2999  ArrayRange range;
3000  if (dim_) {
3001  range.reserve(dim_);
3002  for (unsigned i = 0; i < dim_; ++i)
3003  range.push_back(Interval<unsigned>(0U, shape_[i]));
3004  }
3005  return range;
3006  }
3007 
3008  template <typename Numeric, unsigned Len, unsigned Dim>
3009  bool ArrayND<Numeric, Len, Dim>::isDensity() const {
3010  if (!shapeIsKnown_)
3011  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"isDensity\"");
3012  const Numeric zero = Numeric();
3013  bool hasPositive = false;
3014  if (dim_)
3015  for (unsigned long i = 0; i < len_; ++i) {
3016  // Don't make comparisons whose result can be
3017  // determined in advance by assuming that Numeric
3018  // is an unsigned type. Some compilers will
3019  // complain about it when this template is
3020  // instantiated with such a type.
3021  if (data_[i] == zero)
3022  continue;
3023  if (ComplexComparesFalse<Numeric>::less(zero, data_[i]))
3024  hasPositive = true;
3025  else
3026  return false;
3027  }
3028  else
3029  hasPositive = ComplexComparesFalse<Numeric>::less(zero, localData_[0]);
3030  return hasPositive;
3031  }
3032 
3033  template <typename Numeric, unsigned Len, unsigned Dim>
3034  bool ArrayND<Numeric, Len, Dim>::isZero() const {
3035  if (!shapeIsKnown_)
3036  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"isZero\"");
3037  const Numeric zero = Numeric();
3038  if (dim_) {
3039  for (unsigned long i = 0; i < len_; ++i)
3040  if (data_[i] != zero)
3041  return false;
3042  } else if (localData_[0] != zero)
3043  return false;
3044  return true;
3045  }
3046 
3047  template <typename Numeric, unsigned Len, unsigned Dim>
3048  void ArrayND<Numeric, Len, Dim>::convertLinearIndex(unsigned long l, unsigned* idx, const unsigned idxLen) const {
3049  if (!shapeIsKnown_)
3050  throw npstat::NpstatInvalidArgument(
3051  "Initialize npstat::ArrayND before calling "
3052  "method \"convertLinearIndex\"");
3053  if (!dim_)
3054  throw npstat::NpstatInvalidArgument(
3055  "npstat::ArrayND::convertLinearIndex method "
3056  "can not be used with array of 0 rank");
3057  if (idxLen != dim_)
3058  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::convertLinearIndex: incompatible index length");
3059  if (l >= len_)
3060  throw npstat::NpstatOutOfRange("In npstat::ArrayND::convertLinearIndex: linear index out of range");
3061  assert(idx);
3062 
3063  for (unsigned i = 0; i < dim_; ++i) {
3064  idx[i] = l / strides_[i];
3065  l -= (idx[i] * strides_[i]);
3066  }
3067  }
3068 
3069  template <typename Numeric, unsigned Len, unsigned Dim>
3070  unsigned long ArrayND<Numeric, Len, Dim>::linearIndex(const unsigned* index, unsigned idxLen) const {
3071  if (!shapeIsKnown_)
3072  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"linearIndex\"");
3073  if (!dim_)
3074  throw npstat::NpstatInvalidArgument(
3075  "npstat::ArrayND::linearIndex method "
3076  "can not be used with array of 0 rank");
3077  if (idxLen != dim_)
3078  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::linearIndex: incompatible index length");
3079  assert(index);
3080 
3081  unsigned long idx = 0UL;
3082  for (unsigned i = 0; i < dim_; ++i) {
3083  if (index[i] >= shape_[i])
3084  throw npstat::NpstatOutOfRange("In npstat::ArrayND::linearIndex: index out of range");
3085  idx += index[i] * strides_[i];
3086  }
3087  return idx;
3088  }
3089 
3090  template <typename Numeric, unsigned Len, unsigned Dim>
3091  inline Numeric& ArrayND<Numeric, Len, Dim>::value(const unsigned* index, const unsigned dim) {
3092  if (!shapeIsKnown_)
3093  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"value\"");
3094  if (dim != dim_)
3095  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::value: incompatible index length");
3096  if (dim) {
3097  assert(index);
3098  unsigned long idx = 0UL;
3099  for (unsigned i = 0; i < dim_; ++i)
3100  idx += index[i] * strides_[i];
3101  return data_[idx];
3102  } else
3103  return localData_[0];
3104  }
3105 
3106  template <typename Numeric, unsigned Len, unsigned Dim>
3107  inline const Numeric& ArrayND<Numeric, Len, Dim>::value(const unsigned* index, const unsigned dim) const {
3108  if (!shapeIsKnown_)
3109  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"value\"");
3110  if (dim != dim_)
3111  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::value: incompatible index length");
3112  if (dim) {
3113  assert(index);
3114  unsigned long idx = 0UL;
3115  for (unsigned i = 0; i < dim_; ++i)
3116  idx += index[i] * strides_[i];
3117  return data_[idx];
3118  } else
3119  return localData_[0];
3120  }
3121 
3122  template <typename Numeric, unsigned Len, unsigned Dim>
3123  inline Numeric& ArrayND<Numeric, Len, Dim>::linearValue(const unsigned long index) {
3124  return data_[index];
3125  }
3126 
3127  template <typename Numeric, unsigned Len, unsigned Dim>
3128  inline const Numeric& ArrayND<Numeric, Len, Dim>::linearValue(const unsigned long index) const {
3129  return data_[index];
3130  }
3131 
3132  template <typename Numeric, unsigned Len, unsigned Dim>
3133  inline Numeric& ArrayND<Numeric, Len, Dim>::linearValueAt(const unsigned long index) {
3134  if (index >= len_)
3135  throw npstat::NpstatOutOfRange("In npstat::ArrayND::linearValueAt: linear index out of range");
3136  return data_[index];
3137  }
3138 
3139  template <typename Numeric, unsigned Len, unsigned Dim>
3140  inline const Numeric& ArrayND<Numeric, Len, Dim>::linearValueAt(const unsigned long index) const {
3141  if (index >= len_)
3142  throw npstat::NpstatOutOfRange("In npstat::ArrayND::linearValueAt: linear index out of range");
3143  return data_[index];
3144  }
3145 
3146  template <typename Numeric, unsigned Len, unsigned Dim>
3147  inline unsigned ArrayND<Numeric, Len, Dim>::coordToIndex(const double x, const unsigned idim) const {
3148  if (x <= 0.0)
3149  return 0;
3150  else if (x >= static_cast<double>(shape_[idim] - 1))
3151  return shape_[idim] - 1;
3152  else
3153  return static_cast<unsigned>(std::floor(x + 0.5));
3154  }
3155 
3156  template <typename Numeric, unsigned Len, unsigned Dim>
3157  inline const Numeric& ArrayND<Numeric, Len, Dim>::closest(const double* x, const unsigned dim) const {
3158  if (!shapeIsKnown_)
3159  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"closest\"");
3160  if (dim != dim_)
3161  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::closest: incompatible data length");
3162  if (dim) {
3163  assert(x);
3164  unsigned long idx = 0UL;
3165  for (unsigned i = 0; i < dim_; ++i)
3166  idx += coordToIndex(x[i], i) * strides_[i];
3167  return data_[idx];
3168  } else
3169  return localData_[0];
3170  }
3171 
3172  template <typename Numeric, unsigned Len, unsigned Dim>
3173  inline Numeric& ArrayND<Numeric, Len, Dim>::closest(const double* x, const unsigned dim) {
3174  if (!shapeIsKnown_)
3175  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"closest\"");
3176  if (dim != dim_)
3177  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::closest: incompatible data length");
3178  if (dim) {
3179  assert(x);
3180  unsigned long idx = 0UL;
3181  for (unsigned i = 0; i < dim_; ++i)
3182  idx += coordToIndex(x[i], i) * strides_[i];
3183  return data_[idx];
3184  } else
3185  return localData_[0];
3186  }
3187 
3188  template <typename Numeric, unsigned Len, unsigned Dim>
3189  inline const Numeric& ArrayND<Numeric, Len, Dim>::valueAt(const unsigned* index, const unsigned dim) const {
3190  if (!shapeIsKnown_)
3191  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"valueAt\"");
3192  if (dim != dim_)
3193  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::valueAt: incompatible index length");
3194  if (dim) {
3195  assert(index);
3196  unsigned long idx = 0UL;
3197  for (unsigned i = 0; i < dim_; ++i) {
3198  if (index[i] >= shape_[i])
3199  throw npstat::NpstatOutOfRange("In npstat::ArrayND::valueAt: index out of range");
3200  idx += index[i] * strides_[i];
3201  }
3202  return data_[idx];
3203  } else
3204  return localData_[0];
3205  }
3206 
3207  template <typename Numeric, unsigned Len, unsigned Dim>
3208  inline Numeric& ArrayND<Numeric, Len, Dim>::valueAt(const unsigned* index, const unsigned dim) {
3209  if (!shapeIsKnown_)
3210  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"valueAt\"");
3211  if (dim != dim_)
3212  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::valueAt: incompatible index length");
3213  if (dim) {
3214  assert(index);
3215  unsigned long idx = 0UL;
3216  for (unsigned i = 0; i < dim_; ++i) {
3217  if (index[i] >= shape_[i])
3218  throw npstat::NpstatOutOfRange("In npstat::ArrayND::valueAt: index out of range");
3219  idx += index[i] * strides_[i];
3220  }
3221  return data_[idx];
3222  } else
3223  return localData_[0];
3224  }
3225 
3226  template <typename Numeric, unsigned Len, unsigned Dim>
3227  inline Numeric& ArrayND<Numeric, Len, Dim>::operator()() {
3228  if (!shapeIsKnown_)
3229  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"operator()\"");
3230  if (dim_)
3231  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::operator(): wrong # of args (not rank 0 array)");
3232  return localData_[0];
3233  }
3234 
3235  template <typename Numeric, unsigned Len, unsigned Dim>
3236  inline const Numeric& ArrayND<Numeric, Len, Dim>::operator()() const {
3237  if (!shapeIsKnown_)
3238  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"operator()\"");
3239  if (dim_)
3240  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::operator(): wrong # of args (not rank 0 array)");
3241  return localData_[0];
3242  }
3243 
3244  template <typename Numeric, unsigned Len, unsigned Dim>
3245  inline Numeric& ArrayND<Numeric, Len, Dim>::operator()(const unsigned i) {
3246  if (1U != dim_)
3247  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::operator(): wrong # of args (not rank 1 array)");
3248  return data_[i];
3249  }
3250 
3251  template <typename Numeric, unsigned Len, unsigned Dim>
3252  inline const Numeric& ArrayND<Numeric, Len, Dim>::operator()(const unsigned i) const {
3253  if (1U != dim_)
3254  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::operator(): wrong # of args (not rank 1 array)");
3255  return data_[i];
3256  }
3257 
3258  template <typename Numeric, unsigned Len, unsigned Dim>
3259  const Numeric& ArrayND<Numeric, Len, Dim>::at() const {
3260  if (!shapeIsKnown_)
3261  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"at\"");
3262  if (dim_)
3263  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::at: wrong # of args (not rank 0 array)");
3264  return localData_[0];
3265  }
3266 
3267  template <typename Numeric, unsigned Len, unsigned Dim>
3268  Numeric& ArrayND<Numeric, Len, Dim>::at() {
3269  if (!shapeIsKnown_)
3270  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"at\"");
3271  if (dim_)
3272  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::at: wrong # of args (not rank 0 array)");
3273  return localData_[0];
3274  }
3275 
3276  template <typename Numeric, unsigned Len, unsigned Dim>
3277  const Numeric& ArrayND<Numeric, Len, Dim>::at(const unsigned i0) const {
3278  if (1U != dim_)
3279  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::at: wrong # of args (not rank 1 array)");
3280  if (i0 >= shape_[0])
3281  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 0 out of range (rank 1)");
3282  return data_[i0];
3283  }
3284 
3285  template <typename Numeric, unsigned Len, unsigned Dim>
3286  Numeric& ArrayND<Numeric, Len, Dim>::at(const unsigned i0) {
3287  if (1U != dim_)
3288  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::at: wrong # of args (not rank 1 array)");
3289  if (i0 >= shape_[0])
3290  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 0 out of range (rank 1)");
3291  return data_[i0];
3292  }
3293 
3294  template <typename Numeric, unsigned Len, unsigned Dim>
3295  inline Numeric& ArrayND<Numeric, Len, Dim>::operator()(const unsigned i0, const unsigned i1) {
3296  if (2U != dim_)
3297  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::operator(): wrong # of args (not rank 2 array)");
3298  return data_[i0 * strides_[0] + i1];
3299  }
3300 
3301  template <typename Numeric, unsigned Len, unsigned Dim>
3302  inline const Numeric& ArrayND<Numeric, Len, Dim>::operator()(const unsigned i0, const unsigned i1) const {
3303  if (2U != dim_)
3304  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::operator(): wrong # of args (not rank 2 array)");
3305  return data_[i0 * strides_[0] + i1];
3306  }
3307 
3308  template <typename Numeric, unsigned Len, unsigned Dim>
3309  const Numeric& ArrayND<Numeric, Len, Dim>::at(const unsigned i0, const unsigned i1) const {
3310  if (2U != dim_)
3311  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::at: wrong # of args (not rank 2 array)");
3312  if (i0 >= shape_[0])
3313  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 0 out of range (rank 2)");
3314  if (i1 >= shape_[1])
3315  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 1 out of range (rank 2)");
3316  return data_[i0 * strides_[0] + i1];
3317  }
3318 
3319  template <typename Numeric, unsigned Len, unsigned Dim>
3320  Numeric& ArrayND<Numeric, Len, Dim>::at(const unsigned i0, const unsigned i1) {
3321  if (2U != dim_)
3322  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::at: wrong # of args (not rank 2 array)");
3323  if (i0 >= shape_[0])
3324  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 0 out of range (rank 2)");
3325  if (i1 >= shape_[1])
3326  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 1 out of range (rank 2)");
3327  return data_[i0 * strides_[0] + i1];
3328  }
3329 
3330  template <typename Numeric, unsigned Len, unsigned Dim>
3331  inline const Numeric& ArrayND<Numeric, Len, Dim>::operator()(const unsigned i0,
3332  const unsigned i1,
3333  const unsigned i2) const {
3334  if (3U != dim_)
3335  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::operator(): wrong # of args (not rank 3 array)");
3336  return data_[i0 * strides_[0] + i1 * strides_[1] + i2];
3337  }
3338 
3339  template <typename Numeric, unsigned Len, unsigned Dim>
3340  inline const Numeric& ArrayND<Numeric, Len, Dim>::operator()(const unsigned i0,
3341  const unsigned i1,
3342  const unsigned i2,
3343  const unsigned i3) const {
3344  if (4U != dim_)
3345  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::operator(): wrong # of args (not rank 4 array)");
3346  return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3];
3347  }
3348 
3349  template <typename Numeric, unsigned Len, unsigned Dim>
3350  inline const Numeric& ArrayND<Numeric, Len, Dim>::operator()(
3351  const unsigned i0, const unsigned i1, const unsigned i2, const unsigned i3, const unsigned i4) const {
3352  if (5U != dim_)
3353  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::operator(): wrong # of args (not rank 5 array)");
3354  return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3 * strides_[3] + i4];
3355  }
3356 
3357  template <typename Numeric, unsigned Len, unsigned Dim>
3358  inline const Numeric& ArrayND<Numeric, Len, Dim>::operator()(const unsigned i0,
3359  const unsigned i1,
3360  const unsigned i2,
3361  const unsigned i3,
3362  const unsigned i4,
3363  const unsigned i5) const {
3364  if (6U != dim_)
3365  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::operator(): wrong # of args (not rank 6 array)");
3366  return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3 * strides_[3] + i4 * strides_[4] + i5];
3367  }
3368 
3369  template <typename Numeric, unsigned Len, unsigned Dim>
3370  inline const Numeric& ArrayND<Numeric, Len, Dim>::operator()(const unsigned i0,
3371  const unsigned i1,
3372  const unsigned i2,
3373  const unsigned i3,
3374  const unsigned i4,
3375  const unsigned i5,
3376  const unsigned i6) const {
3377  if (7U != dim_)
3378  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::operator(): wrong # of args (not rank 7 array)");
3379  return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3 * strides_[3] + i4 * strides_[4] +
3380  i5 * strides_[5] + i6];
3381  }
3382 
3383  template <typename Numeric, unsigned Len, unsigned Dim>
3384  inline const Numeric& ArrayND<Numeric, Len, Dim>::operator()(const unsigned i0,
3385  const unsigned i1,
3386  const unsigned i2,
3387  const unsigned i3,
3388  const unsigned i4,
3389  const unsigned i5,
3390  const unsigned i6,
3391  const unsigned i7) const {
3392  if (8U != dim_)
3393  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::operator(): wrong # of args (not rank 8 array)");
3394  return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3 * strides_[3] + i4 * strides_[4] +
3395  i5 * strides_[5] + i6 * strides_[6] + i7];
3396  }
3397 
3398  template <typename Numeric, unsigned Len, unsigned Dim>
3399  inline const Numeric& ArrayND<Numeric, Len, Dim>::operator()(const unsigned i0,
3400  const unsigned i1,
3401  const unsigned i2,
3402  const unsigned i3,
3403  const unsigned i4,
3404  const unsigned i5,
3405  const unsigned i6,
3406  const unsigned i7,
3407  const unsigned i8) const {
3408  if (9U != dim_)
3409  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::operator(): wrong # of args (not rank 9 array)");
3410  return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3 * strides_[3] + i4 * strides_[4] +
3411  i5 * strides_[5] + i6 * strides_[6] + i7 * strides_[7] + i8];
3412  }
3413 
3414  template <typename Numeric, unsigned Len, unsigned Dim>
3415  inline const Numeric& ArrayND<Numeric, Len, Dim>::operator()(const unsigned i0,
3416  const unsigned i1,
3417  const unsigned i2,
3418  const unsigned i3,
3419  const unsigned i4,
3420  const unsigned i5,
3421  const unsigned i6,
3422  const unsigned i7,
3423  const unsigned i8,
3424  const unsigned i9) const {
3425  if (10U != dim_)
3426  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::operator(): wrong # of args (not rank 10 array)");
3427  return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3 * strides_[3] + i4 * strides_[4] +
3428  i5 * strides_[5] + i6 * strides_[6] + i7 * strides_[7] + i8 * strides_[8] + i9];
3429  }
3430 
3431  template <typename Numeric, unsigned Len, unsigned Dim>
3432  inline Numeric& ArrayND<Numeric, Len, Dim>::operator()(const unsigned i0, const unsigned i1, const unsigned i2) {
3433  if (3U != dim_)
3434  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::operator(): wrong # of args (not rank 3 array)");
3435  return data_[i0 * strides_[0] + i1 * strides_[1] + i2];
3436  }
3437 
3438  template <typename Numeric, unsigned Len, unsigned Dim>
3439  inline Numeric& ArrayND<Numeric, Len, Dim>::operator()(const unsigned i0,
3440  const unsigned i1,
3441  const unsigned i2,
3442  const unsigned i3) {
3443  if (4U != dim_)
3444  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::operator(): wrong # of args (not rank 4 array)");
3445  return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3];
3446  }
3447 
3448  template <typename Numeric, unsigned Len, unsigned Dim>
3449  inline Numeric& ArrayND<Numeric, Len, Dim>::operator()(
3450  const unsigned i0, const unsigned i1, const unsigned i2, const unsigned i3, const unsigned i4) {
3451  if (5U != dim_)
3452  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::operator(): wrong # of args (not rank 5 array)");
3453  return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3 * strides_[3] + i4];
3454  }
3455 
3456  template <typename Numeric, unsigned Len, unsigned Dim>
3457  inline Numeric& ArrayND<Numeric, Len, Dim>::operator()(const unsigned i0,
3458  const unsigned i1,
3459  const unsigned i2,
3460  const unsigned i3,
3461  const unsigned i4,
3462  const unsigned i5) {
3463  if (6U != dim_)
3464  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::operator(): wrong # of args (not rank 6 array)");
3465  return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3 * strides_[3] + i4 * strides_[4] + i5];
3466  }
3467 
3468  template <typename Numeric, unsigned Len, unsigned Dim>
3469  inline Numeric& ArrayND<Numeric, Len, Dim>::operator()(const unsigned i0,
3470  const unsigned i1,
3471  const unsigned i2,
3472  const unsigned i3,
3473  const unsigned i4,
3474  const unsigned i5,
3475  const unsigned i6) {
3476  if (7U != dim_)
3477  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::operator(): wrong # of args (not rank 7 array)");
3478  return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3 * strides_[3] + i4 * strides_[4] +
3479  i5 * strides_[5] + i6];
3480  }
3481 
3482  template <typename Numeric, unsigned Len, unsigned Dim>
3483  inline Numeric& ArrayND<Numeric, Len, Dim>::operator()(const unsigned i0,
3484  const unsigned i1,
3485  const unsigned i2,
3486  const unsigned i3,
3487  const unsigned i4,
3488  const unsigned i5,
3489  const unsigned i6,
3490  const unsigned i7) {
3491  if (8U != dim_)
3492  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::operator(): wrong # of args (not rank 8 array)");
3493  return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3 * strides_[3] + i4 * strides_[4] +
3494  i5 * strides_[5] + i6 * strides_[6] + i7];
3495  }
3496 
3497  template <typename Numeric, unsigned Len, unsigned Dim>
3498  inline Numeric& ArrayND<Numeric, Len, Dim>::operator()(const unsigned i0,
3499  const unsigned i1,
3500  const unsigned i2,
3501  const unsigned i3,
3502  const unsigned i4,
3503  const unsigned i5,
3504  const unsigned i6,
3505  const unsigned i7,
3506  const unsigned i8) {
3507  if (9U != dim_)
3508  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::operator(): wrong # of args (not rank 9 array)");
3509  return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3 * strides_[3] + i4 * strides_[4] +
3510  i5 * strides_[5] + i6 * strides_[6] + i7 * strides_[7] + i8];
3511  }
3512 
3513  template <typename Numeric, unsigned Len, unsigned Dim>
3514  inline Numeric& ArrayND<Numeric, Len, Dim>::operator()(const unsigned i0,
3515  const unsigned i1,
3516  const unsigned i2,
3517  const unsigned i3,
3518  const unsigned i4,
3519  const unsigned i5,
3520  const unsigned i6,
3521  const unsigned i7,
3522  const unsigned i8,
3523  const unsigned i9) {
3524  if (10U != dim_)
3525  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::operator(): wrong # of args (not rank 10 array)");
3526  return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3 * strides_[3] + i4 * strides_[4] +
3527  i5 * strides_[5] + i6 * strides_[6] + i7 * strides_[7] + i8 * strides_[8] + i9];
3528  }
3529 
3530  template <typename Numeric, unsigned Len, unsigned Dim>
3531  const Numeric& ArrayND<Numeric, Len, Dim>::at(const unsigned i0, const unsigned i1, const unsigned i2) const {
3532  if (3U != dim_)
3533  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::at: wrong # of args (not rank 3 array)");
3534  if (i0 >= shape_[0])
3535  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 0 out of range (rank 3)");
3536  if (i1 >= shape_[1])
3537  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 1 out of range (rank 3)");
3538  if (i2 >= shape_[2])
3539  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 2 out of range (rank 3)");
3540  return data_[i0 * strides_[0] + i1 * strides_[1] + i2];
3541  }
3542 
3543  template <typename Numeric, unsigned Len, unsigned Dim>
3544  const Numeric& ArrayND<Numeric, Len, Dim>::at(const unsigned i0,
3545  const unsigned i1,
3546  const unsigned i2,
3547  const unsigned i3) const {
3548  if (4U != dim_)
3549  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::at: wrong # of args (not rank 4 array)");
3550  if (i0 >= shape_[0])
3551  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 0 out of range (rank 4)");
3552  if (i1 >= shape_[1])
3553  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 1 out of range (rank 4)");
3554  if (i2 >= shape_[2])
3555  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 2 out of range (rank 4)");
3556  if (i3 >= shape_[3])
3557  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 3 out of range (rank 4)");
3558  return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3];
3559  }
3560 
3561  template <typename Numeric, unsigned Len, unsigned Dim>
3562  const Numeric& ArrayND<Numeric, Len, Dim>::at(
3563  const unsigned i0, const unsigned i1, const unsigned i2, const unsigned i3, const unsigned i4) const {
3564  if (5U != dim_)
3565  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::at: wrong # of args (not rank 5 array)");
3566  if (i0 >= shape_[0])
3567  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 0 out of range (rank 5)");
3568  if (i1 >= shape_[1])
3569  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 1 out of range (rank 5)");
3570  if (i2 >= shape_[2])
3571  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 2 out of range (rank 5)");
3572  if (i3 >= shape_[3])
3573  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 3 out of range (rank 5)");
3574  if (i4 >= shape_[4])
3575  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 4 out of range (rank 5)");
3576  return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3 * strides_[3] + i4];
3577  }
3578 
3579  template <typename Numeric, unsigned Len, unsigned Dim>
3580  const Numeric& ArrayND<Numeric, Len, Dim>::at(const unsigned i0,
3581  const unsigned i1,
3582  const unsigned i2,
3583  const unsigned i3,
3584  const unsigned i4,
3585  const unsigned i5) const {
3586  if (6U != dim_)
3587  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::at: wrong # of args (not rank 6 array)");
3588  if (i0 >= shape_[0])
3589  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 0 out of range (rank 6)");
3590  if (i1 >= shape_[1])
3591  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 1 out of range (rank 6)");
3592  if (i2 >= shape_[2])
3593  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 2 out of range (rank 6)");
3594  if (i3 >= shape_[3])
3595  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 3 out of range (rank 6)");
3596  if (i4 >= shape_[4])
3597  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 4 out of range (rank 6)");
3598  if (i5 >= shape_[5])
3599  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 5 out of range (rank 6)");
3600  return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3 * strides_[3] + i4 * strides_[4] + i5];
3601  }
3602 
3603  template <typename Numeric, unsigned Len, unsigned Dim>
3604  const Numeric& ArrayND<Numeric, Len, Dim>::at(const unsigned i0,
3605  const unsigned i1,
3606  const unsigned i2,
3607  const unsigned i3,
3608  const unsigned i4,
3609  const unsigned i5,
3610  const unsigned i6) const {
3611  if (7U != dim_)
3612  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::at: wrong # of args (not rank 7 array)");
3613  if (i0 >= shape_[0])
3614  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 0 out of range (rank 7)");
3615  if (i1 >= shape_[1])
3616  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 1 out of range (rank 7)");
3617  if (i2 >= shape_[2])
3618  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 2 out of range (rank 7)");
3619  if (i3 >= shape_[3])
3620  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 3 out of range (rank 7)");
3621  if (i4 >= shape_[4])
3622  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 4 out of range (rank 7)");
3623  if (i5 >= shape_[5])
3624  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 5 out of range (rank 7)");
3625  if (i6 >= shape_[6])
3626  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 6 out of range (rank 7)");
3627  return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3 * strides_[3] + i4 * strides_[4] +
3628  i5 * strides_[5] + i6];
3629  }
3630 
3631  template <typename Numeric, unsigned Len, unsigned Dim>
3632  const Numeric& ArrayND<Numeric, Len, Dim>::at(const unsigned i0,
3633  const unsigned i1,
3634  const unsigned i2,
3635  const unsigned i3,
3636  const unsigned i4,
3637  const unsigned i5,
3638  const unsigned i6,
3639  const unsigned i7) const {
3640  if (8U != dim_)
3641  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::at: wrong # of args (not rank 8 array)");
3642  if (i0 >= shape_[0])
3643  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 0 out of range (rank 8)");
3644  if (i1 >= shape_[1])
3645  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 1 out of range (rank 8)");
3646  if (i2 >= shape_[2])
3647  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 2 out of range (rank 8)");
3648  if (i3 >= shape_[3])
3649  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 3 out of range (rank 8)");
3650  if (i4 >= shape_[4])
3651  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 4 out of range (rank 8)");
3652  if (i5 >= shape_[5])
3653  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 5 out of range (rank 8)");
3654  if (i6 >= shape_[6])
3655  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 6 out of range (rank 8)");
3656  if (i7 >= shape_[7])
3657  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 7 out of range (rank 8)");
3658  return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3 * strides_[3] + i4 * strides_[4] +
3659  i5 * strides_[5] + i6 * strides_[6] + i7];
3660  }
3661 
3662  template <typename Numeric, unsigned Len, unsigned Dim>
3663  const Numeric& ArrayND<Numeric, Len, Dim>::at(const unsigned i0,
3664  const unsigned i1,
3665  const unsigned i2,
3666  const unsigned i3,
3667  const unsigned i4,
3668  const unsigned i5,
3669  const unsigned i6,
3670  const unsigned i7,
3671  const unsigned i8) const {
3672  if (9U != dim_)
3673  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::at: wrong # of args (not rank 9 array)");
3674  if (i0 >= shape_[0])
3675  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 0 out of range (rank 9)");
3676  if (i1 >= shape_[1])
3677  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 1 out of range (rank 9)");
3678  if (i2 >= shape_[2])
3679  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 2 out of range (rank 9)");
3680  if (i3 >= shape_[3])
3681  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 3 out of range (rank 9)");
3682  if (i4 >= shape_[4])
3683  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 4 out of range (rank 9)");
3684  if (i5 >= shape_[5])
3685  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 5 out of range (rank 9)");
3686  if (i6 >= shape_[6])
3687  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 6 out of range (rank 9)");
3688  if (i7 >= shape_[7])
3689  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 7 out of range (rank 9)");
3690  if (i8 >= shape_[8])
3691  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 8 out of range (rank 9)");
3692  return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3 * strides_[3] + i4 * strides_[4] +
3693  i5 * strides_[5] + i6 * strides_[6] + i7 * strides_[7] + i8];
3694  }
3695 
3696  template <typename Numeric, unsigned Len, unsigned Dim>
3697  const Numeric& ArrayND<Numeric, Len, Dim>::at(const unsigned i0,
3698  const unsigned i1,
3699  const unsigned i2,
3700  const unsigned i3,
3701  const unsigned i4,
3702  const unsigned i5,
3703  const unsigned i6,
3704  const unsigned i7,
3705  const unsigned i8,
3706  const unsigned i9) const {
3707  if (10U != dim_)
3708  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::at: wrong # of args (not rank 10 array)");
3709  if (i0 >= shape_[0])
3710  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 0 out of range (rank 10)");
3711  if (i1 >= shape_[1])
3712  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 1 out of range (rank 10)");
3713  if (i2 >= shape_[2])
3714  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 2 out of range (rank 10)");
3715  if (i3 >= shape_[3])
3716  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 3 out of range (rank 10)");
3717  if (i4 >= shape_[4])
3718  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 4 out of range (rank 10)");
3719  if (i5 >= shape_[5])
3720  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 5 out of range (rank 10)");
3721  if (i6 >= shape_[6])
3722  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 6 out of range (rank 10)");
3723  if (i7 >= shape_[7])
3724  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 7 out of range (rank 10)");
3725  if (i8 >= shape_[8])
3726  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 8 out of range (rank 10)");
3727  if (i9 >= shape_[9])
3728  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 9 out of range (rank 10)");
3729  return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3 * strides_[3] + i4 * strides_[4] +
3730  i5 * strides_[5] + i6 * strides_[6] + i7 * strides_[7] + i8 * strides_[8] + i9];
3731  }
3732 
3733  template <typename Numeric, unsigned Len, unsigned Dim>
3734  Numeric& ArrayND<Numeric, Len, Dim>::at(const unsigned i0, const unsigned i1, const unsigned i2) {
3735  if (3U != dim_)
3736  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::at: wrong # of args (not rank 3 array)");
3737  if (i0 >= shape_[0])
3738  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 0 out of range (rank 3)");
3739  if (i1 >= shape_[1])
3740  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 1 out of range (rank 3)");
3741  if (i2 >= shape_[2])
3742  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 2 out of range (rank 3)");
3743  return data_[i0 * strides_[0] + i1 * strides_[1] + i2];
3744  }
3745 
3746  template <typename Numeric, unsigned Len, unsigned Dim>
3747  Numeric& ArrayND<Numeric, Len, Dim>::at(const unsigned i0, const unsigned i1, const unsigned i2, const unsigned i3) {
3748  if (4U != dim_)
3749  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::at: wrong # of args (not rank 4 array)");
3750  if (i0 >= shape_[0])
3751  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 0 out of range (rank 4)");
3752  if (i1 >= shape_[1])
3753  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 1 out of range (rank 4)");
3754  if (i2 >= shape_[2])
3755  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 2 out of range (rank 4)");
3756  if (i3 >= shape_[3])
3757  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 3 out of range (rank 4)");
3758  return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3];
3759  }
3760 
3761  template <typename Numeric, unsigned Len, unsigned Dim>
3762  Numeric& ArrayND<Numeric, Len, Dim>::at(
3763  const unsigned i0, const unsigned i1, const unsigned i2, const unsigned i3, const unsigned i4) {
3764  if (5U != dim_)
3765  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::at: wrong # of args (not rank 5 array)");
3766  if (i0 >= shape_[0])
3767  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 0 out of range (rank 5)");
3768  if (i1 >= shape_[1])
3769  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 1 out of range (rank 5)");
3770  if (i2 >= shape_[2])
3771  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 2 out of range (rank 5)");
3772  if (i3 >= shape_[3])
3773  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 3 out of range (rank 5)");
3774  if (i4 >= shape_[4])
3775  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 4 out of range (rank 5)");
3776  return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3 * strides_[3] + i4];
3777  }
3778 
3779  template <typename Numeric, unsigned Len, unsigned Dim>
3780  Numeric& ArrayND<Numeric, Len, Dim>::at(const unsigned i0,
3781  const unsigned i1,
3782  const unsigned i2,
3783  const unsigned i3,
3784  const unsigned i4,
3785  const unsigned i5) {
3786  if (6U != dim_)
3787  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::at: wrong # of args (not rank 6 array)");
3788  if (i0 >= shape_[0])
3789  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 0 out of range (rank 6)");
3790  if (i1 >= shape_[1])
3791  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 1 out of range (rank 6)");
3792  if (i2 >= shape_[2])
3793  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 2 out of range (rank 6)");
3794  if (i3 >= shape_[3])
3795  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 3 out of range (rank 6)");
3796  if (i4 >= shape_[4])
3797  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 4 out of range (rank 6)");
3798  if (i5 >= shape_[5])
3799  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 5 out of range (rank 6)");
3800  return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3 * strides_[3] + i4 * strides_[4] + i5];
3801  }
3802 
3803  template <typename Numeric, unsigned Len, unsigned Dim>
3804  Numeric& ArrayND<Numeric, Len, Dim>::at(const unsigned i0,
3805  const unsigned i1,
3806  const unsigned i2,
3807  const unsigned i3,
3808  const unsigned i4,
3809  const unsigned i5,
3810  const unsigned i6) {
3811  if (7U != dim_)
3812  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::at: wrong # of args (not rank 7 array)");
3813  if (i0 >= shape_[0])
3814  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 0 out of range (rank 7)");
3815  if (i1 >= shape_[1])
3816  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 1 out of range (rank 7)");
3817  if (i2 >= shape_[2])
3818  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 2 out of range (rank 7)");
3819  if (i3 >= shape_[3])
3820  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 3 out of range (rank 7)");
3821  if (i4 >= shape_[4])
3822  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 4 out of range (rank 7)");
3823  if (i5 >= shape_[5])
3824  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 5 out of range (rank 7)");
3825  if (i6 >= shape_[6])
3826  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 6 out of range (rank 7)");
3827  return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3 * strides_[3] + i4 * strides_[4] +
3828  i5 * strides_[5] + i6];
3829  }
3830 
3831  template <typename Numeric, unsigned Len, unsigned Dim>
3832  Numeric& ArrayND<Numeric, Len, Dim>::at(const unsigned i0,
3833  const unsigned i1,
3834  const unsigned i2,
3835  const unsigned i3,
3836  const unsigned i4,
3837  const unsigned i5,
3838  const unsigned i6,
3839  const unsigned i7) {
3840  if (8U != dim_)
3841  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::at: wrong # of args (not rank 8 array)");
3842  if (i0 >= shape_[0])
3843  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 0 out of range (rank 8)");
3844  if (i1 >= shape_[1])
3845  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 1 out of range (rank 8)");
3846  if (i2 >= shape_[2])
3847  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 2 out of range (rank 8)");
3848  if (i3 >= shape_[3])
3849  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 3 out of range (rank 8)");
3850  if (i4 >= shape_[4])
3851  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 4 out of range (rank 8)");
3852  if (i5 >= shape_[5])
3853  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 5 out of range (rank 8)");
3854  if (i6 >= shape_[6])
3855  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 6 out of range (rank 8)");
3856  if (i7 >= shape_[7])
3857  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 7 out of range (rank 8)");
3858  return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3 * strides_[3] + i4 * strides_[4] +
3859  i5 * strides_[5] + i6 * strides_[6] + i7];
3860  }
3861 
3862  template <typename Numeric, unsigned Len, unsigned Dim>
3863  Numeric& ArrayND<Numeric, Len, Dim>::at(const unsigned i0,
3864  const unsigned i1,
3865  const unsigned i2,
3866  const unsigned i3,
3867  const unsigned i4,
3868  const unsigned i5,
3869  const unsigned i6,
3870  const unsigned i7,
3871  const unsigned i8) {
3872  if (9U != dim_)
3873  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::at: wrong # of args (not rank 9 array)");
3874  if (i0 >= shape_[0])
3875  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 0 out of range (rank 9)");
3876  if (i1 >= shape_[1])
3877  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 1 out of range (rank 9)");
3878  if (i2 >= shape_[2])
3879  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 2 out of range (rank 9)");
3880  if (i3 >= shape_[3])
3881  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 3 out of range (rank 9)");
3882  if (i4 >= shape_[4])
3883  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 4 out of range (rank 9)");
3884  if (i5 >= shape_[5])
3885  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 5 out of range (rank 9)");
3886  if (i6 >= shape_[6])
3887  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 6 out of range (rank 9)");
3888  if (i7 >= shape_[7])
3889  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 7 out of range (rank 9)");
3890  if (i8 >= shape_[8])
3891  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 8 out of range (rank 9)");
3892  return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3 * strides_[3] + i4 * strides_[4] +
3893  i5 * strides_[5] + i6 * strides_[6] + i7 * strides_[7] + i8];
3894  }
3895 
3896  template <typename Numeric, unsigned Len, unsigned Dim>
3897  Numeric& ArrayND<Numeric, Len, Dim>::at(const unsigned i0,
3898  const unsigned i1,
3899  const unsigned i2,
3900  const unsigned i3,
3901  const unsigned i4,
3902  const unsigned i5,
3903  const unsigned i6,
3904  const unsigned i7,
3905  const unsigned i8,
3906  const unsigned i9) {
3907  if (10U != dim_)
3908  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::at: wrong # of args (not rank 10 array)");
3909  if (i0 >= shape_[0])
3910  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 0 out of range (rank 10)");
3911  if (i1 >= shape_[1])
3912  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 1 out of range (rank 10)");
3913  if (i2 >= shape_[2])
3914  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 2 out of range (rank 10)");
3915  if (i3 >= shape_[3])
3916  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 3 out of range (rank 10)");
3917  if (i4 >= shape_[4])
3918  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 4 out of range (rank 10)");
3919  if (i5 >= shape_[5])
3920  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 5 out of range (rank 10)");
3921  if (i6 >= shape_[6])
3922  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 6 out of range (rank 10)");
3923  if (i7 >= shape_[7])
3924  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 7 out of range (rank 10)");
3925  if (i8 >= shape_[8])
3926  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 8 out of range (rank 10)");
3927  if (i9 >= shape_[9])
3928  throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 9 out of range (rank 10)");
3929  return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3 * strides_[3] + i4 * strides_[4] +
3930  i5 * strides_[5] + i6 * strides_[6] + i7 * strides_[7] + i8 * strides_[8] + i9];
3931  }
3932 
3933  template <typename Numeric, unsigned Len, unsigned Dim>
3934  template <unsigned Len2, unsigned Dim2>
3935  double ArrayND<Numeric, Len, Dim>::maxAbsDifference(const ArrayND<Numeric, Len2, Dim2>& r) const {
3936  if (!isShapeCompatible(r))
3937  throw npstat::NpstatInvalidArgument(
3938  "In npstat::ArrayND::maxAbsDifference: "
3939  "incompatible argument array shape");
3940  if (dim_) {
3941  double maxd = 0.0;
3942  for (unsigned long i = 0; i < len_; ++i) {
3943  const Numeric rval = r.data_[i];
3944  const double d = absDifference(data_[i], rval);
3945  if (d > maxd)
3946  maxd = d;
3947  }
3948  return maxd;
3949  } else {
3950  const Numeric rval = r.localData_[0];
3951  return absDifference(localData_[0], rval);
3952  }
3953  }
3954 
3955  template <typename Numeric, unsigned Len, unsigned Dim>
3956  template <unsigned Len2, unsigned Dim2>
3957  bool ArrayND<Numeric, Len, Dim>::operator==(const ArrayND<Numeric, Len2, Dim2>& r) const {
3958  if (shapeIsKnown_ != r.shapeIsKnown_)
3959  return false;
3960  if (r.dim_ != dim_)
3961  return false;
3962  if (r.len_ != len_)
3963  return false;
3964  for (unsigned i = 0; i < dim_; ++i)
3965  if (shape_[i] != r.shape_[i])
3966  return false;
3967  for (unsigned i = 0; i < dim_; ++i)
3968  assert(strides_[i] == r.strides_[i]);
3969  for (unsigned long j = 0; j < len_; ++j)
3970  if (data_[j] != r.data_[j])
3971  return false;
3972  return true;
3973  }
3974 
3975  template <typename Numeric, unsigned Len, unsigned Dim>
3976  template <unsigned Len2, unsigned Dim2>
3977  inline bool ArrayND<Numeric, Len, Dim>::operator!=(const ArrayND<Numeric, Len2, Dim2>& r) const {
3978  return !(*this == r);
3979  }
3980 
3981  template <typename Numeric, unsigned Len, unsigned Dim>
3982  template <typename Num2>
3983  ArrayND<Numeric, Len, Dim> ArrayND<Numeric, Len, Dim>::operator*(const Num2& r) const {
3984  if (!shapeIsKnown_)
3985  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"operator*\"");
3986  ArrayND<Numeric, Len, Dim> result(shape_, dim_);
3987  for (unsigned long i = 0; i < len_; ++i)
3988  result.data_[i] = data_[i] * r;
3989  return result;
3990  }
3991 
3992  template <typename Numeric, unsigned Len, unsigned Dim>
3993  template <typename Num2>
3994  ArrayND<Numeric, Len, Dim> ArrayND<Numeric, Len, Dim>::operator/(const Num2& r) const {
3995  if (!shapeIsKnown_)
3996  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"operator/\"");
3997  if (r == Num2())
3998  throw npstat::NpstatRuntimeError("In npstat::ArrayND::operator/: division by zero");
3999  ArrayND<Numeric, Len, Dim> result(shape_, dim_);
4000  for (unsigned long i = 0; i < len_; ++i)
4001  result.data_[i] = data_[i] / r;
4002  return result;
4003  }
4004 
4005  template <typename Numeric, unsigned Len, unsigned Dim>
4006  template <unsigned Len2, unsigned Dim2>
4007  ArrayND<Numeric, Len, Dim> ArrayND<Numeric, Len, Dim>::operator+(const ArrayND<Numeric, Len2, Dim2>& r) const {
4008  if (!isShapeCompatible(r))
4009  throw npstat::NpstatInvalidArgument(
4010  "In npstat::ArrayND::operator+: "
4011  "incompatible argument array shape");
4012  ArrayND<Numeric, Len, Dim> result(shape_, dim_);
4013  for (unsigned long i = 0; i < len_; ++i)
4014  result.data_[i] = data_[i] + r.data_[i];
4015  return result;
4016  }
4017 
4018  template <typename Numeric, unsigned Len, unsigned Dim>
4019  template <unsigned Len2, unsigned Dim2>
4020  ArrayND<Numeric, Len, Dim> ArrayND<Numeric, Len, Dim>::operator-(const ArrayND<Numeric, Len2, Dim2>& r) const {
4021  if (!isShapeCompatible(r))
4022  throw npstat::NpstatInvalidArgument(
4023  "In npstat::ArrayND::operator-: "
4024  "incompatible argument array shape");
4025  ArrayND<Numeric, Len, Dim> result(shape_, dim_);
4026  for (unsigned long i = 0; i < len_; ++i)
4027  result.data_[i] = data_[i] - r.data_[i];
4028  return result;
4029  }
4030 
4031  template <typename Numeric, unsigned Len, unsigned Dim>
4032  inline ArrayND<Numeric, Len, Dim> ArrayND<Numeric, Len, Dim>::operator+() const {
4033  if (!shapeIsKnown_)
4034  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"operator+\"");
4035  return *this;
4036  }
4037 
4038  template <typename Numeric, unsigned Len, unsigned Dim>
4039  ArrayND<Numeric, Len, Dim> ArrayND<Numeric, Len, Dim>::operator-() const {
4040  if (!shapeIsKnown_)
4041  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"operator-\"");
4042  ArrayND<Numeric, Len, Dim> result(shape_, dim_);
4043  for (unsigned long i = 0; i < len_; ++i)
4044  result.data_[i] = -data_[i];
4045  return result;
4046  }
4047 
4048  template <typename Numeric, unsigned Len, unsigned Dim>
4049  template <typename Num2>
4050  ArrayND<Numeric, Len, Dim>& ArrayND<Numeric, Len, Dim>::operator*=(const Num2& r) {
4051  if (!shapeIsKnown_)
4052  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"operator*=\"");
4053  for (unsigned long i = 0; i < len_; ++i)
4054  data_[i] *= r;
4055  return *this;
4056  }
4057 
4058  template <typename Numeric, unsigned Len, unsigned Dim>
4059  ArrayND<Numeric, Len, Dim>& ArrayND<Numeric, Len, Dim>::makeNonNegative() {
4060  if (!shapeIsKnown_)
4061  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"makeNonNegative\"");
4062  const Numeric zero = Numeric();
4063  if (dim_) {
4064  for (unsigned long i = 0; i < len_; ++i)
4065  if (!(ComplexComparesAbs<Numeric>::more(data_[i], zero)))
4066  data_[i] = zero;
4067  } else if (!(ComplexComparesAbs<Numeric>::more(localData_[0], zero)))
4068  localData_[0] = zero;
4069  return *this;
4070  }
4071 
4072  template <typename Numeric, unsigned Len, unsigned Dim>
4073  unsigned ArrayND<Numeric, Len, Dim>::makeCopulaSteps(const double tolerance, const unsigned nCycles) {
4074  if (!shapeIsKnown_)
4075  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"makeCopulaSteps\"");
4076  if (nCycles == 0U)
4077  return 0U;
4078  if (!dim_)
4079  throw npstat::NpstatInvalidArgument(
4080  "npstat::ArrayND::makeCopulaSteps method "
4081  "can not be used with array of 0 rank");
4082 
4083  const Numeric zero = Numeric();
4084  for (unsigned long i = 0; i < len_; ++i)
4085  if (!(ComplexComparesAbs<Numeric>::more(data_[i], zero)))
4086  data_[i] = zero;
4087 
4088  std::vector<Numeric*> axesPtrBuf(dim_);
4089  Numeric** axes = &axesPtrBuf[0];
4090  const Numeric one = static_cast<Numeric>(1);
4091 
4092  // Memory for the axis accumulators
4093  unsigned idxSum = 0;
4094  for (unsigned i = 0; i < dim_; ++i)
4095  idxSum += shape_[i];
4096  std::vector<Numeric> axesBuf(idxSum);
4097  axes[0] = &axesBuf[0];
4098  for (unsigned i = 1; i < dim_; ++i)
4099  axes[i] = axes[i - 1] + shape_[i - 1];
4100 
4101  // Accumulate axis projections
4102  unsigned icycle = 0;
4103  for (; icycle < nCycles; ++icycle) {
4104  for (unsigned i = 0; i < idxSum; ++i)
4105  axesBuf[i] = zero;
4106 
4107  // Accumulate sums for each axis
4108  for (unsigned long idat = 0; idat < len_; ++idat) {
4109  unsigned long l = idat;
4110  for (unsigned i = 0; i < dim_; ++i) {
4111  const unsigned idx = l / strides_[i];
4112  l -= (idx * strides_[i]);
4113  axes[i][idx] += data_[idat];
4114  }
4115  }
4116 
4117  // Make averages out of sums
4118  bool withinTolerance = true;
4119  Numeric totalSum = zero;
4120  for (unsigned i = 0; i < dim_; ++i) {
4121  Numeric axisSum = zero;
4122  const unsigned amax = shape_[i];
4123  for (unsigned a = 0; a < amax; ++a) {
4124  if (axes[i][a] == zero)
4125  throw npstat::NpstatRuntimeError(
4126  "In npstat::ArrayND::makeCopulaSteps: "
4127  "marginal density is zero");
4128  axisSum += axes[i][a];
4129  }
4130  totalSum += axisSum;
4131  const Numeric axisAverage = axisSum / static_cast<Numeric>(amax);
4132  for (unsigned a = 0; a < amax; ++a)
4133  axes[i][a] /= axisAverage;
4134  for (unsigned a = 0; a < amax && withinTolerance; ++a) {
4135  const double adelta = absDifference(axes[i][a], one);
4136  if (adelta > tolerance)
4137  withinTolerance = false;
4138  }
4139  }
4140 
4141  if (withinTolerance)
4142  break;
4143 
4144  const Numeric totalAverage = totalSum / static_cast<Numeric>(len_) / static_cast<Numeric>(dim_);
4145 
4146  // Run over all points again and divide by
4147  // the product of marginals
4148  for (unsigned long idat = 0; idat < len_; ++idat) {
4149  unsigned long l = idat;
4150  for (unsigned i = 0; i < dim_; ++i) {
4151  const unsigned idx = l / strides_[i];
4152  l -= (idx * strides_[i]);
4153  data_[idat] /= axes[i][idx];
4154  }
4155  data_[idat] /= totalAverage;
4156  }
4157  }
4158 
4159  return icycle;
4160  }
4161 
4162  template <typename Numeric, unsigned Len, unsigned Dim>
4163  template <typename Num2>
4164  ArrayND<Numeric, Len, Dim>& ArrayND<Numeric, Len, Dim>::operator/=(const Num2& r) {
4165  if (!shapeIsKnown_)
4166  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"operator/=\"");
4167  if (r == Num2())
4168  throw npstat::NpstatRuntimeError("In npstat::ArrayND::operator/=: division by zero");
4169  for (unsigned long i = 0; i < len_; ++i)
4170  data_[i] /= r;
4171  return *this;
4172  }
4173 
4174  template <typename Numeric, unsigned Len, unsigned Dim>
4175  template <typename Num2, unsigned Len2, unsigned Dim2>
4176  ArrayND<Numeric, Len, Dim>& ArrayND<Numeric, Len, Dim>::operator+=(const ArrayND<Num2, Len2, Dim2>& r) {
4177  if (!isShapeCompatible(r))
4178  throw npstat::NpstatInvalidArgument(
4179  "In npstat::ArrayND::operator+=: "
4180  "incompatible argument array shape");
4181  for (unsigned long i = 0; i < len_; ++i)
4182  data_[i] += r.data_[i];
4183  return *this;
4184  }
4185 
4186  template <typename Numeric, unsigned Len, unsigned Dim>
4187  template <typename Num3, typename Num2, unsigned Len2, unsigned Dim2>
4188  ArrayND<Numeric, Len, Dim>& ArrayND<Numeric, Len, Dim>::addmul(const ArrayND<Num2, Len2, Dim2>& r, const Num3& c) {
4189  if (!isShapeCompatible(r))
4190  throw npstat::NpstatInvalidArgument(
4191  "In npstat::ArrayND::addmul: "
4192  "incompatible argument array shape");
4193  for (unsigned long i = 0; i < len_; ++i)
4194  data_[i] += r.data_[i] * c;
4195  return *this;
4196  }
4197 
4198  template <typename Numeric, unsigned Len, unsigned Dim>
4199  template <typename Num2, unsigned Len2, unsigned Dim2>
4200  ArrayND<Numeric, Len, Dim>& ArrayND<Numeric, Len, Dim>::operator-=(const ArrayND<Num2, Len2, Dim2>& r) {
4201  if (!isShapeCompatible(r))
4202  throw npstat::NpstatInvalidArgument(
4203  "In npstat::ArrayND::operator-=: "
4204  "incompatible argument array shape");
4205  for (unsigned long i = 0; i < len_; ++i)
4206  data_[i] -= r.data_[i];
4207  return *this;
4208  }
4209 
4210  template <typename Numeric, unsigned Len, unsigned Dim>
4211  Numeric ArrayND<Numeric, Len, Dim>::interpolate1(const double* coords, const unsigned dim) const {
4212  if (!shapeIsKnown_)
4213  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"interpolate1\"");
4214  if (dim != dim_)
4215  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::interpolate1: incompatible coordinate length");
4216  if (dim) {
4217  const unsigned maxdim = CHAR_BIT * sizeof(unsigned long);
4218  if (dim_ >= maxdim)
4219  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::interpolate1: array rank is too large");
4220 
4221  double dx[maxdim];
4222  unsigned ix[maxdim];
4223  for (unsigned i = 0; i < dim; ++i) {
4224  const double x = coords[i];
4225  if (x <= 0.0) {
4226  ix[i] = 0;
4227  dx[i] = 0.0;
4228  } else if (x >= static_cast<double>(shape_[i] - 1)) {
4229  ix[i] = shape_[i] - 1;
4230  dx[i] = 0.0;
4231  } else {
4232  ix[i] = static_cast<unsigned>(std::floor(x));
4233  dx[i] = x - ix[i];
4234  }
4235  }
4236 
4237  Numeric sum = Numeric();
4238  const unsigned long maxcycle = 1UL << dim;
4239  for (unsigned long icycle = 0UL; icycle < maxcycle; ++icycle) {
4240  double w = 1.0;
4241  unsigned long icell = 0UL;
4242  for (unsigned i = 0; i < dim; ++i) {
4243  if (icycle & (1UL << i)) {
4244  w *= dx[i];
4245  icell += strides_[i] * (ix[i] + 1U);
4246  } else {
4247  w *= (1.0 - dx[i]);
4248  icell += strides_[i] * ix[i];
4249  }
4250  }
4251  if (w > 0.0)
4252  sum += data_[icell] * static_cast<proper_double>(w);
4253  }
4254  return sum;
4255  } else
4256  return localData_[0];
4257  }
4258 
4259  template <typename Numeric, unsigned Len, unsigned Dim>
4260  Numeric ArrayND<Numeric, Len, Dim>::interpolateLoop(const unsigned level,
4261  const double* coords,
4262  const Numeric* base) const {
4263  const unsigned npoints = shape_[level];
4264  const double x = coords[level];
4265 
4266  unsigned ix, npt = 1;
4267  double dx = 0.0;
4268  if (x < 0.0)
4269  ix = 0;
4270  else if (x > static_cast<double>(npoints - 1))
4271  ix = npoints - 1;
4272  else {
4273  ix = static_cast<unsigned>(std::floor(x));
4274  if (ix)
4275  --ix;
4276  unsigned imax = ix + 3;
4277  while (imax >= npoints) {
4278  if (ix)
4279  --ix;
4280  --imax;
4281  }
4282  dx = x - ix;
4283  npt = imax + 1 - ix;
4284  }
4285  assert(npt >= 1 && npt <= 4);
4286 
4287  Numeric fit[4];
4288  if (level < dim_ - 1)
4289  for (unsigned ipt = 0; ipt < npt; ++ipt)
4290  fit[ipt] = interpolateLoop(level + 1, coords, base + (ix + ipt) * strides_[level]);
4291 
4292  const Numeric* const v = (level == dim_ - 1 ? base + ix : fit);
4293  switch (npt) {
4294  case 1:
4295  return v[0];
4296  case 2:
4297  return interpolate_linear(dx, v[0], v[1]);
4298  case 3:
4299  return interpolate_quadratic(dx, v[0], v[1], v[2]);
4300  case 4:
4301  return interpolate_cubic(dx, v[0], v[1], v[2], v[3]);
4302  default:
4303  assert(0);
4304  return Numeric();
4305  }
4306  }
4307 
4308  template <typename Numeric, unsigned Len, unsigned Dim>
4309  inline Numeric ArrayND<Numeric, Len, Dim>::interpolate3(const double* coords, const unsigned dim) const {
4310  if (!shapeIsKnown_)
4311  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"interpolate3\"");
4312  if (dim != dim_)
4313  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::interpolate3: incompatible coordinate length");
4314  if (dim) {
4315  assert(coords);
4316  return interpolateLoop(0, coords, data_);
4317  } else
4318  return localData_[0];
4319  }
4320 
4321  template <typename Numeric, unsigned Len, unsigned Dim>
4322  template <class Functor>
4323  ArrayND<Numeric, Len, Dim>& ArrayND<Numeric, Len, Dim>::apply(Functor f) {
4324  if (!shapeIsKnown_)
4325  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"apply\"");
4326  for (unsigned long i = 0; i < len_; ++i)
4327  data_[i] = static_cast<Numeric>(f(data_[i]));
4328  return *this;
4329  }
4330 
4331  template <typename Numeric, unsigned Len, unsigned Dim>
4332  template <class Functor>
4333  ArrayND<Numeric, Len, Dim>& ArrayND<Numeric, Len, Dim>::scanInPlace(Functor f) {
4334  if (!shapeIsKnown_)
4335  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"scanInPlace\"");
4336  for (unsigned long i = 0; i < len_; ++i)
4337  f(data_[i]);
4338  return *this;
4339  }
4340 
4341  template <typename Numeric, unsigned Len, unsigned Dim>
4342  ArrayND<Numeric, Len, Dim>& ArrayND<Numeric, Len, Dim>::constFill(const Numeric c) {
4343  if (!shapeIsKnown_)
4344  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"constFill\"");
4345  for (unsigned long i = 0; i < len_; ++i)
4346  data_[i] = c;
4347  return *this;
4348  }
4349 
4350  template <typename Numeric, unsigned Len, unsigned Dim>
4351  inline ArrayND<Numeric, Len, Dim>& ArrayND<Numeric, Len, Dim>::clear() {
4352  return constFill(Numeric());
4353  }
4354 
4355  template <typename Numeric, unsigned Len, unsigned Dim>
4356  ArrayND<Numeric, Len, Dim>& ArrayND<Numeric, Len, Dim>::uninitialize() {
4357  destroyBuffer(data_, localData_);
4358  destroyBuffer(strides_, localStrides_);
4359  destroyBuffer(shape_, localShape_);
4360  localData_[0] = Numeric();
4361  data_ = localData_;
4362  strides_ = nullptr;
4363  shape_ = nullptr;
4364  len_ = 0;
4365  dim_ = 0;
4366  shapeIsKnown_ = false;
4367  return *this;
4368  }
4369 
4370  template <typename Numeric, unsigned Len, unsigned Dim>
4371  ArrayND<Numeric, Len, Dim>& ArrayND<Numeric, Len, Dim>::makeUnit() {
4372  if (!shapeIsKnown_)
4373  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"makeUnit\"");
4374  if (dim_ < 2)
4375  throw npstat::NpstatInvalidArgument(
4376  "npstat::ArrayND::makeUnit method "
4377  "can not be used with arrays of rank less than 2");
4378  constFill(Numeric());
4379  unsigned long stride = 0UL;
4380  const unsigned dimlen = shape_[0];
4381  for (unsigned i = 0; i < dim_; ++i) {
4382  if (shape_[i] != dimlen)
4383  throw npstat::NpstatInvalidArgument(
4384  "npstat::ArrayND::makeUnit method needs "
4385  "the array span to be the same in ech dimension");
4386  stride += strides_[i];
4387  }
4388  const Numeric one(static_cast<Numeric>(1));
4389  for (unsigned i = 0; i < dimlen; ++i)
4390  data_[i * stride] = one;
4391  return *this;
4392  }
4393 
4394  template <typename Numeric, unsigned Len, unsigned Dim>
4395  void ArrayND<Numeric, Len, Dim>::linearFillLoop(
4396  const unsigned level, const double s0, const unsigned long idx, const double shift, const double* coeffs) {
4397  const unsigned imax = shape_[level];
4398  const double c = coeffs[level];
4399  if (level == dim_ - 1) {
4400  Numeric* d = &data_[idx];
4401  for (unsigned i = 0; i < imax; ++i) {
4402  // Note that we want to add "shift" only at the
4403  // very end. This might improve the numerical
4404  // precision of the result.
4405  const double sum = s0 + c * i + shift;
4406  d[i] = static_cast<Numeric>(sum);
4407  }
4408  } else {
4409  const unsigned long stride = strides_[level];
4410  for (unsigned i = 0; i < imax; ++i)
4411  linearFillLoop(level + 1, s0 + c * i, idx + i * stride, shift, coeffs);
4412  }
4413  }
4414 
4415  template <typename Numeric, unsigned Len, unsigned Dim>
4416  ArrayND<Numeric, Len, Dim>& ArrayND<Numeric, Len, Dim>::linearFill(const double* coeffs,
4417  const unsigned dimCoeffs,
4418  const double shift) {
4419  // Make sure the object has been initialized
4420  if (!shapeIsKnown_)
4421  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"linearFill\"");
4422  if (dim_ != dimCoeffs)
4423  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::linearFill: incompatible number of coefficients");
4424  if (dim_) {
4425  assert(coeffs);
4426  linearFillLoop(0U, 0.0, 0UL, shift, coeffs);
4427  } else
4428  localData_[0] = static_cast<Numeric>(shift);
4429  return *this;
4430  }
4431 
4432  template <typename Numeric, unsigned Len, unsigned Dim>
4433  template <class Functor>
4434  void ArrayND<Numeric, Len, Dim>::functorFillLoop(const unsigned level,
4435  const unsigned long idx,
4436  Functor f,
4437  unsigned* farg) {
4438  const unsigned imax = shape_[level];
4439  if (level == dim_ - 1) {
4440  Numeric* d = &data_[idx];
4441  const unsigned* myarg = farg;
4442  for (unsigned i = 0; i < imax; ++i) {
4443  farg[level] = i;
4444  d[i] = static_cast<Numeric>(f(myarg, dim_));
4445  }
4446  } else {
4447  const unsigned long stride = strides_[level];
4448  for (unsigned i = 0; i < imax; ++i) {
4449  farg[level] = i;
4450  functorFillLoop(level + 1, idx + i * stride, f, farg);
4451  }
4452  }
4453  }
4454 
4455  template <typename Numeric, unsigned Len, unsigned Dim>
4456  template <typename Accumulator>
4457  void ArrayND<Numeric, Len, Dim>::convertToLastDimCdfLoop(
4458  ArrayND* sumSlice, const unsigned level, unsigned long idx0, unsigned long idxSlice, const bool useTrapezoids) {
4459  static const proper_double half = 0.5;
4460  const unsigned imax = shape_[level];
4461  if (level == dim_ - 1) {
4462  Accumulator acc = Accumulator();
4463  Numeric* data = data_ + idx0;
4464  if (useTrapezoids) {
4465  Numeric oldval = Numeric();
4466  for (unsigned i = 0; i < imax; ++i) {
4467  acc += (data[i] + oldval) * half;
4468  oldval = data[i];
4469  data[i] = static_cast<Numeric>(acc);
4470  }
4471  acc += oldval * half;
4472  } else
4473  for (unsigned i = 0; i < imax; ++i) {
4474  acc += data[i];
4475  data[i] = static_cast<Numeric>(acc);
4476  }
4477  if (sumSlice->dim_)
4478  sumSlice->data_[idxSlice] = static_cast<Numeric>(acc);
4479  else
4480  sumSlice->localData_[0] = static_cast<Numeric>(acc);
4481  } else {
4482  const unsigned long stride = strides_[level];
4483  unsigned long sumStride = 0UL;
4484  if (sumSlice->dim_)
4485  sumStride = sumSlice->strides_[level];
4486  for (unsigned i = 0; i < imax; ++i) {
4487  convertToLastDimCdfLoop<Accumulator>(sumSlice, level + 1, idx0, idxSlice, useTrapezoids);
4488  idx0 += stride;
4489  idxSlice += sumStride;
4490  }
4491  }
4492  }
4493 
4494  template <typename Numeric, unsigned Len, unsigned Dim>
4495  template <typename Accumulator>
4496  inline void ArrayND<Numeric, Len, Dim>::convertToLastDimCdf(ArrayND* sumSlice, const bool useTrapezoids) {
4497  if (!shapeIsKnown_)
4498  throw npstat::NpstatInvalidArgument(
4499  "Initialize npstat::ArrayND before calling "
4500  "method \"convertToLastDimCdf\"");
4501  if (!dim_)
4502  throw npstat::NpstatInvalidArgument(
4503  "npstat::ArrayND::convertToLastDimCdf method "
4504  "can not be used with array of 0 rank");
4505  assert(sumSlice);
4506  if (!sumSlice->shapeIsKnown_)
4507  throw npstat::NpstatInvalidArgument(
4508  "In npstat::ArrayND::convertToLastDimCdf: "
4509  "uninitialized argument array");
4510  convertToLastDimCdfLoop<Accumulator>(sumSlice, 0U, 0UL, 0UL, useTrapezoids);
4511  }
4512 
4513  template <typename Numeric, unsigned Len, unsigned Dim>
4514  template <class Functor>
4515  ArrayND<Numeric, Len, Dim>& ArrayND<Numeric, Len, Dim>::functorFill(Functor f) {
4516  if (!shapeIsKnown_)
4517  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"functorFill\"");
4518  if (dim_) {
4519  unsigned localIndex[Dim];
4520  unsigned* index = makeBuffer(dim_, localIndex, Dim);
4521  functorFillLoop(0U, 0UL, f, index);
4522  destroyBuffer(index, localIndex);
4523  } else
4524  localData_[0] = static_cast<Numeric>(f(static_cast<unsigned*>(nullptr), 0U));
4525  return *this;
4526  }
4527 
4528  template <typename Numeric, unsigned Len, unsigned Dim>
4529  template <typename Num2, unsigned Len2, unsigned Dim2>
4530  bool ArrayND<Numeric, Len, Dim>::isClose(const ArrayND<Num2, Len2, Dim2>& r, const double eps) const {
4531  if (eps < 0.0)
4532  throw npstat::NpstatDomainError("In npstat::ArrayND::isClose: tolerance must not be negative");
4533  if (!isShapeCompatible(r))
4534  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::isClose: incompatible argument array shape");
4535  if (dim_) {
4536  for (unsigned long i = 0; i < len_; ++i) {
4537  const Numeric rval = r.data_[i];
4538  if (static_cast<double>(absDifference(data_[i], rval)) > eps)
4539  return false;
4540  }
4541  } else {
4542  const Numeric rval = r.localData_[0];
4543  if (static_cast<double>(absDifference(localData_[0], rval)) > eps)
4544  return false;
4545  }
4546  return true;
4547  }
4548 
4549  template <typename Numeric, unsigned Len, unsigned Dim>
4550  template <typename Num2, unsigned Len2, unsigned Dim2>
4551  ArrayND<Numeric, Len, Dim> ArrayND<Numeric, Len, Dim>::outer(const ArrayND<Num2, Len2, Dim2>& r) const {
4552  return ArrayND<Numeric, Len, Dim>(*this, r);
4553  }
4554 
4555  template <typename Numeric, unsigned Len, unsigned Dim>
4556  void ArrayND<Numeric, Len, Dim>::contractLoop(unsigned thisLevel,
4557  const unsigned resLevel,
4558  const unsigned pos1,
4559  const unsigned pos2,
4560  unsigned long idxThis,
4561  unsigned long idxRes,
4562  ArrayND& result) const {
4563  while (thisLevel == pos1 || thisLevel == pos2)
4564  ++thisLevel;
4565  assert(thisLevel < dim_);
4566 
4567  if (resLevel == result.dim_ - 1) {
4568  const unsigned ncontract = shape_[pos1];
4569  const unsigned imax = result.shape_[resLevel];
4570  const unsigned long stride = strides_[pos1] + strides_[pos2];
4571  for (unsigned i = 0; i < imax; ++i) {
4572  const Numeric* tmp = data_ + (idxThis + i * strides_[thisLevel]);
4573  Numeric sum = Numeric();
4574  for (unsigned j = 0; j < ncontract; ++j)
4575  sum += tmp[j * stride];
4576  result.data_[idxRes + i] = sum;
4577  }
4578  } else {
4579  const unsigned imax = result.shape_[resLevel];
4580  assert(imax == shape_[thisLevel]);
4581  for (unsigned i = 0; i < imax; ++i) {
4582  contractLoop(thisLevel + 1, resLevel + 1, pos1, pos2, idxThis, idxRes, result);
4583  idxThis += strides_[thisLevel];
4584  idxRes += result.strides_[resLevel];
4585  }
4586  }
4587  }
4588 
4589  template <typename Numeric, unsigned Len, unsigned Dim>
4590  ArrayND<Numeric, Len, Dim> ArrayND<Numeric, Len, Dim>::contract(const unsigned pos1, const unsigned pos2) const {
4591  if (!shapeIsKnown_)
4592  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"contract\"");
4593  if (!(pos1 < dim_ && pos2 < dim_ && pos1 != pos2))
4594  throw npstat::NpstatInvalidArgument(
4595  "In npstat::ArrayND::contract: "
4596  "incompatible contraction indices");
4597  if (shape_[pos1] != shape_[pos2])
4598  throw npstat::NpstatInvalidArgument(
4599  "In npstat::ArrayND::contract: incompatible "
4600  "length of contracted dimensions");
4601 
4602  // Construct the new shape
4603  unsigned newshapeBuf[Dim];
4604  unsigned* newshape = makeBuffer(dim_ - 2, newshapeBuf, Dim);
4605  unsigned ishap = 0;
4606  for (unsigned i = 0; i < dim_; ++i)
4607  if (i != pos1 && i != pos2)
4608  newshape[ishap++] = shape_[i];
4609 
4610  // Form the result array
4611  ArrayND<Numeric, Len, Dim> result(newshape, ishap);
4612  if (ishap)
4613  contractLoop(0, 0, pos1, pos2, 0UL, 0UL, result);
4614  else {
4615  // We are just calculating the trace
4616  Numeric sum = Numeric();
4617  const unsigned imax = shape_[0];
4618  const unsigned long stride = strides_[0] + strides_[1];
4619  for (unsigned i = 0; i < imax; ++i)
4620  sum += data_[i * stride];
4621  result() = sum;
4622  }
4623 
4624  destroyBuffer(newshape, newshapeBuf);
4625  return result;
4626  }
4627 
4628  template <typename Numeric, unsigned Len, unsigned Dim>
4629  void ArrayND<Numeric, Len, Dim>::transposeLoop(const unsigned level,
4630  const unsigned pos1,
4631  const unsigned pos2,
4632  unsigned long idxThis,
4633  unsigned long idxRes,
4634  ArrayND& result) const {
4635  const unsigned imax = shape_[level];
4636  const unsigned long mystride = strides_[level];
4637  const unsigned relevel = level == pos1 ? pos2 : (level == pos2 ? pos1 : level);
4638  const unsigned long restride = result.strides_[relevel];
4639  const bool ready = (level == dim_ - 1);
4640  for (unsigned i = 0; i < imax; ++i) {
4641  if (ready)
4642  result.data_[idxRes] = data_[idxThis];
4643  else
4644  transposeLoop(level + 1, pos1, pos2, idxThis, idxRes, result);
4645  idxThis += mystride;
4646  idxRes += restride;
4647  }
4648  }
4649 
4650  template <typename Numeric, unsigned Len, unsigned Dim>
4651  template <typename Num2>
4652  Num2 ArrayND<Numeric, Len, Dim>::sum() const {
4653  if (!shapeIsKnown_)
4654  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"sum\"");
4655  Num2 sum = Num2();
4656  for (unsigned long i = 0; i < len_; ++i)
4657  sum += data_[i];
4658  return sum;
4659  }
4660 
4661  template <typename Numeric, unsigned Len, unsigned Dim>
4662  template <typename Num2>
4663  Num2 ArrayND<Numeric, Len, Dim>::sumsq() const {
4664  if (!shapeIsKnown_)
4665  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"sumsq\"");
4666  Num2 sum = Num2();
4667  for (unsigned long i = 0; i < len_; ++i) {
4668  const Num2 absval = absValue(data_[i]);
4669  sum += absval * absval;
4670  }
4671  return sum;
4672  }
4673 
4674  template <typename Numeric, unsigned Len, unsigned Dim>
4675  Numeric ArrayND<Numeric, Len, Dim>::min() const {
4676  if (!shapeIsKnown_)
4677  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"min\"");
4678  if (dim_) {
4679  Numeric minval(data_[0]);
4680  for (unsigned long i = 1UL; i < len_; ++i)
4681  if (ComplexComparesAbs<Numeric>::less(data_[i], minval))
4682  minval = data_[i];
4683  return minval;
4684  } else
4685  return localData_[0];
4686  }
4687 
4688  template <typename Numeric, unsigned Len, unsigned Dim>
4689  Numeric ArrayND<Numeric, Len, Dim>::min(unsigned* index, const unsigned indexLen) const {
4690  if (!shapeIsKnown_)
4691  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"min\"");
4692  if (indexLen != dim_)
4693  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::min: incompatible index length");
4694  if (dim_) {
4695  unsigned long minind = 0UL;
4696  Numeric minval(data_[0]);
4697  for (unsigned long i = 1UL; i < len_; ++i)
4698  if (ComplexComparesAbs<Numeric>::less(data_[i], minval)) {
4699  minval = data_[i];
4700  minind = i;
4701  }
4702  convertLinearIndex(minind, index, indexLen);
4703  return minval;
4704  } else
4705  return localData_[0];
4706  }
4707 
4708  template <typename Numeric, unsigned Len, unsigned Dim>
4709  Numeric ArrayND<Numeric, Len, Dim>::max() const {
4710  if (!shapeIsKnown_)
4711  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"max\"");
4712  if (dim_) {
4713  Numeric maxval(data_[0]);
4714  for (unsigned long i = 1UL; i < len_; ++i)
4715  if (ComplexComparesAbs<Numeric>::less(maxval, data_[i]))
4716  maxval = data_[i];
4717  return maxval;
4718  } else
4719  return localData_[0];
4720  }
4721 
4722  template <typename Numeric, unsigned Len, unsigned Dim>
4723  Numeric ArrayND<Numeric, Len, Dim>::max(unsigned* index, const unsigned indexLen) const {
4724  if (!shapeIsKnown_)
4725  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"max\"");
4726  if (indexLen != dim_)
4727  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::max: incompatible index length");
4728  if (dim_) {
4729  unsigned long maxind = 0UL;
4730  Numeric maxval(data_[0]);
4731  for (unsigned long i = 1UL; i < len_; ++i)
4732  if (ComplexComparesAbs<Numeric>::less(maxval, data_[i])) {
4733  maxval = data_[i];
4734  maxind = i;
4735  }
4736  convertLinearIndex(maxind, index, indexLen);
4737  return maxval;
4738  } else
4739  return localData_[0];
4740  }
4741 
4742  // Faster function for 2d transpose
4743  template <typename Numeric, unsigned Len, unsigned Dim>
4744  ArrayND<Numeric, Len, Dim> ArrayND<Numeric, Len, Dim>::transpose() const {
4745  if (dim_ != 2)
4746  throw npstat::NpstatInvalidArgument(
4747  "npstat::ArrayND::transpose method "
4748  "can not be used with arrays of rank other than 2");
4749  unsigned newshape[2];
4750  newshape[0] = shape_[1];
4751  newshape[1] = shape_[0];
4752  ArrayND<Numeric, Len, Dim> result(newshape, dim_);
4753  const unsigned imax = shape_[0];
4754  const unsigned jmax = shape_[1];
4755  for (unsigned i = 0; i < imax; ++i)
4756  for (unsigned j = 0; j < jmax; ++j)
4757  result.data_[j * imax + i] = data_[i * jmax + j];
4758  return result;
4759  }
4760 
4761  template <typename Numeric, unsigned Len, unsigned Dim>
4762  template <typename Accumulator>
4763  ArrayND<Numeric, Len, Dim> ArrayND<Numeric, Len, Dim>::derivative(const double inscale) const {
4764  if (!shapeIsKnown_)
4765  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"derivative\"");
4766  if (!dim_)
4767  throw npstat::NpstatInvalidArgument(
4768  "npstat::ArrayND::derivative method "
4769  "can not be used with array of 0 rank");
4770 
4771  const typename ProperDblFromCmpl<Accumulator>::type scale = inscale;
4772  const unsigned maxdim = CHAR_BIT * sizeof(unsigned long);
4773  if (dim_ >= maxdim)
4774  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::derivative: array rank is too large");
4775  const unsigned long maxcycle = 1UL << dim_;
4776 
4777  ArrayShape sh;
4778  sh.reserve(dim_);
4779  for (unsigned i = 0; i < dim_; ++i) {
4780  if (shape_[i] <= 1U)
4781  throw npstat::NpstatInvalidArgument(
4782  "In npstat::ArrayND::derivative: in some dimendions "
4783  "array size is too small");
4784  sh.push_back(shape_[i] - 1U);
4785  }
4786 
4787  ArrayND result(sh);
4788  const unsigned long rLen = result.length();
4789  for (unsigned long ilin = 0; ilin < rLen; ++ilin) {
4790  result.convertLinearIndex(ilin, &sh[0], dim_);
4791 
4792  Accumulator deriv = Accumulator();
4793  for (unsigned long icycle = 0UL; icycle < maxcycle; ++icycle) {
4794  unsigned long icell = 0UL;
4795  unsigned n1 = 0U;
4796  for (unsigned i = 0; i < dim_; ++i) {
4797  if (icycle & (1UL << i)) {
4798  ++n1;
4799  icell += strides_[i] * (sh[i] + 1);
4800  } else
4801  icell += strides_[i] * sh[i];
4802  }
4803  if ((dim_ - n1) % 2U)
4804  deriv -= data_[icell];
4805  else
4806  deriv += data_[icell];
4807  }
4808  result.data_[ilin] = static_cast<Numeric>(deriv * scale);
4809  }
4810 
4811  return result;
4812  }
4813 
4814  template <typename Numeric, unsigned Len, unsigned Dim>
4815  template <typename Accumulator>
4816  Accumulator ArrayND<Numeric, Len, Dim>::sumBelowLoop(const unsigned level,
4817  unsigned long idx0,
4818  const unsigned* limit) const {
4819  Accumulator cdf = Accumulator();
4820  const unsigned imax = limit[level] + 1U;
4821  if (level == dim_ - 1) {
4822  Numeric* base = data_ + idx0;
4823  for (unsigned i = 0; i < imax; ++i)
4824  cdf += base[i];
4825  } else {
4826  const unsigned long stride = strides_[level];
4827  for (unsigned i = 0; i < imax; ++i, idx0 += stride)
4828  cdf += sumBelowLoop<Accumulator>(level + 1, idx0, limit);
4829  }
4830  return cdf;
4831  }
4832 
4833  template <typename Numeric, unsigned Len, unsigned Dim>
4834  template <typename Accumulator>
4835  Accumulator ArrayND<Numeric, Len, Dim>::cdfValue(const unsigned* index, const unsigned indexLen) const {
4836  if (!shapeIsKnown_)
4837  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"cdfValue\"");
4838  if (!dim_)
4839  throw npstat::NpstatInvalidArgument(
4840  "npstat::ArrayND::cdfValue method "
4841  "can not be used with array of 0 rank");
4842  if (indexLen != dim_)
4843  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::cdfValue: incompatible index length");
4844  for (unsigned i = 0; i < indexLen; ++i)
4845  if (index[i] >= shape_[i])
4846  throw npstat::NpstatOutOfRange("In npstat::ArrayND::cdfValue: index out of range");
4847  return sumBelowLoop<Accumulator>(0, 0U, index);
4848  }
4849 
4850  template <typename Numeric, unsigned Len, unsigned Dim>
4851  template <typename Accumulator>
4852  ArrayND<Numeric, Len, Dim> ArrayND<Numeric, Len, Dim>::cdfArray(const double inscale) const {
4853  if (!shapeIsKnown_)
4854  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"cdfArray\"");
4855  if (!dim_)
4856  throw npstat::NpstatInvalidArgument(
4857  "npstat::ArrayND::cdfArray method "
4858  "can not be used with array of 0 rank");
4859 
4860  const proper_double scale = inscale;
4861  const unsigned maxdim = CHAR_BIT * sizeof(unsigned long);
4862  if (dim_ >= maxdim)
4863  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::cdfArray: array rank is too large");
4864  const unsigned long maxcycle = 1UL << dim_;
4865 
4866  ArrayShape sh;
4867  sh.reserve(dim_);
4868  for (unsigned i = 0; i < dim_; ++i)
4869  sh.push_back(shape_[i] + 1U);
4870 
4871  ArrayND<Accumulator> result(sh);
4872 
4873  unsigned* psh = &sh[0];
4874  const unsigned long len = result.length();
4875  for (unsigned long ipre = 0; ipre < len; ++ipre) {
4876  result.convertLinearIndex(ipre, psh, dim_);
4877  Accumulator deriv = Accumulator();
4878  bool has0 = false;
4879  for (unsigned i = 0; i < dim_; ++i)
4880  if (psh[i]-- == 0U) {
4881  has0 = true;
4882  break;
4883  }
4884  if (!has0) {
4885  for (unsigned long icycle = 0UL; icycle < maxcycle; ++icycle) {
4886  unsigned long icell = 0UL;
4887  unsigned n1 = 0U;
4888  for (unsigned i = 0; i < dim_; ++i) {
4889  if (icycle & (1UL << i)) {
4890  ++n1;
4891  icell += result.strides_[i] * (psh[i] + 1);
4892  } else
4893  icell += result.strides_[i] * psh[i];
4894  }
4895  if (n1 < dim_) {
4896  if ((dim_ - n1) % 2U)
4897  deriv += result.data_[icell];
4898  else
4899  deriv -= result.data_[icell];
4900  }
4901  }
4902  deriv += static_cast<Accumulator>(value(psh, dim_) * scale);
4903  }
4904  result.data_[ipre] = deriv;
4905  }
4906 
4907  // The "return" will convert Accumulator type into Numeric
4908  return result;
4909  }
4910 
4911  template <typename Numeric, unsigned Len, unsigned Dim>
4912  ArrayND<Numeric, Len, Dim> ArrayND<Numeric, Len, Dim>::transpose(const unsigned pos1, const unsigned pos2) const {
4913  if (!shapeIsKnown_)
4914  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"transpose\"");
4915  if (!(pos1 < dim_ && pos2 < dim_ && pos1 != pos2))
4916  throw npstat::NpstatInvalidArgument(
4917  "In npstat::ArrayND::transpose: "
4918  "incompatible transposition indices");
4919  if (dim_ == 2)
4920  return transpose();
4921  else {
4922  // Construct the new shape
4923  unsigned newshapeBuf[Dim];
4924  unsigned* newshape = makeBuffer(dim_, newshapeBuf, Dim);
4925  copyBuffer(newshape, shape_, dim_);
4926  std::swap(newshape[pos1], newshape[pos2]);
4927 
4928  // Form the result array
4929  ArrayND<Numeric, Len, Dim> result(newshape, dim_);
4930 
4931  // Fill the result array
4932  transposeLoop(0, pos1, pos2, 0UL, 0UL, result);
4933 
4934  destroyBuffer(newshape, newshapeBuf);
4935  return result;
4936  }
4937  }
4938 
4939  template <typename Numeric, unsigned Len, unsigned Dim>
4940  template <typename Num2, unsigned Len2, unsigned Dim2>
4941  void ArrayND<Numeric, Len, Dim>::multiMirror(ArrayND<Num2, Len2, Dim2>* out) const {
4942  assert(out);
4943  if (!shapeIsKnown_)
4944  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"multiMirror\"");
4945  if (!out->shapeIsKnown_)
4946  *out = ArrayND<Num2, Len2, Dim2>(doubleShape(shape()));
4947  if (dim_ != out->dim_)
4948  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::multiMirror: incompatible argument array rank");
4949 
4950  if (dim_) {
4951  const unsigned* dshape = out->shape_;
4952  for (unsigned i = 0; i < dim_; ++i)
4953  if (dshape[i] != shape_[i] * 2U)
4954  throw npstat::NpstatInvalidArgument(
4955  "In npstat::ArrayND::multiMirror: "
4956  "incompatible argument array shape");
4957 
4958  if (dim_ >= CHAR_BIT * sizeof(unsigned long))
4959  throw npstat::NpstatInvalidArgument(
4960  "In npstat::ArrayND::multiMirror: "
4961  "array rank is too large");
4962  const unsigned long maxcycle = 1UL << dim_;
4963  std::vector<unsigned> indexbuf(dim_ * 2U);
4964  unsigned* idx = &indexbuf[0];
4965  unsigned* mirror = idx + dim_;
4966 
4967  for (unsigned long ipt = 0; ipt < len_; ++ipt) {
4968  unsigned long l = ipt;
4969  for (unsigned i = 0; i < dim_; ++i) {
4970  idx[i] = l / strides_[i];
4971  l -= (idx[i] * strides_[i]);
4972  }
4973  for (unsigned long icycle = 0UL; icycle < maxcycle; ++icycle) {
4974  for (unsigned i = 0; i < dim_; ++i) {
4975  if (icycle & (1UL << i))
4976  mirror[i] = dshape[i] - idx[i] - 1U;
4977  else
4978  mirror[i] = idx[i];
4979  }
4980  out->value(mirror, dim_) = data_[ipt];
4981  }
4982  }
4983  } else
4984  out->localData_[0] = static_cast<Num2>(localData_[0]);
4985  }
4986 
4987  template <typename Numeric, unsigned Len, unsigned Dim>
4988  template <typename Num2, unsigned Len2, unsigned Dim2>
4989  void ArrayND<Numeric, Len, Dim>::rotate(const unsigned* shifts,
4990  const unsigned lenShifts,
4991  ArrayND<Num2, Len2, Dim2>* rotated) const {
4992  assert(rotated);
4993  if (!shapeIsKnown_)
4994  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"rotate\"");
4995  // Can't rotate into itself -- it will be a mess
4996  if ((void*)rotated == (void*)this)
4997  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::rotate: can not rotate array into itself");
4998  if (!rotated->shapeIsKnown_)
4999  *rotated = *this;
5000  if (dim_ != rotated->dim_)
5001  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::rotate: incompatible argument array rank");
5002  if (lenShifts != dim_)
5003  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::rotate: incompatible dimensionality of shifts");
5004 
5005  if (dim_) {
5006  assert(shifts);
5007  if (dim_ > CHAR_BIT * sizeof(unsigned long))
5008  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::rotate: array rank is too large");
5009  unsigned buf[CHAR_BIT * sizeof(unsigned long)];
5010  clearBuffer(buf, dim_);
5011  (const_cast<ArrayND*>(this))
5012  ->flatCircularLoop(0U, 0UL, 0UL, buf, shape_, shifts, *rotated, scast_assign_right<Numeric, Num2>());
5013  } else
5014  rotated->localData_[0] = static_cast<Num2>(localData_[0]);
5015  }
5016 
5017  template <typename Numeric, unsigned Len, unsigned Dim>
5018  template <typename Num2, unsigned Len2, unsigned Dim2>
5019  void ArrayND<Numeric, Len, Dim>::dotProductLoop(const unsigned level,
5020  unsigned long idx0,
5021  unsigned long idx1,
5022  unsigned long idx2,
5023  const ArrayND<Num2, Len2, Dim2>& r,
5024  ArrayND& result) const {
5025  // idx0 -- this object
5026  // idx1 -- dot product argument
5027  // idx2 -- result
5028  if (level == result.dim_) {
5029  Numeric sum = Numeric();
5030  const unsigned imax = r.shape_[0];
5031  const unsigned rstride = r.strides_[0];
5032  const Numeric* l = data_ + idx0;
5033  const Num2* ri = r.data_ + idx1;
5034  for (unsigned i = 0; i < imax; ++i)
5035  sum += l[i] * ri[i * rstride];
5036  result.data_[idx2] = sum;
5037  } else {
5038  const unsigned imax = result.shape_[level];
5039  for (unsigned i = 0; i < imax; ++i) {
5040  dotProductLoop(level + 1, idx0, idx1, idx2, r, result);
5041  idx2 += result.strides_[level];
5042  if (level < dim_ - 1)
5043  idx0 += strides_[level];
5044  else
5045  idx1 += r.strides_[level + 2 - dim_];
5046  }
5047  }
5048  }
5049 
5050  template <typename Numeric, unsigned Len, unsigned Dim>
5051  template <typename Num2, unsigned Len2, unsigned Dim2>
5052  ArrayND<Numeric, Len, Dim> ArrayND<Numeric, Len, Dim>::dot(const ArrayND<Num2, Len2, Dim2>& r) const {
5053  if (!dim_)
5054  throw npstat::NpstatInvalidArgument(
5055  "npstat::ArrayND::dot method "
5056  "can not be used with array of 0 rank");
5057  if (!r.dim_)
5058  throw npstat::NpstatInvalidArgument(
5059  "npstat::ArrayND::dot method "
5060  "can not be used with argument array of 0 rank");
5061  if (shape_[dim_ - 1] != r.shape_[0])
5062  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::dot: incompatible argument array shape");
5063 
5064  if (dim_ == 1 && r.dim_ == 1) {
5065  // Special case: the result is of 0 rank
5066  ArrayND<Numeric, Len, Dim> result(static_cast<unsigned*>(nullptr), 0U);
5067  Numeric sum = Numeric();
5068  const unsigned imax = shape_[0];
5069  for (unsigned i = 0; i < imax; ++i)
5070  sum += data_[i] * r.data_[i];
5071  result() = sum;
5072  return result;
5073  } else {
5074  unsigned newshapeBuf[2 * Dim];
5075  unsigned* newshape = makeBuffer(dim_ + r.dim_ - 2, newshapeBuf, 2 * Dim);
5076  copyBuffer(newshape, shape_, dim_ - 1);
5077  copyBuffer(newshape + (dim_ - 1), r.shape_ + 1, r.dim_ - 1);
5078  ArrayND<Numeric, Len, Dim> result(newshape, dim_ + r.dim_ - 2);
5079 
5080  dotProductLoop(0U, 0UL, 0UL, 0UL, r, result);
5081 
5082  destroyBuffer(newshape, newshapeBuf);
5083  return result;
5084  }
5085  }
5086 
5087  template <typename Numeric, unsigned Len, unsigned Dim>
5088  inline unsigned ArrayND<Numeric, Len, Dim>::span(const unsigned dim) const {
5089  if (dim >= dim_)
5090  throw npstat::NpstatOutOfRange("In npstat::ArrayND::span: dimension number is out of range");
5091  return shape_[dim];
5092  }
5093 
5094  template <typename Numeric, unsigned Len, unsigned Dim>
5095  unsigned ArrayND<Numeric, Len, Dim>::maximumSpan() const {
5096  unsigned maxspan = 0;
5097  for (unsigned i = 0; i < dim_; ++i)
5098  if (shape_[i] > maxspan)
5099  maxspan = shape_[i];
5100  return maxspan;
5101  }
5102 
5103  template <typename Numeric, unsigned Len, unsigned Dim>
5104  unsigned ArrayND<Numeric, Len, Dim>::minimumSpan() const {
5105  if (dim_ == 0)
5106  return 0U;
5107  unsigned minspan = shape_[0];
5108  for (unsigned i = 1; i < dim_; ++i)
5109  if (shape_[i] < minspan)
5110  minspan = shape_[i];
5111  return minspan;
5112  }
5113 
5114  template <typename Numeric, unsigned Len, unsigned Dim>
5115  inline const Numeric& ArrayND<Numeric, Len, Dim>::cl() const {
5116  if (!shapeIsKnown_)
5117  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"cl\"");
5118  if (dim_)
5119  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::cl: wrong # of args (not rank 0 array)");
5120  return localData_[0];
5121  }
5122 
5123  template <typename Numeric, unsigned Len, unsigned Dim>
5124  inline const Numeric& ArrayND<Numeric, Len, Dim>::cl(const double i0) const {
5125  if (1U != dim_)
5126  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::cl: wrong # of args (not rank 1 array)");
5127  return data_[coordToIndex(i0, 0)];
5128  }
5129 
5130  template <typename Numeric, unsigned Len, unsigned Dim>
5131  inline const Numeric& ArrayND<Numeric, Len, Dim>::cl(const double i0, const double i1) const {
5132  if (2U != dim_)
5133  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::cl: wrong # of args (not rank 2 array)");
5134  return data_[coordToIndex(i0, 0) * strides_[0] + coordToIndex(i1, 1)];
5135  }
5136 
5137  template <typename Numeric, unsigned Len, unsigned Dim>
5138  inline const Numeric& ArrayND<Numeric, Len, Dim>::cl(const double i0, const double i1, const double i2) const {
5139  if (3U != dim_)
5140  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::cl: wrong # of args (not rank 3 array)");
5141  return data_[coordToIndex(i0, 0) * strides_[0] + coordToIndex(i1, 1) * strides_[1] + coordToIndex(i2, 2)];
5142  }
5143 
5144  template <typename Numeric, unsigned Len, unsigned Dim>
5145  inline const Numeric& ArrayND<Numeric, Len, Dim>::cl(const double i0,
5146  const double i1,
5147  const double i2,
5148  const double i3) const {
5149  if (4U != dim_)
5150  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::cl: wrong # of args (not rank 4 array)");
5151  return data_[coordToIndex(i0, 0) * strides_[0] + coordToIndex(i1, 1) * strides_[1] +
5152  coordToIndex(i2, 2) * strides_[2] + coordToIndex(i3, 3)];
5153  }
5154 
5155  template <typename Numeric, unsigned Len, unsigned Dim>
5156  inline const Numeric& ArrayND<Numeric, Len, Dim>::cl(
5157  const double i0, const double i1, const double i2, const double i3, const double i4) const {
5158  if (5U != dim_)
5159  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::cl: wrong # of args (not rank 5 array)");
5160  return data_[coordToIndex(i0, 0) * strides_[0] + coordToIndex(i1, 1) * strides_[1] +
5161  coordToIndex(i2, 2) * strides_[2] + coordToIndex(i3, 3) * strides_[3] + coordToIndex(i4, 4)];
5162  }
5163 
5164  template <typename Numeric, unsigned Len, unsigned Dim>
5165  inline const Numeric& ArrayND<Numeric, Len, Dim>::cl(
5166  const double i0, const double i1, const double i2, const double i3, const double i4, const double i5) const {
5167  if (6U != dim_)
5168  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::cl: wrong # of args (not rank 6 array)");
5169  return data_[coordToIndex(i0, 0) * strides_[0] + coordToIndex(i1, 1) * strides_[1] +
5170  coordToIndex(i2, 2) * strides_[2] + coordToIndex(i3, 3) * strides_[3] +
5171  coordToIndex(i4, 4) * strides_[4] + coordToIndex(i5, 5)];
5172  }
5173 
5174  template <typename Numeric, unsigned Len, unsigned Dim>
5175  inline const Numeric& ArrayND<Numeric, Len, Dim>::cl(const double i0,
5176  const double i1,
5177  const double i2,
5178  const double i3,
5179  const double i4,
5180  const double i5,
5181  const double i6) const {
5182  if (7U != dim_)
5183  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::cl: wrong # of args (not rank 7 array)");
5184  return data_[coordToIndex(i0, 0) * strides_[0] + coordToIndex(i1, 1) * strides_[1] +
5185  coordToIndex(i2, 2) * strides_[2] + coordToIndex(i3, 3) * strides_[3] +
5186  coordToIndex(i4, 4) * strides_[4] + coordToIndex(i5, 5) * strides_[5] + coordToIndex(i6, 6)];
5187  }
5188 
5189  template <typename Numeric, unsigned Len, unsigned Dim>
5190  inline const Numeric& ArrayND<Numeric, Len, Dim>::cl(const double i0,
5191  const double i1,
5192  const double i2,
5193  const double i3,
5194  const double i4,
5195  const double i5,
5196  const double i6,
5197  const double i7) const {
5198  if (8U != dim_)
5199  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::cl: wrong # of args (not rank 8 array)");
5200  return data_[coordToIndex(i0, 0) * strides_[0] + coordToIndex(i1, 1) * strides_[1] +
5201  coordToIndex(i2, 2) * strides_[2] + coordToIndex(i3, 3) * strides_[3] +
5202  coordToIndex(i4, 4) * strides_[4] + coordToIndex(i5, 5) * strides_[5] +
5203  coordToIndex(i6, 6) * strides_[6] + coordToIndex(i7, 7)];
5204  }
5205 
5206  template <typename Numeric, unsigned Len, unsigned Dim>
5207  inline const Numeric& ArrayND<Numeric, Len, Dim>::cl(const double i0,
5208  const double i1,
5209  const double i2,
5210  const double i3,
5211  const double i4,
5212  const double i5,
5213  const double i6,
5214  const double i7,
5215  const double i8) const {
5216  if (9U != dim_)
5217  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::cl: wrong # of args (not rank 9 array)");
5218  return data_[coordToIndex(i0, 0) * strides_[0] + coordToIndex(i1, 1) * strides_[1] +
5219  coordToIndex(i2, 2) * strides_[2] + coordToIndex(i3, 3) * strides_[3] +
5220  coordToIndex(i4, 4) * strides_[4] + coordToIndex(i5, 5) * strides_[5] +
5221  coordToIndex(i6, 6) * strides_[6] + coordToIndex(i7, 7) * strides_[7] + coordToIndex(i8, 8)];
5222  }
5223 
5224  template <typename Numeric, unsigned Len, unsigned Dim>
5225  inline const Numeric& ArrayND<Numeric, Len, Dim>::cl(const double i0,
5226  const double i1,
5227  const double i2,
5228  const double i3,
5229  const double i4,
5230  const double i5,
5231  const double i6,
5232  const double i7,
5233  const double i8,
5234  const double i9) const {
5235  if (10U != dim_)
5236  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::cl: wrong # of args (not rank 10 array)");
5237  return data_[coordToIndex(i0, 0) * strides_[0] + coordToIndex(i1, 1) * strides_[1] +
5238  coordToIndex(i2, 2) * strides_[2] + coordToIndex(i3, 3) * strides_[3] +
5239  coordToIndex(i4, 4) * strides_[4] + coordToIndex(i5, 5) * strides_[5] +
5240  coordToIndex(i6, 6) * strides_[6] + coordToIndex(i7, 7) * strides_[7] +
5241  coordToIndex(i8, 8) * strides_[8] + coordToIndex(i9, 9)];
5242  }
5243 
5244  template <typename Numeric, unsigned Len, unsigned Dim>
5245  inline Numeric& ArrayND<Numeric, Len, Dim>::cl() {
5246  if (!shapeIsKnown_)
5247  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"cl\"");
5248  if (dim_)
5249  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::cl: wrong # of args (not rank 0 array)");
5250  return localData_[0];
5251  }
5252 
5253  template <typename Numeric, unsigned Len, unsigned Dim>
5254  inline Numeric& ArrayND<Numeric, Len, Dim>::cl(const double i0) {
5255  if (1U != dim_)
5256  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::cl: wrong # of args (not rank 1 array)");
5257  return data_[coordToIndex(i0, 0)];
5258  }
5259 
5260  template <typename Numeric, unsigned Len, unsigned Dim>
5261  inline Numeric& ArrayND<Numeric, Len, Dim>::cl(const double i0, const double i1) {
5262  if (2U != dim_)
5263  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::cl: wrong # of args (not rank 2 array)");
5264  return data_[coordToIndex(i0, 0) * strides_[0] + coordToIndex(i1, 1)];
5265  }
5266 
5267  template <typename Numeric, unsigned Len, unsigned Dim>
5268  inline Numeric& ArrayND<Numeric, Len, Dim>::cl(const double i0, const double i1, const double i2) {
5269  if (3U != dim_)
5270  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::cl: wrong # of args (not rank 3 array)");
5271  return data_[coordToIndex(i0, 0) * strides_[0] + coordToIndex(i1, 1) * strides_[1] + coordToIndex(i2, 2)];
5272  }
5273 
5274  template <typename Numeric, unsigned Len, unsigned Dim>
5275  inline Numeric& ArrayND<Numeric, Len, Dim>::cl(const double i0, const double i1, const double i2, const double i3) {
5276  if (4U != dim_)
5277  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::cl: wrong # of args (not rank 4 array)");
5278  return data_[coordToIndex(i0, 0) * strides_[0] + coordToIndex(i1, 1) * strides_[1] +
5279  coordToIndex(i2, 2) * strides_[2] + coordToIndex(i3, 3)];
5280  }
5281 
5282  template <typename Numeric, unsigned Len, unsigned Dim>
5283  inline Numeric& ArrayND<Numeric, Len, Dim>::cl(
5284  const double i0, const double i1, const double i2, const double i3, const double i4) {
5285  if (5U != dim_)
5286  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::cl: wrong # of args (not rank 5 array)");
5287  return data_[coordToIndex(i0, 0) * strides_[0] + coordToIndex(i1, 1) * strides_[1] +
5288  coordToIndex(i2, 2) * strides_[2] + coordToIndex(i3, 3) * strides_[3] + coordToIndex(i4, 4)];
5289  }
5290 
5291  template <typename Numeric, unsigned Len, unsigned Dim>
5292  inline Numeric& ArrayND<Numeric, Len, Dim>::cl(
5293  const double i0, const double i1, const double i2, const double i3, const double i4, const double i5) {
5294  if (6U != dim_)
5295  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::cl: wrong # of args (not rank 6 array)");
5296  return data_[coordToIndex(i0, 0) * strides_[0] + coordToIndex(i1, 1) * strides_[1] +
5297  coordToIndex(i2, 2) * strides_[2] + coordToIndex(i3, 3) * strides_[3] +
5298  coordToIndex(i4, 4) * strides_[4] + coordToIndex(i5, 5)];
5299  }
5300 
5301  template <typename Numeric, unsigned Len, unsigned Dim>
5302  inline Numeric& ArrayND<Numeric, Len, Dim>::cl(const double i0,
5303  const double i1,
5304  const double i2,
5305  const double i3,
5306  const double i4,
5307  const double i5,
5308  const double i6) {
5309  if (7U != dim_)
5310  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::cl: wrong # of args (not rank 7 array)");
5311  return data_[coordToIndex(i0, 0) * strides_[0] + coordToIndex(i1, 1) * strides_[1] +
5312  coordToIndex(i2, 2) * strides_[2] + coordToIndex(i3, 3) * strides_[3] +
5313  coordToIndex(i4, 4) * strides_[4] + coordToIndex(i5, 5) * strides_[5] + coordToIndex(i6, 6)];
5314  }
5315 
5316  template <typename Numeric, unsigned Len, unsigned Dim>
5317  inline Numeric& ArrayND<Numeric, Len, Dim>::cl(const double i0,
5318  const double i1,
5319  const double i2,
5320  const double i3,
5321  const double i4,
5322  const double i5,
5323  const double i6,
5324  const double i7) {
5325  if (8U != dim_)
5326  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::cl: wrong # of args (not rank 8 array)");
5327  return data_[coordToIndex(i0, 0) * strides_[0] + coordToIndex(i1, 1) * strides_[1] +
5328  coordToIndex(i2, 2) * strides_[2] + coordToIndex(i3, 3) * strides_[3] +
5329  coordToIndex(i4, 4) * strides_[4] + coordToIndex(i5, 5) * strides_[5] +
5330  coordToIndex(i6, 6) * strides_[6] + coordToIndex(i7, 7)];
5331  }
5332 
5333  template <typename Numeric, unsigned Len, unsigned Dim>
5334  inline Numeric& ArrayND<Numeric, Len, Dim>::cl(const double i0,
5335  const double i1,
5336  const double i2,
5337  const double i3,
5338  const double i4,
5339  const double i5,
5340  const double i6,
5341  const double i7,
5342  const double i8) {
5343  if (9U != dim_)
5344  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::cl: wrong # of args (not rank 9 array)");
5345  return data_[coordToIndex(i0, 0) * strides_[0] + coordToIndex(i1, 1) * strides_[1] +
5346  coordToIndex(i2, 2) * strides_[2] + coordToIndex(i3, 3) * strides_[3] +
5347  coordToIndex(i4, 4) * strides_[4] + coordToIndex(i5, 5) * strides_[5] +
5348  coordToIndex(i6, 6) * strides_[6] + coordToIndex(i7, 7) * strides_[7] + coordToIndex(i8, 8)];
5349  }
5350 
5351  template <typename Numeric, unsigned Len, unsigned Dim>
5352  inline Numeric& ArrayND<Numeric, Len, Dim>::cl(const double i0,
5353  const double i1,
5354  const double i2,
5355  const double i3,
5356  const double i4,
5357  const double i5,
5358  const double i6,
5359  const double i7,
5360  const double i8,
5361  const double i9) {
5362  if (10U != dim_)
5363  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::cl: wrong # of args (not rank 10 array)");
5364  return data_[coordToIndex(i0, 0) * strides_[0] + coordToIndex(i1, 1) * strides_[1] +
5365  coordToIndex(i2, 2) * strides_[2] + coordToIndex(i3, 3) * strides_[3] +
5366  coordToIndex(i4, 4) * strides_[4] + coordToIndex(i5, 5) * strides_[5] +
5367  coordToIndex(i6, 6) * strides_[6] + coordToIndex(i7, 7) * strides_[7] +
5368  coordToIndex(i8, 8) * strides_[8] + coordToIndex(i9, 9)];
5369  }
5370 
5371  template <typename Numeric, unsigned StackLen, unsigned StackDim>
5372  const char* ArrayND<Numeric, StackLen, StackDim>::classname() {
5373  static const std::string name(gs::template_class_name<Numeric>("npstat::ArrayND"));
5374  return name.c_str();
5375  }
5376 
5377  template <typename Numeric, unsigned StackLen, unsigned StackDim>
5378  bool ArrayND<Numeric, StackLen, StackDim>::write(std::ostream& os) const {
5379  if (!shapeIsKnown_)
5380  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"write\"");
5381  gs::write_pod_vector(os, shape());
5382  return !os.fail() && (dim_ ? gs::write_array(os, data_, len_) : gs::write_item(os, localData_[0], false));
5383  }
5384 
5385  template <typename Numeric, unsigned Len, unsigned Dim>
5386  void ArrayND<Numeric, Len, Dim>::restore(const gs::ClassId& id, std::istream& in, ArrayND* array) {
5387  static const gs::ClassId current(gs::ClassId::makeId<ArrayND<Numeric, Len, Dim> >());
5388  current.ensureSameId(id);
5389 
5390  ArrayShape rshape;
5391  gs::read_pod_vector(in, &rshape);
5392  if (in.fail())
5393  throw gs::IOReadFailure("In npstat::ArrayND::restore: input stream failure (checkpoint 0)");
5394 
5395  assert(array);
5396  array->uninitialize();
5397  array->dim_ = rshape.size();
5398  array->shapeIsKnown_ = true;
5399  array->len_ = 1UL;
5400  if (array->dim_) {
5401  array->shape_ = makeBuffer(array->dim_, array->localShape_, Dim);
5402  for (unsigned i = 0; i < array->dim_; ++i) {
5403  array->shape_[i] = rshape[i];
5404  assert(array->shape_[i]);
5405  array->len_ *= array->shape_[i];
5406  }
5407  array->buildStrides();
5408  array->data_ = makeBuffer(array->len_, array->localData_, Len);
5409  gs::read_array(in, array->data_, array->len_);
5410  } else
5411  gs::restore_item(in, array->localData_, false);
5412  if (in.fail())
5413  throw gs::IOReadFailure("In npstat::ArrayND::restore: input stream failure (checkpoint 1)");
5414  }
5415 
5416  template <typename Numeric, unsigned Len, unsigned StackDim>
5417  template <typename Num2, unsigned Len2, unsigned Dim2>
5418  void ArrayND<Numeric, Len, StackDim>::exportSubrange(const unsigned* corner,
5419  const unsigned lenCorner,
5420  ArrayND<Num2, Len2, Dim2>* out) const {
5421  if (!shapeIsKnown_)
5422  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"exportSubrange\"");
5423  if (dim_ != lenCorner)
5424  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::exportSubrange: incompatible corner index length");
5425  assert(out);
5426  if (!out->shapeIsKnown_)
5427  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::exportSubrange: uninitialized argument array");
5428  if (out->dim_ != dim_)
5429  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::exportSubrange: incompatible argument array rank");
5430 
5431  if (dim_) {
5432  assert(corner);
5433  if (dim_ > CHAR_BIT * sizeof(unsigned long))
5434  throw npstat::NpstatInvalidArgument(
5435  "In npstat::ArrayND::exportSubrange: "
5436  "array rank is too large");
5437  unsigned toBuf[CHAR_BIT * sizeof(unsigned long)];
5438  clearBuffer(toBuf, dim_);
5439  (const_cast<ArrayND*>(this))
5440  ->commonSubrangeLoop(0U, 0UL, 0UL, corner, out->shape_, toBuf, *out, scast_assign_right<Numeric, Num2>());
5441  } else
5442  out->localData_[0] = static_cast<Num2>(localData_[0]);
5443  }
5444 
5445  template <typename Numeric, unsigned Len, unsigned StackDim>
5446  template <typename Num2, unsigned Len2, unsigned Dim2>
5447  void ArrayND<Numeric, Len, StackDim>::importSubrange(const unsigned* corner,
5448  const unsigned lenCorner,
5449  const ArrayND<Num2, Len2, Dim2>& from) {
5450  if (!shapeIsKnown_)
5451  throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"importSubrange\"");
5452  if (dim_ != lenCorner)
5453  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::importSubrange: incompatible corner index length");
5454  if (!from.shapeIsKnown_)
5455  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::importSubrange: uninitialized argument array");
5456  if (from.dim_ != dim_)
5457  throw npstat::NpstatInvalidArgument("In npstat::ArrayND::importSubrange: incompatible argument array rank");
5458 
5459  if (dim_) {
5460  assert(corner);
5461  if (dim_ > CHAR_BIT * sizeof(unsigned long))
5462  throw npstat::NpstatInvalidArgument(
5463  "In npstat::ArrayND::importSubrange: "
5464  "array rank is too large");
5465  unsigned toBuf[CHAR_BIT * sizeof(unsigned long)];
5466  clearBuffer(toBuf, dim_);
5467  commonSubrangeLoop(0U,
5468  0UL,
5469  0UL,
5470  corner,
5471  from.shape_,
5472  toBuf,
5473  const_cast<ArrayND<Num2, Len2, Dim2>&>(from),
5474  scast_assign_left<Numeric, Num2>());
5475  } else
5476  localData_[0] = static_cast<Numeric>(from.localData_[0]);
5477  }
5478 } // namespace npstat
5479 
5480 #endif // NPSTAT_ARRAYND_HH_
dot
T dot(const Basic3DVector &v) const
Scalar product, or "dot" product, with a vector of same type.
Definition: Basic3DVectorLD.h:212
personalPlayback.level
level
Definition: personalPlayback.py:22
npstat::ArrayND::localData_
Numeric localData_[StackLen]
Definition: ArrayND.h:1048
npstat::ArrayND::jointSliceScan
void jointSliceScan(ArrayND< Num2, Len2, Dim2 > &slice, const unsigned *fixedIndices, const unsigned *fixedIndexValues, unsigned nFixedIndices, Functor binaryFunct)
Definition: ArrayND.h:1933
npstat::scast_pluseq_left
Definition: SimpleFunctors.h:357
bTagCombinedSVVariables_cff.indices
indices
Definition: bTagCombinedSVVariables_cff.py:67
npstat::ArrayND::dualCircularLoop
void dualCircularLoop(unsigned level, unsigned long idx0, unsigned long idx1, const unsigned *thisCorner, const unsigned *range, const unsigned *otherCorner, ArrayND< Num2, Len2, Dim2 > &other, Functor binaryFunct)
Definition: ArrayND.h:1392
npoints
static const int npoints
Definition: NuclearInteractionFTFSimulator.h:38
npstat::ArrayND::operator!=
bool operator!=(const ArrayND< Numeric, Len2, Dim2 > &) const
Definition: ArrayND.h:3977
geometryDiff.transpose
def transpose(a)
Definition: geometryDiff.py:39
FastTimerService_cff.range
range
Definition: FastTimerService_cff.py:34
npstat::ArrayND::cdfArray
ArrayND cdfArray(double scale=1.0) const
npstat::ArrayND::interpolate3
Numeric interpolate3(const double *x, unsigned xDim) const
Definition: ArrayND.h:4309
photonAnalyzer_cfi.rMax
rMax
Definition: photonAnalyzer_cfi.py:91
npstat::copyBuffer
void copyBuffer(T1 *dest, const T2 *source, const unsigned long len)
Definition: allocators.h:38
testProducerWithPsetDescEmpty_cfi.i2
i2
Definition: testProducerWithPsetDescEmpty_cfi.py:46
npstat::ArrayND::circularFlatLoop
void circularFlatLoop(unsigned level, unsigned long idx0, unsigned long idx1, const unsigned *thisCorner, const unsigned *range, const unsigned *otherCorner, ArrayND< Num2, Len2, Dim2 > &other, Functor binaryFunct)
Definition: ArrayND.h:1467
alignBH_cfg.fixed
fixed
Definition: alignBH_cfg.py:54
npstat::ArrayND::operator/=
ArrayND & operator/=(const Num2 &r)
testProducerWithPsetDescEmpty_cfi.i3
i3
Definition: testProducerWithPsetDescEmpty_cfi.py:47
npstat::ArrayND::shapeData
const unsigned * shapeData() const
Definition: ArrayND.h:248
absDifference.h
Calculate absolute value of a difference between two numbers for an extended set of types.
mps_fire.i
i
Definition: mps_fire.py:428
SiPixelPI::one
Definition: SiPixelPayloadInspectorHelper.h:39
npstat::ArrayND::jointSliceLoop
void jointSliceLoop(unsigned level, unsigned long idx0, unsigned level1, unsigned long idx1, Num2 *sliceData, const unsigned long *sliceStrides, const unsigned *fixedIndices, const unsigned *fixedIndexValues, unsigned nFixedIndices, Functor binaryFunctor)
npstat::ArrayND::convertToLastDimCdfLoop
void convertToLastDimCdfLoop(ArrayND *sumSlice, unsigned level, unsigned long idx0, unsigned long idxSlice, bool useTrapezoids)
Definition: ArrayND.h:4457
npstat::ArrayND::applySlice
void applySlice(ArrayND< Num2, Len2, Dim2 > &slice, const unsigned *fixedIndices, unsigned nFixedIndices, Functor binaryFunct)
Definition: ArrayND.h:2354
npstat::AbsVisitor
Definition: AbsVisitor.h:19
funct::false
false
Definition: Factorize.h:29
npstat::ArrayND::operator=
ArrayND & operator=(const ArrayND &)
Definition: ArrayND.h:2914
operator*
MatrixMeschach operator*(const MatrixMeschach &mat1, const MatrixMeschach &mat2)
Definition: MatrixMeschach.cc:117
npstat::ArrayND::write
bool write(std::ostream &of) const
Definition: ArrayND.h:5378
L1EGammaCrystalsEmulatorProducer_cfi.scale
scale
Definition: L1EGammaCrystalsEmulatorProducer_cfi.py:10
npstat::ArrayND::restore
static void restore(const gs::ClassId &id, std::istream &in, ArrayND *array)
Definition: ArrayND.h:5386
npstat::ArrayND::functorFill
ArrayND & functorFill(Functor f)
npstat::ArrayND::isClose
bool isClose(const ArrayND< Num2, Len2, Dim2 > &r, double eps) const
Definition: ArrayND.h:4530
npstat::ArrayND::exportSubrange
void exportSubrange(const unsigned *fromCorner, unsigned lenCorner, ArrayND< Num2, Len2, Dim2 > *dest) const
Definition: ArrayND.h:5418
npstat::ArrayND::uninitialize
ArrayND & uninitialize()
Definition: ArrayND.h:4356
npstat::interpolate_cubic
T interpolate_cubic(const double x, const T &f0, const T &f1, const T &f2, const T &f3)
Definition: interpolate.h:44
f
double f[11][100]
Definition: MuScleFitUtils.cc:78
npstat::ArrayND::operator+=
ArrayND & operator+=(const ArrayND< Num2, Len2, Dim2 > &r)
npstat::ArrayND::classname
static const char * classname()
Definition: ArrayND.h:5372
npstat::Interval::min
const Numeric min() const
Definition: Interval.h:78
npstat::ArrayND::importSlice
void importSlice(const ArrayND< Num2, Len2, Dim2 > &slice, const unsigned *fixedIndices, const unsigned *fixedIndexValues, unsigned nFixedIndices)
Definition: ArrayND.h:722
npstat::ArrayND::outerProductLoop
void outerProductLoop(unsigned level, unsigned long idx0, unsigned long idx1, unsigned long idx2, const ArrayND< Num1, Len1, Dim1 > &a1, const ArrayND< Num2, Len2, Dim2 > &a2)
Definition: ArrayND.h:2845
min
T min(T a, T b)
Definition: MathUtil.h:58
npstat::ArrayRange
Definition: ArrayRange.h:21
npstat::ArrayND::jointScan
void jointScan(ArrayND< Num2, Len2, Dim2 > &other, Functor binaryFunct)
Definition: ArrayND.h:2342
npstat::ArrayND::transpose
ArrayND transpose() const
Definition: ArrayND.h:4744
testProducerWithPsetDescEmpty_cfi.x2
x2
Definition: testProducerWithPsetDescEmpty_cfi.py:28
npstat::ArrayND::importMemSlice
void importMemSlice(const Num2 *buffer, unsigned long bufLen, const unsigned *fixedIndices, const unsigned *fixedIndexValues, unsigned nFixedIndices)
Definition: ArrayND.h:738
assign
void assign(const std::vector< float > &vec, float &a, float &b, float &c, float &d)
Definition: PileupJetIdAlgo.cc:139
npstat::ArrayND::shape_
unsigned * shape_
Definition: ArrayND.h:1055
gpuVertexFinder::eps
WorkSpace int float eps
Definition: gpuClusterTracksDBSCAN.h:18
npstat::ArrayND::exportMemSlice
void exportMemSlice(Num2 *buffer, unsigned long bufLen, const unsigned *fixedIndices, const unsigned *fixedIndexValues, unsigned nFixedIndices) const
Definition: ArrayND.h:710
npstat::ArrayND::isDensity
bool isDensity() const
Definition: ArrayND.h:3009
testProducerWithPsetDescEmpty_cfi.i1
i1
Definition: testProducerWithPsetDescEmpty_cfi.py:45
n0
int n0
Definition: AMPTWrapper.h:44
npstat::ArrayND::operator==
bool operator==(const ArrayND< Numeric, Len2, Dim2 > &) const
Definition: ArrayND.h:3957
npstat::ArrayND::sum
Num2 sum() const
Definition: ArrayND.h:4652
ProperDblFromCmpl.h
Compile-time deduction of the underlying floating point type from the given complex type.
npstat::ArrayND::minimumSpan
unsigned minimumSpan() const
Definition: ArrayND.h:5104
npstat::ArrayND::projectLoop2
void projectLoop2(unsigned level, unsigned long idx0, unsigned level1, unsigned long idx1, ArrayND< Num2, Len2, Dim2 > *projection, AbsVisitor< Numeric, Num3 > &projector, const unsigned *projectedIndices, unsigned nProjectedIndices, Op fcn) const
Definition: ArrayND.h:2199
to
cms::cuda::assert
assert(be >=bs)
funct::derivative
Derivative< X, A >::type derivative(const A &_)
Definition: Derivative.h:18
npstat::ArrayND::multiplyBySlice
ArrayND & multiplyBySlice(const ArrayND< Num2, Len2, Dim2 > &slice, const unsigned *fixedIndices, unsigned nFixedIndices)
Definition: ArrayND.h:768
npstat::AbsArrayProjector::clear
virtual void clear()=0
npstat::ArrayND::contract
ArrayND contract(unsigned pos1, unsigned pos2) const
Definition: ArrayND.h:4590
npstat::ArrayND::operator*
ArrayND operator*(const Num2 &r) const
npstat::ArrayND::verifyProjectionCompatibility
void verifyProjectionCompatibility(const ArrayND< Num2, Len2, Dim2 > &projection, const unsigned *projectedIndices, unsigned nProjectedIndices) const
Definition: ArrayND.h:2055
npstat::ArrayND::circularFlatScan
void circularFlatScan(ArrayND< Num2, Len2, Dim2 > &other, const unsigned *thisCorner, const unsigned *range, const unsigned *otherCorner, unsigned arrLen, Functor binaryFunct)
Definition: ArrayND.h:1508
npstat::ArrayND::contractLoop
void contractLoop(unsigned thisLevel, unsigned resLevel, unsigned pos1, unsigned pos2, unsigned long idxThis, unsigned long idxRes, ArrayND &result) const
Definition: ArrayND.h:4556
mps_check.array
array
Definition: mps_check.py:216
npstat::ArrayND::linearFillLoop
void linearFillLoop(unsigned level, double s0, unsigned long idx, double shift, const double *coeffs)
Definition: ArrayND.h:4395
L1TowerCalibrationProducer_cfi.fit
fit
Definition: L1TowerCalibrationProducer_cfi.py:36
apply
Vec apply(Vec v, F f)
Definition: ExtVec.h:81
DDAxes::x
SiPixelPI::zero
Definition: SiPixelPayloadInspectorHelper.h:39
npstat::ArrayND::importSubrange
void importSubrange(const unsigned *fromCorner, unsigned lenCorner, const ArrayND< Num2, Len2, Dim2 > &from)
Definition: ArrayND.h:5447
npstat::ArrayND::marginalize
ArrayND marginalize(const ArrayND< Num2, Len2, Dim2 > &prior, const unsigned *indexMap, unsigned mapLen) const
fcn
void fcn(int &, double *, double &, double *, int)
Definition: LASBarrelAlgorithm.cc:388
npstat::ArrayND::operator/
ArrayND operator/(const Num2 &r) const
npstat::ArrayND::rank
unsigned rank() const
Definition: ArrayND.h:242
findQualityFiles.v
v
Definition: findQualityFiles.py:179
npstat::ArrayND::jointMemSliceScan
void jointMemSliceScan(Num2 *buffer, unsigned long bufLen, const unsigned *fixedIndices, const unsigned *fixedIndexValues, unsigned nFixedIndices, Functor binaryFunct)
Definition: ArrayND.h:1948
createJobs.tmp
tmp
align.sh
Definition: createJobs.py:716
npstat::ArrayND::makeNonNegative
ArrayND & makeNonNegative()
Definition: ArrayND.h:4059
watchdog.const
const
Definition: watchdog.py:83
npstat::ArrayND::data_
Numeric * data_
Definition: ArrayND.h:1049
AbsVisitor.h
Interface for piecemeal processing of a data collection.
npstat::ArrayND::inPlaceMul
ArrayND & inPlaceMul(const ArrayND< Num2, Len2, Dim2 > &r)
Definition: ArrayND.h:551
npstat::ArrayND::verifySliceCompatibility
unsigned long verifySliceCompatibility(const ArrayND< Num2, Len2, Dim2 > &slice, const unsigned *fixedIndices, const unsigned *fixedIndexValues, unsigned nFixedIndices) const
Definition: ArrayND.h:1765
svgfig.rotate
def rotate(angle, cx=0, cy=0)
Definition: svgfig.py:705
npstat::interpolate_quadratic
T interpolate_quadratic(const double x, const T &f0, const T &f1, const T &f2)
Definition: interpolate.h:33
testProducerWithPsetDescEmpty_cfi.a2
a2
Definition: testProducerWithPsetDescEmpty_cfi.py:35
heavyIonCSV_trainingSettings.idx
idx
Definition: heavyIonCSV_trainingSettings.py:5
npstat::ArrayShape
std::vector< unsigned > ArrayShape
Definition: ArrayShape.h:21
npstat::ArrayND::transposeLoop
void transposeLoop(unsigned level, unsigned pos1, unsigned pos2, unsigned long idxThis, unsigned long idxRes, ArrayND &result) const
Definition: ArrayND.h:4629
npstat::ArrayND::constFill
ArrayND & constFill(Numeric c)
Definition: ArrayND.h:4342
edmScanValgrind.buffer
buffer
Definition: edmScanValgrind.py:171
npstat::clearBuffer
void clearBuffer(T *buf, const unsigned long len)
Definition: allocators.h:86
npstat::ArrayND::marginalizeLoop
void marginalizeLoop(unsigned level, unsigned long idx, unsigned levelRes, unsigned long idxRes, const ArrayND< Num2, Len2, Dim2 > &prior, const unsigned *indexMap, ArrayND &res) const
Definition: ArrayND.h:1537
data-class-funcs.classname
classname
Definition: data-class-funcs.py:33
GetRecoTauVFromDQM_MC_cff.cl
cl
Definition: GetRecoTauVFromDQM_MC_cff.py:38
npstat::ArrayND::derivative
ArrayND derivative(double scale=1.0) const
npstat::ArrayND::isCompatible
bool isCompatible(const ArrayShape &shape) const
Definition: ArrayND.h:2370
npstat::ArrayND::ArrayND
ArrayND()
Definition: ArrayND.h:2488
npstat
Definition: AbsArrayProjector.h:14
testProducerWithPsetDescEmpty_cfi.x1
x1
Definition: testProducerWithPsetDescEmpty_cfi.py:33
npstat::ArrayND::localShape_
unsigned localShape_[StackDim]
Definition: ArrayND.h:1054
dqmdumpme.last
last
Definition: dqmdumpme.py:56
npstat::ArrayND::jointSubrangeScan
void jointSubrangeScan(ArrayND< Num2, Len2, Dim2 > &other, const unsigned *thisCorner, const unsigned *range, const unsigned *otherCorner, unsigned arrLen, Functor binaryFunct)
Definition: ArrayND.h:1504
operator-
MatrixMeschach operator-(const MatrixMeschach &mat1, const MatrixMeschach &mat2)
Definition: MatrixMeschach.cc:167
npstat::absDifference
Private::AbsReturnType< T >::type absDifference(const T &v1, const T &v2)
Definition: absDifference.h:73
npstat::NpstatRuntimeError
Definition: NpstatException.h:46
npstat::ArrayND::closest
Numeric & closest(const double *x, unsigned xDim)
Definition: ArrayND.h:3173
npstat::ArrayND::sumBelowLoop
Accumulator sumBelowLoop(unsigned level, unsigned long idx0, const unsigned *limit) const
Definition: ArrayND.h:4816
std::swap
void swap(edm::DataFrameContainer &lhs, edm::DataFrameContainer &rhs)
Definition: DataFrameContainer.h:209
npstat::ComplexComparesFalse
Definition: ComplexComparesFalse.h:23
npstat::ArrayND::version
static unsigned version()
Definition: ArrayND.h:1044
ComplexComparesAbs.h
Ordering extended to complex numbers by comparing their magnitudes.
npstat::ArrayND::scanInPlace
ArrayND & scanInPlace(Functor f)
npstat::ArrayND::data
const Numeric * data() const
Definition: ArrayND.h:236
w
const double w
Definition: UKUtility.cc:23
npstat::ArrayND::length
unsigned long length() const
Definition: ArrayND.h:233
npstat::ArrayND::maximumSpan
unsigned maximumSpan() const
Definition: ArrayND.h:5095
SimpleFunctors.h
Interface definitions and concrete simple functors for a variety of functor-based calculations.
npstat::ArrayND::functorFillLoop
void functorFillLoop(unsigned level, unsigned long idx, Functor f, unsigned *farg)
Definition: ArrayND.h:4434
npstat::ArrayND::verifyBufferSliceCompatibility
unsigned long verifyBufferSliceCompatibility(unsigned long bufLen, const unsigned *fixedIndices, const unsigned *fixedIndexValues, unsigned nFixedIndices, unsigned long *sliceStrides) const
Definition: ArrayND.h:1822
npstat::ArrayND::len_
unsigned long len_
Definition: ArrayND.h:1057
npstat::ArrayND::value
Numeric & value(const unsigned *index, unsigned indexLen)
Definition: ArrayND.h:3091
npstat::Interval
Definition: Interval.h:25
npstat::ArrayND::~ArrayND
~ArrayND()
Definition: ArrayND.h:2907
npstat::ArrayND::at
Numeric & at()
Definition: ArrayND.h:3268
AbsArrayProjector.h
Interface definition for functors used to make array projections.
trackingPlots.other
other
Definition: trackingPlots.py:1460
npstat::ArrayND::operator*=
ArrayND & operator*=(const Num2 &r)
npstat::ArrayND::project
void project(ArrayND< Num2, Len2, Dim2 > *projection, AbsArrayProjector< Numeric, Num3 > &projector, const unsigned *projectedIndices, unsigned nProjectedIndices) const
Definition: ArrayND.h:2104
npstat::ArrayND::scaleBySliceLoop
void scaleBySliceLoop(unsigned level, unsigned long idx0, unsigned level1, unsigned long idx1, ArrayND< Num2, Len2, Dim2 > &slice, const unsigned *fixedIndices, unsigned nFixedIndices, Functor binaryFunct)
Definition: ArrayND.h:2302
npstat::NpstatInvalidArgument
Definition: NpstatException.h:38
npstat::ArrayND::max
Numeric max() const
Definition: ArrayND.h:4709
npstat::ArrayND::operator+
ArrayND operator+() const
Definition: ArrayND.h:4032
npstat::AbsVisitor::process
virtual void process(const Input &value)=0
npstat::ArrayND::marginalizeInnerLoop
Numeric marginalizeInnerLoop(unsigned long idx, unsigned levelPr, unsigned long idxPr, const ArrayND< Num2, Len2, Dim2 > &prior, const unsigned *indexMap) const
Definition: ArrayND.h:1512
npstat::ArrayND::dualCircularScan
void dualCircularScan(ArrayND< Num2, Len2, Dim2 > &other, const unsigned *thisCorner, const unsigned *range, const unsigned *otherCorner, unsigned arrLen, Functor binaryFunct)
Definition: ArrayND.h:1505
npstat::ArrayND::convertToLastDimCdf
void convertToLastDimCdf(ArrayND *sumSlice, bool useTrapezoids)
Definition: ArrayND.h:4496
npstat::ArrayND::isShapeKnown
bool isShapeKnown() const
Definition: ArrayND.h:239
npstat::ArrayND::subtractFromProjection
void subtractFromProjection(ArrayND< Num2, Len2, Dim2 > *projection, AbsArrayProjector< Numeric, Num3 > &projector, const unsigned *projectedIndices, unsigned nProjectedIndices) const
Definition: ArrayND.h:2154
npstat::ArrayND::buildStrides
void buildStrides()
Definition: ArrayND.h:2478
mitigatedMETSequence_cff.U
U
Definition: mitigatedMETSequence_cff.py:36
npstat::ArrayND::scaleBySliceInnerLoop
void scaleBySliceInnerLoop(unsigned level, unsigned long idx0, Num2 &scale, const unsigned *projectedIndices, unsigned nProjectedIndices, Functor binaryFunct)
Definition: ArrayND.h:2281
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
npstat::ArrayND::coordToIndex
unsigned coordToIndex(double coord, unsigned idim) const
Definition: ArrayND.h:3147
me_macro_check_loop_prerequisites
#define me_macro_check_loop_prerequisites(METHOD, INNERLOOP)
Definition: ArrayND.h:1330
npstat::ArrayND::copyRangeLoopFunct
void copyRangeLoopFunct(unsigned level, unsigned long idx0, unsigned long idx1, const ArrayND< Num2, Len2, Dim2 > &r, const ArrayRange &range, Functor f)
Definition: ArrayND.h:2565
npstat::ArrayND::min
Numeric min() const
Definition: ArrayND.h:4675
allocators.h
Utilities related to memory management.
operator==
bool operator==(const QGLikelihoodParameters &lhs, const QGLikelihoodCategory &rhs)
Test if parameters are compatible with category.
Definition: QGLikelihoodObject.h:49
funct::true
true
Definition: Factorize.h:173
operator!=
bool operator!=(DTCELinkId const &lhs, DTCELinkId const &rhs)
Definition: DTCELinkId.h:81
operator*=
Basic3DVector & operator*=(T t)
Scaling by a scalar value (multiplication)
Definition: Basic3DVectorLD.h:195
npstat::ArrayND::interpolateLoop
Numeric interpolateLoop(unsigned level, const double *x, const Numeric *base) const
Definition: ArrayND.h:4260
operator-=
Basic3DVector & operator-=(const Basic3DVector< U > &p)
Definition: Basic3DVectorLD.h:184
npstat::ArrayND::multiMirror
void multiMirror(ArrayND< Num2, Len2, Dim2 > *out) const
Definition: ArrayND.h:4941
npstat::doubleShape
ArrayShape doubleShape(const ArrayShape &inputShape)
Definition: ArrayShape.cc:148
a
double a
Definition: hdecay.h:119
HLT_FULL_cff.maxIterations
maxIterations
Definition: HLT_FULL_cff.py:13241
clear
void clear(HadCaloObj &c)
Definition: data.h:124
npstat::ArrayND::addToProjection
void addToProjection(ArrayND< Num2, Len2, Dim2 > *projection, AbsArrayProjector< Numeric, Num3 > &projector, const unsigned *projectedIndices, unsigned nProjectedIndices) const
Definition: ArrayND.h:2129
npstat::ComplexComparesAbs
Definition: ComplexComparesAbs.h:24
npstat::ArrayND::exportSlice
void exportSlice(ArrayND< Num2, Len2, Dim2 > *slice, const unsigned *fixedIndices, const unsigned *fixedIndexValues, unsigned nFixedIndices) const
Definition: ArrayND.h:696
SiStripPI::max
Definition: SiStripPayloadInspectorHelper.h:169
runTheMatrix.mapped
mapped
Definition: runTheMatrix.py:352
npstat::BoxND
Definition: BoxND.h:24
npstat::ArrayND::classId
gs::ClassId classId() const
Definition: ArrayND.h:1039
npstat::ArrayND::sliceShape
ArrayShape sliceShape(const unsigned *fixedIndices, unsigned nFixedIndices) const
Definition: ArrayND.h:1733
npstat::ArrayND::projectInnerLoop2
void projectInnerLoop2(unsigned level, unsigned long idx0, AbsVisitor< Numeric, Num2 > &projector, const unsigned *projectedIndices, unsigned nProjectedIndices) const
Definition: ArrayND.h:2179
npstat::ArrayND::clear
ArrayND & clear()
Definition: ArrayND.h:4351
npstat::ArrayND::operator()
Numeric & operator()()
Definition: ArrayND.h:3227
recoMuon::in
Definition: RecoMuonEnumerators.h:6
npstat::destroyBuffer
void destroyBuffer(T *thisBuffer, const T *stackBuffer)
Definition: allocators.h:31
npstat::ArrayND::isShapeCompatible
bool isShapeCompatible(const ArrayND< Num2, Len2, Dim2 > &r) const
Definition: ArrayND.h:2386
npstat::ArrayND::value_type
Numeric value_type
Definition: ArrayND.h:53
npstat::ArrayND::commonSubrangeLoop
void commonSubrangeLoop(unsigned level, unsigned long idx0, unsigned long idx1, const unsigned *thisCorner, const unsigned *range, const unsigned *otherCorner, ArrayND< Num2, Len2, Dim2 > &other, Functor binaryFunct)
Definition: ArrayND.h:1359
npstat::ArrayND::span
unsigned span(unsigned dim) const
Definition: ArrayND.h:5088
npstat::ArrayND::linearFill
ArrayND & linearFill(const double *coeff, unsigned coeffLen, double c)
Definition: ArrayND.h:4416
npstat::ArrayND::interpolate1
Numeric interpolate1(const double *x, unsigned xDim) const
Definition: ArrayND.h:4211
npstat::scast_minuseq_left
Definition: SimpleFunctors.h:369
EmbeddingLHEProducer_cfi.mirror
mirror
Definition: EmbeddingLHEProducer_cfi.py:9
tolerance
const double tolerance
Definition: HGCalGeomParameters.cc:29
PreciseType.h
Compile-time deduction of an appropriate precise numeric type.
npstat::ArrayND::operator-
ArrayND operator-() const
Definition: ArrayND.h:4039
npstat::ArrayND::proper_double
ProperDblFromCmpl< Numeric >::type proper_double
Definition: ArrayND.h:54
npstat::absValue
Private::AbsReturnType< T >::type absValue(const T &v1)
Definition: absDifference.h:82
npstat::ArrayND
Definition: ArrayND.h:48
callgraph.farg
farg
Definition: callgraph.py:9
npstat::ArrayND::valueAt
Numeric & valueAt(const unsigned *index, unsigned indexLen)
Definition: ArrayND.h:3208
svgfig.template
def template(fileName, svg, replaceme="REPLACEME")
Definition: svgfig.py:521
npstat::ArrayND::dim_
unsigned dim_
Definition: ArrayND.h:1058
npstat::ProperDblFromCmpl::type
double type
Definition: ProperDblFromCmpl.h:20
npstat::ArrayND::convertLinearIndex
void convertLinearIndex(unsigned long l, unsigned *index, unsigned indexLen) const
Definition: ArrayND.h:3048
npstat::ArrayND::dotProductLoop
void dotProductLoop(unsigned level, unsigned long idx0, unsigned long idx1, unsigned long idx2, const ArrayND< Num2, Len2, Dim2 > &r, ArrayND &result) const
Definition: ArrayND.h:5019
npstat::ArrayND::strides
const unsigned long * strides() const
Definition: ArrayND.h:263
npstat::ArrayND::shapeIsKnown_
bool shapeIsKnown_
Definition: ArrayND.h:1060
res
Definition: Electron.h:6
npstat::ArrayND::maxAbsDifference
double maxAbsDifference(const ArrayND< Numeric, Len2, Dim2 > &) const
Definition: ArrayND.h:3935
npstat::ArrayND::makeCopulaSteps
unsigned makeCopulaSteps(double tolerance, unsigned maxIterations)
Definition: ArrayND.h:4073
npstat::ArrayND::strides_
unsigned long * strides_
Definition: ArrayND.h:1052
visDQMUpload.buf
buf
Definition: visDQMUpload.py:154
npstat::NpstatDomainError
Definition: NpstatException.h:54
npstat::Interval::max
const Numeric max() const
Definition: Interval.h:81
npstat::ArrayND::isZero
bool isZero() const
Definition: ArrayND.h:3034
npstat::multeq_left
Definition: SimpleFunctors.h:279
cmsLHEtoEOSManager.l
l
Definition: cmsLHEtoEOSManager.py:204
npstat::ArrayND::cl
Numeric & cl()
Definition: ArrayND.h:5245
alignCSCRings.r
r
Definition: alignCSCRings.py:93
writeEcalDQMStatus.write
write
Definition: writeEcalDQMStatus.py:48
npstat::ArrayND::linearValueAt
Numeric & linearValueAt(unsigned long index)
Definition: ArrayND.h:3133
L1DTConfigBti_cff.LL
LL
Definition: L1DTConfigBti_cff.py:25
npstat::NpstatOutOfRange
Definition: NpstatException.h:30
ArrayRange.h
Multidimensional range of array indices.
npstat::ArrayND::shape
ArrayShape shape() const
Definition: ArrayND.h:2989
npstat::ArrayND::flatCircularScan
void flatCircularScan(ArrayND< Num2, Len2, Dim2 > &other, const unsigned *thisCorner, const unsigned *range, const unsigned *otherCorner, unsigned arrLen, Functor binaryFunct)
Definition: ArrayND.h:1507
npstat::ArrayND::cdfValue
Num2 cdfValue(const unsigned *index, unsigned indexLen) const
npstat::AbsArrayProjector
Definition: AbsArrayProjector.h:20
npstat::ArrayND::setData
ArrayND & setData(const Num2 *data, unsigned long dataLength)
npstat::ArrayND::sumsq
Num2 sumsq() const
Definition: ArrayND.h:4663
remoteMonitoring_LED_IterMethod_cfg.limit
limit
Definition: remoteMonitoring_LED_IterMethod_cfg.py:427
npstat::AbsArrayProjector::result
virtual Result result()=0
edm::shift
static unsigned const int shift
Definition: LuminosityBlockID.cc:7
npstat::AbsVisitor::clear
virtual void clear()=0
Geom::operator/
T operator/(const OnePiRange< T > &a, const OnePiRange< T > &b)
Division.
Definition: OnePiRange.h:122
npstat::BoxND::dim
unsigned long dim() const
Definition: BoxND.h:54
npstat::ArrayND::projectLoop
void projectLoop(unsigned level, unsigned long idx0, unsigned level1, unsigned long idx1, unsigned *currentIndex, ArrayND< Num2, Len2, Dim2 > *projection, AbsArrayProjector< Numeric, Num3 > &projector, const unsigned *projectedIndices, unsigned nProjectedIndices, Op fcn) const
Definition: ArrayND.h:1994
relativeConstraints.value
value
Definition: relativeConstraints.py:53
genVertex_cff.x
x
Definition: genVertex_cff.py:12
multiplicitycorr_cfi.scaleFactor
scaleFactor
Definition: multiplicitycorr_cfi.py:7
npstat::scast_assign_left
Definition: SimpleFunctors.h:345
npstat::makeBuffer
T * makeBuffer(unsigned sizeNeeded, T *stackBuffer, unsigned sizeofStackBuffer)
Definition: allocators.h:22
operator/=
Basic3DVector & operator/=(T t)
Scaling by a scalar value (division)
Definition: Basic3DVectorLD.h:203
Skims_PA_cff.name
name
Definition: Skims_PA_cff.py:17
TCMET_cfi.corner
corner
Definition: TCMET_cfi.py:34
data
char data[epos_bytes_allocation]
Definition: EPOS_Wrapper.h:79
npstat::AbsVisitor::result
virtual Result result()=0
bookConverter.prior
prior
Definition: bookConverter.py:146
npstat::ArrayND::addmul
ArrayND & addmul(const ArrayND< Num2, Len2, Dim2 > &r, const Num3 &c)
AlignmentPI::index
index
Definition: AlignmentPayloadInspectorHelper.h:46
ComplexComparesFalse.h
Ordering extended to complex numbers by always returning "false".
npstat::ArrayND::flatCircularLoop
void flatCircularLoop(unsigned level, unsigned long idx0, unsigned long idx1, const unsigned *thisCorner, const unsigned *range, const unsigned *otherCorner, ArrayND< Num2, Len2, Dim2 > &other, Functor binaryFunct)
Definition: ArrayND.h:1428
npstat::ArrayND::fullRange
ArrayRange fullRange() const
Definition: ArrayND.h:2996
MillePedeFileConverter_cfg.out
out
Definition: MillePedeFileConverter_cfg.py:31
ztail.d
d
Definition: ztail.py:151
npstat::ArrayND::assign
ArrayND & assign(const ArrayND< Num2, Len2, Dim2 > &, Functor f)
mps_fire.result
result
Definition: mps_fire.py:311
interpolate.h
Low-order polynomial interpolation (up to cubic) for equidistant coordinates.
npstat::AbsArrayProjector::process
virtual void process(const unsigned *index, unsigned indexLen, unsigned long linearIndex, const Input &value)=0
npstat::interpolate_linear
T interpolate_linear(const double x, const T &f0, const T &f1)
Definition: interpolate.h:23
npstat::ArrayND::linearValue
Numeric & linearValue(unsigned long index)
Definition: ArrayND.h:3123
npstat::ArrayND::apply
ArrayND & apply(Functor f)
npstat::ArrayND::dot
ArrayND dot(const ArrayND< Num2, Len2, Dim2 > &r) const
SurfaceOrientation::outer
Definition: Surface.h:19
c
auto & c
Definition: CAHitNtupletGeneratorKernelsImpl.h:46
dqmiolumiharvest.j
j
Definition: dqmiolumiharvest.py:66
operator+=
Basic3DVector & operator+=(const Basic3DVector< U > &p)
Definition: Basic3DVectorLD.h:174
npstat::ArrayND::operator-=
ArrayND & operator-=(const ArrayND< Num2, Len2, Dim2 > &r)
npstat::ArrayND::projectInnerLoop
void projectInnerLoop(unsigned level, unsigned long idx0, unsigned *currentIndex, AbsArrayProjector< Numeric, Num2 > &projector, const unsigned *projectedIndices, unsigned nProjectedIndices) const
Definition: ArrayND.h:1970
NpstatException.h
Exceptions for the npstat namespace.
newFWLiteAna.base
base
Definition: newFWLiteAna.py:92
npstat::ArrayND::processSubrangeLoop
void processSubrangeLoop(unsigned level, unsigned long idx0, unsigned *currentIndex, AbsArrayProjector< Numeric, Num2 > &f, const BoxND< Integer > &subrange) const
Definition: ArrayND.h:2408
mergeAndRegister.slice
slice
Definition: mergeAndRegister.py:121
npstat::ArrayND::rotate
void rotate(const unsigned *shifts, unsigned lenShifts, ArrayND< Num2, Len2, Dim2 > *rotated) const
Definition: ArrayND.h:4989
PVValHelper::dx
Definition: PVValidationHelpers.h:49
fit
Definition: CombinedChiSquaredLikelihood.h:6
operator+
MatrixMeschach operator+(const MatrixMeschach &mat1, const MatrixMeschach &mat2)
Definition: MatrixMeschach.cc:154
npstat::ArrayND::linearIndex
unsigned long linearIndex(const unsigned *idx, unsigned idxLen) const
Definition: ArrayND.h:3070
npstat::ArrayND::processSubrange
void processSubrange(AbsArrayProjector< Numeric, Num2 > &f, const BoxND< Integer > &subrange) const
Definition: ArrayND.h:2447
gpuPixelDoublets::stride
auto stride
Definition: gpuPixelDoubletsAlgos.h:80
mps_fire.dest
dest
Definition: mps_fire.py:179
npstat::scast_assign_right
Definition: SimpleFunctors.h:351
npstat::ArrayND::buildFromShapePtr
void buildFromShapePtr(const unsigned *, unsigned)
Definition: ArrayND.h:1641
npstat::ArrayND::localStrides_
unsigned long localStrides_[StackDim]
Definition: ArrayND.h:1051
npstat::ArrayND::makeUnit
ArrayND & makeUnit()
Definition: ArrayND.h:4371
npstat::ArrayND::outer
ArrayND outer(const ArrayND< Num2, Len2, Dim2 > &r) const
Accumulator
Definition: Accumulator.h:12
Generated for CMSSW Reference Manual by doxygen 1.8.16