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