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