Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:21:29

0001 #ifndef NPSTAT_ARRAYND_HH_
0002 #define NPSTAT_ARRAYND_HH_
0003 
0004 /*!
0005 // \file ArrayND.h
0006 //
0007 // \brief Arbitrary-dimensional array template
0008 //
0009 // Author: I. Volobouev
0010 //
0011 // October 2009
0012 */
0013 
0014 #include <cassert>
0015 
0016 #include "Alignment/Geners/interface/ClassId.hh"
0017 
0018 #include "JetMETCorrections/InterpolationTables/interface/SimpleFunctors.h"
0019 #include "JetMETCorrections/InterpolationTables/interface/ArrayRange.h"
0020 #include "JetMETCorrections/InterpolationTables/interface/AbsArrayProjector.h"
0021 #include "JetMETCorrections/InterpolationTables/interface/AbsVisitor.h"
0022 #include "JetMETCorrections/InterpolationTables/interface/PreciseType.h"
0023 #include "JetMETCorrections/InterpolationTables/interface/ProperDblFromCmpl.h"
0024 
0025 namespace npstat {
0026   /**
0027     // A class for multidimensional array manipulation. A number of methods
0028     // of this class will work only if dimensionality is limited by
0029     // CHAR_BIT*sizeof(unsigned long)-1 (which is 31 and 63 on 32- and 64-bit
0030     // architectures, respectively).
0031     //
0032     // Depending on how much space is provided with the "StackLen" template
0033     // parameter, the array data will be placed either on the stack or on the
0034     // heap. By default, array data leaves on the heap unless the array has
0035     // rank 0.
0036     //
0037     // Depending on how much space is provided with the "StackDim" template
0038     // parameter, the array strides will be placed either on the stack or
0039     // on the heap. By default, strides will be placed on the stack in case
0040     // the array dimensionality is ten or less.
0041     //
0042     // The "Numeric" type must have a default constructor (of course,
0043     // pointers to arbitrary types can be used as well).
0044     //
0045     // Both StackLen and StackDim parameters must be positive.
0046     */
0047   template <typename Numeric, unsigned StackLen = 1U, unsigned StackDim = 10U>
0048   class ArrayND {
0049     template <typename Num2, unsigned Len2, unsigned Dim2>
0050     friend class ArrayND;
0051 
0052   public:
0053     typedef Numeric value_type;
0054     typedef typename ProperDblFromCmpl<Numeric>::type proper_double;
0055 
0056     /**
0057         // Default constructor creates an uninitialized array.
0058         // Only three things can be done safely with such an array:
0059         //
0060         // 1) Assigning it from another array (initialized or not).
0061         //
0062         // 2) Passing it as an argument to the class static method "restore".
0063         //
0064         // 3) Calling the "uninitialize" method.
0065         //
0066         // Any other operation results in an undefined behavior (often,
0067         // an exception is thrown). Note that initialized array can not
0068         // be assigned from uninitialized one.
0069         */
0070     ArrayND();
0071 
0072     //@{
0073     /**
0074         // Constructor which creates arrays with the given shape.
0075         // The array data remains undefined. Simple inilitalization
0076         // of the data can be performed using methods clear() or
0077         // constFill(SomeValue). More complicated initialization
0078         // can be done by "linearFill", "functorFill", or by setting
0079         // every array element to a desired value.
0080         */
0081     explicit ArrayND(const ArrayShape& shape);
0082     ArrayND(const unsigned* shape, unsigned dim);
0083     //@}
0084 
0085     /** The copy constructor */
0086     ArrayND(const ArrayND&);
0087 
0088     /**
0089         // Converting constructor. It looks more general than the copy
0090         // constructor, but the actual copy constructor has to be created
0091         // anyway -- otherwise the compiler will generate an incorrect
0092         // default copy constructor. Note that existence of this
0093         // constructor essentially disables data type safety for copying
0094         // arrays -- but the code significantly gains in convenience.
0095         */
0096     template <typename Num2, unsigned Len2, unsigned Dim2>
0097     ArrayND(const ArrayND<Num2, Len2, Dim2>&);
0098 
0099     /**
0100         // Converting constructor where the array values are filled
0101         // by a functor using values of another array as arguments
0102         */
0103     template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
0104     ArrayND(const ArrayND<Num2, Len2, Dim2>&, Functor f);
0105 
0106     /** Constructor from a subrange of another array */
0107     template <typename Num2, unsigned Len2, unsigned Dim2>
0108     ArrayND(const ArrayND<Num2, Len2, Dim2>& from, const ArrayRange& fromRange);
0109 
0110     /** Similar constructor with a transforming functor */
0111     template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
0112     ArrayND(const ArrayND<Num2, Len2, Dim2>& from, const ArrayRange& fromRange, Functor f);
0113 
0114     /**
0115         // Constructor from a slice of another array. The data of the
0116         // constructed array remains undefined. The argument "indices"
0117         // lists either the array indices whose numbers will be fixed
0118         // when slicing is performed or the indices which will be iterated
0119         // over during projections (for example, array values may be
0120         // summed over these indices). These indices will be excluded
0121         // from the constructed array. The created array can be subsequently
0122         // used with methods "exportSlice", "importSlice", "project", etc.
0123         // of the parent array "slicedArray".
0124         */
0125     template <typename Num2, unsigned Len2, unsigned Dim2>
0126     ArrayND(const ArrayND<Num2, Len2, Dim2>& slicedArray, const unsigned* indices, unsigned nIndices);
0127 
0128     /** Outer product constructor */
0129     template <typename Num1, unsigned Len1, unsigned Dim1, typename Num2, unsigned Len2, unsigned Dim2>
0130     ArrayND(const ArrayND<Num1, Len1, Dim1>& a1, const ArrayND<Num2, Len2, Dim2>& a2);
0131 
0132     //@{
0133     /** 
0134         // Constructor in which the spans are explicitly provided
0135         // for each dimension. The array data remains undefined.
0136         */
0137     explicit ArrayND(unsigned n0);
0138     ArrayND(unsigned n0, unsigned n1);
0139     ArrayND(unsigned n0, unsigned n1, unsigned n2);
0140     ArrayND(unsigned n0, unsigned n1, unsigned n2, unsigned n3);
0141     ArrayND(unsigned n0, unsigned n1, unsigned n2, unsigned n3, unsigned n4);
0142     ArrayND(unsigned n0, unsigned n1, unsigned n2, unsigned n3, unsigned n4, unsigned n5);
0143     ArrayND(unsigned n0, unsigned n1, unsigned n2, unsigned n3, unsigned n4, unsigned n5, unsigned n6);
0144     ArrayND(unsigned n0, unsigned n1, unsigned n2, unsigned n3, unsigned n4, unsigned n5, unsigned n6, unsigned n7);
0145     ArrayND(unsigned n0,
0146             unsigned n1,
0147             unsigned n2,
0148             unsigned n3,
0149             unsigned n4,
0150             unsigned n5,
0151             unsigned n6,
0152             unsigned n7,
0153             unsigned n8);
0154     ArrayND(unsigned n0,
0155             unsigned n1,
0156             unsigned n2,
0157             unsigned n3,
0158             unsigned n4,
0159             unsigned n5,
0160             unsigned n6,
0161             unsigned n7,
0162             unsigned n8,
0163             unsigned n9);
0164     //@}
0165 
0166     /** Destructor */
0167     ~ArrayND();
0168 
0169     /**
0170         // Assignment operator. The shape of the array on the right
0171         // must be compatible with the shape of the array on the left.
0172         // The only exception is when the array on the left has no shape
0173         // at all (i.e., it was created by the default constructor or
0174         // its "uninitialize" method was called). In this case the array
0175         // on the left will assume the shape of the array on the right.
0176         */
0177     ArrayND& operator=(const ArrayND&);
0178 
0179     /** Converting assignment operator */
0180     template <typename Num2, unsigned Len2, unsigned Dim2>
0181     ArrayND& operator=(const ArrayND<Num2, Len2, Dim2>&);
0182 
0183     /** Converting assignment method with a transforming functor */
0184     template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
0185     ArrayND& assign(const ArrayND<Num2, Len2, Dim2>&, Functor f);
0186 
0187     /**
0188         // The function which can "uninitialize" the array to the same
0189         // state as produced by the default constructor. Can be applied
0190         // in order to force the assignment operators to work.
0191         */
0192     ArrayND& uninitialize();
0193 
0194     //@{
0195     /**
0196         // Element access using multidimensional array index
0197         // (no bounds checking). The length of the index array
0198         // must be equal to the rank of this object.
0199         */
0200     Numeric& value(const unsigned* index, unsigned indexLen);
0201     const Numeric& value(const unsigned* index, unsigned indexLen) const;
0202     //@}
0203 
0204     //@{
0205     /**
0206         // Element access using multidimensional array index
0207         // (with bounds checking)
0208         */
0209     Numeric& valueAt(const unsigned* index, unsigned indexLen);
0210     const Numeric& valueAt(const unsigned* index, unsigned indexLen) const;
0211     //@}
0212 
0213     //@{
0214     /** Element access using linear index (no bounds checking) */
0215     Numeric& linearValue(unsigned long index);
0216     const Numeric& linearValue(unsigned long index) const;
0217     //@}
0218 
0219     //@{
0220     /** Element access using linear index (with bounds checking) */
0221     Numeric& linearValueAt(unsigned long index);
0222     const Numeric& linearValueAt(unsigned long index) const;
0223     //@}
0224 
0225     /** Convert linear index into multidimensional index */
0226     void convertLinearIndex(unsigned long l, unsigned* index, unsigned indexLen) const;
0227 
0228     /** Convert multidimensional index into linear index */
0229     unsigned long linearIndex(const unsigned* idx, unsigned idxLen) const;
0230 
0231     // Some inspectors
0232     /** Total number of data array elements */
0233     inline unsigned long length() const { return len_; }
0234 
0235     /** Linearized data */
0236     inline const Numeric* data() const { return data_; }
0237 
0238     /** Check whether the array has been initialized */
0239     inline bool isShapeKnown() const { return shapeIsKnown_; }
0240 
0241     /** The number of array dimensions */
0242     inline unsigned rank() const { return dim_; }
0243 
0244     /** Get the complete shape */
0245     ArrayShape shape() const;
0246 
0247     /** Shape data as a C-style array */
0248     inline const unsigned* shapeData() const { return shape_; }
0249 
0250     /** Get the complete range */
0251     ArrayRange fullRange() const;
0252 
0253     /** Get the number of elements in some particular dimension */
0254     unsigned span(unsigned dim) const;
0255 
0256     /** Maximum span among all dimensions */
0257     unsigned maximumSpan() const;
0258 
0259     /** Minimum span among all dimensions */
0260     unsigned minimumSpan() const;
0261 
0262     /** Get the strides */
0263     inline const unsigned long* strides() const { return strides_; }
0264 
0265     /** Check if all array elements are zero */
0266     bool isZero() const;
0267 
0268     /**
0269         // This method checks whether all array elements are
0270         // non-negative and, in addition, there is at least
0271         // one positive element
0272         */
0273     bool isDensity() const;
0274 
0275     /** This method modifies all the data in one statement */
0276     template <typename Num2>
0277     ArrayND& setData(const Num2* data, unsigned long dataLength);
0278 
0279     /** Compare two arrays for equality */
0280     template <unsigned Len2, unsigned Dim2>
0281     bool operator==(const ArrayND<Numeric, Len2, Dim2>&) const;
0282 
0283     /** Logical negation of operator== */
0284     template <unsigned Len2, unsigned Dim2>
0285     bool operator!=(const ArrayND<Numeric, Len2, Dim2>&) const;
0286 
0287     /** Largest absolute difference with another bin-compatible array */
0288     template <unsigned Len2, unsigned Dim2>
0289     double maxAbsDifference(const ArrayND<Numeric, Len2, Dim2>&) const;
0290 
0291     /** operator+ returns a copy of this array */
0292     ArrayND operator+() const;
0293 
0294     /** operator- applies the unary minus operator to every element */
0295     ArrayND operator-() const;
0296 
0297     /** addition of two arrays */
0298     template <unsigned Len2, unsigned Dim2>
0299     ArrayND operator+(const ArrayND<Numeric, Len2, Dim2>& r) const;
0300 
0301     /** subtraction of two arrays */
0302     template <unsigned Len2, unsigned Dim2>
0303     ArrayND operator-(const ArrayND<Numeric, Len2, Dim2>& r) const;
0304 
0305     /** multiplication by a scalar */
0306     template <typename Num2>
0307     ArrayND operator*(const Num2& r) const;
0308 
0309     /** division by a scalar */
0310     template <typename Num2>
0311     ArrayND operator/(const Num2& r) const;
0312 
0313     //@{
0314     /**
0315         // In-place operators. Note that these work faster than the binary
0316         // versions, i.e., A += B is much faster than A = A + B.
0317         */
0318     template <typename Num2>
0319     ArrayND& operator*=(const Num2& r);
0320 
0321     template <typename Num2>
0322     ArrayND& operator/=(const Num2& r);
0323 
0324     template <typename Num2, unsigned Len2, unsigned Dim2>
0325     ArrayND& operator+=(const ArrayND<Num2, Len2, Dim2>& r);
0326 
0327     template <typename Num2, unsigned Len2, unsigned Dim2>
0328     ArrayND& operator-=(const ArrayND<Num2, Len2, Dim2>& r);
0329     //@}
0330 
0331     /** This method is equivalent to (but faster than) += r*c */
0332     template <typename Num3, typename Num2, unsigned Len2, unsigned Dim2>
0333     ArrayND& addmul(const ArrayND<Num2, Len2, Dim2>& r, const Num3& c);
0334 
0335     /** Outer product as a method (see also the outer product constructor) */
0336     template <typename Num2, unsigned Len2, unsigned Dim2>
0337     ArrayND outer(const ArrayND<Num2, Len2, Dim2>& r) const;
0338 
0339     /**
0340         // Contraction of a pair of indices. Note that the array length
0341         // must be the same in both dimensions.
0342         */
0343     ArrayND contract(unsigned pos1, unsigned pos2) const;
0344 
0345     /**
0346         // Here, dot product corresponds to outer product followed
0347         // by the contraction over two indices -- the last index
0348         // of this object and the first index of the argument.
0349         */
0350     template <typename Num2, unsigned Len2, unsigned Dim2>
0351     ArrayND dot(const ArrayND<Num2, Len2, Dim2>& r) const;
0352 
0353     /**
0354         // The intent of this method is to marginalize
0355         // over a set of indices with a prior. Essentially, we are
0356         // calculating integrals akin to p(y) = Integral f(y|x) g(x) dx
0357         // in which all functions are represented on an equidistant grid.
0358         // If needed, multiplication of the result by the grid cell size
0359         // should be performed after this function. "indexMap" specifies
0360         // how the indices of the prior array (which is like g(x)) are
0361         // mapped into the indices of this array (which is like f(y|x)).
0362         // The number of elements in the map, "mapLen", must be equal to
0363         // the rank of the prior. Dimension 0 of the prior corresponds
0364         // to the dimension indexMap[0] of this array, dimension 1
0365         // corresponds to indexMap[1], etc.
0366         */
0367     template <typename Num2, unsigned Len2, unsigned Dim2>
0368     ArrayND marginalize(const ArrayND<Num2, Len2, Dim2>& prior, const unsigned* indexMap, unsigned mapLen) const;
0369 
0370     /** Transposed array */
0371     ArrayND transpose(unsigned pos1, unsigned pos2) const;
0372 
0373     /** Transpose without arguments can be invoked for 2-d arrays only */
0374     ArrayND transpose() const;
0375 
0376     /**
0377         // Sum of all array elements which uses Num2 type as accumulator.
0378         // Typically, the precision and dynamic range of Num2 should be
0379         // suitably larger than the precision and dynamic range of Numeric.
0380         // For example, if Numeric is float then Num2 should be double, etc.
0381         */
0382     template <typename Num2>
0383     Num2 sum() const;
0384 
0385     /**
0386         // Sum of absolute values squared which uses Num2 as accumulator.
0387         // Function std::abs(Numeric) must exist.
0388         */
0389     template <typename Num2>
0390     Num2 sumsq() const;
0391 
0392     /**
0393         // Mixed derivative over all directions. Useful for generating
0394         // densities from distribution functions. The resulting array
0395         // will have one less point in each dimension. Class Num2 is
0396         // used as accumulator for calculations. static_cast from
0397         // Num2 to Numeric must exist. The result is multiplied by the
0398         // scale factor provided.
0399         */
0400     template <typename Num2>
0401     ArrayND derivative(double scale = 1.0) const;
0402 
0403     /**
0404         // The operation inverse to "derivative". Constructs multivariate
0405         // cumulative density function.
0406         */
0407     template <typename Num2>
0408     ArrayND cdfArray(double scale = 1.0) const;
0409 
0410     /**
0411         // Calculate just one multivariate cumulative density function
0412         // value. Point with given index will be included in the sum.
0413         */
0414     template <typename Num2>
0415     Num2 cdfValue(const unsigned* index, unsigned indexLen) const;
0416 
0417     /**
0418         // The next function turns the array data into the conditional
0419         // cumulative density function for the last dimension. "Num2"
0420         // is the type of accumulator class used. The cdf is stored
0421         // in such a way that the cdf value of 0 is skipped (the first
0422         // stored value is the sum which includes the 0th bin). The slice
0423         // is filled with the sum of values. The "useTrapezoids" parameter
0424         // specifies whether trapezoidal integration formula should be
0425         // utilized (rectangular integration is used in case
0426         // "useTrapezoids" value is "false").
0427         */
0428     template <typename Num2>
0429     void convertToLastDimCdf(ArrayND* sumSlice, bool useTrapezoids);
0430 
0431     /** Minimum array element */
0432     Numeric min() const;
0433 
0434     /** Minimum array element and its index */
0435     Numeric min(unsigned* index, unsigned indexLen) const;
0436 
0437     /** Maximum array element */
0438     Numeric max() const;
0439 
0440     /** Maximum array element and its index */
0441     Numeric max(unsigned* index, unsigned indexLen) const;
0442 
0443     //@{
0444     /**
0445         // Closest value accessor (works as if the array allows access
0446         // with non-integer indices). For example, the second point
0447         // in some dimension will be accessed in case the coordinate
0448         // in that dimension is between 0.5 and 1.5. This function can be
0449         // used, for example, for implementing simple N-D histogramming
0450         // or for closest value interpolation and extrapolation.
0451         */
0452     Numeric& closest(const double* x, unsigned xDim);
0453     const Numeric& closest(const double* x, unsigned xDim) const;
0454     //@}
0455 
0456     /**
0457         // Multilinear interpolation. Closest value extrapolation is used
0458         // in case some index is outside of the array bounds. Note that
0459         // this function works only if the array dimensionality is less
0460         // than CHAR_BIT*sizeof(unsigned long). x is the "coordinate"
0461         // which coincides with array index for x equal to unsigned
0462         // integers.
0463         */
0464     Numeric interpolate1(const double* x, unsigned xDim) const;
0465 
0466     /**
0467         // Multicubic interpolation. Closest value extrapolation is used
0468         // in case some index is outside of the array bounds. This
0469         // function is much slower than "interpolate1" (in the current
0470         // implementation, a recursive algorithm is used).
0471         */
0472     Numeric interpolate3(const double* x, unsigned xDim) const;
0473 
0474     /**
0475         // This method applies a single-argument functor to each
0476         // element of the array (in-place). The result returned
0477         // by the functor becomes the new value of the element. There
0478         // must be a conversion (static cast) from the functor result to
0479         // the "Numeric" type. The method returns *this which allows
0480         // for chaining of such methods. Use the transforming constructor
0481         // if you want a new array instead.
0482         */
0483     template <class Functor>
0484     ArrayND& apply(Functor f);
0485 
0486     /**
0487         // This method applies a single-argument functor to each
0488         // element of the array. The result returned by the functor
0489         // is ignored inside the scan. Depending on what the functor does,
0490         // the array values may or may not be modified (they can be modified
0491         // if the functor takes its argument via a non-const reference).
0492         */
0493     template <class Functor>
0494     ArrayND& scanInPlace(Functor f);
0495 
0496     /** This method fills the array data with a constant value */
0497     ArrayND& constFill(Numeric c);
0498 
0499     /** Zero the array out (every datum becomes Numeric()) */
0500     ArrayND& clear();
0501 
0502     /**
0503         // This method fills the array with a linear combination
0504         // of the index values. For example, a 2-d array element with indices
0505         // i, k will be set to (coeff[0]*i + coeff[1]*k + c). There must be
0506         // a conversion (static cast) from double into "Numeric".
0507         */
0508     ArrayND& linearFill(const double* coeff, unsigned coeffLen, double c);
0509 
0510     /**
0511         // This method fills the array from a functor
0512         // which takes (const unsigned* index, unsigned indexLen)
0513         // arguments. There must be a conversion (static cast) from
0514         // the functor result to the "Numeric" type.
0515         */
0516     template <class Functor>
0517     ArrayND& functorFill(Functor f);
0518 
0519     /**
0520         // This method can be used for arrays with rank
0521         // of at least 2 whose length is the same in all dimensions.
0522         // It puts static_cast<Numeric>(1) on the main diagonal and
0523         // Numeric() everywhere else.
0524         */
0525     ArrayND& makeUnit();
0526 
0527     /** This method turns all negative elements into zeros */
0528     ArrayND& makeNonNegative();
0529 
0530     /**
0531         // This method accumulates marginals and divides
0532         // the array (treated as a distribution) by the product of the
0533         // marginals. Several iterations like this turn the distribution
0534         // into a copula. If the array contains negative elements, they
0535         // are turned into zeros before the iterations are performed.
0536         // The function returns the actual number of iteration performed
0537         // when the given tolerance was reached for all marginals.
0538         */
0539     unsigned makeCopulaSteps(double tolerance, unsigned maxIterations);
0540 
0541     /**
0542         // Loop over all elements of two compatible arrays and apply
0543         // a binary functor
0544         */
0545     template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
0546     void jointScan(ArrayND<Num2, Len2, Dim2>& other, Functor binaryFunct);
0547 
0548     /** Convenience method for element-by-element in-place multiplication
0549         */
0550     template <typename Num2, unsigned Len2, unsigned Dim2>
0551     inline ArrayND& inPlaceMul(const ArrayND<Num2, Len2, Dim2>& r) {
0552       jointScan(const_cast<ArrayND<Num2, Len2, Dim2>&>(r), multeq_left<Numeric, Num2>());
0553       return *this;
0554     }
0555 
0556     /**
0557         // Loop over subranges in two arrays in such a way that the functor
0558         // is called only if the indices on both sides are valid. The topology
0559         // of both arrays is assumed to be box-like (flat). The starting
0560         // corner in this object (where cycling begins) is provided by the
0561         // argument "thisCorner". The "range" argument specifies the width
0562         // of the processed patch in each dimension. The corner of the "other"
0563         // array where cycling begins is provided by the "otherCorner"
0564         // argument. The "arrLen" argument specifies the number of elements
0565         // in "thisCorner", "range", and "otherCorner" arrays. It should be
0566         // equal to the rank of either of the two ArrayND arrays.
0567         //
0568         // Note that there is no good way for this method to assume constness
0569         // of this or "other" array: this becomes apparent only after the 
0570         // functor has been specified. Apply const_cast judiciously as needed,
0571         // other solutions of this problem are not any better.
0572         */
0573     template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
0574     void jointSubrangeScan(ArrayND<Num2, Len2, Dim2>& other,
0575                            const unsigned* thisCorner,
0576                            const unsigned* range,
0577                            const unsigned* otherCorner,
0578                            unsigned arrLen,
0579                            Functor binaryFunct);
0580 
0581     /**
0582         // Method similar to "jointSubrangeScan" in which the topology of
0583         // both arrays is assumed to be hypertoroidal (circular buffer in
0584         // every dimension)
0585         */
0586     template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
0587     void dualCircularScan(ArrayND<Num2, Len2, Dim2>& other,
0588                           const unsigned* thisCorner,
0589                           const unsigned* range,
0590                           const unsigned* otherCorner,
0591                           unsigned arrLen,
0592                           Functor binaryFunct);
0593 
0594     /**
0595         // Method similar to "jointSubrangeScan" in which the topology of
0596         // this array is assumed to be flat and the other array hypertoroidal
0597         */
0598     template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
0599     void flatCircularScan(ArrayND<Num2, Len2, Dim2>& other,
0600                           const unsigned* thisCorner,
0601                           const unsigned* range,
0602                           const unsigned* otherCorner,
0603                           unsigned arrLen,
0604                           Functor binaryFunct);
0605 
0606     /**
0607         // Method similar to "jointSubrangeScan" in which the topology of
0608         // this array is assumed to be hypertoroidal and the other array flat
0609         */
0610     template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
0611     void circularFlatScan(ArrayND<Num2, Len2, Dim2>& other,
0612                           const unsigned* thisCorner,
0613                           const unsigned* range,
0614                           const unsigned* otherCorner,
0615                           unsigned arrLen,
0616                           Functor binaryFunct);
0617 
0618     /**
0619         // This method runs over a subrange of the array
0620         // and calls the argument functor on every point. This
0621         // method will not call "clear" or "result" functions of
0622         // the argument functor.
0623         */
0624     template <typename Num2, typename Integer>
0625     void processSubrange(AbsArrayProjector<Numeric, Num2>& f, const BoxND<Integer>& subrange) const;
0626 
0627     /**
0628         // Copy a hyperrectangular subrange of this array potentially
0629         // completely overwriting the destination array. The starting
0630         // corner in this object where copying begins is provided by
0631         // the first two arguments. The subrange size is defined by
0632         // the shape of the destination array.
0633         */
0634     template <typename Num2, unsigned Len2, unsigned Dim2>
0635     void exportSubrange(const unsigned* fromCorner, unsigned lenCorner, ArrayND<Num2, Len2, Dim2>* dest) const;
0636 
0637     /** The inverse operation to "exportSubrange" */
0638     template <typename Num2, unsigned Len2, unsigned Dim2>
0639     void importSubrange(const unsigned* fromCorner, unsigned lenCorner, const ArrayND<Num2, Len2, Dim2>& from);
0640 
0641     /**
0642         // Check that all elements of this array differ from the
0643         // corresponding elements of another array by at most "eps".
0644         // Equivalent to maxAbsDifference(r) <= eps (but usually faster).
0645         */
0646     template <typename Num2, unsigned Len2, unsigned Dim2>
0647     bool isClose(const ArrayND<Num2, Len2, Dim2>& r, double eps) const;
0648 
0649     /** Check compatibility with another shape */
0650     bool isCompatible(const ArrayShape& shape) const;
0651 
0652     /**
0653         // Check shape compatibility with another array. Equivalent to
0654         // but faster than isCompatible(r.shape()).
0655         */
0656     template <typename Num2, unsigned Len2, unsigned Dim2>
0657     bool isShapeCompatible(const ArrayND<Num2, Len2, Dim2>& r) const;
0658 
0659     /**
0660         // Joint cycle over the data of this array and the slice.
0661         // The array to which the "slice" argument refers should normally
0662         // be created by the slicing constructor using this array as
0663         // the argument. The "fixedIndices" argument should be the same
0664         // as the "indices" argument in that constructor. This method
0665         // is to be used for import/export of slice data and in-place
0666         // operations (addition, multiplication, etc).
0667         */
0668     template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
0669     void jointSliceScan(ArrayND<Num2, Len2, Dim2>& slice,
0670                         const unsigned* fixedIndices,
0671                         const unsigned* fixedIndexValues,
0672                         unsigned nFixedIndices,
0673                         Functor binaryFunct);
0674 
0675     /**
0676         // Joint cycle over a slice this array and some memory buffer.
0677         // The length of the buffer must be equal to the length of the
0678         // slice, as in "jointSliceScan". This method is to be used
0679         // for import/export of slice data and in-place operations
0680         // (addition, multiplication, etc) with memory managed not
0681         // by ArrayND but in some other manner.
0682         */
0683     template <typename Num2, class Functor>
0684     void jointMemSliceScan(Num2* buffer,
0685                            unsigned long bufLen,
0686                            const unsigned* fixedIndices,
0687                            const unsigned* fixedIndexValues,
0688                            unsigned nFixedIndices,
0689                            Functor binaryFunct);
0690 
0691     /** Figure out the slice shape without actually making the slice */
0692     ArrayShape sliceShape(const unsigned* fixedIndices, unsigned nFixedIndices) const;
0693 
0694     /** Convenience method for exporting a slice of this array */
0695     template <typename Num2, unsigned Len2, unsigned Dim2>
0696     inline void exportSlice(ArrayND<Num2, Len2, Dim2>* slice,
0697                             const unsigned* fixedIndices,
0698                             const unsigned* fixedIndexValues,
0699                             unsigned nFixedIndices) const {
0700       assert(slice);
0701       (const_cast<ArrayND*>(this))
0702           ->jointSliceScan(*slice, fixedIndices, fixedIndexValues, nFixedIndices, scast_assign_right<Numeric, Num2>());
0703     }
0704 
0705     /** 
0706         // Convenience method for exporting a slice of this array
0707         // into a memory buffer
0708         */
0709     template <typename Num2>
0710     inline void exportMemSlice(Num2* buffer,
0711                                unsigned long bufLen,
0712                                const unsigned* fixedIndices,
0713                                const unsigned* fixedIndexValues,
0714                                unsigned nFixedIndices) const {
0715       (const_cast<ArrayND*>(this))
0716           ->jointMemSliceScan(
0717               buffer, bufLen, fixedIndices, fixedIndexValues, nFixedIndices, scast_assign_right<Numeric, Num2>());
0718     }
0719 
0720     /** Convenience method for importing a slice into this array */
0721     template <typename Num2, unsigned Len2, unsigned Dim2>
0722     inline void importSlice(const ArrayND<Num2, Len2, Dim2>& slice,
0723                             const unsigned* fixedIndices,
0724                             const unsigned* fixedIndexValues,
0725                             unsigned nFixedIndices) {
0726       jointSliceScan(const_cast<ArrayND<Num2, Len2, Dim2>&>(slice),
0727                      fixedIndices,
0728                      fixedIndexValues,
0729                      nFixedIndices,
0730                      scast_assign_left<Numeric, Num2>());
0731     }
0732 
0733     /**
0734         // Convenience method for importing a slice into this array
0735         // from a memory buffer
0736         */
0737     template <typename Num2>
0738     inline void importMemSlice(const Num2* buffer,
0739                                unsigned long bufLen,
0740                                const unsigned* fixedIndices,
0741                                const unsigned* fixedIndexValues,
0742                                unsigned nFixedIndices) {
0743       jointMemSliceScan(const_cast<Num2*>(buffer),
0744                         bufLen,
0745                         fixedIndices,
0746                         fixedIndexValues,
0747                         nFixedIndices,
0748                         scast_assign_left<Numeric, Num2>());
0749     }
0750 
0751     /**
0752         // This method applies the values in the slice to all other
0753         // coresponding values in the array. This can be used, for example,
0754         // to multiply/divide by some factor which varies across the slice.
0755         // The slice values will be used as the right functor argument.
0756         */
0757     template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
0758     void applySlice(ArrayND<Num2, Len2, Dim2>& slice,
0759                     const unsigned* fixedIndices,
0760                     unsigned nFixedIndices,
0761                     Functor binaryFunct);
0762 
0763     /**
0764         // Convenience method which multiplies the array by a scale factor
0765         // which varies across the slice
0766         */
0767     template <typename Num2, unsigned Len2, unsigned Dim2>
0768     inline ArrayND& multiplyBySlice(const ArrayND<Num2, Len2, Dim2>& slice,
0769                                     const unsigned* fixedIndices,
0770                                     unsigned nFixedIndices) {
0771       applySlice(
0772           const_cast<ArrayND<Num2, Len2, Dim2>&>(slice), fixedIndices, nFixedIndices, multeq_left<Numeric, Num2>());
0773       return *this;
0774     }
0775 
0776     //@{
0777     /**
0778         // This method fills a projection. The array to which
0779         // "projection" argument points should normally be created by
0780         // the slicing constructor using this array as an argument.
0781         // "projectedIndices" should be the same as "indices" specified
0782         // during the slice creation.
0783         */
0784     template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
0785     void project(ArrayND<Num2, Len2, Dim2>* projection,
0786                  AbsArrayProjector<Numeric, Num3>& projector,
0787                  const unsigned* projectedIndices,
0788                  unsigned nProjectedIndices) const;
0789 
0790     template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
0791     void project(ArrayND<Num2, Len2, Dim2>* projection,
0792                  AbsVisitor<Numeric, Num3>& projector,
0793                  const unsigned* projectedIndices,
0794                  unsigned nProjectedIndices) const;
0795     //@}
0796 
0797     //@{
0798     /**
0799         // Similar method to "project", but projections are added to
0800         // (or subtracted from) the existing projection data instead of
0801         // replacing them
0802         */
0803     template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
0804     void addToProjection(ArrayND<Num2, Len2, Dim2>* projection,
0805                          AbsArrayProjector<Numeric, Num3>& projector,
0806                          const unsigned* projectedIndices,
0807                          unsigned nProjectedIndices) const;
0808 
0809     template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
0810     void subtractFromProjection(ArrayND<Num2, Len2, Dim2>* projection,
0811                                 AbsArrayProjector<Numeric, Num3>& projector,
0812                                 const unsigned* projectedIndices,
0813                                 unsigned nProjectedIndices) const;
0814 
0815     template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
0816     void addToProjection(ArrayND<Num2, Len2, Dim2>* projection,
0817                          AbsVisitor<Numeric, Num3>& projector,
0818                          const unsigned* projectedIndices,
0819                          unsigned nProjectedIndices) const;
0820 
0821     template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
0822     void subtractFromProjection(ArrayND<Num2, Len2, Dim2>* projection,
0823                                 AbsVisitor<Numeric, Num3>& projector,
0824                                 const unsigned* projectedIndices,
0825                                 unsigned nProjectedIndices) const;
0826     //@}
0827 
0828     /**
0829         // Rotation. Place the result into another array. The elements
0830         // with indices 0 in the current array will become elements with
0831         // indices "shifts" in the rotated array.
0832         */
0833     template <typename Num2, unsigned Len2, unsigned Dim2>
0834     void rotate(const unsigned* shifts, unsigned lenShifts, ArrayND<Num2, Len2, Dim2>* rotated) const;
0835 
0836     /**
0837         // Fill another array with all possible mirror images
0838         // of this one. This other array must have twice the span
0839         // in each dimension.
0840         */
0841     template <typename Num2, unsigned Len2, unsigned Dim2>
0842     void multiMirror(ArrayND<Num2, Len2, Dim2>* out) const;
0843 
0844     //@{
0845     /**
0846         // Fortran-style subscripting without bounds checking (of course,
0847         // with indices starting at 0).
0848         */
0849     Numeric& operator()();
0850     const Numeric& operator()() const;
0851 
0852     Numeric& operator()(unsigned i0);
0853     const Numeric& operator()(unsigned i0) const;
0854 
0855     Numeric& operator()(unsigned i0, unsigned i1);
0856     const Numeric& operator()(unsigned i0, unsigned i1) const;
0857 
0858     Numeric& operator()(unsigned i0, unsigned i1, unsigned i2);
0859     const Numeric& operator()(unsigned i0, unsigned i1, unsigned i2) const;
0860 
0861     Numeric& operator()(unsigned i0, unsigned i1, unsigned i2, unsigned i3);
0862     const Numeric& operator()(unsigned i0, unsigned i1, unsigned i2, unsigned i3) const;
0863 
0864     Numeric& operator()(unsigned i0, unsigned i1, unsigned i2, unsigned i3, unsigned i4);
0865     const Numeric& operator()(unsigned i0, unsigned i1, unsigned i2, unsigned i3, unsigned i4) const;
0866 
0867     Numeric& operator()(unsigned i0, unsigned i1, unsigned i2, unsigned i3, unsigned i4, unsigned i5);
0868     const Numeric& operator()(unsigned i0, unsigned i1, unsigned i2, unsigned i3, unsigned i4, unsigned i5) const;
0869 
0870     Numeric& operator()(unsigned i0, unsigned i1, unsigned i2, unsigned i3, unsigned i4, unsigned i5, unsigned i6);
0871     const Numeric& operator()(
0872         unsigned i0, unsigned i1, unsigned i2, unsigned i3, unsigned i4, unsigned i5, unsigned i6) const;
0873 
0874     Numeric& operator()(
0875         unsigned i0, unsigned i1, unsigned i2, unsigned i3, unsigned i4, unsigned i5, unsigned i6, unsigned i7);
0876     const Numeric& operator()(
0877         unsigned i0, unsigned i1, unsigned i2, unsigned i3, unsigned i4, unsigned i5, unsigned i6, unsigned i7) const;
0878 
0879     Numeric& operator()(unsigned i0,
0880                         unsigned i1,
0881                         unsigned i2,
0882                         unsigned i3,
0883                         unsigned i4,
0884                         unsigned i5,
0885                         unsigned i6,
0886                         unsigned i7,
0887                         unsigned i8);
0888     const Numeric& operator()(unsigned i0,
0889                               unsigned i1,
0890                               unsigned i2,
0891                               unsigned i3,
0892                               unsigned i4,
0893                               unsigned i5,
0894                               unsigned i6,
0895                               unsigned i7,
0896                               unsigned i8) const;
0897 
0898     Numeric& operator()(unsigned i0,
0899                         unsigned i1,
0900                         unsigned i2,
0901                         unsigned i3,
0902                         unsigned i4,
0903                         unsigned i5,
0904                         unsigned i6,
0905                         unsigned i7,
0906                         unsigned i8,
0907                         unsigned i9);
0908     const Numeric& operator()(unsigned i0,
0909                               unsigned i1,
0910                               unsigned i2,
0911                               unsigned i3,
0912                               unsigned i4,
0913                               unsigned i5,
0914                               unsigned i6,
0915                               unsigned i7,
0916                               unsigned i8,
0917                               unsigned i9) const;
0918     //@}
0919 
0920     //@{
0921     /**
0922         // Fortran-style subscripting with bounds checking (of course,
0923         // with indices starting at 0).
0924         */
0925     Numeric& at();
0926     const Numeric& at() const;
0927 
0928     Numeric& at(unsigned i0);
0929     const Numeric& at(unsigned i0) const;
0930 
0931     Numeric& at(unsigned i0, unsigned i1);
0932     const Numeric& at(unsigned i0, unsigned i1) const;
0933 
0934     Numeric& at(unsigned i0, unsigned i1, unsigned i2);
0935     const Numeric& at(unsigned i0, unsigned i1, unsigned i2) const;
0936 
0937     Numeric& at(unsigned i0, unsigned i1, unsigned i2, unsigned i3);
0938     const Numeric& at(unsigned i0, unsigned i1, unsigned i2, unsigned i3) const;
0939 
0940     Numeric& at(unsigned i0, unsigned i1, unsigned i2, unsigned i3, unsigned i4);
0941     const Numeric& at(unsigned i0, unsigned i1, unsigned i2, unsigned i3, unsigned i4) const;
0942 
0943     Numeric& at(unsigned i0, unsigned i1, unsigned i2, unsigned i3, unsigned i4, unsigned i5);
0944     const Numeric& at(unsigned i0, unsigned i1, unsigned i2, unsigned i3, unsigned i4, unsigned i5) const;
0945 
0946     Numeric& at(unsigned i0, unsigned i1, unsigned i2, unsigned i3, unsigned i4, unsigned i5, unsigned i6);
0947     const Numeric& at(unsigned i0, unsigned i1, unsigned i2, unsigned i3, unsigned i4, unsigned i5, unsigned i6) const;
0948 
0949     Numeric& at(unsigned i0, unsigned i1, unsigned i2, unsigned i3, unsigned i4, unsigned i5, unsigned i6, unsigned i7);
0950     const Numeric& at(
0951         unsigned i0, unsigned i1, unsigned i2, unsigned i3, unsigned i4, unsigned i5, unsigned i6, unsigned i7) const;
0952 
0953     Numeric& at(unsigned i0,
0954                 unsigned i1,
0955                 unsigned i2,
0956                 unsigned i3,
0957                 unsigned i4,
0958                 unsigned i5,
0959                 unsigned i6,
0960                 unsigned i7,
0961                 unsigned i8);
0962     const Numeric& at(unsigned i0,
0963                       unsigned i1,
0964                       unsigned i2,
0965                       unsigned i3,
0966                       unsigned i4,
0967                       unsigned i5,
0968                       unsigned i6,
0969                       unsigned i7,
0970                       unsigned i8) const;
0971 
0972     Numeric& at(unsigned i0,
0973                 unsigned i1,
0974                 unsigned i2,
0975                 unsigned i3,
0976                 unsigned i4,
0977                 unsigned i5,
0978                 unsigned i6,
0979                 unsigned i7,
0980                 unsigned i8,
0981                 unsigned i9);
0982     const Numeric& at(unsigned i0,
0983                       unsigned i1,
0984                       unsigned i2,
0985                       unsigned i3,
0986                       unsigned i4,
0987                       unsigned i5,
0988                       unsigned i6,
0989                       unsigned i7,
0990                       unsigned i8,
0991                       unsigned i9) const;
0992     //@}
0993 
0994     //@{
0995     /**
0996         // Subscripting by continuous coordinate.
0997         // Works similar to the "closest" method.
0998         */
0999     Numeric& cl();
1000     const Numeric& cl() const;
1001 
1002     Numeric& cl(double x0);
1003     const Numeric& cl(double x0) const;
1004 
1005     Numeric& cl(double x0, double x1);
1006     const Numeric& cl(double x0, double x1) const;
1007 
1008     Numeric& cl(double x0, double x1, double x2);
1009     const Numeric& cl(double x0, double x1, double x2) const;
1010 
1011     Numeric& cl(double x0, double x1, double x2, double x3);
1012     const Numeric& cl(double x0, double x1, double x2, double x3) const;
1013 
1014     Numeric& cl(double x0, double x1, double x2, double x3, double x4);
1015     const Numeric& cl(double x0, double x1, double x2, double x3, double x4) const;
1016 
1017     Numeric& cl(double x0, double x1, double x2, double x3, double x4, double x5);
1018     const Numeric& cl(double x0, double x1, double x2, double x3, double x4, double x5) const;
1019 
1020     Numeric& cl(double x0, double x1, double x2, double x3, double x4, double x5, double x6);
1021     const Numeric& cl(double x0, double x1, double x2, double x3, double x4, double x5, double x6) const;
1022 
1023     Numeric& cl(double x0, double x1, double x2, double x3, double x4, double x5, double x6, double x7);
1024     const Numeric& cl(double x0, double x1, double x2, double x3, double x4, double x5, double x6, double x7) const;
1025 
1026     Numeric& cl(double x0, double x1, double x2, double x3, double x4, double x5, double x6, double x7, double x8);
1027     const Numeric& cl(
1028         double x0, double x1, double x2, double x3, double x4, double x5, double x6, double x7, double x8) const;
1029 
1030     Numeric& cl(
1031         double x0, double x1, double x2, double x3, double x4, double x5, double x6, double x7, double x8, double x9);
1032     const Numeric& cl(
1033         double x0, double x1, double x2, double x3, double x4, double x5, double x6, double x7, double x8, double x9)
1034         const;
1035     //@}
1036 
1037     //@{
1038     /** Methods related to "geners" I/O */
1039     inline gs::ClassId classId() const { return gs::ClassId(*this); }
1040     bool write(std::ostream& of) const;
1041     //@}
1042 
1043     static const char* classname();
1044     static inline unsigned version() { return 1; }
1045     static void restore(const gs::ClassId& id, std::istream& in, ArrayND* array);
1046 
1047   private:
1048     Numeric localData_[StackLen];
1049     Numeric* data_;
1050 
1051     unsigned long localStrides_[StackDim];
1052     unsigned long* strides_;
1053 
1054     unsigned localShape_[StackDim];
1055     unsigned* shape_;
1056 
1057     unsigned long len_;
1058     unsigned dim_;
1059 
1060     bool shapeIsKnown_;
1061 
1062     // Basic initialization from unsigned* shape and dimensionality
1063     void buildFromShapePtr(const unsigned*, unsigned);
1064 
1065     // Build strides_ array out of the shape_ array
1066     void buildStrides();
1067 
1068     // Recursive implementation of nested loops for "linearFill"
1069     void linearFillLoop(unsigned level, double s0, unsigned long idx, double shift, const double* coeffs);
1070 
1071     // Recursive implementation of nested loops for "functorFill"
1072     template <class Functor>
1073     void functorFillLoop(unsigned level, unsigned long idx, Functor f, unsigned* farg);
1074 
1075     // Recursive implementation of nested loops for "interpolate3"
1076     Numeric interpolateLoop(unsigned level, const double* x, const Numeric* base) const;
1077 
1078     // Recursive implementation of nested loops for the outer product
1079     template <typename Num1, unsigned Len1, unsigned Dim1, typename Num2, unsigned Len2, unsigned Dim2>
1080     void outerProductLoop(unsigned level,
1081                           unsigned long idx0,
1082                           unsigned long idx1,
1083                           unsigned long idx2,
1084                           const ArrayND<Num1, Len1, Dim1>& a1,
1085                           const ArrayND<Num2, Len2, Dim2>& a2);
1086 
1087     // Recursive implementation of nested loops for contraction
1088     void contractLoop(unsigned thisLevel,
1089                       unsigned resLevel,
1090                       unsigned pos1,
1091                       unsigned pos2,
1092                       unsigned long idxThis,
1093                       unsigned long idxRes,
1094                       ArrayND& result) const;
1095 
1096     // Recursive implementation of nested loops for transposition
1097     void transposeLoop(unsigned level,
1098                        unsigned pos1,
1099                        unsigned pos2,
1100                        unsigned long idxThis,
1101                        unsigned long idxRes,
1102                        ArrayND& result) const;
1103 
1104     // Recursive implementation of nested loops for the dot product
1105     template <typename Num2, unsigned Len2, unsigned Dim2>
1106     void dotProductLoop(unsigned level,
1107                         unsigned long idx0,
1108                         unsigned long idx1,
1109                         unsigned long idx2,
1110                         const ArrayND<Num2, Len2, Dim2>& r,
1111                         ArrayND& result) const;
1112 
1113     // Recursive implementation of nested loops for marginalization
1114     template <typename Num2, unsigned Len2, unsigned Dim2>
1115     Numeric marginalizeInnerLoop(unsigned long idx,
1116                                  unsigned levelPr,
1117                                  unsigned long idxPr,
1118                                  const ArrayND<Num2, Len2, Dim2>& prior,
1119                                  const unsigned* indexMap) const;
1120     template <typename Num2, unsigned Len2, unsigned Dim2>
1121     void marginalizeLoop(unsigned level,
1122                          unsigned long idx,
1123                          unsigned levelRes,
1124                          unsigned long idxRes,
1125                          const ArrayND<Num2, Len2, Dim2>& prior,
1126                          const unsigned* indexMap,
1127                          ArrayND& res) const;
1128 
1129     // Recursive implementation of nested loops for range copy
1130     // with functor modification of elements
1131     template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
1132     void copyRangeLoopFunct(unsigned level,
1133                             unsigned long idx0,
1134                             unsigned long idx1,
1135                             const ArrayND<Num2, Len2, Dim2>& r,
1136                             const ArrayRange& range,
1137                             Functor f);
1138 
1139     // Loop over subrange in such a way that the functor is called
1140     // only if indices on both sides are valid. The topology of both
1141     // arrays is that of the hyperplane (flat).
1142     template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
1143     void commonSubrangeLoop(unsigned level,
1144                             unsigned long idx0,
1145                             unsigned long idx1,
1146                             const unsigned* thisCorner,
1147                             const unsigned* range,
1148                             const unsigned* otherCorner,
1149                             ArrayND<Num2, Len2, Dim2>& other,
1150                             Functor binaryFunct);
1151 
1152     // Similar loop with the topology of the hypertorus for both
1153     // arrays (all indices of both arrays are wrapped around)
1154     template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
1155     void dualCircularLoop(unsigned level,
1156                           unsigned long idx0,
1157                           unsigned long idx1,
1158                           const unsigned* thisCorner,
1159                           const unsigned* range,
1160                           const unsigned* otherCorner,
1161                           ArrayND<Num2, Len2, Dim2>& other,
1162                           Functor binaryFunct);
1163 
1164     // Similar loop in which the topology of this array is assumed
1165     // to be flat and the topology of the destination array is that
1166     // of the hypertorus. Due to the asymmetry of the topologies,
1167     // "const" function is not provided (use const_cast as appropriate).
1168     template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
1169     void flatCircularLoop(unsigned level,
1170                           unsigned long idx0,
1171                           unsigned long idx1,
1172                           const unsigned* thisCorner,
1173                           const unsigned* range,
1174                           const unsigned* otherCorner,
1175                           ArrayND<Num2, Len2, Dim2>& other,
1176                           Functor binaryFunct);
1177 
1178     // Similar loop in which the topology of this array is assumed
1179     // to be hypertoroidal and the topology of the destination array
1180     // is flat.
1181     template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
1182     void circularFlatLoop(unsigned level,
1183                           unsigned long idx0,
1184                           unsigned long idx1,
1185                           const unsigned* thisCorner,
1186                           const unsigned* range,
1187                           const unsigned* otherCorner,
1188                           ArrayND<Num2, Len2, Dim2>& other,
1189                           Functor binaryFunct);
1190 
1191     // Slice compatibility verification
1192     template <typename Num2, unsigned Len2, unsigned Dim2>
1193     unsigned long verifySliceCompatibility(const ArrayND<Num2, Len2, Dim2>& slice,
1194                                            const unsigned* fixedIndices,
1195                                            const unsigned* fixedIndexValues,
1196                                            unsigned nFixedIndices) const;
1197 
1198     // Buffer compatibility verification with a slice
1199     unsigned long verifyBufferSliceCompatibility(unsigned long bufLen,
1200                                                  const unsigned* fixedIndices,
1201                                                  const unsigned* fixedIndexValues,
1202                                                  unsigned nFixedIndices,
1203                                                  unsigned long* sliceStrides) const;
1204 
1205     // Recursive implementation of nested loops for slice operations
1206     template <typename Num2, class Functor>
1207     void jointSliceLoop(unsigned level,
1208                         unsigned long idx0,
1209                         unsigned level1,
1210                         unsigned long idx1,
1211                         Num2* sliceData,
1212                         const unsigned long* sliceStrides,
1213                         const unsigned* fixedIndices,
1214                         const unsigned* fixedIndexValues,
1215                         unsigned nFixedIndices,
1216                         Functor binaryFunctor);
1217 
1218     // Recursive implementation of nested loops for "applySlice"
1219     template <typename Num2, class Functor>
1220     void scaleBySliceInnerLoop(unsigned level,
1221                                unsigned long idx0,
1222                                Num2& scale,
1223                                const unsigned* projectedIndices,
1224                                unsigned nProjectedIndices,
1225                                Functor binaryFunct);
1226 
1227     template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
1228     void scaleBySliceLoop(unsigned level,
1229                           unsigned long idx0,
1230                           unsigned level1,
1231                           unsigned long idx1,
1232                           ArrayND<Num2, Len2, Dim2>& slice,
1233                           const unsigned* fixedIndices,
1234                           unsigned nFixedIndices,
1235                           Functor binaryFunct);
1236 
1237     // Recursive implementation of nested loops for projections
1238     template <typename Num2>
1239     void projectInnerLoop(unsigned level,
1240                           unsigned long idx0,
1241                           unsigned* currentIndex,
1242                           AbsArrayProjector<Numeric, Num2>& projector,
1243                           const unsigned* projectedIndices,
1244                           unsigned nProjectedIndices) const;
1245 
1246     template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3, class Op>
1247     void projectLoop(unsigned level,
1248                      unsigned long idx0,
1249                      unsigned level1,
1250                      unsigned long idx1,
1251                      unsigned* currentIndex,
1252                      ArrayND<Num2, Len2, Dim2>* projection,
1253                      AbsArrayProjector<Numeric, Num3>& projector,
1254                      const unsigned* projectedIndices,
1255                      unsigned nProjectedIndices,
1256                      Op fcn) const;
1257 
1258     // Note that "projectLoop2" is almost identical to "projectLoop"
1259     // while "projectInnerLoop2" is almost identical to "projectInnerLoop".
1260     // It would make a lot of sense to combine these functions into
1261     // the same code and then partially specialize the little piece
1262     // where the "AbsVisitor" or "AbsArrayProjector" is actually called.
1263     // Unfortunately, "AbsVisitor" and "AbsArrayProjector" are
1264     // templates themselves, and it is not possible in C++ to partially
1265     // specialize a function template (that is, even if we can specialize
1266     // on "AbsVisitor" vs. "AbsArrayProjector", there is no way to
1267     // specialize on their parameter types).
1268     template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3, class Op>
1269     void projectLoop2(unsigned level,
1270                       unsigned long idx0,
1271                       unsigned level1,
1272                       unsigned long idx1,
1273                       ArrayND<Num2, Len2, Dim2>* projection,
1274                       AbsVisitor<Numeric, Num3>& projector,
1275                       const unsigned* projectedIndices,
1276                       unsigned nProjectedIndices,
1277                       Op fcn) const;
1278 
1279     template <typename Num2>
1280     void projectInnerLoop2(unsigned level,
1281                            unsigned long idx0,
1282                            AbsVisitor<Numeric, Num2>& projector,
1283                            const unsigned* projectedIndices,
1284                            unsigned nProjectedIndices) const;
1285 
1286     template <typename Num2, typename Integer>
1287     void processSubrangeLoop(unsigned level,
1288                              unsigned long idx0,
1289                              unsigned* currentIndex,
1290                              AbsArrayProjector<Numeric, Num2>& f,
1291                              const BoxND<Integer>& subrange) const;
1292 
1293     // Sum of all points with the given index and below
1294     template <typename Accumulator>
1295     Accumulator sumBelowLoop(unsigned level, unsigned long idx0, const unsigned* limit) const;
1296 
1297     // Loop for "convertToLastDimCdf"
1298     template <typename Accumulator>
1299     void convertToLastDimCdfLoop(
1300         ArrayND* sumSlice, unsigned level, unsigned long idx0, unsigned long idxSlice, bool useTrapezoids);
1301 
1302     // Convert a coordinate into index.
1303     // No checking whether "idim" is within limits.
1304     unsigned coordToIndex(double coord, unsigned idim) const;
1305 
1306     // Verify that projection array is compatible with this one
1307     template <typename Num2, unsigned Len2, unsigned Dim2>
1308     void verifyProjectionCompatibility(const ArrayND<Num2, Len2, Dim2>& projection,
1309                                        const unsigned* projectedIndices,
1310                                        unsigned nProjectedIndices) const;
1311   };
1312 }  // namespace npstat
1313 
1314 #include <cmath>
1315 #include <climits>
1316 #include <algorithm>
1317 #include <sstream>
1318 #include "JetMETCorrections/InterpolationTables/interface/NpstatException.h"
1319 
1320 #include "Alignment/Geners/interface/GenericIO.hh"
1321 #include "Alignment/Geners/interface/IOIsUnsigned.hh"
1322 
1323 #include "JetMETCorrections/InterpolationTables/interface/allocators.h"
1324 
1325 #include "JetMETCorrections/InterpolationTables/interface/interpolate.h"
1326 #include "JetMETCorrections/InterpolationTables/interface/absDifference.h"
1327 #include "JetMETCorrections/InterpolationTables/interface/ComplexComparesFalse.h"
1328 #include "JetMETCorrections/InterpolationTables/interface/ComplexComparesAbs.h"
1329 
1330 #define me_macro_check_loop_prerequisites(METHOD, INNERLOOP) /**/                                               \
1331   template <typename Numeric, unsigned Len, unsigned Dim>                                                       \
1332   template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>                                         \
1333   void ArrayND<Numeric, Len, Dim>::METHOD(ArrayND<Num2, Len2, Dim2>& other,                                     \
1334                                           const unsigned* thisCorner,                                           \
1335                                           const unsigned* range,                                                \
1336                                           const unsigned* otherCorner,                                          \
1337                                           const unsigned arrLen,                                                \
1338                                           Functor binaryFunct) {                                                \
1339     if (!shapeIsKnown_)                                                                                         \
1340       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"" #METHOD "\"");  \
1341     if (!other.shapeIsKnown_)                                                                                   \
1342       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::" #METHOD ": uninitialized argument array");     \
1343     if (dim_ != other.dim_)                                                                                     \
1344       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::" #METHOD ": incompatible argument array rank"); \
1345     if (arrLen != dim_)                                                                                         \
1346       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::" #METHOD ": incompatible index length");        \
1347     if (dim_) {                                                                                                 \
1348       assert(thisCorner);                                                                                       \
1349       assert(range);                                                                                            \
1350       assert(otherCorner);                                                                                      \
1351       INNERLOOP(0U, 0UL, 0UL, thisCorner, range, otherCorner, other, binaryFunct);                              \
1352     } else                                                                                                      \
1353       binaryFunct(localData_[0], other.localData_[0]);                                                          \
1354   }
1355 
1356 namespace npstat {
1357   template <typename Numeric, unsigned Len, unsigned Dim>
1358   template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
1359   void ArrayND<Numeric, Len, Dim>::commonSubrangeLoop(unsigned level,
1360                                                       unsigned long idx0,
1361                                                       unsigned long idx1,
1362                                                       const unsigned* thisCorner,
1363                                                       const unsigned* range,
1364                                                       const unsigned* otherCorner,
1365                                                       ArrayND<Num2, Len2, Dim2>& r,
1366                                                       Functor binaryFunct) {
1367     const unsigned imax = range[level];
1368 
1369     if (level == dim_ - 1) {
1370       Numeric* left = data_ + (idx0 + thisCorner[level]);
1371       Numeric* const lMax = data_ + (idx0 + shape_[level]);
1372       Num2* right = r.data_ + (idx1 + otherCorner[level]);
1373       Num2* const rMax = r.data_ + (idx1 + r.shape_[level]);
1374 
1375       for (unsigned i = 0; i < imax && left < lMax && right < rMax; ++i)
1376         binaryFunct(*left++, *right++);
1377     } else {
1378       const unsigned long leftStride = strides_[level];
1379       const unsigned long leftMax = idx0 + shape_[level] * leftStride;
1380       idx0 += thisCorner[level] * leftStride;
1381       const unsigned long rightStride = r.strides_[level];
1382       const unsigned long rightMax = idx1 + r.shape_[level] * rightStride;
1383       idx1 += otherCorner[level] * rightStride;
1384 
1385       for (unsigned i = 0; i < imax && idx0 < leftMax && idx1 < rightMax; ++i, idx0 += leftStride, idx1 += rightStride)
1386         commonSubrangeLoop(level + 1, idx0, idx1, thisCorner, range, otherCorner, r, binaryFunct);
1387     }
1388   }
1389 
1390   template <typename Numeric, unsigned Len, unsigned Dim>
1391   template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
1392   void ArrayND<Numeric, Len, Dim>::dualCircularLoop(unsigned level,
1393                                                     unsigned long idx0,
1394                                                     unsigned long idx1,
1395                                                     const unsigned* thisCorner,
1396                                                     const unsigned* range,
1397                                                     const unsigned* otherCorner,
1398                                                     ArrayND<Num2, Len2, Dim2>& r,
1399                                                     Functor binaryFunct) {
1400     const unsigned imax = range[level];
1401     const unsigned leftShift = thisCorner[level];
1402     const unsigned leftPeriod = shape_[level];
1403     const unsigned rightShift = otherCorner[level];
1404     const unsigned rightPeriod = r.shape_[level];
1405 
1406     if (level == dim_ - 1) {
1407       Numeric* left = data_ + idx0;
1408       Num2* right = r.data_ + idx1;
1409       for (unsigned i = 0; i < imax; ++i)
1410         binaryFunct(left[(i + leftShift) % leftPeriod], right[(i + rightShift) % rightPeriod]);
1411     } else {
1412       const unsigned long leftStride = strides_[level];
1413       const unsigned long rightStride = r.strides_[level];
1414       for (unsigned i = 0; i < imax; ++i)
1415         dualCircularLoop(level + 1,
1416                          idx0 + ((i + leftShift) % leftPeriod) * leftStride,
1417                          idx1 + ((i + rightShift) % rightPeriod) * rightStride,
1418                          thisCorner,
1419                          range,
1420                          otherCorner,
1421                          r,
1422                          binaryFunct);
1423     }
1424   }
1425 
1426   template <typename Numeric, unsigned Len, unsigned Dim>
1427   template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
1428   void ArrayND<Numeric, Len, Dim>::flatCircularLoop(unsigned level,
1429                                                     unsigned long idx0,
1430                                                     unsigned long idx1,
1431                                                     const unsigned* thisCorner,
1432                                                     const unsigned* range,
1433                                                     const unsigned* otherCorner,
1434                                                     ArrayND<Num2, Len2, Dim2>& r,
1435                                                     Functor binaryFunct) {
1436     const unsigned imax = range[level];
1437     const unsigned rightShift = otherCorner[level];
1438     const unsigned rightPeriod = r.shape_[level];
1439 
1440     if (level == dim_ - 1) {
1441       Numeric* left = data_ + (idx0 + thisCorner[level]);
1442       Numeric* const lMax = data_ + (idx0 + shape_[level]);
1443       Num2* right = r.data_ + idx1;
1444 
1445       for (unsigned i = 0; i < imax && left < lMax; ++i)
1446         binaryFunct(*left++, right[(i + rightShift) % rightPeriod]);
1447     } else {
1448       const unsigned long leftStride = strides_[level];
1449       const unsigned long leftMax = idx0 + shape_[level] * leftStride;
1450       idx0 += thisCorner[level] * leftStride;
1451       const unsigned long rightStride = r.strides_[level];
1452 
1453       for (unsigned i = 0; i < imax && idx0 < leftMax; ++i, idx0 += leftStride)
1454         flatCircularLoop(level + 1,
1455                          idx0,
1456                          idx1 + ((i + rightShift) % rightPeriod) * rightStride,
1457                          thisCorner,
1458                          range,
1459                          otherCorner,
1460                          r,
1461                          binaryFunct);
1462     }
1463   }
1464 
1465   template <typename Numeric, unsigned Len, unsigned Dim>
1466   template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
1467   void ArrayND<Numeric, Len, Dim>::circularFlatLoop(unsigned level,
1468                                                     unsigned long idx0,
1469                                                     unsigned long idx1,
1470                                                     const unsigned* thisCorner,
1471                                                     const unsigned* range,
1472                                                     const unsigned* otherCorner,
1473                                                     ArrayND<Num2, Len2, Dim2>& r,
1474                                                     Functor binaryFunct) {
1475     const unsigned imax = range[level];
1476     const unsigned leftShift = thisCorner[level];
1477     const unsigned leftPeriod = shape_[level];
1478 
1479     if (level == dim_ - 1) {
1480       Numeric* left = data_ + idx0;
1481       Num2* right = r.data_ + (idx1 + otherCorner[level]);
1482       Num2* const rMax = r.data_ + (idx1 + r.shape_[level]);
1483 
1484       for (unsigned i = 0; i < imax && right < rMax; ++i)
1485         binaryFunct(left[(i + leftShift) % leftPeriod], *right++);
1486     } else {
1487       const unsigned long leftStride = strides_[level];
1488       const unsigned long rightStride = r.strides_[level];
1489       const unsigned long rightMax = idx1 + r.shape_[level] * rightStride;
1490       idx1 += otherCorner[level] * rightStride;
1491 
1492       for (unsigned i = 0; i < imax && idx1 < rightMax; ++i, idx1 += rightStride)
1493         circularFlatLoop(level + 1,
1494                          idx0 + ((i + leftShift) % leftPeriod) * leftStride,
1495                          idx1,
1496                          thisCorner,
1497                          range,
1498                          otherCorner,
1499                          r,
1500                          binaryFunct);
1501     }
1502   }
1503 
1504   me_macro_check_loop_prerequisites(jointSubrangeScan, commonSubrangeLoop)
1505       me_macro_check_loop_prerequisites(dualCircularScan, dualCircularLoop)
1506           me_macro_check_loop_prerequisites(flatCircularScan,
1507                                             flatCircularLoop) me_macro_check_loop_prerequisites(circularFlatScan,
1508                                                                                                 circularFlatLoop)
1509 
1510               template <typename Numeric, unsigned StackLen, unsigned StackDim>
1511               template <typename Num2, unsigned Len2, unsigned Dim2>
1512               Numeric ArrayND<Numeric, StackLen, StackDim>::marginalizeInnerLoop(unsigned long idx,
1513                                                                                  const unsigned levelPr,
1514                                                                                  unsigned long idxPr,
1515                                                                                  const ArrayND<Num2, Len2, Dim2>& prior,
1516                                                                                  const unsigned* indexMap) const {
1517     Numeric sum = Numeric();
1518     const unsigned long myStride = strides_[indexMap[levelPr]];
1519     const unsigned imax = prior.shape_[levelPr];
1520     assert(imax == shape_[indexMap[levelPr]]);
1521     if (levelPr == prior.dim_ - 1) {
1522       for (unsigned i = 0; i < imax; ++i)
1523         sum += data_[idx + i * myStride] * prior.data_[idxPr++];
1524     } else {
1525       const unsigned long priorStride = prior.strides_[levelPr];
1526       for (unsigned i = 0; i < imax; ++i) {
1527         sum += marginalizeInnerLoop(idx, levelPr + 1U, idxPr, prior, indexMap);
1528         idx += myStride;
1529         idxPr += priorStride;
1530       }
1531     }
1532     return sum;
1533   }
1534 
1535   template <typename Numeric, unsigned StackLen, unsigned StackDim>
1536   template <typename Num2, unsigned Len2, unsigned Dim2>
1537   void ArrayND<Numeric, StackLen, StackDim>::marginalizeLoop(const unsigned level,
1538                                                              unsigned long idx,
1539                                                              const unsigned levelRes,
1540                                                              unsigned long idxRes,
1541                                                              const ArrayND<Num2, Len2, Dim2>& prior,
1542                                                              const unsigned* indexMap,
1543                                                              ArrayND& result) const {
1544     if (level == dim_) {
1545       const Numeric res = marginalizeInnerLoop(idx, 0U, 0UL, prior, indexMap);
1546       if (result.dim_)
1547         result.data_[idxRes] = res;
1548       else
1549         result.localData_[0] = res;
1550     } else {
1551       // Check if this level is mapped or not
1552       bool mapped = false;
1553       for (unsigned i = 0; i < prior.dim_; ++i)
1554         if (level == indexMap[i]) {
1555           mapped = true;
1556           break;
1557         }
1558       if (mapped)
1559         marginalizeLoop(level + 1U, idx, levelRes, idxRes, prior, indexMap, result);
1560       else {
1561         const unsigned imax = shape_[level];
1562         const unsigned long myStride = strides_[level];
1563         const unsigned long resStride = result.strides_[levelRes];
1564         for (unsigned i = 0; i < imax; ++i) {
1565           marginalizeLoop(level + 1U, idx, levelRes + 1U, idxRes, prior, indexMap, result);
1566           idx += myStride;
1567           idxRes += resStride;
1568         }
1569       }
1570     }
1571   }
1572 
1573   template <typename Numeric, unsigned StackLen, unsigned StackDim>
1574   template <typename Num2, unsigned Len2, unsigned Dim2>
1575   ArrayND<Numeric, StackLen, StackDim> ArrayND<Numeric, StackLen, StackDim>::marginalize(
1576       const ArrayND<Num2, Len2, Dim2>& prior, const unsigned* indexMap, const unsigned mapLen) const {
1577     if (!shapeIsKnown_)
1578       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"marginalize\"");
1579     if (!(prior.dim_ && prior.dim_ <= dim_))
1580       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::marginalize: incompatible argument array rank");
1581     const unsigned resultDim = dim_ - prior.dim_;
1582 
1583     // Check that the index map is reasonable
1584     if (mapLen != prior.dim_)
1585       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::marginalize: incompatible index map length");
1586     assert(indexMap);
1587     for (unsigned i = 0; i < mapLen; ++i) {
1588       const unsigned thisInd = indexMap[i];
1589       if (shape_[thisInd] != prior.shape_[i])
1590         throw npstat::NpstatInvalidArgument(
1591             "In npstat::ArrayND::marginalize: "
1592             "incompatible argument array dimensions");
1593       if (thisInd >= dim_)
1594         throw npstat::NpstatOutOfRange("In npstat::ArrayND::marginalize: index map entry out of range");
1595       for (unsigned j = 0; j < i; ++j)
1596         if (indexMap[j] == thisInd)
1597           throw npstat::NpstatInvalidArgument(
1598               "In npstat::ArrayND::marginalize: "
1599               "duplicate entry in the index map");
1600     }
1601 
1602     // Build the shape for the array of results
1603     ArrayShape newShape;
1604     newShape.reserve(resultDim);
1605     for (unsigned i = 0; i < dim_; ++i) {
1606       bool mapped = false;
1607       for (unsigned j = 0; j < mapLen; ++j)
1608         if (indexMap[j] == i) {
1609           mapped = true;
1610           break;
1611         }
1612       if (!mapped)
1613         newShape.push_back(shape_[i]);
1614     }
1615 
1616     ArrayND result(newShape);
1617     assert(result.dim_ == resultDim);
1618     bool calculated = false;
1619     if (resultDim == 0) {
1620       calculated = true;
1621       for (unsigned i = 0; i < dim_; ++i)
1622         if (indexMap[i] != i) {
1623           calculated = false;
1624           break;
1625         }
1626       if (calculated) {
1627         Numeric sum = Numeric();
1628         for (unsigned long i = 0; i < len_; ++i)
1629           sum += data_[i] * prior.data_[i];
1630         result.localData_[0] = sum;
1631       }
1632     }
1633 
1634     if (!calculated)
1635       marginalizeLoop(0U, 0UL, 0U, 0UL, prior, indexMap, result);
1636 
1637     return result;
1638   }
1639 
1640   template <typename Numeric, unsigned Len, unsigned Dim>
1641   void ArrayND<Numeric, Len, Dim>::buildFromShapePtr(const unsigned* sizes, const unsigned dim) {
1642     dim_ = dim;
1643     if (dim_) {
1644       assert(sizes);
1645       for (unsigned i = 0; i < dim_; ++i)
1646         if (sizes[i] == 0)
1647           throw npstat::NpstatInvalidArgument(
1648               "In npstat::ArrayND::buildFromShapePtr: "
1649               "detected span of zero");
1650 
1651       // Copy the array shape and figure out the array length
1652       shape_ = makeBuffer(dim_, localShape_, Dim);
1653       for (unsigned i = 0; i < dim_; ++i) {
1654         shape_[i] = sizes[i];
1655         len_ *= shape_[i];
1656       }
1657 
1658       // Figure out the array strides
1659       buildStrides();
1660 
1661       // Allocate the data array
1662       data_ = makeBuffer(len_, localData_, Len);
1663     } else {
1664       localData_[0] = Numeric();
1665       data_ = localData_;
1666     }
1667   }
1668 
1669   template <typename Numeric, unsigned StackLen, unsigned StackDim>
1670   template <typename Num2, unsigned Len2, unsigned Dim2>
1671   ArrayND<Numeric, StackLen, StackDim>::ArrayND(const ArrayND<Num2, Len2, Dim2>& slicedArray,
1672                                                 const unsigned* fixedIndices,
1673                                                 const unsigned nFixedIndices)
1674       : data_(0),
1675         strides_(nullptr),
1676         shape_(nullptr),
1677         len_(1UL),
1678         dim_(slicedArray.dim_ - nFixedIndices),
1679         shapeIsKnown_(true) {
1680     if (nFixedIndices) {
1681       assert(fixedIndices);
1682       if (nFixedIndices > slicedArray.dim_)
1683         throw npstat::NpstatInvalidArgument("In npstat::ArrayND slicing constructor: too many fixed indices");
1684       if (!slicedArray.shapeIsKnown_)
1685         throw npstat::NpstatInvalidArgument(
1686             "In npstat::ArrayND slicing constructor: "
1687             "uninitialized argument array");
1688 
1689       // Check that the fixed indices are within range
1690       for (unsigned j = 0; j < nFixedIndices; ++j)
1691         if (fixedIndices[j] >= slicedArray.dim_)
1692           throw npstat::NpstatOutOfRange(
1693               "In npstat::ArrayND slicing "
1694               "constructor: fixed index out of range");
1695 
1696       // Build the shape for the slice
1697       shape_ = makeBuffer(dim_, localShape_, StackDim);
1698       unsigned idim = 0;
1699       for (unsigned i = 0; i < slicedArray.dim_; ++i) {
1700         bool fixed = false;
1701         for (unsigned j = 0; j < nFixedIndices; ++j)
1702           if (fixedIndices[j] == i) {
1703             fixed = true;
1704             break;
1705           }
1706         if (!fixed) {
1707           assert(idim < dim_);
1708           shape_[idim++] = slicedArray.shape_[i];
1709         }
1710       }
1711       assert(idim == dim_);
1712 
1713       if (dim_) {
1714         // Copy the array shape and figure out the array length
1715         for (unsigned i = 0; i < dim_; ++i)
1716           len_ *= shape_[i];
1717 
1718         // Figure out the array strides
1719         buildStrides();
1720 
1721         // Allocate the data array
1722         data_ = makeBuffer(len_, localData_, StackLen);
1723       } else {
1724         localData_[0] = Numeric();
1725         data_ = localData_;
1726       }
1727     } else {
1728       new (this) ArrayND(slicedArray);
1729     }
1730   }
1731 
1732   template <typename Numeric, unsigned StackLen, unsigned StackDim>
1733   ArrayShape ArrayND<Numeric, StackLen, StackDim>::sliceShape(const unsigned* fixedIndices,
1734                                                               const unsigned nFixedIndices) const {
1735     if (!shapeIsKnown_)
1736       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"sliceShape\"");
1737     if (nFixedIndices) {
1738       assert(fixedIndices);
1739       if (nFixedIndices > dim_)
1740         throw npstat::NpstatInvalidArgument("In npstat::ArrayND::sliceShape: too many fixed indices");
1741       for (unsigned j = 0; j < nFixedIndices; ++j)
1742         if (fixedIndices[j] >= dim_)
1743           throw npstat::NpstatOutOfRange(
1744               "In npstat::ArrayND::sliceShape: "
1745               "fixed index out of range");
1746       ArrayShape sh;
1747       sh.reserve(dim_ - nFixedIndices);
1748       for (unsigned i = 0; i < dim_; ++i) {
1749         bool fixed = false;
1750         for (unsigned j = 0; j < nFixedIndices; ++j)
1751           if (fixedIndices[j] == i) {
1752             fixed = true;
1753             break;
1754           }
1755         if (!fixed)
1756           sh.push_back(shape_[i]);
1757       }
1758       return sh;
1759     } else
1760       return ArrayShape(shape_, shape_ + dim_);
1761   }
1762 
1763   template <typename Numeric, unsigned StackLen, unsigned StackDim>
1764   template <typename Num2, unsigned Len2, unsigned Dim2>
1765   unsigned long ArrayND<Numeric, StackLen, StackDim>::verifySliceCompatibility(const ArrayND<Num2, Len2, Dim2>& slice,
1766                                                                                const unsigned* fixedIndices,
1767                                                                                const unsigned* fixedIndexValues,
1768                                                                                const unsigned nFixedIndices) const {
1769     if (!(nFixedIndices && nFixedIndices <= dim_))
1770       throw npstat::NpstatInvalidArgument(
1771           "In npstat::ArrayND::verifySliceCompatibility: "
1772           "invalid number of fixed indices");
1773     if (!shapeIsKnown_)
1774       throw npstat::NpstatInvalidArgument(
1775           "Initialize npstat::ArrayND before calling "
1776           "method \"verifySliceCompatibility\"");
1777     if (!slice.shapeIsKnown_)
1778       throw npstat::NpstatInvalidArgument(
1779           "In npstat::ArrayND::verifySliceCompatibility: "
1780           "uninitialized argument array");
1781     if (slice.dim_ != dim_ - nFixedIndices)
1782       throw npstat::NpstatInvalidArgument(
1783           "In npstat::ArrayND::verifySliceCompatibility: "
1784           "incompatible argument array rank");
1785     assert(fixedIndices);
1786     assert(fixedIndexValues);
1787 
1788     for (unsigned j = 0; j < nFixedIndices; ++j)
1789       if (fixedIndices[j] >= dim_)
1790         throw npstat::NpstatOutOfRange(
1791             "In npstat::ArrayND::verifySliceCompatibility: "
1792             "fixed index out of range");
1793 
1794     // Check slice shape compatibility
1795     unsigned long idx = 0UL;
1796     unsigned sliceDim = 0U;
1797     for (unsigned i = 0; i < dim_; ++i) {
1798       bool fixed = false;
1799       for (unsigned j = 0; j < nFixedIndices; ++j)
1800         if (fixedIndices[j] == i) {
1801           fixed = true;
1802           if (fixedIndexValues[j] >= shape_[i])
1803             throw npstat::NpstatOutOfRange(
1804                 "In npstat::ArrayND::verifySliceCompatibility: "
1805                 "fixed index value out of range");
1806           idx += fixedIndexValues[j] * strides_[i];
1807           break;
1808         }
1809       if (!fixed) {
1810         if (shape_[i] != slice.shape_[sliceDim])
1811           throw npstat::NpstatInvalidArgument(
1812               "In npstat::ArrayND::verifySliceCompatibility: "
1813               "incompatible argument array dimensions");
1814         ++sliceDim;
1815       }
1816     }
1817     assert(sliceDim == slice.dim_);
1818     return idx;
1819   }
1820 
1821   template <typename Numeric, unsigned StackLen, unsigned StackDim>
1822   unsigned long ArrayND<Numeric, StackLen, StackDim>::verifyBufferSliceCompatibility(const unsigned long bufLen,
1823                                                                                      const unsigned* fixedIndices,
1824                                                                                      const unsigned* fixedIndexValues,
1825                                                                                      const unsigned nFixedIndices,
1826                                                                                      unsigned long* sliceStrides) const {
1827     if (!(nFixedIndices && nFixedIndices <= dim_))
1828       throw npstat::NpstatInvalidArgument(
1829           "In npstat::ArrayND::verifyBufferSliceCompatibility: "
1830           "invalid number of fixed indices");
1831     if (!shapeIsKnown_)
1832       throw npstat::NpstatInvalidArgument(
1833           "Initialize npstat::ArrayND before calling "
1834           "method \"verifyBufferSliceCompatibility\"");
1835     assert(fixedIndices);
1836     assert(fixedIndexValues);
1837 
1838     for (unsigned j = 0; j < nFixedIndices; ++j)
1839       if (fixedIndices[j] >= dim_)
1840         throw npstat::NpstatOutOfRange(
1841             "In npstat::ArrayND::verifyBufferSliceCompatibility: "
1842             "fixed index out of range");
1843 
1844     // Figure out slice shape. Store it, temporarily, in sliceStrides.
1845     unsigned long idx = 0UL;
1846     unsigned sliceDim = 0U;
1847     for (unsigned i = 0; i < dim_; ++i) {
1848       bool fixed = false;
1849       for (unsigned j = 0; j < nFixedIndices; ++j)
1850         if (fixedIndices[j] == i) {
1851           fixed = true;
1852           if (fixedIndexValues[j] >= shape_[i])
1853             throw npstat::NpstatOutOfRange(
1854                 "In npstat::ArrayND::verifyBufferSliceCompatibility:"
1855                 " fixed index value out of range");
1856           idx += fixedIndexValues[j] * strides_[i];
1857           break;
1858         }
1859       if (!fixed)
1860         sliceStrides[sliceDim++] = shape_[i];
1861     }
1862     assert(sliceDim + nFixedIndices == dim_);
1863 
1864     // Convert the slice shape into slice strides
1865     unsigned long expectedBufLen = 1UL;
1866     if (sliceDim) {
1867       unsigned long shapeJ = sliceStrides[sliceDim - 1];
1868       sliceStrides[sliceDim - 1] = 1UL;
1869       for (unsigned j = sliceDim - 1; j > 0; --j) {
1870         const unsigned long nextStride = sliceStrides[j] * shapeJ;
1871         shapeJ = sliceStrides[j - 1];
1872         sliceStrides[j - 1] = nextStride;
1873       }
1874       expectedBufLen = sliceStrides[0] * shapeJ;
1875     }
1876     if (expectedBufLen != bufLen)
1877       throw npstat::NpstatInvalidArgument(
1878           "In npstat::ArrayND::verifyBufferSliceCompatibility: "
1879           "invalid memory buffer length");
1880 
1881     return idx;
1882   }
1883 
1884   template <typename Numeric, unsigned StackLen, unsigned StackDim>
1885   template <typename Num2, class Op>
1886   void ArrayND<Numeric, StackLen, StackDim>::jointSliceLoop(const unsigned level,
1887                                                             const unsigned long idx0,
1888                                                             const unsigned level1,
1889                                                             const unsigned long idx1,
1890                                                             Num2* sliceData,
1891                                                             const unsigned long* sliceStrides,
1892                                                             const unsigned* fixedIndices,
1893                                                             const unsigned* fixedIndexValues,
1894                                                             const unsigned nFixedIndices,
1895                                                             Op fcn) {
1896     bool fixed = false;
1897     for (unsigned j = 0; j < nFixedIndices; ++j)
1898       if (fixedIndices[j] == level) {
1899         fixed = true;
1900         break;
1901       }
1902     if (fixed)
1903       jointSliceLoop(
1904           level + 1, idx0, level1, idx1, sliceData, sliceStrides, fixedIndices, fixedIndexValues, nFixedIndices, fcn);
1905     else {
1906       const unsigned imax = shape_[level];
1907       const unsigned long stride = strides_[level];
1908 
1909       if (level1 == dim_ - nFixedIndices - 1) {
1910         sliceData += idx1;
1911         Numeric* localData = data_ + idx0;
1912         for (unsigned i = 0; i < imax; ++i)
1913           fcn(localData[i * stride], sliceData[i]);
1914       } else {
1915         const unsigned long stride2 = sliceStrides[level1];
1916         for (unsigned i = 0; i < imax; ++i)
1917           jointSliceLoop(level + 1,
1918                          idx0 + i * stride,
1919                          level1 + 1,
1920                          idx1 + i * stride2,
1921                          sliceData,
1922                          sliceStrides,
1923                          fixedIndices,
1924                          fixedIndexValues,
1925                          nFixedIndices,
1926                          fcn);
1927       }
1928     }
1929   }
1930 
1931   template <typename Numeric, unsigned StackLen, unsigned StackDim>
1932   template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
1933   void ArrayND<Numeric, StackLen, StackDim>::jointSliceScan(ArrayND<Num2, Len2, Dim2>& slice,
1934                                                             const unsigned* fixedIndices,
1935                                                             const unsigned* fixedIndexValues,
1936                                                             const unsigned nFixedIndices,
1937                                                             Functor binaryFunct) {
1938     const unsigned long idx = verifySliceCompatibility(slice, fixedIndices, fixedIndexValues, nFixedIndices);
1939     if (slice.dim_)
1940       jointSliceLoop(
1941           0U, idx, 0U, 0UL, slice.data_, slice.strides_, fixedIndices, fixedIndexValues, nFixedIndices, binaryFunct);
1942     else
1943       binaryFunct(data_[idx], slice.localData_[0]);
1944   }
1945 
1946   template <typename Numeric, unsigned StackLen, unsigned StackDim>
1947   template <typename Num2, class Functor>
1948   void ArrayND<Numeric, StackLen, StackDim>::jointMemSliceScan(Num2* slice,
1949                                                                const unsigned long len,
1950                                                                const unsigned* fixedIndices,
1951                                                                const unsigned* fixedIndexValues,
1952                                                                unsigned nFixedIndices,
1953                                                                Functor binaryFunct) {
1954     assert(slice);
1955     if (dim_ > CHAR_BIT * sizeof(unsigned long))
1956       throw npstat::NpstatInvalidArgument(
1957           "In npstat::ArrayND::jointMemSliceScan: "
1958           "rank of this array is too large");
1959     unsigned long sliceStrides[CHAR_BIT * sizeof(unsigned long)];
1960     const unsigned long idx =
1961         verifyBufferSliceCompatibility(len, fixedIndices, fixedIndexValues, nFixedIndices, sliceStrides);
1962     if (dim_ > nFixedIndices)
1963       jointSliceLoop(0U, idx, 0U, 0UL, slice, sliceStrides, fixedIndices, fixedIndexValues, nFixedIndices, binaryFunct);
1964     else
1965       binaryFunct(data_[idx], *slice);
1966   }
1967 
1968   template <typename Numeric, unsigned StackLen, unsigned StackDim>
1969   template <typename Num2>
1970   void ArrayND<Numeric, StackLen, StackDim>::projectInnerLoop(const unsigned level,
1971                                                               unsigned long idx0,
1972                                                               unsigned* currentIndex,
1973                                                               AbsArrayProjector<Numeric, Num2>& projector,
1974                                                               const unsigned* projectedIndices,
1975                                                               const unsigned nProjectedIndices) const {
1976     // level :  dimension number among indices which are being iterated
1977     const unsigned idx = projectedIndices[level];
1978     const unsigned imax = shape_[idx];
1979     const unsigned long stride = strides_[idx];
1980     const bool last = (level == nProjectedIndices - 1);
1981 
1982     for (unsigned i = 0; i < imax; ++i) {
1983       currentIndex[idx] = i;
1984       if (last)
1985         projector.process(currentIndex, dim_, idx0, data_[idx0]);
1986       else
1987         projectInnerLoop(level + 1, idx0, currentIndex, projector, projectedIndices, nProjectedIndices);
1988       idx0 += stride;
1989     }
1990   }
1991 
1992   template <typename Numeric, unsigned StackLen, unsigned StackDim>
1993   template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3, class Op>
1994   void ArrayND<Numeric, StackLen, StackDim>::projectLoop(const unsigned level,
1995                                                          const unsigned long idx0,
1996                                                          const unsigned level1,
1997                                                          const unsigned long idx1,
1998                                                          unsigned* currentIndex,
1999                                                          ArrayND<Num2, Len2, Dim2>* projection,
2000                                                          AbsArrayProjector<Numeric, Num3>& projector,
2001                                                          const unsigned* projectedIndices,
2002                                                          const unsigned nProjectedIndices,
2003                                                          Op fcn) const {
2004     // level        : dimension number in this array
2005     // level1       : dimension number in the projection
2006     // idx0         : linear index in this array
2007     // idx1         : linear index in the projection
2008     // currentIndex : cycled over in this loop with the exception of the
2009     //                dimensions which are iterated over to build the
2010     //                projection
2011     if (level == dim_) {
2012       assert(level1 == projection->dim_);
2013       projector.clear();
2014       projectInnerLoop(0U, idx0, currentIndex, projector, projectedIndices, nProjectedIndices);
2015       if (projection->dim_)
2016         fcn(projection->data_[idx1], projector.result());
2017       else
2018         fcn(projection->localData_[0], projector.result());
2019     } else {
2020       bool iterated = false;
2021       for (unsigned j = 0; j < nProjectedIndices; ++j)
2022         if (projectedIndices[j] == level) {
2023           iterated = true;
2024           break;
2025         }
2026       if (iterated) {
2027         // This index will be iterated over inside "projectInnerLoop"
2028         projectLoop(
2029             level + 1, idx0, level1, idx1, currentIndex, projection, projector, projectedIndices, nProjectedIndices, fcn);
2030       } else {
2031         const unsigned imax = shape_[level];
2032         const unsigned long stride = strides_[level];
2033         // We will not be able to get here if projection->dim_ is 0.
2034         // Therefore, it is safe to access projection->strides_.
2035         const unsigned long stride2 = projection->strides_[level1];
2036         for (unsigned i = 0; i < imax; ++i) {
2037           currentIndex[level] = i;
2038           projectLoop(level + 1,
2039                       idx0 + i * stride,
2040                       level1 + 1,
2041                       idx1 + i * stride2,
2042                       currentIndex,
2043                       projection,
2044                       projector,
2045                       projectedIndices,
2046                       nProjectedIndices,
2047                       fcn);
2048         }
2049       }
2050     }
2051   }
2052 
2053   template <typename Numeric, unsigned StackLen, unsigned StackDim>
2054   template <typename Num2, unsigned Len2, unsigned Dim2>
2055   void ArrayND<Numeric, StackLen, StackDim>::verifyProjectionCompatibility(const ArrayND<Num2, Len2, Dim2>& projection,
2056                                                                            const unsigned* projectedIndices,
2057                                                                            const unsigned nProjectedIndices) const {
2058     if (!(nProjectedIndices && nProjectedIndices <= dim_))
2059       throw npstat::NpstatInvalidArgument(
2060           "In npstat::ArrayND::verifyProjectionCompatibility: "
2061           "invalid number of projected indices");
2062     if (!shapeIsKnown_)
2063       throw npstat::NpstatInvalidArgument(
2064           "Initialize npstat::ArrayND before calling "
2065           "method \"verifyProjectionCompatibility\"");
2066     if (!projection.shapeIsKnown_)
2067       throw npstat::NpstatInvalidArgument(
2068           "In npstat::ArrayND::verifyProjectionCompatibility: "
2069           "uninitialized argument array");
2070     if (projection.dim_ != dim_ - nProjectedIndices)
2071       throw npstat::NpstatInvalidArgument(
2072           "In npstat::ArrayND::verifyProjectionCompatibility: "
2073           "incompatible argument array rank");
2074     assert(projectedIndices);
2075 
2076     for (unsigned j = 0; j < nProjectedIndices; ++j)
2077       if (projectedIndices[j] >= dim_)
2078         throw npstat::NpstatOutOfRange(
2079             "In npstat::ArrayND::verifyProjectionCompatibility: "
2080             "projected index out of range");
2081 
2082     // Check projection shape compatibility
2083     unsigned sliceDim = 0U;
2084     for (unsigned i = 0; i < dim_; ++i) {
2085       bool fixed = false;
2086       for (unsigned j = 0; j < nProjectedIndices; ++j)
2087         if (projectedIndices[j] == i) {
2088           fixed = true;
2089           break;
2090         }
2091       if (!fixed) {
2092         if (shape_[i] != projection.shape_[sliceDim])
2093           throw npstat::NpstatInvalidArgument(
2094               "In npstat::ArrayND::verifyProjectionCompatibility: "
2095               "incompatible argument array dimensions");
2096         ++sliceDim;
2097       }
2098     }
2099     assert(sliceDim == projection.dim_);
2100   }
2101 
2102   template <typename Numeric, unsigned StackLen, unsigned StackDim>
2103   template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
2104   void ArrayND<Numeric, StackLen, StackDim>::project(ArrayND<Num2, Len2, Dim2>* projection,
2105                                                      AbsArrayProjector<Numeric, Num3>& projector,
2106                                                      const unsigned* projectedIndices,
2107                                                      const unsigned nProjectedIndices) const {
2108     assert(projection);
2109     verifyProjectionCompatibility(*projection, projectedIndices, nProjectedIndices);
2110     unsigned ibuf[StackDim];
2111     unsigned* buf = makeBuffer(dim_, ibuf, StackDim);
2112     for (unsigned i = 0; i < dim_; ++i)
2113       buf[i] = 0U;
2114     projectLoop(0U,
2115                 0UL,
2116                 0U,
2117                 0UL,
2118                 buf,
2119                 projection,
2120                 projector,
2121                 projectedIndices,
2122                 nProjectedIndices,
2123                 scast_assign_left<Num2, Num3>());
2124     destroyBuffer(buf, ibuf);
2125   }
2126 
2127   template <typename Numeric, unsigned StackLen, unsigned StackDim>
2128   template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
2129   void ArrayND<Numeric, StackLen, StackDim>::addToProjection(ArrayND<Num2, Len2, Dim2>* projection,
2130                                                              AbsArrayProjector<Numeric, Num3>& projector,
2131                                                              const unsigned* projectedIndices,
2132                                                              const unsigned nProjectedIndices) const {
2133     assert(projection);
2134     verifyProjectionCompatibility(*projection, projectedIndices, nProjectedIndices);
2135     unsigned ibuf[StackDim];
2136     unsigned* buf = makeBuffer(dim_, ibuf, StackDim);
2137     for (unsigned i = 0; i < dim_; ++i)
2138       buf[i] = 0U;
2139     projectLoop(0U,
2140                 0UL,
2141                 0U,
2142                 0UL,
2143                 buf,
2144                 projection,
2145                 projector,
2146                 projectedIndices,
2147                 nProjectedIndices,
2148                 scast_pluseq_left<Num2, Num3>());
2149     destroyBuffer(buf, ibuf);
2150   }
2151 
2152   template <typename Numeric, unsigned StackLen, unsigned StackDim>
2153   template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
2154   void ArrayND<Numeric, StackLen, StackDim>::subtractFromProjection(ArrayND<Num2, Len2, Dim2>* projection,
2155                                                                     AbsArrayProjector<Numeric, Num3>& projector,
2156                                                                     const unsigned* projectedIndices,
2157                                                                     const unsigned nProjectedIndices) const {
2158     assert(projection);
2159     verifyProjectionCompatibility(*projection, projectedIndices, nProjectedIndices);
2160     unsigned ibuf[StackDim];
2161     unsigned* buf = makeBuffer(dim_, ibuf, StackDim);
2162     for (unsigned i = 0; i < dim_; ++i)
2163       buf[i] = 0U;
2164     projectLoop(0U,
2165                 0UL,
2166                 0U,
2167                 0UL,
2168                 buf,
2169                 projection,
2170                 projector,
2171                 projectedIndices,
2172                 nProjectedIndices,
2173                 scast_minuseq_left<Num2, Num3>());
2174     destroyBuffer(buf, ibuf);
2175   }
2176 
2177   template <typename Numeric, unsigned StackLen, unsigned StackDim>
2178   template <typename Num2>
2179   void ArrayND<Numeric, StackLen, StackDim>::projectInnerLoop2(const unsigned level,
2180                                                                const unsigned long idx0,
2181                                                                AbsVisitor<Numeric, Num2>& projector,
2182                                                                const unsigned* projectedIndices,
2183                                                                const unsigned nProjectedIndices) const {
2184     const unsigned idx = projectedIndices[level];
2185     const unsigned imax = shape_[idx];
2186     const unsigned long stride = strides_[idx];
2187     const bool last = (level == nProjectedIndices - 1);
2188 
2189     for (unsigned i = 0; i < imax; ++i) {
2190       if (last)
2191         projector.process(data_[idx0 + i * stride]);
2192       else
2193         projectInnerLoop2(level + 1, idx0 + i * stride, projector, projectedIndices, nProjectedIndices);
2194     }
2195   }
2196 
2197   template <typename Numeric, unsigned StackLen, unsigned StackDim>
2198   template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3, class Op>
2199   void ArrayND<Numeric, StackLen, StackDim>::projectLoop2(const unsigned level,
2200                                                           const unsigned long idx0,
2201                                                           const unsigned level1,
2202                                                           const unsigned long idx1,
2203                                                           ArrayND<Num2, Len2, Dim2>* projection,
2204                                                           AbsVisitor<Numeric, Num3>& projector,
2205                                                           const unsigned* projectedIndices,
2206                                                           const unsigned nProjectedIndices,
2207                                                           Op fcn) const {
2208     if (level == dim_) {
2209       assert(level1 == projection->dim_);
2210       projector.clear();
2211       projectInnerLoop2(0U, idx0, projector, projectedIndices, nProjectedIndices);
2212       if (projection->dim_)
2213         fcn(projection->data_[idx1], projector.result());
2214       else
2215         fcn(projection->localData_[0], projector.result());
2216     } else {
2217       bool fixed = false;
2218       for (unsigned j = 0; j < nProjectedIndices; ++j)
2219         if (projectedIndices[j] == level) {
2220           fixed = true;
2221           break;
2222         }
2223       if (fixed) {
2224         projectLoop2(level + 1, idx0, level1, idx1, projection, projector, projectedIndices, nProjectedIndices, fcn);
2225       } else {
2226         const unsigned imax = shape_[level];
2227         const unsigned long stride = strides_[level];
2228         const unsigned long stride2 = projection->strides_[level1];
2229         for (unsigned i = 0; i < imax; ++i)
2230           projectLoop2(level + 1,
2231                        idx0 + i * stride,
2232                        level1 + 1,
2233                        idx1 + i * stride2,
2234                        projection,
2235                        projector,
2236                        projectedIndices,
2237                        nProjectedIndices,
2238                        fcn);
2239       }
2240     }
2241   }
2242 
2243   template <typename Numeric, unsigned StackLen, unsigned StackDim>
2244   template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
2245   void ArrayND<Numeric, StackLen, StackDim>::project(ArrayND<Num2, Len2, Dim2>* projection,
2246                                                      AbsVisitor<Numeric, Num3>& projector,
2247                                                      const unsigned* projectedIndices,
2248                                                      const unsigned nProjectedIndices) const {
2249     assert(projection);
2250     verifyProjectionCompatibility(*projection, projectedIndices, nProjectedIndices);
2251     projectLoop2(
2252         0U, 0UL, 0U, 0UL, projection, projector, projectedIndices, nProjectedIndices, scast_assign_left<Num2, Num3>());
2253   }
2254 
2255   template <typename Numeric, unsigned StackLen, unsigned StackDim>
2256   template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
2257   void ArrayND<Numeric, StackLen, StackDim>::addToProjection(ArrayND<Num2, Len2, Dim2>* projection,
2258                                                              AbsVisitor<Numeric, Num3>& projector,
2259                                                              const unsigned* projectedIndices,
2260                                                              const unsigned nProjectedIndices) const {
2261     assert(projection);
2262     verifyProjectionCompatibility(*projection, projectedIndices, nProjectedIndices);
2263     projectLoop2(
2264         0U, 0UL, 0U, 0UL, projection, projector, projectedIndices, nProjectedIndices, scast_pluseq_left<Num2, Num3>());
2265   }
2266 
2267   template <typename Numeric, unsigned StackLen, unsigned StackDim>
2268   template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
2269   void ArrayND<Numeric, StackLen, StackDim>::subtractFromProjection(ArrayND<Num2, Len2, Dim2>* projection,
2270                                                                     AbsVisitor<Numeric, Num3>& projector,
2271                                                                     const unsigned* projectedIndices,
2272                                                                     const unsigned nProjectedIndices) const {
2273     assert(projection);
2274     verifyProjectionCompatibility(*projection, projectedIndices, nProjectedIndices);
2275     projectLoop2(
2276         0U, 0UL, 0U, 0UL, projection, projector, projectedIndices, nProjectedIndices, scast_minuseq_left<Num2, Num3>());
2277   }
2278 
2279   template <typename Numeric, unsigned StackLen, unsigned StackDim>
2280   template <typename Num2, class Functor>
2281   void ArrayND<Numeric, StackLen, StackDim>::scaleBySliceInnerLoop(const unsigned level,
2282                                                                    const unsigned long idx0,
2283                                                                    Num2& scale,
2284                                                                    const unsigned* projectedIndices,
2285                                                                    const unsigned nProjectedIndices,
2286                                                                    Functor binaryFunct) {
2287     const unsigned idx = projectedIndices[level];
2288     const unsigned imax = shape_[idx];
2289     const unsigned long stride = strides_[idx];
2290 
2291     if (level == nProjectedIndices - 1) {
2292       Numeric* data = data_ + idx0;
2293       for (unsigned i = 0; i < imax; ++i)
2294         binaryFunct(data[i * stride], scale);
2295     } else
2296       for (unsigned i = 0; i < imax; ++i)
2297         scaleBySliceInnerLoop(level + 1, idx0 + i * stride, scale, projectedIndices, nProjectedIndices, binaryFunct);
2298   }
2299 
2300   template <typename Numeric, unsigned StackLen, unsigned StackDim>
2301   template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
2302   void ArrayND<Numeric, StackLen, StackDim>::scaleBySliceLoop(unsigned level,
2303                                                               unsigned long idx0,
2304                                                               unsigned level1,
2305                                                               unsigned long idx1,
2306                                                               ArrayND<Num2, Len2, Dim2>& slice,
2307                                                               const unsigned* projectedIndices,
2308                                                               const unsigned nProjectedIndices,
2309                                                               Functor binaryFunct) {
2310     if (level == dim_) {
2311       assert(level1 == slice.dim_);
2312       Num2& scaleFactor = slice.dim_ ? slice.data_[idx1] : slice.localData_[0];
2313       scaleBySliceInnerLoop(0U, idx0, scaleFactor, projectedIndices, nProjectedIndices, binaryFunct);
2314     } else {
2315       bool fixed = false;
2316       for (unsigned j = 0; j < nProjectedIndices; ++j)
2317         if (projectedIndices[j] == level) {
2318           fixed = true;
2319           break;
2320         }
2321       if (fixed) {
2322         scaleBySliceLoop(level + 1, idx0, level1, idx1, slice, projectedIndices, nProjectedIndices, binaryFunct);
2323       } else {
2324         const unsigned imax = shape_[level];
2325         const unsigned long stride = strides_[level];
2326         const unsigned long stride2 = slice.strides_[level1];
2327         for (unsigned i = 0; i < imax; ++i)
2328           scaleBySliceLoop(level + 1,
2329                            idx0 + i * stride,
2330                            level1 + 1,
2331                            idx1 + i * stride2,
2332                            slice,
2333                            projectedIndices,
2334                            nProjectedIndices,
2335                            binaryFunct);
2336       }
2337     }
2338   }
2339 
2340   template <typename Numeric, unsigned StackLen, unsigned StackDim>
2341   template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
2342   void ArrayND<Numeric, StackLen, StackDim>::jointScan(ArrayND<Num2, Len2, Dim2>& r, Functor binaryFunct) {
2343     if (!isShapeCompatible(r))
2344       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::jointScan: incompatible argument array shape");
2345     if (dim_)
2346       for (unsigned long i = 0; i < len_; ++i)
2347         binaryFunct(data_[i], r.data_[i]);
2348     else
2349       binaryFunct(localData_[0], r.localData_[0]);
2350   }
2351 
2352   template <typename Numeric, unsigned StackLen, unsigned StackDim>
2353   template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
2354   void ArrayND<Numeric, StackLen, StackDim>::applySlice(ArrayND<Num2, Len2, Dim2>& slice,
2355                                                         const unsigned* fixedIndices,
2356                                                         const unsigned nFixedIndices,
2357                                                         Functor binaryFunct) {
2358     if (nFixedIndices) {
2359       verifyProjectionCompatibility(slice, fixedIndices, nFixedIndices);
2360       if (slice.dim_ == 0U)
2361         for (unsigned long i = 0; i < len_; ++i)
2362           binaryFunct(data_[i], slice.localData_[0]);
2363       else
2364         scaleBySliceLoop(0U, 0UL, 0U, 0UL, slice, fixedIndices, nFixedIndices, binaryFunct);
2365     } else
2366       jointScan(slice, binaryFunct);
2367   }
2368 
2369   template <typename Numeric, unsigned StackLen, unsigned StackDim>
2370   inline bool ArrayND<Numeric, StackLen, StackDim>::isCompatible(const ArrayShape& shape) const {
2371     if (!shapeIsKnown_)
2372       return false;
2373     if (dim_ != shape.size())
2374       return false;
2375     if (dim_) {
2376       for (unsigned i = 0; i < dim_; ++i)
2377         if (shape_[i] != shape[i])
2378           return false;
2379     } else
2380       assert(len_ == 1UL);
2381     return true;
2382   }
2383 
2384   template <typename Numeric, unsigned StackLen, unsigned StackDim>
2385   template <typename Num2, unsigned Len2, unsigned Dim2>
2386   inline bool ArrayND<Numeric, StackLen, StackDim>::isShapeCompatible(const ArrayND<Num2, Len2, Dim2>& r) const {
2387     if (!shapeIsKnown_)
2388       return false;
2389     if (!r.shapeIsKnown_)
2390       return false;
2391     if (dim_ != r.dim_)
2392       return false;
2393     if (len_ != r.len_)
2394       return false;
2395     if (dim_) {
2396       assert(shape_);
2397       assert(r.shape_);
2398       for (unsigned i = 0; i < dim_; ++i)
2399         if (shape_[i] != r.shape_[i])
2400           return false;
2401     } else
2402       assert(len_ == 1UL);
2403     return true;
2404   }
2405 
2406   template <typename Numeric, unsigned Len, unsigned Dim>
2407   template <typename Num2, typename Integer>
2408   void ArrayND<Numeric, Len, Dim>::processSubrangeLoop(const unsigned level,
2409                                                        unsigned long idx0,
2410                                                        unsigned* currentIndex,
2411                                                        AbsArrayProjector<Numeric, Num2>& f,
2412                                                        const BoxND<Integer>& subrange) const {
2413     // Deal with possible negative limits first
2414     const Interval<Integer>& levelRange(subrange[level]);
2415     long long int iminl = static_cast<long long int>(levelRange.min());
2416     if (iminl < 0LL)
2417       iminl = 0LL;
2418     long long int imaxl = static_cast<long long int>(levelRange.max());
2419     if (imaxl < 0LL)
2420       imaxl = 0LL;
2421 
2422     // Now deal with possible out-of-range limits
2423     const unsigned imin = static_cast<unsigned>(iminl);
2424     unsigned imax = static_cast<unsigned>(imaxl);
2425     if (imax > shape_[level])
2426       imax = shape_[level];
2427 
2428     if (level == dim_ - 1) {
2429       idx0 += imin;
2430       for (unsigned i = imin; i < imax; ++i, ++idx0) {
2431         currentIndex[level] = i;
2432         f.process(currentIndex, dim_, idx0, data_[idx0]);
2433       }
2434     } else {
2435       const unsigned long stride = strides_[level];
2436       idx0 += imin * stride;
2437       for (unsigned i = imin; i < imax; ++i) {
2438         currentIndex[level] = i;
2439         processSubrangeLoop(level + 1U, idx0, currentIndex, f, subrange);
2440         idx0 += stride;
2441       }
2442     }
2443   }
2444 
2445   template <typename Numeric, unsigned Len, unsigned StackDim>
2446   template <typename Num2, typename Integer>
2447   void ArrayND<Numeric, Len, StackDim>::processSubrange(AbsArrayProjector<Numeric, Num2>& f,
2448                                                         const BoxND<Integer>& subrange) const {
2449     if (!shapeIsKnown_)
2450       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"processSubrange\"");
2451     if (!dim_)
2452       throw npstat::NpstatInvalidArgument(
2453           "npstat::ArrayND::processSubrange method "
2454           "can not be used with array of 0 rank");
2455     if (dim_ != subrange.dim())
2456       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::processSubrange: incompatible subrange rank");
2457     unsigned ibuf[StackDim];
2458     unsigned* buf = makeBuffer(dim_, ibuf, StackDim);
2459     for (unsigned i = 0; i < dim_; ++i)
2460       buf[i] = 0U;
2461     processSubrangeLoop(0U, 0UL, buf, f, subrange);
2462     destroyBuffer(buf, ibuf);
2463   }
2464 
2465   template <typename Numeric, unsigned Len, unsigned Dim>
2466   template <typename Num2>
2467   inline ArrayND<Numeric, Len, Dim>& ArrayND<Numeric, Len, Dim>::setData(const Num2* data,
2468                                                                          const unsigned long dataLength) {
2469     if (!shapeIsKnown_)
2470       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"setData\"");
2471     if (dataLength != len_)
2472       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::setData: incompatible input data length");
2473     copyBuffer(data_, data, dataLength);
2474     return *this;
2475   }
2476 
2477   template <typename Numeric, unsigned Len, unsigned Dim>
2478   void ArrayND<Numeric, Len, Dim>::buildStrides() {
2479     assert(dim_);
2480     if (strides_ == nullptr)
2481       strides_ = makeBuffer(dim_, localStrides_, Dim);
2482     strides_[dim_ - 1] = 1UL;
2483     for (unsigned j = dim_ - 1; j > 0; --j)
2484       strides_[j - 1] = strides_[j] * shape_[j];
2485   }
2486 
2487   template <typename Numeric, unsigned StackLen, unsigned StackDim>
2488   inline ArrayND<Numeric, StackLen, StackDim>::ArrayND()
2489       : data_(nullptr), strides_(nullptr), shape_(nullptr), len_(0UL), dim_(0U), shapeIsKnown_(false) {
2490     localData_[0] = Numeric();
2491     data_ = localData_;
2492   }
2493 
2494   template <typename Numeric, unsigned Len, unsigned Dim>
2495   ArrayND<Numeric, Len, Dim>::ArrayND(const ArrayND& r)
2496       : data_(nullptr), strides_(nullptr), shape_(nullptr), len_(r.len_), dim_(r.dim_), shapeIsKnown_(r.shapeIsKnown_) {
2497     if (dim_) {
2498       // Copy the shape
2499       shape_ = makeBuffer(dim_, localShape_, Dim);
2500       copyBuffer(shape_, r.shape_, dim_);
2501 
2502       // Copy the strides
2503       strides_ = makeBuffer(dim_, localStrides_, Dim);
2504       copyBuffer(strides_, r.strides_, dim_);
2505 
2506       // Copy the data
2507       data_ = makeBuffer(len_, localData_, Len);
2508       copyBuffer(data_, r.data_, len_);
2509     } else {
2510       assert(len_ == 1UL);
2511       localData_[0] = r.localData_[0];
2512       data_ = localData_;
2513     }
2514   }
2515 
2516   template <typename Numeric, unsigned Len, unsigned Dim>
2517   template <typename Num2, unsigned Len2, unsigned Dim2>
2518   ArrayND<Numeric, Len, Dim>::ArrayND(const ArrayND<Num2, Len2, Dim2>& r)
2519       : data_(0), strides_(nullptr), shape_(nullptr), len_(r.len_), dim_(r.dim_), shapeIsKnown_(r.shapeIsKnown_) {
2520     if (dim_) {
2521       // Copy the shape
2522       shape_ = makeBuffer(dim_, localShape_, Dim);
2523       copyBuffer(shape_, r.shape_, dim_);
2524 
2525       // Copy the strides
2526       strides_ = makeBuffer(dim_, localStrides_, Dim);
2527       copyBuffer(strides_, r.strides_, dim_);
2528 
2529       // Copy the data
2530       data_ = makeBuffer(len_, localData_, Len);
2531       copyBuffer(data_, r.data_, len_);
2532     } else {
2533       assert(len_ == 1UL);
2534       localData_[0] = static_cast<Numeric>(r.localData_[0]);
2535       data_ = localData_;
2536     }
2537   }
2538 
2539   template <typename Numeric, unsigned Len, unsigned Dim>
2540   template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
2541   ArrayND<Numeric, Len, Dim>::ArrayND(const ArrayND<Num2, Len2, Dim2>& r, Functor f)
2542       : data_(0), strides_(nullptr), shape_(nullptr), len_(r.len_), dim_(r.dim_), shapeIsKnown_(r.shapeIsKnown_) {
2543     if (dim_) {
2544       // Copy the shape
2545       shape_ = makeBuffer(dim_, localShape_, Dim);
2546       copyBuffer(shape_, r.shape_, dim_);
2547 
2548       // Copy the strides
2549       strides_ = makeBuffer(dim_, localStrides_, Dim);
2550       copyBuffer(strides_, r.strides_, dim_);
2551 
2552       // Copy the data
2553       data_ = makeBuffer(len_, localData_, Len);
2554       for (unsigned long i = 0; i < len_; ++i)
2555         data_[i] = static_cast<Numeric>(f(r.data_[i]));
2556     } else {
2557       assert(len_ == 1UL);
2558       localData_[0] = static_cast<Numeric>(f(r.localData_[0]));
2559       data_ = localData_;
2560     }
2561   }
2562 
2563   template <typename Numeric, unsigned Len, unsigned Dim>
2564   template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
2565   void ArrayND<Numeric, Len, Dim>::copyRangeLoopFunct(const unsigned level,
2566                                                       unsigned long idx0,
2567                                                       unsigned long idx1,
2568                                                       const ArrayND<Num2, Len2, Dim2>& r,
2569                                                       const ArrayRange& range,
2570                                                       Functor f) {
2571     const unsigned imax = shape_[level];
2572     if (level == dim_ - 1) {
2573       Numeric* to = data_ + idx0;
2574       const Num2* from = r.data_ + (idx1 + range[level].min());
2575       for (unsigned i = 0; i < imax; ++i)
2576         *to++ = static_cast<Numeric>(f(*from++));
2577     } else {
2578       const unsigned long fromstride = r.strides_[level];
2579       const unsigned long tostride = strides_[level];
2580       idx1 += range[level].min() * fromstride;
2581       for (unsigned i = 0; i < imax; ++i) {
2582         copyRangeLoopFunct(level + 1, idx0, idx1, r, range, f);
2583         idx0 += tostride;
2584         idx1 += fromstride;
2585       }
2586     }
2587   }
2588 
2589   template <typename Numeric, unsigned Len, unsigned Dim>
2590   template <typename Num2, unsigned Len2, unsigned Dim2>
2591   ArrayND<Numeric, Len, Dim>::ArrayND(const ArrayND<Num2, Len2, Dim2>& r, const ArrayRange& range)
2592       : data_(0), strides_(nullptr), shape_(nullptr), len_(r.len_), dim_(r.dim_), shapeIsKnown_(r.shapeIsKnown_) {
2593     if (!range.isCompatible(r.shape_, r.dim_))
2594       throw npstat::NpstatInvalidArgument("In npstat::ArrayND subrange constructor: invalid subrange");
2595     if (dim_) {
2596       len_ = range.rangeSize();
2597       if (!len_)
2598         throw npstat::NpstatInvalidArgument("In npstat::ArrayND subrange constructor: empty subrange");
2599 
2600       // Figure out the shape
2601       shape_ = makeBuffer(dim_, localShape_, Dim);
2602       range.rangeLength(shape_, dim_);
2603 
2604       // Figure out the strides
2605       buildStrides();
2606 
2607       // Allocate the data array
2608       data_ = makeBuffer(len_, localData_, Len);
2609 
2610       // Copy the data
2611       if (dim_ > CHAR_BIT * sizeof(unsigned long))
2612         throw npstat::NpstatInvalidArgument(
2613             "In npstat::ArrayND subrange constructor: "
2614             "input array rank is too large");
2615       unsigned lolim[CHAR_BIT * sizeof(unsigned long)];
2616       range.lowerLimits(lolim, dim_);
2617       unsigned toBuf[CHAR_BIT * sizeof(unsigned long)];
2618       clearBuffer(toBuf, dim_);
2619       (const_cast<ArrayND<Num2, Len2, Dim2>&>(r))
2620           .commonSubrangeLoop(0U, 0UL, 0UL, lolim, shape_, toBuf, *this, scast_assign_right<Num2, Numeric>());
2621     } else {
2622       assert(len_ == 1UL);
2623       localData_[0] = static_cast<Numeric>(r.localData_[0]);
2624       data_ = localData_;
2625     }
2626   }
2627 
2628   template <typename Numeric, unsigned Len, unsigned Dim>
2629   template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
2630   ArrayND<Numeric, Len, Dim>::ArrayND(const ArrayND<Num2, Len2, Dim2>& r, const ArrayRange& range, Functor f)
2631       : data_(0), strides_(nullptr), shape_(nullptr), len_(r.len_), dim_(r.dim_), shapeIsKnown_(r.shapeIsKnown_) {
2632     if (!range.isCompatible(r.shape_, r.dim_))
2633       throw npstat::NpstatInvalidArgument(
2634           "In npstat::ArrayND transforming subrange constructor: "
2635           "incompatible subrange");
2636     if (dim_) {
2637       len_ = range.rangeSize();
2638       if (!len_)
2639         throw npstat::NpstatInvalidArgument(
2640             "In npstat::ArrayND transforming subrange constructor: "
2641             "empty subrange");
2642 
2643       // Figure out the shape
2644       shape_ = makeBuffer(dim_, localShape_, Dim);
2645       for (unsigned i = 0; i < dim_; ++i)
2646         shape_[i] = range[i].length();
2647 
2648       // Figure out the strides
2649       buildStrides();
2650 
2651       // Allocate the data array
2652       data_ = makeBuffer(len_, localData_, Len);
2653 
2654       // Transform the data
2655       copyRangeLoopFunct(0U, 0UL, 0UL, r, range, f);
2656     } else {
2657       assert(len_ == 1UL);
2658       localData_[0] = static_cast<Numeric>(f(r.localData_[0]));
2659       data_ = localData_;
2660     }
2661   }
2662 
2663   template <typename Numeric, unsigned Len, unsigned Dim>
2664   ArrayND<Numeric, Len, Dim>::ArrayND(const ArrayShape& sh)
2665       : data_(nullptr), strides_(nullptr), shape_(nullptr), len_(1UL), shapeIsKnown_(true) {
2666     const unsigned sz = sh.size();
2667     buildFromShapePtr(sz ? &sh[0] : nullptr, sz);
2668   }
2669 
2670   template <typename Numeric, unsigned Len, unsigned Dim>
2671   ArrayND<Numeric, Len, Dim>::ArrayND(const unsigned* sizes, const unsigned dim)
2672       : data_(0), strides_(nullptr), shape_(nullptr), len_(1UL), shapeIsKnown_(true) {
2673     buildFromShapePtr(sizes, dim);
2674   }
2675 
2676   template <typename Numeric, unsigned Len, unsigned Dim>
2677   ArrayND<Numeric, Len, Dim>::ArrayND(const unsigned n0)
2678       : data_(0), strides_(nullptr), shape_(nullptr), len_(1UL), shapeIsKnown_(true) {
2679     const unsigned dim = 1U;
2680     unsigned sizes[dim];
2681     sizes[0] = n0;
2682     buildFromShapePtr(sizes, dim);
2683   }
2684 
2685   template <typename Numeric, unsigned Len, unsigned Dim>
2686   ArrayND<Numeric, Len, Dim>::ArrayND(const unsigned n0, const unsigned n1)
2687       : data_(0), strides_(nullptr), shape_(nullptr), len_(1UL), shapeIsKnown_(true) {
2688     const unsigned dim = 2U;
2689     unsigned sizes[dim];
2690     sizes[0] = n0;
2691     sizes[1] = n1;
2692     buildFromShapePtr(sizes, dim);
2693   }
2694 
2695   template <typename Numeric, unsigned Len, unsigned Dim>
2696   ArrayND<Numeric, Len, Dim>::ArrayND(const unsigned n0, const unsigned n1, const unsigned n2)
2697       : data_(0), strides_(nullptr), shape_(nullptr), len_(1UL), shapeIsKnown_(true) {
2698     const unsigned dim = 3U;
2699     unsigned sizes[dim];
2700     sizes[0] = n0;
2701     sizes[1] = n1;
2702     sizes[2] = n2;
2703     buildFromShapePtr(sizes, dim);
2704   }
2705 
2706   template <typename Numeric, unsigned Len, unsigned Dim>
2707   ArrayND<Numeric, Len, Dim>::ArrayND(const unsigned n0, const unsigned n1, const unsigned n2, const unsigned n3)
2708       : data_(0), strides_(nullptr), shape_(nullptr), len_(1UL), shapeIsKnown_(true) {
2709     const unsigned dim = 4U;
2710     unsigned sizes[dim];
2711     sizes[0] = n0;
2712     sizes[1] = n1;
2713     sizes[2] = n2;
2714     sizes[3] = n3;
2715     buildFromShapePtr(sizes, dim);
2716   }
2717 
2718   template <typename Numeric, unsigned Len, unsigned Dim>
2719   ArrayND<Numeric, Len, Dim>::ArrayND(
2720       const unsigned n0, const unsigned n1, const unsigned n2, const unsigned n3, const unsigned n4)
2721       : data_(0), strides_(nullptr), shape_(nullptr), len_(1UL), shapeIsKnown_(true) {
2722     const unsigned dim = 5U;
2723     unsigned sizes[dim];
2724     sizes[0] = n0;
2725     sizes[1] = n1;
2726     sizes[2] = n2;
2727     sizes[3] = n3;
2728     sizes[4] = n4;
2729     buildFromShapePtr(sizes, dim);
2730   }
2731 
2732   template <typename Numeric, unsigned Len, unsigned Dim>
2733   ArrayND<Numeric, Len, Dim>::ArrayND(
2734       const unsigned n0, const unsigned n1, const unsigned n2, const unsigned n3, const unsigned n4, const unsigned n5)
2735       : data_(0), strides_(nullptr), shape_(nullptr), len_(1UL), shapeIsKnown_(true) {
2736     const unsigned dim = 6U;
2737     unsigned sizes[dim];
2738     sizes[0] = n0;
2739     sizes[1] = n1;
2740     sizes[2] = n2;
2741     sizes[3] = n3;
2742     sizes[4] = n4;
2743     sizes[5] = n5;
2744     buildFromShapePtr(sizes, dim);
2745   }
2746 
2747   template <typename Numeric, unsigned Len, unsigned Dim>
2748   ArrayND<Numeric, Len, Dim>::ArrayND(const unsigned n0,
2749                                       const unsigned n1,
2750                                       const unsigned n2,
2751                                       const unsigned n3,
2752                                       const unsigned n4,
2753                                       const unsigned n5,
2754                                       const unsigned n6)
2755       : data_(0), strides_(nullptr), shape_(nullptr), len_(1UL), shapeIsKnown_(true) {
2756     const unsigned dim = 7U;
2757     unsigned sizes[dim];
2758     sizes[0] = n0;
2759     sizes[1] = n1;
2760     sizes[2] = n2;
2761     sizes[3] = n3;
2762     sizes[4] = n4;
2763     sizes[5] = n5;
2764     sizes[6] = n6;
2765     buildFromShapePtr(sizes, dim);
2766   }
2767 
2768   template <typename Numeric, unsigned Len, unsigned Dim>
2769   ArrayND<Numeric, Len, Dim>::ArrayND(const unsigned n0,
2770                                       const unsigned n1,
2771                                       const unsigned n2,
2772                                       const unsigned n3,
2773                                       const unsigned n4,
2774                                       const unsigned n5,
2775                                       const unsigned n6,
2776                                       const unsigned n7)
2777       : data_(0), strides_(nullptr), shape_(nullptr), len_(1UL), shapeIsKnown_(true) {
2778     const unsigned dim = 8U;
2779     unsigned sizes[dim];
2780     sizes[0] = n0;
2781     sizes[1] = n1;
2782     sizes[2] = n2;
2783     sizes[3] = n3;
2784     sizes[4] = n4;
2785     sizes[5] = n5;
2786     sizes[6] = n6;
2787     sizes[7] = n7;
2788     buildFromShapePtr(sizes, dim);
2789   }
2790 
2791   template <typename Numeric, unsigned Len, unsigned Dim>
2792   ArrayND<Numeric, Len, Dim>::ArrayND(const unsigned n0,
2793                                       const unsigned n1,
2794                                       const unsigned n2,
2795                                       const unsigned n3,
2796                                       const unsigned n4,
2797                                       const unsigned n5,
2798                                       const unsigned n6,
2799                                       const unsigned n7,
2800                                       const unsigned n8)
2801       : data_(0), strides_(nullptr), shape_(nullptr), len_(1UL), shapeIsKnown_(true) {
2802     const unsigned dim = 9U;
2803     unsigned sizes[dim];
2804     sizes[0] = n0;
2805     sizes[1] = n1;
2806     sizes[2] = n2;
2807     sizes[3] = n3;
2808     sizes[4] = n4;
2809     sizes[5] = n5;
2810     sizes[6] = n6;
2811     sizes[7] = n7;
2812     sizes[8] = n8;
2813     buildFromShapePtr(sizes, dim);
2814   }
2815 
2816   template <typename Numeric, unsigned Len, unsigned Dim>
2817   ArrayND<Numeric, Len, Dim>::ArrayND(const unsigned n0,
2818                                       const unsigned n1,
2819                                       const unsigned n2,
2820                                       const unsigned n3,
2821                                       const unsigned n4,
2822                                       const unsigned n5,
2823                                       const unsigned n6,
2824                                       const unsigned n7,
2825                                       const unsigned n8,
2826                                       const unsigned n9)
2827       : data_(0), strides_(nullptr), shape_(nullptr), len_(1UL), shapeIsKnown_(true) {
2828     const unsigned dim = 10U;
2829     unsigned sizes[dim];
2830     sizes[0] = n0;
2831     sizes[1] = n1;
2832     sizes[2] = n2;
2833     sizes[3] = n3;
2834     sizes[4] = n4;
2835     sizes[5] = n5;
2836     sizes[6] = n6;
2837     sizes[7] = n7;
2838     sizes[8] = n8;
2839     sizes[9] = n9;
2840     buildFromShapePtr(sizes, dim);
2841   }
2842 
2843   template <typename Numeric, unsigned Len, unsigned Dim>
2844   template <typename Num1, unsigned Len1, unsigned Dim1, typename Num2, unsigned Len2, unsigned Dim2>
2845   void ArrayND<Numeric, Len, Dim>::outerProductLoop(const unsigned level,
2846                                                     unsigned long idx0,
2847                                                     unsigned long idx1,
2848                                                     unsigned long idx2,
2849                                                     const ArrayND<Num1, Len1, Dim1>& a1,
2850                                                     const ArrayND<Num2, Len2, Dim2>& a2) {
2851     const unsigned imax = shape_[level];
2852     if (level == dim_ - 1) {
2853       for (unsigned i = 0; i < imax; ++i)
2854         data_[idx0 + i] = a1.data_[idx1] * a2.data_[idx2 + i];
2855     } else {
2856       for (unsigned i = 0; i < imax; ++i) {
2857         outerProductLoop(level + 1, idx0, idx1, idx2, a1, a2);
2858         idx0 += strides_[level];
2859         if (level < a1.dim_)
2860           idx1 += a1.strides_[level];
2861         else
2862           idx2 += a2.strides_[level - a1.dim_];
2863       }
2864     }
2865   }
2866 
2867   template <typename Numeric, unsigned Len, unsigned Dim>
2868   template <typename Num1, unsigned Len1, unsigned Dim1, typename Num2, unsigned Len2, unsigned Dim2>
2869   ArrayND<Numeric, Len, Dim>::ArrayND(const ArrayND<Num1, Len1, Dim1>& a1, const ArrayND<Num2, Len2, Dim2>& a2)
2870       : data_(0), strides_(nullptr), shape_(nullptr), len_(1UL), dim_(a1.dim_ + a2.dim_), shapeIsKnown_(true) {
2871     if (!(a1.shapeIsKnown_ && a2.shapeIsKnown_))
2872       throw npstat::NpstatInvalidArgument(
2873           "In npstat::ArrayND outer product constructor: "
2874           "uninitialized argument array");
2875     if (dim_) {
2876       shape_ = makeBuffer(dim_, localShape_, Dim);
2877       copyBuffer(shape_, a1.shape_, a1.dim_);
2878       copyBuffer(shape_ + a1.dim_, a2.shape_, a2.dim_);
2879 
2880       for (unsigned i = 0; i < dim_; ++i) {
2881         assert(shape_[i]);
2882         len_ *= shape_[i];
2883       }
2884 
2885       // Figure out the array strides
2886       buildStrides();
2887 
2888       // Allocate the data array
2889       data_ = makeBuffer(len_, localData_, Len);
2890 
2891       // Fill the data array
2892       if (a1.dim_ == 0) {
2893         for (unsigned long i = 0; i < len_; ++i)
2894           data_[i] = a1.localData_[0] * a2.data_[i];
2895       } else if (a2.dim_ == 0) {
2896         for (unsigned long i = 0; i < len_; ++i)
2897           data_[i] = a1.data_[i] * a2.localData_[0];
2898       } else
2899         outerProductLoop(0U, 0UL, 0UL, 0UL, a1, a2);
2900     } else {
2901       localData_[0] = a1.localData_[0] * a2.localData_[0];
2902       data_ = localData_;
2903     }
2904   }
2905 
2906   template <typename Numeric, unsigned Len, unsigned Dim>
2907   inline ArrayND<Numeric, Len, Dim>::~ArrayND() {
2908     destroyBuffer(data_, localData_);
2909     destroyBuffer(strides_, localStrides_);
2910     destroyBuffer(shape_, localShape_);
2911   }
2912 
2913   template <typename Numeric, unsigned Len, unsigned Dim>
2914   ArrayND<Numeric, Len, Dim>& ArrayND<Numeric, Len, Dim>::operator=(const ArrayND& r) {
2915     if (this == &r)
2916       return *this;
2917     if (shapeIsKnown_) {
2918       if (!r.shapeIsKnown_)
2919         throw npstat::NpstatInvalidArgument(
2920             "In npstat::ArrayND assignment operator: "
2921             "uninitialized argument array");
2922       if (!isShapeCompatible(r))
2923         throw npstat::NpstatInvalidArgument(
2924             "In npstat::ArrayND assignment operator: "
2925             "incompatible argument array shape");
2926       if (dim_)
2927         copyBuffer(data_, r.data_, len_);
2928       else
2929         localData_[0] = r.localData_[0];
2930     } else {
2931       // This object is uninitialized. If the object on the
2932       // right is itself initialized, make an in-place copy.
2933       if (r.shapeIsKnown_)
2934         new (this) ArrayND(r);
2935     }
2936     return *this;
2937   }
2938 
2939   template <typename Numeric, unsigned Len, unsigned Dim>
2940   template <typename Num2, unsigned Len2, unsigned Dim2>
2941   ArrayND<Numeric, Len, Dim>& ArrayND<Numeric, Len, Dim>::operator=(const ArrayND<Num2, Len2, Dim2>& r) {
2942     if ((void*)this == (void*)(&r))
2943       return *this;
2944     if (shapeIsKnown_) {
2945       if (!r.shapeIsKnown_)
2946         throw npstat::NpstatInvalidArgument(
2947             "In npstat::ArrayND assignment operator: "
2948             "uninitialized argument array");
2949       if (!isShapeCompatible(r))
2950         throw npstat::NpstatInvalidArgument(
2951             "In npstat::ArrayND assignment operator: "
2952             "incompatible argument array shape");
2953       if (dim_)
2954         copyBuffer(data_, r.data_, len_);
2955       else
2956         localData_[0] = static_cast<Numeric>(r.localData_[0]);
2957     } else {
2958       // This object is uninitialized. If the object on the
2959       // right is itself initialized, make an in-place copy.
2960       if (r.shapeIsKnown_)
2961         new (this) ArrayND(r);
2962     }
2963     return *this;
2964   }
2965 
2966   template <typename Numeric, unsigned Len, unsigned Dim>
2967   template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
2968   ArrayND<Numeric, Len, Dim>& ArrayND<Numeric, Len, Dim>::assign(const ArrayND<Num2, Len2, Dim2>& r, Functor f) {
2969     if (shapeIsKnown_) {
2970       if (!r.shapeIsKnown_)
2971         throw npstat::NpstatInvalidArgument("In npstat::ArrayND::assign: uninitialized argument array");
2972       if (!isShapeCompatible(r))
2973         throw npstat::NpstatInvalidArgument("In npstat::ArrayND::assign: incompatible argument array shape");
2974       if (dim_)
2975         for (unsigned long i = 0; i < len_; ++i)
2976           data_[i] = static_cast<Numeric>(f(r.data_[i]));
2977       else
2978         localData_[0] = static_cast<Numeric>(f(r.localData_[0]));
2979     } else {
2980       // This object is uninitialized. If the object on the
2981       // right is itself initialized, build new array in place.
2982       if (r.shapeIsKnown_)
2983         new (this) ArrayND(r, f);
2984     }
2985     return *this;
2986   }
2987 
2988   template <typename Numeric, unsigned Len, unsigned Dim>
2989   inline ArrayShape ArrayND<Numeric, Len, Dim>::shape() const {
2990     if (!shapeIsKnown_)
2991       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"shape\"");
2992     return ArrayShape(shape_, shape_ + dim_);
2993   }
2994 
2995   template <typename Numeric, unsigned Len, unsigned Dim>
2996   inline ArrayRange ArrayND<Numeric, Len, Dim>::fullRange() const {
2997     if (!shapeIsKnown_)
2998       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"fullRange\"");
2999     ArrayRange range;
3000     if (dim_) {
3001       range.reserve(dim_);
3002       for (unsigned i = 0; i < dim_; ++i)
3003         range.push_back(Interval<unsigned>(0U, shape_[i]));
3004     }
3005     return range;
3006   }
3007 
3008   template <typename Numeric, unsigned Len, unsigned Dim>
3009   bool ArrayND<Numeric, Len, Dim>::isDensity() const {
3010     if (!shapeIsKnown_)
3011       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"isDensity\"");
3012     const Numeric zero = Numeric();
3013     bool hasPositive = false;
3014     if (dim_)
3015       for (unsigned long i = 0; i < len_; ++i) {
3016         // Don't make comparisons whose result can be
3017         // determined in advance by assuming that Numeric
3018         // is an unsigned type. Some compilers will
3019         // complain about it when this template is
3020         // instantiated with such a type.
3021         if (data_[i] == zero)
3022           continue;
3023         if (ComplexComparesFalse<Numeric>::less(zero, data_[i]))
3024           hasPositive = true;
3025         else
3026           return false;
3027       }
3028     else
3029       hasPositive = ComplexComparesFalse<Numeric>::less(zero, localData_[0]);
3030     return hasPositive;
3031   }
3032 
3033   template <typename Numeric, unsigned Len, unsigned Dim>
3034   bool ArrayND<Numeric, Len, Dim>::isZero() const {
3035     if (!shapeIsKnown_)
3036       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"isZero\"");
3037     const Numeric zero = Numeric();
3038     if (dim_) {
3039       for (unsigned long i = 0; i < len_; ++i)
3040         if (data_[i] != zero)
3041           return false;
3042     } else if (localData_[0] != zero)
3043       return false;
3044     return true;
3045   }
3046 
3047   template <typename Numeric, unsigned Len, unsigned Dim>
3048   void ArrayND<Numeric, Len, Dim>::convertLinearIndex(unsigned long l, unsigned* idx, const unsigned idxLen) const {
3049     if (!shapeIsKnown_)
3050       throw npstat::NpstatInvalidArgument(
3051           "Initialize npstat::ArrayND before calling "
3052           "method \"convertLinearIndex\"");
3053     if (!dim_)
3054       throw npstat::NpstatInvalidArgument(
3055           "npstat::ArrayND::convertLinearIndex method "
3056           "can not be used with array of 0 rank");
3057     if (idxLen != dim_)
3058       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::convertLinearIndex: incompatible index length");
3059     if (l >= len_)
3060       throw npstat::NpstatOutOfRange("In npstat::ArrayND::convertLinearIndex: linear index out of range");
3061     assert(idx);
3062 
3063     for (unsigned i = 0; i < dim_; ++i) {
3064       idx[i] = l / strides_[i];
3065       l -= (idx[i] * strides_[i]);
3066     }
3067   }
3068 
3069   template <typename Numeric, unsigned Len, unsigned Dim>
3070   unsigned long ArrayND<Numeric, Len, Dim>::linearIndex(const unsigned* index, unsigned idxLen) const {
3071     if (!shapeIsKnown_)
3072       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"linearIndex\"");
3073     if (!dim_)
3074       throw npstat::NpstatInvalidArgument(
3075           "npstat::ArrayND::linearIndex method "
3076           "can not be used with array of 0 rank");
3077     if (idxLen != dim_)
3078       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::linearIndex: incompatible index length");
3079     assert(index);
3080 
3081     unsigned long idx = 0UL;
3082     for (unsigned i = 0; i < dim_; ++i) {
3083       if (index[i] >= shape_[i])
3084         throw npstat::NpstatOutOfRange("In npstat::ArrayND::linearIndex: index out of range");
3085       idx += index[i] * strides_[i];
3086     }
3087     return idx;
3088   }
3089 
3090   template <typename Numeric, unsigned Len, unsigned Dim>
3091   inline Numeric& ArrayND<Numeric, Len, Dim>::value(const unsigned* index, const unsigned dim) {
3092     if (!shapeIsKnown_)
3093       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"value\"");
3094     if (dim != dim_)
3095       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::value: incompatible index length");
3096     if (dim) {
3097       assert(index);
3098       unsigned long idx = 0UL;
3099       for (unsigned i = 0; i < dim_; ++i)
3100         idx += index[i] * strides_[i];
3101       return data_[idx];
3102     } else
3103       return localData_[0];
3104   }
3105 
3106   template <typename Numeric, unsigned Len, unsigned Dim>
3107   inline const Numeric& ArrayND<Numeric, Len, Dim>::value(const unsigned* index, const unsigned dim) const {
3108     if (!shapeIsKnown_)
3109       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"value\"");
3110     if (dim != dim_)
3111       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::value: incompatible index length");
3112     if (dim) {
3113       assert(index);
3114       unsigned long idx = 0UL;
3115       for (unsigned i = 0; i < dim_; ++i)
3116         idx += index[i] * strides_[i];
3117       return data_[idx];
3118     } else
3119       return localData_[0];
3120   }
3121 
3122   template <typename Numeric, unsigned Len, unsigned Dim>
3123   inline Numeric& ArrayND<Numeric, Len, Dim>::linearValue(const unsigned long index) {
3124     return data_[index];
3125   }
3126 
3127   template <typename Numeric, unsigned Len, unsigned Dim>
3128   inline const Numeric& ArrayND<Numeric, Len, Dim>::linearValue(const unsigned long index) const {
3129     return data_[index];
3130   }
3131 
3132   template <typename Numeric, unsigned Len, unsigned Dim>
3133   inline Numeric& ArrayND<Numeric, Len, Dim>::linearValueAt(const unsigned long index) {
3134     if (index >= len_)
3135       throw npstat::NpstatOutOfRange("In npstat::ArrayND::linearValueAt: linear index out of range");
3136     return data_[index];
3137   }
3138 
3139   template <typename Numeric, unsigned Len, unsigned Dim>
3140   inline const Numeric& ArrayND<Numeric, Len, Dim>::linearValueAt(const unsigned long index) const {
3141     if (index >= len_)
3142       throw npstat::NpstatOutOfRange("In npstat::ArrayND::linearValueAt: linear index out of range");
3143     return data_[index];
3144   }
3145 
3146   template <typename Numeric, unsigned Len, unsigned Dim>
3147   inline unsigned ArrayND<Numeric, Len, Dim>::coordToIndex(const double x, const unsigned idim) const {
3148     if (x <= 0.0)
3149       return 0;
3150     else if (x >= static_cast<double>(shape_[idim] - 1))
3151       return shape_[idim] - 1;
3152     else
3153       return static_cast<unsigned>(std::floor(x + 0.5));
3154   }
3155 
3156   template <typename Numeric, unsigned Len, unsigned Dim>
3157   inline const Numeric& ArrayND<Numeric, Len, Dim>::closest(const double* x, const unsigned dim) const {
3158     if (!shapeIsKnown_)
3159       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"closest\"");
3160     if (dim != dim_)
3161       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::closest: incompatible data length");
3162     if (dim) {
3163       assert(x);
3164       unsigned long idx = 0UL;
3165       for (unsigned i = 0; i < dim_; ++i)
3166         idx += coordToIndex(x[i], i) * strides_[i];
3167       return data_[idx];
3168     } else
3169       return localData_[0];
3170   }
3171 
3172   template <typename Numeric, unsigned Len, unsigned Dim>
3173   inline Numeric& ArrayND<Numeric, Len, Dim>::closest(const double* x, const unsigned dim) {
3174     if (!shapeIsKnown_)
3175       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"closest\"");
3176     if (dim != dim_)
3177       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::closest: incompatible data length");
3178     if (dim) {
3179       assert(x);
3180       unsigned long idx = 0UL;
3181       for (unsigned i = 0; i < dim_; ++i)
3182         idx += coordToIndex(x[i], i) * strides_[i];
3183       return data_[idx];
3184     } else
3185       return localData_[0];
3186   }
3187 
3188   template <typename Numeric, unsigned Len, unsigned Dim>
3189   inline const Numeric& ArrayND<Numeric, Len, Dim>::valueAt(const unsigned* index, const unsigned dim) const {
3190     if (!shapeIsKnown_)
3191       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"valueAt\"");
3192     if (dim != dim_)
3193       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::valueAt: incompatible index length");
3194     if (dim) {
3195       assert(index);
3196       unsigned long idx = 0UL;
3197       for (unsigned i = 0; i < dim_; ++i) {
3198         if (index[i] >= shape_[i])
3199           throw npstat::NpstatOutOfRange("In npstat::ArrayND::valueAt: index out of range");
3200         idx += index[i] * strides_[i];
3201       }
3202       return data_[idx];
3203     } else
3204       return localData_[0];
3205   }
3206 
3207   template <typename Numeric, unsigned Len, unsigned Dim>
3208   inline Numeric& ArrayND<Numeric, Len, Dim>::valueAt(const unsigned* index, const unsigned dim) {
3209     if (!shapeIsKnown_)
3210       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"valueAt\"");
3211     if (dim != dim_)
3212       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::valueAt: incompatible index length");
3213     if (dim) {
3214       assert(index);
3215       unsigned long idx = 0UL;
3216       for (unsigned i = 0; i < dim_; ++i) {
3217         if (index[i] >= shape_[i])
3218           throw npstat::NpstatOutOfRange("In npstat::ArrayND::valueAt: index out of range");
3219         idx += index[i] * strides_[i];
3220       }
3221       return data_[idx];
3222     } else
3223       return localData_[0];
3224   }
3225 
3226   template <typename Numeric, unsigned Len, unsigned Dim>
3227   inline Numeric& ArrayND<Numeric, Len, Dim>::operator()() {
3228     if (!shapeIsKnown_)
3229       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"operator()\"");
3230     if (dim_)
3231       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::operator(): wrong # of args (not rank 0 array)");
3232     return localData_[0];
3233   }
3234 
3235   template <typename Numeric, unsigned Len, unsigned Dim>
3236   inline const Numeric& ArrayND<Numeric, Len, Dim>::operator()() const {
3237     if (!shapeIsKnown_)
3238       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"operator()\"");
3239     if (dim_)
3240       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::operator(): wrong # of args (not rank 0 array)");
3241     return localData_[0];
3242   }
3243 
3244   template <typename Numeric, unsigned Len, unsigned Dim>
3245   inline Numeric& ArrayND<Numeric, Len, Dim>::operator()(const unsigned i) {
3246     if (1U != dim_)
3247       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::operator(): wrong # of args (not rank 1 array)");
3248     return data_[i];
3249   }
3250 
3251   template <typename Numeric, unsigned Len, unsigned Dim>
3252   inline const Numeric& ArrayND<Numeric, Len, Dim>::operator()(const unsigned i) const {
3253     if (1U != dim_)
3254       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::operator(): wrong # of args (not rank 1 array)");
3255     return data_[i];
3256   }
3257 
3258   template <typename Numeric, unsigned Len, unsigned Dim>
3259   const Numeric& ArrayND<Numeric, Len, Dim>::at() const {
3260     if (!shapeIsKnown_)
3261       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"at\"");
3262     if (dim_)
3263       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::at: wrong # of args (not rank 0 array)");
3264     return localData_[0];
3265   }
3266 
3267   template <typename Numeric, unsigned Len, unsigned Dim>
3268   Numeric& ArrayND<Numeric, Len, Dim>::at() {
3269     if (!shapeIsKnown_)
3270       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"at\"");
3271     if (dim_)
3272       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::at: wrong # of args (not rank 0 array)");
3273     return localData_[0];
3274   }
3275 
3276   template <typename Numeric, unsigned Len, unsigned Dim>
3277   const Numeric& ArrayND<Numeric, Len, Dim>::at(const unsigned i0) const {
3278     if (1U != dim_)
3279       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::at: wrong # of args (not rank 1 array)");
3280     if (i0 >= shape_[0])
3281       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 0 out of range (rank 1)");
3282     return data_[i0];
3283   }
3284 
3285   template <typename Numeric, unsigned Len, unsigned Dim>
3286   Numeric& ArrayND<Numeric, Len, Dim>::at(const unsigned i0) {
3287     if (1U != dim_)
3288       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::at: wrong # of args (not rank 1 array)");
3289     if (i0 >= shape_[0])
3290       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 0 out of range (rank 1)");
3291     return data_[i0];
3292   }
3293 
3294   template <typename Numeric, unsigned Len, unsigned Dim>
3295   inline Numeric& ArrayND<Numeric, Len, Dim>::operator()(const unsigned i0, const unsigned i1) {
3296     if (2U != dim_)
3297       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::operator(): wrong # of args (not rank 2 array)");
3298     return data_[i0 * strides_[0] + i1];
3299   }
3300 
3301   template <typename Numeric, unsigned Len, unsigned Dim>
3302   inline const Numeric& ArrayND<Numeric, Len, Dim>::operator()(const unsigned i0, const unsigned i1) const {
3303     if (2U != dim_)
3304       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::operator(): wrong # of args (not rank 2 array)");
3305     return data_[i0 * strides_[0] + i1];
3306   }
3307 
3308   template <typename Numeric, unsigned Len, unsigned Dim>
3309   const Numeric& ArrayND<Numeric, Len, Dim>::at(const unsigned i0, const unsigned i1) const {
3310     if (2U != dim_)
3311       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::at: wrong # of args (not rank 2 array)");
3312     if (i0 >= shape_[0])
3313       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 0 out of range (rank 2)");
3314     if (i1 >= shape_[1])
3315       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 1 out of range (rank 2)");
3316     return data_[i0 * strides_[0] + i1];
3317   }
3318 
3319   template <typename Numeric, unsigned Len, unsigned Dim>
3320   Numeric& ArrayND<Numeric, Len, Dim>::at(const unsigned i0, const unsigned i1) {
3321     if (2U != dim_)
3322       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::at: wrong # of args (not rank 2 array)");
3323     if (i0 >= shape_[0])
3324       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 0 out of range (rank 2)");
3325     if (i1 >= shape_[1])
3326       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 1 out of range (rank 2)");
3327     return data_[i0 * strides_[0] + i1];
3328   }
3329 
3330   template <typename Numeric, unsigned Len, unsigned Dim>
3331   inline const Numeric& ArrayND<Numeric, Len, Dim>::operator()(const unsigned i0,
3332                                                                const unsigned i1,
3333                                                                const unsigned i2) const {
3334     if (3U != dim_)
3335       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::operator(): wrong # of args (not rank 3 array)");
3336     return data_[i0 * strides_[0] + i1 * strides_[1] + i2];
3337   }
3338 
3339   template <typename Numeric, unsigned Len, unsigned Dim>
3340   inline const Numeric& ArrayND<Numeric, Len, Dim>::operator()(const unsigned i0,
3341                                                                const unsigned i1,
3342                                                                const unsigned i2,
3343                                                                const unsigned i3) const {
3344     if (4U != dim_)
3345       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::operator(): wrong # of args (not rank 4 array)");
3346     return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3];
3347   }
3348 
3349   template <typename Numeric, unsigned Len, unsigned Dim>
3350   inline const Numeric& ArrayND<Numeric, Len, Dim>::operator()(
3351       const unsigned i0, const unsigned i1, const unsigned i2, const unsigned i3, const unsigned i4) const {
3352     if (5U != dim_)
3353       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::operator(): wrong # of args (not rank 5 array)");
3354     return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3 * strides_[3] + i4];
3355   }
3356 
3357   template <typename Numeric, unsigned Len, unsigned Dim>
3358   inline const Numeric& ArrayND<Numeric, Len, Dim>::operator()(const unsigned i0,
3359                                                                const unsigned i1,
3360                                                                const unsigned i2,
3361                                                                const unsigned i3,
3362                                                                const unsigned i4,
3363                                                                const unsigned i5) const {
3364     if (6U != dim_)
3365       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::operator(): wrong # of args (not rank 6 array)");
3366     return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3 * strides_[3] + i4 * strides_[4] + i5];
3367   }
3368 
3369   template <typename Numeric, unsigned Len, unsigned Dim>
3370   inline const Numeric& ArrayND<Numeric, Len, Dim>::operator()(const unsigned i0,
3371                                                                const unsigned i1,
3372                                                                const unsigned i2,
3373                                                                const unsigned i3,
3374                                                                const unsigned i4,
3375                                                                const unsigned i5,
3376                                                                const unsigned i6) const {
3377     if (7U != dim_)
3378       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::operator(): wrong # of args (not rank 7 array)");
3379     return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3 * strides_[3] + i4 * strides_[4] +
3380                  i5 * strides_[5] + i6];
3381   }
3382 
3383   template <typename Numeric, unsigned Len, unsigned Dim>
3384   inline const Numeric& ArrayND<Numeric, Len, Dim>::operator()(const unsigned i0,
3385                                                                const unsigned i1,
3386                                                                const unsigned i2,
3387                                                                const unsigned i3,
3388                                                                const unsigned i4,
3389                                                                const unsigned i5,
3390                                                                const unsigned i6,
3391                                                                const unsigned i7) const {
3392     if (8U != dim_)
3393       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::operator(): wrong # of args (not rank 8 array)");
3394     return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3 * strides_[3] + i4 * strides_[4] +
3395                  i5 * strides_[5] + i6 * strides_[6] + i7];
3396   }
3397 
3398   template <typename Numeric, unsigned Len, unsigned Dim>
3399   inline const Numeric& ArrayND<Numeric, Len, Dim>::operator()(const unsigned i0,
3400                                                                const unsigned i1,
3401                                                                const unsigned i2,
3402                                                                const unsigned i3,
3403                                                                const unsigned i4,
3404                                                                const unsigned i5,
3405                                                                const unsigned i6,
3406                                                                const unsigned i7,
3407                                                                const unsigned i8) const {
3408     if (9U != dim_)
3409       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::operator(): wrong # of args (not rank 9 array)");
3410     return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3 * strides_[3] + i4 * strides_[4] +
3411                  i5 * strides_[5] + i6 * strides_[6] + i7 * strides_[7] + i8];
3412   }
3413 
3414   template <typename Numeric, unsigned Len, unsigned Dim>
3415   inline const Numeric& ArrayND<Numeric, Len, Dim>::operator()(const unsigned i0,
3416                                                                const unsigned i1,
3417                                                                const unsigned i2,
3418                                                                const unsigned i3,
3419                                                                const unsigned i4,
3420                                                                const unsigned i5,
3421                                                                const unsigned i6,
3422                                                                const unsigned i7,
3423                                                                const unsigned i8,
3424                                                                const unsigned i9) const {
3425     if (10U != dim_)
3426       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::operator(): wrong # of args (not rank 10 array)");
3427     return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3 * strides_[3] + i4 * strides_[4] +
3428                  i5 * strides_[5] + i6 * strides_[6] + i7 * strides_[7] + i8 * strides_[8] + i9];
3429   }
3430 
3431   template <typename Numeric, unsigned Len, unsigned Dim>
3432   inline Numeric& ArrayND<Numeric, Len, Dim>::operator()(const unsigned i0, const unsigned i1, const unsigned i2) {
3433     if (3U != dim_)
3434       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::operator(): wrong # of args (not rank 3 array)");
3435     return data_[i0 * strides_[0] + i1 * strides_[1] + i2];
3436   }
3437 
3438   template <typename Numeric, unsigned Len, unsigned Dim>
3439   inline Numeric& ArrayND<Numeric, Len, Dim>::operator()(const unsigned i0,
3440                                                          const unsigned i1,
3441                                                          const unsigned i2,
3442                                                          const unsigned i3) {
3443     if (4U != dim_)
3444       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::operator(): wrong # of args (not rank 4 array)");
3445     return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3];
3446   }
3447 
3448   template <typename Numeric, unsigned Len, unsigned Dim>
3449   inline Numeric& ArrayND<Numeric, Len, Dim>::operator()(
3450       const unsigned i0, const unsigned i1, const unsigned i2, const unsigned i3, const unsigned i4) {
3451     if (5U != dim_)
3452       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::operator(): wrong # of args (not rank 5 array)");
3453     return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3 * strides_[3] + i4];
3454   }
3455 
3456   template <typename Numeric, unsigned Len, unsigned Dim>
3457   inline Numeric& ArrayND<Numeric, Len, Dim>::operator()(const unsigned i0,
3458                                                          const unsigned i1,
3459                                                          const unsigned i2,
3460                                                          const unsigned i3,
3461                                                          const unsigned i4,
3462                                                          const unsigned i5) {
3463     if (6U != dim_)
3464       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::operator(): wrong # of args (not rank 6 array)");
3465     return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3 * strides_[3] + i4 * strides_[4] + i5];
3466   }
3467 
3468   template <typename Numeric, unsigned Len, unsigned Dim>
3469   inline Numeric& ArrayND<Numeric, Len, Dim>::operator()(const unsigned i0,
3470                                                          const unsigned i1,
3471                                                          const unsigned i2,
3472                                                          const unsigned i3,
3473                                                          const unsigned i4,
3474                                                          const unsigned i5,
3475                                                          const unsigned i6) {
3476     if (7U != dim_)
3477       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::operator(): wrong # of args (not rank 7 array)");
3478     return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3 * strides_[3] + i4 * strides_[4] +
3479                  i5 * strides_[5] + i6];
3480   }
3481 
3482   template <typename Numeric, unsigned Len, unsigned Dim>
3483   inline Numeric& ArrayND<Numeric, Len, Dim>::operator()(const unsigned i0,
3484                                                          const unsigned i1,
3485                                                          const unsigned i2,
3486                                                          const unsigned i3,
3487                                                          const unsigned i4,
3488                                                          const unsigned i5,
3489                                                          const unsigned i6,
3490                                                          const unsigned i7) {
3491     if (8U != dim_)
3492       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::operator(): wrong # of args (not rank 8 array)");
3493     return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3 * strides_[3] + i4 * strides_[4] +
3494                  i5 * strides_[5] + i6 * strides_[6] + i7];
3495   }
3496 
3497   template <typename Numeric, unsigned Len, unsigned Dim>
3498   inline Numeric& ArrayND<Numeric, Len, Dim>::operator()(const unsigned i0,
3499                                                          const unsigned i1,
3500                                                          const unsigned i2,
3501                                                          const unsigned i3,
3502                                                          const unsigned i4,
3503                                                          const unsigned i5,
3504                                                          const unsigned i6,
3505                                                          const unsigned i7,
3506                                                          const unsigned i8) {
3507     if (9U != dim_)
3508       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::operator(): wrong # of args (not rank 9 array)");
3509     return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3 * strides_[3] + i4 * strides_[4] +
3510                  i5 * strides_[5] + i6 * strides_[6] + i7 * strides_[7] + i8];
3511   }
3512 
3513   template <typename Numeric, unsigned Len, unsigned Dim>
3514   inline Numeric& ArrayND<Numeric, Len, Dim>::operator()(const unsigned i0,
3515                                                          const unsigned i1,
3516                                                          const unsigned i2,
3517                                                          const unsigned i3,
3518                                                          const unsigned i4,
3519                                                          const unsigned i5,
3520                                                          const unsigned i6,
3521                                                          const unsigned i7,
3522                                                          const unsigned i8,
3523                                                          const unsigned i9) {
3524     if (10U != dim_)
3525       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::operator(): wrong # of args (not rank 10 array)");
3526     return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3 * strides_[3] + i4 * strides_[4] +
3527                  i5 * strides_[5] + i6 * strides_[6] + i7 * strides_[7] + i8 * strides_[8] + i9];
3528   }
3529 
3530   template <typename Numeric, unsigned Len, unsigned Dim>
3531   const Numeric& ArrayND<Numeric, Len, Dim>::at(const unsigned i0, const unsigned i1, const unsigned i2) const {
3532     if (3U != dim_)
3533       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::at: wrong # of args (not rank 3 array)");
3534     if (i0 >= shape_[0])
3535       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 0 out of range (rank 3)");
3536     if (i1 >= shape_[1])
3537       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 1 out of range (rank 3)");
3538     if (i2 >= shape_[2])
3539       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 2 out of range (rank 3)");
3540     return data_[i0 * strides_[0] + i1 * strides_[1] + i2];
3541   }
3542 
3543   template <typename Numeric, unsigned Len, unsigned Dim>
3544   const Numeric& ArrayND<Numeric, Len, Dim>::at(const unsigned i0,
3545                                                 const unsigned i1,
3546                                                 const unsigned i2,
3547                                                 const unsigned i3) const {
3548     if (4U != dim_)
3549       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::at: wrong # of args (not rank 4 array)");
3550     if (i0 >= shape_[0])
3551       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 0 out of range (rank 4)");
3552     if (i1 >= shape_[1])
3553       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 1 out of range (rank 4)");
3554     if (i2 >= shape_[2])
3555       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 2 out of range (rank 4)");
3556     if (i3 >= shape_[3])
3557       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 3 out of range (rank 4)");
3558     return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3];
3559   }
3560 
3561   template <typename Numeric, unsigned Len, unsigned Dim>
3562   const Numeric& ArrayND<Numeric, Len, Dim>::at(
3563       const unsigned i0, const unsigned i1, const unsigned i2, const unsigned i3, const unsigned i4) const {
3564     if (5U != dim_)
3565       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::at: wrong # of args (not rank 5 array)");
3566     if (i0 >= shape_[0])
3567       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 0 out of range (rank 5)");
3568     if (i1 >= shape_[1])
3569       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 1 out of range (rank 5)");
3570     if (i2 >= shape_[2])
3571       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 2 out of range (rank 5)");
3572     if (i3 >= shape_[3])
3573       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 3 out of range (rank 5)");
3574     if (i4 >= shape_[4])
3575       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 4 out of range (rank 5)");
3576     return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3 * strides_[3] + i4];
3577   }
3578 
3579   template <typename Numeric, unsigned Len, unsigned Dim>
3580   const Numeric& ArrayND<Numeric, Len, Dim>::at(const unsigned i0,
3581                                                 const unsigned i1,
3582                                                 const unsigned i2,
3583                                                 const unsigned i3,
3584                                                 const unsigned i4,
3585                                                 const unsigned i5) const {
3586     if (6U != dim_)
3587       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::at: wrong # of args (not rank 6 array)");
3588     if (i0 >= shape_[0])
3589       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 0 out of range (rank 6)");
3590     if (i1 >= shape_[1])
3591       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 1 out of range (rank 6)");
3592     if (i2 >= shape_[2])
3593       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 2 out of range (rank 6)");
3594     if (i3 >= shape_[3])
3595       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 3 out of range (rank 6)");
3596     if (i4 >= shape_[4])
3597       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 4 out of range (rank 6)");
3598     if (i5 >= shape_[5])
3599       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 5 out of range (rank 6)");
3600     return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3 * strides_[3] + i4 * strides_[4] + i5];
3601   }
3602 
3603   template <typename Numeric, unsigned Len, unsigned Dim>
3604   const Numeric& ArrayND<Numeric, Len, Dim>::at(const unsigned i0,
3605                                                 const unsigned i1,
3606                                                 const unsigned i2,
3607                                                 const unsigned i3,
3608                                                 const unsigned i4,
3609                                                 const unsigned i5,
3610                                                 const unsigned i6) const {
3611     if (7U != dim_)
3612       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::at: wrong # of args (not rank 7 array)");
3613     if (i0 >= shape_[0])
3614       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 0 out of range (rank 7)");
3615     if (i1 >= shape_[1])
3616       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 1 out of range (rank 7)");
3617     if (i2 >= shape_[2])
3618       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 2 out of range (rank 7)");
3619     if (i3 >= shape_[3])
3620       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 3 out of range (rank 7)");
3621     if (i4 >= shape_[4])
3622       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 4 out of range (rank 7)");
3623     if (i5 >= shape_[5])
3624       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 5 out of range (rank 7)");
3625     if (i6 >= shape_[6])
3626       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 6 out of range (rank 7)");
3627     return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3 * strides_[3] + i4 * strides_[4] +
3628                  i5 * strides_[5] + i6];
3629   }
3630 
3631   template <typename Numeric, unsigned Len, unsigned Dim>
3632   const Numeric& ArrayND<Numeric, Len, Dim>::at(const unsigned i0,
3633                                                 const unsigned i1,
3634                                                 const unsigned i2,
3635                                                 const unsigned i3,
3636                                                 const unsigned i4,
3637                                                 const unsigned i5,
3638                                                 const unsigned i6,
3639                                                 const unsigned i7) const {
3640     if (8U != dim_)
3641       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::at: wrong # of args (not rank 8 array)");
3642     if (i0 >= shape_[0])
3643       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 0 out of range (rank 8)");
3644     if (i1 >= shape_[1])
3645       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 1 out of range (rank 8)");
3646     if (i2 >= shape_[2])
3647       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 2 out of range (rank 8)");
3648     if (i3 >= shape_[3])
3649       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 3 out of range (rank 8)");
3650     if (i4 >= shape_[4])
3651       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 4 out of range (rank 8)");
3652     if (i5 >= shape_[5])
3653       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 5 out of range (rank 8)");
3654     if (i6 >= shape_[6])
3655       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 6 out of range (rank 8)");
3656     if (i7 >= shape_[7])
3657       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 7 out of range (rank 8)");
3658     return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3 * strides_[3] + i4 * strides_[4] +
3659                  i5 * strides_[5] + i6 * strides_[6] + i7];
3660   }
3661 
3662   template <typename Numeric, unsigned Len, unsigned Dim>
3663   const Numeric& ArrayND<Numeric, Len, Dim>::at(const unsigned i0,
3664                                                 const unsigned i1,
3665                                                 const unsigned i2,
3666                                                 const unsigned i3,
3667                                                 const unsigned i4,
3668                                                 const unsigned i5,
3669                                                 const unsigned i6,
3670                                                 const unsigned i7,
3671                                                 const unsigned i8) const {
3672     if (9U != dim_)
3673       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::at: wrong # of args (not rank 9 array)");
3674     if (i0 >= shape_[0])
3675       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 0 out of range (rank 9)");
3676     if (i1 >= shape_[1])
3677       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 1 out of range (rank 9)");
3678     if (i2 >= shape_[2])
3679       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 2 out of range (rank 9)");
3680     if (i3 >= shape_[3])
3681       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 3 out of range (rank 9)");
3682     if (i4 >= shape_[4])
3683       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 4 out of range (rank 9)");
3684     if (i5 >= shape_[5])
3685       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 5 out of range (rank 9)");
3686     if (i6 >= shape_[6])
3687       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 6 out of range (rank 9)");
3688     if (i7 >= shape_[7])
3689       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 7 out of range (rank 9)");
3690     if (i8 >= shape_[8])
3691       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 8 out of range (rank 9)");
3692     return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3 * strides_[3] + i4 * strides_[4] +
3693                  i5 * strides_[5] + i6 * strides_[6] + i7 * strides_[7] + i8];
3694   }
3695 
3696   template <typename Numeric, unsigned Len, unsigned Dim>
3697   const Numeric& ArrayND<Numeric, Len, Dim>::at(const unsigned i0,
3698                                                 const unsigned i1,
3699                                                 const unsigned i2,
3700                                                 const unsigned i3,
3701                                                 const unsigned i4,
3702                                                 const unsigned i5,
3703                                                 const unsigned i6,
3704                                                 const unsigned i7,
3705                                                 const unsigned i8,
3706                                                 const unsigned i9) const {
3707     if (10U != dim_)
3708       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::at: wrong # of args (not rank 10 array)");
3709     if (i0 >= shape_[0])
3710       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 0 out of range (rank 10)");
3711     if (i1 >= shape_[1])
3712       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 1 out of range (rank 10)");
3713     if (i2 >= shape_[2])
3714       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 2 out of range (rank 10)");
3715     if (i3 >= shape_[3])
3716       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 3 out of range (rank 10)");
3717     if (i4 >= shape_[4])
3718       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 4 out of range (rank 10)");
3719     if (i5 >= shape_[5])
3720       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 5 out of range (rank 10)");
3721     if (i6 >= shape_[6])
3722       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 6 out of range (rank 10)");
3723     if (i7 >= shape_[7])
3724       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 7 out of range (rank 10)");
3725     if (i8 >= shape_[8])
3726       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 8 out of range (rank 10)");
3727     if (i9 >= shape_[9])
3728       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 9 out of range (rank 10)");
3729     return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3 * strides_[3] + i4 * strides_[4] +
3730                  i5 * strides_[5] + i6 * strides_[6] + i7 * strides_[7] + i8 * strides_[8] + i9];
3731   }
3732 
3733   template <typename Numeric, unsigned Len, unsigned Dim>
3734   Numeric& ArrayND<Numeric, Len, Dim>::at(const unsigned i0, const unsigned i1, const unsigned i2) {
3735     if (3U != dim_)
3736       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::at: wrong # of args (not rank 3 array)");
3737     if (i0 >= shape_[0])
3738       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 0 out of range (rank 3)");
3739     if (i1 >= shape_[1])
3740       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 1 out of range (rank 3)");
3741     if (i2 >= shape_[2])
3742       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 2 out of range (rank 3)");
3743     return data_[i0 * strides_[0] + i1 * strides_[1] + i2];
3744   }
3745 
3746   template <typename Numeric, unsigned Len, unsigned Dim>
3747   Numeric& ArrayND<Numeric, Len, Dim>::at(const unsigned i0, const unsigned i1, const unsigned i2, const unsigned i3) {
3748     if (4U != dim_)
3749       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::at: wrong # of args (not rank 4 array)");
3750     if (i0 >= shape_[0])
3751       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 0 out of range (rank 4)");
3752     if (i1 >= shape_[1])
3753       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 1 out of range (rank 4)");
3754     if (i2 >= shape_[2])
3755       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 2 out of range (rank 4)");
3756     if (i3 >= shape_[3])
3757       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 3 out of range (rank 4)");
3758     return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3];
3759   }
3760 
3761   template <typename Numeric, unsigned Len, unsigned Dim>
3762   Numeric& ArrayND<Numeric, Len, Dim>::at(
3763       const unsigned i0, const unsigned i1, const unsigned i2, const unsigned i3, const unsigned i4) {
3764     if (5U != dim_)
3765       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::at: wrong # of args (not rank 5 array)");
3766     if (i0 >= shape_[0])
3767       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 0 out of range (rank 5)");
3768     if (i1 >= shape_[1])
3769       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 1 out of range (rank 5)");
3770     if (i2 >= shape_[2])
3771       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 2 out of range (rank 5)");
3772     if (i3 >= shape_[3])
3773       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 3 out of range (rank 5)");
3774     if (i4 >= shape_[4])
3775       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 4 out of range (rank 5)");
3776     return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3 * strides_[3] + i4];
3777   }
3778 
3779   template <typename Numeric, unsigned Len, unsigned Dim>
3780   Numeric& ArrayND<Numeric, Len, Dim>::at(const unsigned i0,
3781                                           const unsigned i1,
3782                                           const unsigned i2,
3783                                           const unsigned i3,
3784                                           const unsigned i4,
3785                                           const unsigned i5) {
3786     if (6U != dim_)
3787       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::at: wrong # of args (not rank 6 array)");
3788     if (i0 >= shape_[0])
3789       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 0 out of range (rank 6)");
3790     if (i1 >= shape_[1])
3791       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 1 out of range (rank 6)");
3792     if (i2 >= shape_[2])
3793       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 2 out of range (rank 6)");
3794     if (i3 >= shape_[3])
3795       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 3 out of range (rank 6)");
3796     if (i4 >= shape_[4])
3797       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 4 out of range (rank 6)");
3798     if (i5 >= shape_[5])
3799       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 5 out of range (rank 6)");
3800     return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3 * strides_[3] + i4 * strides_[4] + i5];
3801   }
3802 
3803   template <typename Numeric, unsigned Len, unsigned Dim>
3804   Numeric& ArrayND<Numeric, Len, Dim>::at(const unsigned i0,
3805                                           const unsigned i1,
3806                                           const unsigned i2,
3807                                           const unsigned i3,
3808                                           const unsigned i4,
3809                                           const unsigned i5,
3810                                           const unsigned i6) {
3811     if (7U != dim_)
3812       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::at: wrong # of args (not rank 7 array)");
3813     if (i0 >= shape_[0])
3814       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 0 out of range (rank 7)");
3815     if (i1 >= shape_[1])
3816       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 1 out of range (rank 7)");
3817     if (i2 >= shape_[2])
3818       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 2 out of range (rank 7)");
3819     if (i3 >= shape_[3])
3820       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 3 out of range (rank 7)");
3821     if (i4 >= shape_[4])
3822       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 4 out of range (rank 7)");
3823     if (i5 >= shape_[5])
3824       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 5 out of range (rank 7)");
3825     if (i6 >= shape_[6])
3826       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 6 out of range (rank 7)");
3827     return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3 * strides_[3] + i4 * strides_[4] +
3828                  i5 * strides_[5] + i6];
3829   }
3830 
3831   template <typename Numeric, unsigned Len, unsigned Dim>
3832   Numeric& ArrayND<Numeric, Len, Dim>::at(const unsigned i0,
3833                                           const unsigned i1,
3834                                           const unsigned i2,
3835                                           const unsigned i3,
3836                                           const unsigned i4,
3837                                           const unsigned i5,
3838                                           const unsigned i6,
3839                                           const unsigned i7) {
3840     if (8U != dim_)
3841       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::at: wrong # of args (not rank 8 array)");
3842     if (i0 >= shape_[0])
3843       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 0 out of range (rank 8)");
3844     if (i1 >= shape_[1])
3845       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 1 out of range (rank 8)");
3846     if (i2 >= shape_[2])
3847       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 2 out of range (rank 8)");
3848     if (i3 >= shape_[3])
3849       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 3 out of range (rank 8)");
3850     if (i4 >= shape_[4])
3851       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 4 out of range (rank 8)");
3852     if (i5 >= shape_[5])
3853       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 5 out of range (rank 8)");
3854     if (i6 >= shape_[6])
3855       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 6 out of range (rank 8)");
3856     if (i7 >= shape_[7])
3857       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 7 out of range (rank 8)");
3858     return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3 * strides_[3] + i4 * strides_[4] +
3859                  i5 * strides_[5] + i6 * strides_[6] + i7];
3860   }
3861 
3862   template <typename Numeric, unsigned Len, unsigned Dim>
3863   Numeric& ArrayND<Numeric, Len, Dim>::at(const unsigned i0,
3864                                           const unsigned i1,
3865                                           const unsigned i2,
3866                                           const unsigned i3,
3867                                           const unsigned i4,
3868                                           const unsigned i5,
3869                                           const unsigned i6,
3870                                           const unsigned i7,
3871                                           const unsigned i8) {
3872     if (9U != dim_)
3873       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::at: wrong # of args (not rank 9 array)");
3874     if (i0 >= shape_[0])
3875       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 0 out of range (rank 9)");
3876     if (i1 >= shape_[1])
3877       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 1 out of range (rank 9)");
3878     if (i2 >= shape_[2])
3879       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 2 out of range (rank 9)");
3880     if (i3 >= shape_[3])
3881       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 3 out of range (rank 9)");
3882     if (i4 >= shape_[4])
3883       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 4 out of range (rank 9)");
3884     if (i5 >= shape_[5])
3885       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 5 out of range (rank 9)");
3886     if (i6 >= shape_[6])
3887       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 6 out of range (rank 9)");
3888     if (i7 >= shape_[7])
3889       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 7 out of range (rank 9)");
3890     if (i8 >= shape_[8])
3891       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 8 out of range (rank 9)");
3892     return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3 * strides_[3] + i4 * strides_[4] +
3893                  i5 * strides_[5] + i6 * strides_[6] + i7 * strides_[7] + i8];
3894   }
3895 
3896   template <typename Numeric, unsigned Len, unsigned Dim>
3897   Numeric& ArrayND<Numeric, Len, Dim>::at(const unsigned i0,
3898                                           const unsigned i1,
3899                                           const unsigned i2,
3900                                           const unsigned i3,
3901                                           const unsigned i4,
3902                                           const unsigned i5,
3903                                           const unsigned i6,
3904                                           const unsigned i7,
3905                                           const unsigned i8,
3906                                           const unsigned i9) {
3907     if (10U != dim_)
3908       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::at: wrong # of args (not rank 10 array)");
3909     if (i0 >= shape_[0])
3910       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 0 out of range (rank 10)");
3911     if (i1 >= shape_[1])
3912       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 1 out of range (rank 10)");
3913     if (i2 >= shape_[2])
3914       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 2 out of range (rank 10)");
3915     if (i3 >= shape_[3])
3916       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 3 out of range (rank 10)");
3917     if (i4 >= shape_[4])
3918       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 4 out of range (rank 10)");
3919     if (i5 >= shape_[5])
3920       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 5 out of range (rank 10)");
3921     if (i6 >= shape_[6])
3922       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 6 out of range (rank 10)");
3923     if (i7 >= shape_[7])
3924       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 7 out of range (rank 10)");
3925     if (i8 >= shape_[8])
3926       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 8 out of range (rank 10)");
3927     if (i9 >= shape_[9])
3928       throw npstat::NpstatOutOfRange("In npstat::ArrayND::at: index 9 out of range (rank 10)");
3929     return data_[i0 * strides_[0] + i1 * strides_[1] + i2 * strides_[2] + i3 * strides_[3] + i4 * strides_[4] +
3930                  i5 * strides_[5] + i6 * strides_[6] + i7 * strides_[7] + i8 * strides_[8] + i9];
3931   }
3932 
3933   template <typename Numeric, unsigned Len, unsigned Dim>
3934   template <unsigned Len2, unsigned Dim2>
3935   double ArrayND<Numeric, Len, Dim>::maxAbsDifference(const ArrayND<Numeric, Len2, Dim2>& r) const {
3936     if (!isShapeCompatible(r))
3937       throw npstat::NpstatInvalidArgument(
3938           "In npstat::ArrayND::maxAbsDifference: "
3939           "incompatible argument array shape");
3940     if (dim_) {
3941       double maxd = 0.0;
3942       for (unsigned long i = 0; i < len_; ++i) {
3943         const Numeric rval = r.data_[i];
3944         const double d = absDifference(data_[i], rval);
3945         if (d > maxd)
3946           maxd = d;
3947       }
3948       return maxd;
3949     } else {
3950       const Numeric rval = r.localData_[0];
3951       return absDifference(localData_[0], rval);
3952     }
3953   }
3954 
3955   template <typename Numeric, unsigned Len, unsigned Dim>
3956   template <unsigned Len2, unsigned Dim2>
3957   bool ArrayND<Numeric, Len, Dim>::operator==(const ArrayND<Numeric, Len2, Dim2>& r) const {
3958     if (shapeIsKnown_ != r.shapeIsKnown_)
3959       return false;
3960     if (r.dim_ != dim_)
3961       return false;
3962     if (r.len_ != len_)
3963       return false;
3964     for (unsigned i = 0; i < dim_; ++i)
3965       if (shape_[i] != r.shape_[i])
3966         return false;
3967     for (unsigned i = 0; i < dim_; ++i)
3968       assert(strides_[i] == r.strides_[i]);
3969     for (unsigned long j = 0; j < len_; ++j)
3970       if (data_[j] != r.data_[j])
3971         return false;
3972     return true;
3973   }
3974 
3975   template <typename Numeric, unsigned Len, unsigned Dim>
3976   template <unsigned Len2, unsigned Dim2>
3977   inline bool ArrayND<Numeric, Len, Dim>::operator!=(const ArrayND<Numeric, Len2, Dim2>& r) const {
3978     return !(*this == r);
3979   }
3980 
3981   template <typename Numeric, unsigned Len, unsigned Dim>
3982   template <typename Num2>
3983   ArrayND<Numeric, Len, Dim> ArrayND<Numeric, Len, Dim>::operator*(const Num2& r) const {
3984     if (!shapeIsKnown_)
3985       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"operator*\"");
3986     ArrayND<Numeric, Len, Dim> result(shape_, dim_);
3987     for (unsigned long i = 0; i < len_; ++i)
3988       result.data_[i] = data_[i] * r;
3989     return result;
3990   }
3991 
3992   template <typename Numeric, unsigned Len, unsigned Dim>
3993   template <typename Num2>
3994   ArrayND<Numeric, Len, Dim> ArrayND<Numeric, Len, Dim>::operator/(const Num2& r) const {
3995     if (!shapeIsKnown_)
3996       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"operator/\"");
3997     if (r == Num2())
3998       throw npstat::NpstatRuntimeError("In npstat::ArrayND::operator/: division by zero");
3999     ArrayND<Numeric, Len, Dim> result(shape_, dim_);
4000     for (unsigned long i = 0; i < len_; ++i)
4001       result.data_[i] = data_[i] / r;
4002     return result;
4003   }
4004 
4005   template <typename Numeric, unsigned Len, unsigned Dim>
4006   template <unsigned Len2, unsigned Dim2>
4007   ArrayND<Numeric, Len, Dim> ArrayND<Numeric, Len, Dim>::operator+(const ArrayND<Numeric, Len2, Dim2>& r) const {
4008     if (!isShapeCompatible(r))
4009       throw npstat::NpstatInvalidArgument(
4010           "In npstat::ArrayND::operator+: "
4011           "incompatible argument array shape");
4012     ArrayND<Numeric, Len, Dim> result(shape_, dim_);
4013     for (unsigned long i = 0; i < len_; ++i)
4014       result.data_[i] = data_[i] + r.data_[i];
4015     return result;
4016   }
4017 
4018   template <typename Numeric, unsigned Len, unsigned Dim>
4019   template <unsigned Len2, unsigned Dim2>
4020   ArrayND<Numeric, Len, Dim> ArrayND<Numeric, Len, Dim>::operator-(const ArrayND<Numeric, Len2, Dim2>& r) const {
4021     if (!isShapeCompatible(r))
4022       throw npstat::NpstatInvalidArgument(
4023           "In npstat::ArrayND::operator-: "
4024           "incompatible argument array shape");
4025     ArrayND<Numeric, Len, Dim> result(shape_, dim_);
4026     for (unsigned long i = 0; i < len_; ++i)
4027       result.data_[i] = data_[i] - r.data_[i];
4028     return result;
4029   }
4030 
4031   template <typename Numeric, unsigned Len, unsigned Dim>
4032   inline ArrayND<Numeric, Len, Dim> ArrayND<Numeric, Len, Dim>::operator+() const {
4033     if (!shapeIsKnown_)
4034       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"operator+\"");
4035     return *this;
4036   }
4037 
4038   template <typename Numeric, unsigned Len, unsigned Dim>
4039   ArrayND<Numeric, Len, Dim> ArrayND<Numeric, Len, Dim>::operator-() const {
4040     if (!shapeIsKnown_)
4041       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"operator-\"");
4042     ArrayND<Numeric, Len, Dim> result(shape_, dim_);
4043     for (unsigned long i = 0; i < len_; ++i)
4044       result.data_[i] = -data_[i];
4045     return result;
4046   }
4047 
4048   template <typename Numeric, unsigned Len, unsigned Dim>
4049   template <typename Num2>
4050   ArrayND<Numeric, Len, Dim>& ArrayND<Numeric, Len, Dim>::operator*=(const Num2& r) {
4051     if (!shapeIsKnown_)
4052       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"operator*=\"");
4053     for (unsigned long i = 0; i < len_; ++i)
4054       data_[i] *= r;
4055     return *this;
4056   }
4057 
4058   template <typename Numeric, unsigned Len, unsigned Dim>
4059   ArrayND<Numeric, Len, Dim>& ArrayND<Numeric, Len, Dim>::makeNonNegative() {
4060     if (!shapeIsKnown_)
4061       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"makeNonNegative\"");
4062     const Numeric zero = Numeric();
4063     if (dim_) {
4064       for (unsigned long i = 0; i < len_; ++i)
4065         if (!(ComplexComparesAbs<Numeric>::more(data_[i], zero)))
4066           data_[i] = zero;
4067     } else if (!(ComplexComparesAbs<Numeric>::more(localData_[0], zero)))
4068       localData_[0] = zero;
4069     return *this;
4070   }
4071 
4072   template <typename Numeric, unsigned Len, unsigned Dim>
4073   unsigned ArrayND<Numeric, Len, Dim>::makeCopulaSteps(const double tolerance, const unsigned nCycles) {
4074     if (!shapeIsKnown_)
4075       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"makeCopulaSteps\"");
4076     if (nCycles == 0U)
4077       return 0U;
4078     if (!dim_)
4079       throw npstat::NpstatInvalidArgument(
4080           "npstat::ArrayND::makeCopulaSteps method "
4081           "can not be used with array of 0 rank");
4082 
4083     const Numeric zero = Numeric();
4084     for (unsigned long i = 0; i < len_; ++i)
4085       if (!(ComplexComparesAbs<Numeric>::more(data_[i], zero)))
4086         data_[i] = zero;
4087 
4088     std::vector<Numeric*> axesPtrBuf(dim_);
4089     Numeric** axes = &axesPtrBuf[0];
4090     const Numeric one = static_cast<Numeric>(1);
4091 
4092     // Memory for the axis accumulators
4093     unsigned idxSum = 0;
4094     for (unsigned i = 0; i < dim_; ++i)
4095       idxSum += shape_[i];
4096     std::vector<Numeric> axesBuf(idxSum);
4097     axes[0] = &axesBuf[0];
4098     for (unsigned i = 1; i < dim_; ++i)
4099       axes[i] = axes[i - 1] + shape_[i - 1];
4100 
4101     // Accumulate axis projections
4102     unsigned icycle = 0;
4103     for (; icycle < nCycles; ++icycle) {
4104       for (unsigned i = 0; i < idxSum; ++i)
4105         axesBuf[i] = zero;
4106 
4107       // Accumulate sums for each axis
4108       for (unsigned long idat = 0; idat < len_; ++idat) {
4109         unsigned long l = idat;
4110         for (unsigned i = 0; i < dim_; ++i) {
4111           const unsigned idx = l / strides_[i];
4112           l -= (idx * strides_[i]);
4113           axes[i][idx] += data_[idat];
4114         }
4115       }
4116 
4117       // Make averages out of sums
4118       bool withinTolerance = true;
4119       Numeric totalSum = zero;
4120       for (unsigned i = 0; i < dim_; ++i) {
4121         Numeric axisSum = zero;
4122         const unsigned amax = shape_[i];
4123         for (unsigned a = 0; a < amax; ++a) {
4124           if (axes[i][a] == zero)
4125             throw npstat::NpstatRuntimeError(
4126                 "In npstat::ArrayND::makeCopulaSteps: "
4127                 "marginal density is zero");
4128           axisSum += axes[i][a];
4129         }
4130         totalSum += axisSum;
4131         const Numeric axisAverage = axisSum / static_cast<Numeric>(amax);
4132         for (unsigned a = 0; a < amax; ++a)
4133           axes[i][a] /= axisAverage;
4134         for (unsigned a = 0; a < amax && withinTolerance; ++a) {
4135           const double adelta = absDifference(axes[i][a], one);
4136           if (adelta > tolerance)
4137             withinTolerance = false;
4138         }
4139       }
4140 
4141       if (withinTolerance)
4142         break;
4143 
4144       const Numeric totalAverage = totalSum / static_cast<Numeric>(len_) / static_cast<Numeric>(dim_);
4145 
4146       // Run over all points again and divide by
4147       // the product of marginals
4148       for (unsigned long idat = 0; idat < len_; ++idat) {
4149         unsigned long l = idat;
4150         for (unsigned i = 0; i < dim_; ++i) {
4151           const unsigned idx = l / strides_[i];
4152           l -= (idx * strides_[i]);
4153           data_[idat] /= axes[i][idx];
4154         }
4155         data_[idat] /= totalAverage;
4156       }
4157     }
4158 
4159     return icycle;
4160   }
4161 
4162   template <typename Numeric, unsigned Len, unsigned Dim>
4163   template <typename Num2>
4164   ArrayND<Numeric, Len, Dim>& ArrayND<Numeric, Len, Dim>::operator/=(const Num2& r) {
4165     if (!shapeIsKnown_)
4166       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"operator/=\"");
4167     if (r == Num2())
4168       throw npstat::NpstatRuntimeError("In npstat::ArrayND::operator/=: division by zero");
4169     for (unsigned long i = 0; i < len_; ++i)
4170       data_[i] /= r;
4171     return *this;
4172   }
4173 
4174   template <typename Numeric, unsigned Len, unsigned Dim>
4175   template <typename Num2, unsigned Len2, unsigned Dim2>
4176   ArrayND<Numeric, Len, Dim>& ArrayND<Numeric, Len, Dim>::operator+=(const ArrayND<Num2, Len2, Dim2>& r) {
4177     if (!isShapeCompatible(r))
4178       throw npstat::NpstatInvalidArgument(
4179           "In npstat::ArrayND::operator+=: "
4180           "incompatible argument array shape");
4181     for (unsigned long i = 0; i < len_; ++i)
4182       data_[i] += r.data_[i];
4183     return *this;
4184   }
4185 
4186   template <typename Numeric, unsigned Len, unsigned Dim>
4187   template <typename Num3, typename Num2, unsigned Len2, unsigned Dim2>
4188   ArrayND<Numeric, Len, Dim>& ArrayND<Numeric, Len, Dim>::addmul(const ArrayND<Num2, Len2, Dim2>& r, const Num3& c) {
4189     if (!isShapeCompatible(r))
4190       throw npstat::NpstatInvalidArgument(
4191           "In npstat::ArrayND::addmul: "
4192           "incompatible argument array shape");
4193     for (unsigned long i = 0; i < len_; ++i)
4194       data_[i] += r.data_[i] * c;
4195     return *this;
4196   }
4197 
4198   template <typename Numeric, unsigned Len, unsigned Dim>
4199   template <typename Num2, unsigned Len2, unsigned Dim2>
4200   ArrayND<Numeric, Len, Dim>& ArrayND<Numeric, Len, Dim>::operator-=(const ArrayND<Num2, Len2, Dim2>& r) {
4201     if (!isShapeCompatible(r))
4202       throw npstat::NpstatInvalidArgument(
4203           "In npstat::ArrayND::operator-=: "
4204           "incompatible argument array shape");
4205     for (unsigned long i = 0; i < len_; ++i)
4206       data_[i] -= r.data_[i];
4207     return *this;
4208   }
4209 
4210   template <typename Numeric, unsigned Len, unsigned Dim>
4211   Numeric ArrayND<Numeric, Len, Dim>::interpolate1(const double* coords, const unsigned dim) const {
4212     if (!shapeIsKnown_)
4213       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"interpolate1\"");
4214     if (dim != dim_)
4215       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::interpolate1: incompatible coordinate length");
4216     if (dim) {
4217       const unsigned maxdim = CHAR_BIT * sizeof(unsigned long);
4218       if (dim_ >= maxdim)
4219         throw npstat::NpstatInvalidArgument("In npstat::ArrayND::interpolate1: array rank is too large");
4220 
4221       double dx[maxdim];
4222       unsigned ix[maxdim];
4223       for (unsigned i = 0; i < dim; ++i) {
4224         const double x = coords[i];
4225         if (x <= 0.0) {
4226           ix[i] = 0;
4227           dx[i] = 0.0;
4228         } else if (x >= static_cast<double>(shape_[i] - 1)) {
4229           ix[i] = shape_[i] - 1;
4230           dx[i] = 0.0;
4231         } else {
4232           ix[i] = static_cast<unsigned>(std::floor(x));
4233           dx[i] = x - ix[i];
4234         }
4235       }
4236 
4237       Numeric sum = Numeric();
4238       const unsigned long maxcycle = 1UL << dim;
4239       for (unsigned long icycle = 0UL; icycle < maxcycle; ++icycle) {
4240         double w = 1.0;
4241         unsigned long icell = 0UL;
4242         for (unsigned i = 0; i < dim; ++i) {
4243           if (icycle & (1UL << i)) {
4244             w *= dx[i];
4245             icell += strides_[i] * (ix[i] + 1U);
4246           } else {
4247             w *= (1.0 - dx[i]);
4248             icell += strides_[i] * ix[i];
4249           }
4250         }
4251         if (w > 0.0)
4252           sum += data_[icell] * static_cast<proper_double>(w);
4253       }
4254       return sum;
4255     } else
4256       return localData_[0];
4257   }
4258 
4259   template <typename Numeric, unsigned Len, unsigned Dim>
4260   Numeric ArrayND<Numeric, Len, Dim>::interpolateLoop(const unsigned level,
4261                                                       const double* coords,
4262                                                       const Numeric* base) const {
4263     const unsigned npoints = shape_[level];
4264     const double x = coords[level];
4265 
4266     unsigned ix, npt = 1;
4267     double dx = 0.0;
4268     if (x < 0.0)
4269       ix = 0;
4270     else if (x > static_cast<double>(npoints - 1))
4271       ix = npoints - 1;
4272     else {
4273       ix = static_cast<unsigned>(std::floor(x));
4274       if (ix)
4275         --ix;
4276       unsigned imax = ix + 3;
4277       while (imax >= npoints) {
4278         if (ix)
4279           --ix;
4280         --imax;
4281       }
4282       dx = x - ix;
4283       npt = imax + 1 - ix;
4284     }
4285     assert(npt >= 1 && npt <= 4);
4286 
4287     Numeric fit[4];
4288     if (level < dim_ - 1)
4289       for (unsigned ipt = 0; ipt < npt; ++ipt)
4290         fit[ipt] = interpolateLoop(level + 1, coords, base + (ix + ipt) * strides_[level]);
4291 
4292     const Numeric* const v = (level == dim_ - 1 ? base + ix : fit);
4293     switch (npt) {
4294       case 1:
4295         return v[0];
4296       case 2:
4297         return interpolate_linear(dx, v[0], v[1]);
4298       case 3:
4299         return interpolate_quadratic(dx, v[0], v[1], v[2]);
4300       case 4:
4301         return interpolate_cubic(dx, v[0], v[1], v[2], v[3]);
4302       default:
4303         assert(0);
4304         return Numeric();
4305     }
4306   }
4307 
4308   template <typename Numeric, unsigned Len, unsigned Dim>
4309   inline Numeric ArrayND<Numeric, Len, Dim>::interpolate3(const double* coords, const unsigned dim) const {
4310     if (!shapeIsKnown_)
4311       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"interpolate3\"");
4312     if (dim != dim_)
4313       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::interpolate3: incompatible coordinate length");
4314     if (dim) {
4315       assert(coords);
4316       return interpolateLoop(0, coords, data_);
4317     } else
4318       return localData_[0];
4319   }
4320 
4321   template <typename Numeric, unsigned Len, unsigned Dim>
4322   template <class Functor>
4323   ArrayND<Numeric, Len, Dim>& ArrayND<Numeric, Len, Dim>::apply(Functor f) {
4324     if (!shapeIsKnown_)
4325       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"apply\"");
4326     for (unsigned long i = 0; i < len_; ++i)
4327       data_[i] = static_cast<Numeric>(f(data_[i]));
4328     return *this;
4329   }
4330 
4331   template <typename Numeric, unsigned Len, unsigned Dim>
4332   template <class Functor>
4333   ArrayND<Numeric, Len, Dim>& ArrayND<Numeric, Len, Dim>::scanInPlace(Functor f) {
4334     if (!shapeIsKnown_)
4335       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"scanInPlace\"");
4336     for (unsigned long i = 0; i < len_; ++i)
4337       f(data_[i]);
4338     return *this;
4339   }
4340 
4341   template <typename Numeric, unsigned Len, unsigned Dim>
4342   ArrayND<Numeric, Len, Dim>& ArrayND<Numeric, Len, Dim>::constFill(const Numeric c) {
4343     if (!shapeIsKnown_)
4344       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"constFill\"");
4345     for (unsigned long i = 0; i < len_; ++i)
4346       data_[i] = c;
4347     return *this;
4348   }
4349 
4350   template <typename Numeric, unsigned Len, unsigned Dim>
4351   inline ArrayND<Numeric, Len, Dim>& ArrayND<Numeric, Len, Dim>::clear() {
4352     return constFill(Numeric());
4353   }
4354 
4355   template <typename Numeric, unsigned Len, unsigned Dim>
4356   ArrayND<Numeric, Len, Dim>& ArrayND<Numeric, Len, Dim>::uninitialize() {
4357     destroyBuffer(data_, localData_);
4358     destroyBuffer(strides_, localStrides_);
4359     destroyBuffer(shape_, localShape_);
4360     localData_[0] = Numeric();
4361     data_ = localData_;
4362     strides_ = nullptr;
4363     shape_ = nullptr;
4364     len_ = 0;
4365     dim_ = 0;
4366     shapeIsKnown_ = false;
4367     return *this;
4368   }
4369 
4370   template <typename Numeric, unsigned Len, unsigned Dim>
4371   ArrayND<Numeric, Len, Dim>& ArrayND<Numeric, Len, Dim>::makeUnit() {
4372     if (!shapeIsKnown_)
4373       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"makeUnit\"");
4374     if (dim_ < 2)
4375       throw npstat::NpstatInvalidArgument(
4376           "npstat::ArrayND::makeUnit method "
4377           "can not be used with arrays of rank less than 2");
4378     constFill(Numeric());
4379     unsigned long stride = 0UL;
4380     const unsigned dimlen = shape_[0];
4381     for (unsigned i = 0; i < dim_; ++i) {
4382       if (shape_[i] != dimlen)
4383         throw npstat::NpstatInvalidArgument(
4384             "npstat::ArrayND::makeUnit method needs "
4385             "the array span to be the same in ech dimension");
4386       stride += strides_[i];
4387     }
4388     const Numeric one(static_cast<Numeric>(1));
4389     for (unsigned i = 0; i < dimlen; ++i)
4390       data_[i * stride] = one;
4391     return *this;
4392   }
4393 
4394   template <typename Numeric, unsigned Len, unsigned Dim>
4395   void ArrayND<Numeric, Len, Dim>::linearFillLoop(
4396       const unsigned level, const double s0, const unsigned long idx, const double shift, const double* coeffs) {
4397     const unsigned imax = shape_[level];
4398     const double c = coeffs[level];
4399     if (level == dim_ - 1) {
4400       Numeric* d = &data_[idx];
4401       for (unsigned i = 0; i < imax; ++i) {
4402         // Note that we want to add "shift" only at the
4403         // very end. This might improve the numerical
4404         // precision of the result.
4405         const double sum = s0 + c * i + shift;
4406         d[i] = static_cast<Numeric>(sum);
4407       }
4408     } else {
4409       const unsigned long stride = strides_[level];
4410       for (unsigned i = 0; i < imax; ++i)
4411         linearFillLoop(level + 1, s0 + c * i, idx + i * stride, shift, coeffs);
4412     }
4413   }
4414 
4415   template <typename Numeric, unsigned Len, unsigned Dim>
4416   ArrayND<Numeric, Len, Dim>& ArrayND<Numeric, Len, Dim>::linearFill(const double* coeffs,
4417                                                                      const unsigned dimCoeffs,
4418                                                                      const double shift) {
4419     // Make sure the object has been initialized
4420     if (!shapeIsKnown_)
4421       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"linearFill\"");
4422     if (dim_ != dimCoeffs)
4423       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::linearFill: incompatible number of coefficients");
4424     if (dim_) {
4425       assert(coeffs);
4426       linearFillLoop(0U, 0.0, 0UL, shift, coeffs);
4427     } else
4428       localData_[0] = static_cast<Numeric>(shift);
4429     return *this;
4430   }
4431 
4432   template <typename Numeric, unsigned Len, unsigned Dim>
4433   template <class Functor>
4434   void ArrayND<Numeric, Len, Dim>::functorFillLoop(const unsigned level,
4435                                                    const unsigned long idx,
4436                                                    Functor f,
4437                                                    unsigned* farg) {
4438     const unsigned imax = shape_[level];
4439     if (level == dim_ - 1) {
4440       Numeric* d = &data_[idx];
4441       const unsigned* myarg = farg;
4442       for (unsigned i = 0; i < imax; ++i) {
4443         farg[level] = i;
4444         d[i] = static_cast<Numeric>(f(myarg, dim_));
4445       }
4446     } else {
4447       const unsigned long stride = strides_[level];
4448       for (unsigned i = 0; i < imax; ++i) {
4449         farg[level] = i;
4450         functorFillLoop(level + 1, idx + i * stride, f, farg);
4451       }
4452     }
4453   }
4454 
4455   template <typename Numeric, unsigned Len, unsigned Dim>
4456   template <typename Accumulator>
4457   void ArrayND<Numeric, Len, Dim>::convertToLastDimCdfLoop(
4458       ArrayND* sumSlice, const unsigned level, unsigned long idx0, unsigned long idxSlice, const bool useTrapezoids) {
4459     static const proper_double half = 0.5;
4460     const unsigned imax = shape_[level];
4461     if (level == dim_ - 1) {
4462       Accumulator acc = Accumulator();
4463       Numeric* data = data_ + idx0;
4464       if (useTrapezoids) {
4465         Numeric oldval = Numeric();
4466         for (unsigned i = 0; i < imax; ++i) {
4467           acc += (data[i] + oldval) * half;
4468           oldval = data[i];
4469           data[i] = static_cast<Numeric>(acc);
4470         }
4471         acc += oldval * half;
4472       } else
4473         for (unsigned i = 0; i < imax; ++i) {
4474           acc += data[i];
4475           data[i] = static_cast<Numeric>(acc);
4476         }
4477       if (sumSlice->dim_)
4478         sumSlice->data_[idxSlice] = static_cast<Numeric>(acc);
4479       else
4480         sumSlice->localData_[0] = static_cast<Numeric>(acc);
4481     } else {
4482       const unsigned long stride = strides_[level];
4483       unsigned long sumStride = 0UL;
4484       if (sumSlice->dim_)
4485         sumStride = sumSlice->strides_[level];
4486       for (unsigned i = 0; i < imax; ++i) {
4487         convertToLastDimCdfLoop<Accumulator>(sumSlice, level + 1, idx0, idxSlice, useTrapezoids);
4488         idx0 += stride;
4489         idxSlice += sumStride;
4490       }
4491     }
4492   }
4493 
4494   template <typename Numeric, unsigned Len, unsigned Dim>
4495   template <typename Accumulator>
4496   inline void ArrayND<Numeric, Len, Dim>::convertToLastDimCdf(ArrayND* sumSlice, const bool useTrapezoids) {
4497     if (!shapeIsKnown_)
4498       throw npstat::NpstatInvalidArgument(
4499           "Initialize npstat::ArrayND before calling "
4500           "method \"convertToLastDimCdf\"");
4501     if (!dim_)
4502       throw npstat::NpstatInvalidArgument(
4503           "npstat::ArrayND::convertToLastDimCdf method "
4504           "can not be used with array of 0 rank");
4505     assert(sumSlice);
4506     if (!sumSlice->shapeIsKnown_)
4507       throw npstat::NpstatInvalidArgument(
4508           "In npstat::ArrayND::convertToLastDimCdf: "
4509           "uninitialized argument array");
4510     convertToLastDimCdfLoop<Accumulator>(sumSlice, 0U, 0UL, 0UL, useTrapezoids);
4511   }
4512 
4513   template <typename Numeric, unsigned Len, unsigned Dim>
4514   template <class Functor>
4515   ArrayND<Numeric, Len, Dim>& ArrayND<Numeric, Len, Dim>::functorFill(Functor f) {
4516     if (!shapeIsKnown_)
4517       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"functorFill\"");
4518     if (dim_) {
4519       unsigned localIndex[Dim];
4520       unsigned* index = makeBuffer(dim_, localIndex, Dim);
4521       functorFillLoop(0U, 0UL, f, index);
4522       destroyBuffer(index, localIndex);
4523     } else
4524       localData_[0] = static_cast<Numeric>(f(static_cast<unsigned*>(nullptr), 0U));
4525     return *this;
4526   }
4527 
4528   template <typename Numeric, unsigned Len, unsigned Dim>
4529   template <typename Num2, unsigned Len2, unsigned Dim2>
4530   bool ArrayND<Numeric, Len, Dim>::isClose(const ArrayND<Num2, Len2, Dim2>& r, const double eps) const {
4531     if (eps < 0.0)
4532       throw npstat::NpstatDomainError("In npstat::ArrayND::isClose: tolerance must not be negative");
4533     if (!isShapeCompatible(r))
4534       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::isClose: incompatible argument array shape");
4535     if (dim_) {
4536       for (unsigned long i = 0; i < len_; ++i) {
4537         const Numeric rval = r.data_[i];
4538         if (static_cast<double>(absDifference(data_[i], rval)) > eps)
4539           return false;
4540       }
4541     } else {
4542       const Numeric rval = r.localData_[0];
4543       if (static_cast<double>(absDifference(localData_[0], rval)) > eps)
4544         return false;
4545     }
4546     return true;
4547   }
4548 
4549   template <typename Numeric, unsigned Len, unsigned Dim>
4550   template <typename Num2, unsigned Len2, unsigned Dim2>
4551   ArrayND<Numeric, Len, Dim> ArrayND<Numeric, Len, Dim>::outer(const ArrayND<Num2, Len2, Dim2>& r) const {
4552     return ArrayND<Numeric, Len, Dim>(*this, r);
4553   }
4554 
4555   template <typename Numeric, unsigned Len, unsigned Dim>
4556   void ArrayND<Numeric, Len, Dim>::contractLoop(unsigned thisLevel,
4557                                                 const unsigned resLevel,
4558                                                 const unsigned pos1,
4559                                                 const unsigned pos2,
4560                                                 unsigned long idxThis,
4561                                                 unsigned long idxRes,
4562                                                 ArrayND& result) const {
4563     while (thisLevel == pos1 || thisLevel == pos2)
4564       ++thisLevel;
4565     assert(thisLevel < dim_);
4566 
4567     if (resLevel == result.dim_ - 1) {
4568       const unsigned ncontract = shape_[pos1];
4569       const unsigned imax = result.shape_[resLevel];
4570       const unsigned long stride = strides_[pos1] + strides_[pos2];
4571       for (unsigned i = 0; i < imax; ++i) {
4572         const Numeric* tmp = data_ + (idxThis + i * strides_[thisLevel]);
4573         Numeric sum = Numeric();
4574         for (unsigned j = 0; j < ncontract; ++j)
4575           sum += tmp[j * stride];
4576         result.data_[idxRes + i] = sum;
4577       }
4578     } else {
4579       const unsigned imax = result.shape_[resLevel];
4580       assert(imax == shape_[thisLevel]);
4581       for (unsigned i = 0; i < imax; ++i) {
4582         contractLoop(thisLevel + 1, resLevel + 1, pos1, pos2, idxThis, idxRes, result);
4583         idxThis += strides_[thisLevel];
4584         idxRes += result.strides_[resLevel];
4585       }
4586     }
4587   }
4588 
4589   template <typename Numeric, unsigned Len, unsigned Dim>
4590   ArrayND<Numeric, Len, Dim> ArrayND<Numeric, Len, Dim>::contract(const unsigned pos1, const unsigned pos2) const {
4591     if (!shapeIsKnown_)
4592       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"contract\"");
4593     if (!(pos1 < dim_ && pos2 < dim_ && pos1 != pos2))
4594       throw npstat::NpstatInvalidArgument(
4595           "In npstat::ArrayND::contract: "
4596           "incompatible contraction indices");
4597     if (shape_[pos1] != shape_[pos2])
4598       throw npstat::NpstatInvalidArgument(
4599           "In npstat::ArrayND::contract: incompatible "
4600           "length of contracted dimensions");
4601 
4602     // Construct the new shape
4603     unsigned newshapeBuf[Dim];
4604     unsigned* newshape = makeBuffer(dim_ - 2, newshapeBuf, Dim);
4605     unsigned ishap = 0;
4606     for (unsigned i = 0; i < dim_; ++i)
4607       if (i != pos1 && i != pos2)
4608         newshape[ishap++] = shape_[i];
4609 
4610     // Form the result array
4611     ArrayND<Numeric, Len, Dim> result(newshape, ishap);
4612     if (ishap)
4613       contractLoop(0, 0, pos1, pos2, 0UL, 0UL, result);
4614     else {
4615       // We are just calculating the trace
4616       Numeric sum = Numeric();
4617       const unsigned imax = shape_[0];
4618       const unsigned long stride = strides_[0] + strides_[1];
4619       for (unsigned i = 0; i < imax; ++i)
4620         sum += data_[i * stride];
4621       result() = sum;
4622     }
4623 
4624     destroyBuffer(newshape, newshapeBuf);
4625     return result;
4626   }
4627 
4628   template <typename Numeric, unsigned Len, unsigned Dim>
4629   void ArrayND<Numeric, Len, Dim>::transposeLoop(const unsigned level,
4630                                                  const unsigned pos1,
4631                                                  const unsigned pos2,
4632                                                  unsigned long idxThis,
4633                                                  unsigned long idxRes,
4634                                                  ArrayND& result) const {
4635     const unsigned imax = shape_[level];
4636     const unsigned long mystride = strides_[level];
4637     const unsigned relevel = level == pos1 ? pos2 : (level == pos2 ? pos1 : level);
4638     const unsigned long restride = result.strides_[relevel];
4639     const bool ready = (level == dim_ - 1);
4640     for (unsigned i = 0; i < imax; ++i) {
4641       if (ready)
4642         result.data_[idxRes] = data_[idxThis];
4643       else
4644         transposeLoop(level + 1, pos1, pos2, idxThis, idxRes, result);
4645       idxThis += mystride;
4646       idxRes += restride;
4647     }
4648   }
4649 
4650   template <typename Numeric, unsigned Len, unsigned Dim>
4651   template <typename Num2>
4652   Num2 ArrayND<Numeric, Len, Dim>::sum() const {
4653     if (!shapeIsKnown_)
4654       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"sum\"");
4655     Num2 sum = Num2();
4656     for (unsigned long i = 0; i < len_; ++i)
4657       sum += data_[i];
4658     return sum;
4659   }
4660 
4661   template <typename Numeric, unsigned Len, unsigned Dim>
4662   template <typename Num2>
4663   Num2 ArrayND<Numeric, Len, Dim>::sumsq() const {
4664     if (!shapeIsKnown_)
4665       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"sumsq\"");
4666     Num2 sum = Num2();
4667     for (unsigned long i = 0; i < len_; ++i) {
4668       const Num2 absval = absValue(data_[i]);
4669       sum += absval * absval;
4670     }
4671     return sum;
4672   }
4673 
4674   template <typename Numeric, unsigned Len, unsigned Dim>
4675   Numeric ArrayND<Numeric, Len, Dim>::min() const {
4676     if (!shapeIsKnown_)
4677       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"min\"");
4678     if (dim_) {
4679       Numeric minval(data_[0]);
4680       for (unsigned long i = 1UL; i < len_; ++i)
4681         if (ComplexComparesAbs<Numeric>::less(data_[i], minval))
4682           minval = data_[i];
4683       return minval;
4684     } else
4685       return localData_[0];
4686   }
4687 
4688   template <typename Numeric, unsigned Len, unsigned Dim>
4689   Numeric ArrayND<Numeric, Len, Dim>::min(unsigned* index, const unsigned indexLen) const {
4690     if (!shapeIsKnown_)
4691       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"min\"");
4692     if (indexLen != dim_)
4693       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::min: incompatible index length");
4694     if (dim_) {
4695       unsigned long minind = 0UL;
4696       Numeric minval(data_[0]);
4697       for (unsigned long i = 1UL; i < len_; ++i)
4698         if (ComplexComparesAbs<Numeric>::less(data_[i], minval)) {
4699           minval = data_[i];
4700           minind = i;
4701         }
4702       convertLinearIndex(minind, index, indexLen);
4703       return minval;
4704     } else
4705       return localData_[0];
4706   }
4707 
4708   template <typename Numeric, unsigned Len, unsigned Dim>
4709   Numeric ArrayND<Numeric, Len, Dim>::max() const {
4710     if (!shapeIsKnown_)
4711       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"max\"");
4712     if (dim_) {
4713       Numeric maxval(data_[0]);
4714       for (unsigned long i = 1UL; i < len_; ++i)
4715         if (ComplexComparesAbs<Numeric>::less(maxval, data_[i]))
4716           maxval = data_[i];
4717       return maxval;
4718     } else
4719       return localData_[0];
4720   }
4721 
4722   template <typename Numeric, unsigned Len, unsigned Dim>
4723   Numeric ArrayND<Numeric, Len, Dim>::max(unsigned* index, const unsigned indexLen) const {
4724     if (!shapeIsKnown_)
4725       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"max\"");
4726     if (indexLen != dim_)
4727       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::max: incompatible index length");
4728     if (dim_) {
4729       unsigned long maxind = 0UL;
4730       Numeric maxval(data_[0]);
4731       for (unsigned long i = 1UL; i < len_; ++i)
4732         if (ComplexComparesAbs<Numeric>::less(maxval, data_[i])) {
4733           maxval = data_[i];
4734           maxind = i;
4735         }
4736       convertLinearIndex(maxind, index, indexLen);
4737       return maxval;
4738     } else
4739       return localData_[0];
4740   }
4741 
4742   // Faster function for 2d transpose
4743   template <typename Numeric, unsigned Len, unsigned Dim>
4744   ArrayND<Numeric, Len, Dim> ArrayND<Numeric, Len, Dim>::transpose() const {
4745     if (dim_ != 2)
4746       throw npstat::NpstatInvalidArgument(
4747           "npstat::ArrayND::transpose method "
4748           "can not be used with arrays of rank other than 2");
4749     unsigned newshape[2];
4750     newshape[0] = shape_[1];
4751     newshape[1] = shape_[0];
4752     ArrayND<Numeric, Len, Dim> result(newshape, dim_);
4753     const unsigned imax = shape_[0];
4754     const unsigned jmax = shape_[1];
4755     for (unsigned i = 0; i < imax; ++i)
4756       for (unsigned j = 0; j < jmax; ++j)
4757         result.data_[j * imax + i] = data_[i * jmax + j];
4758     return result;
4759   }
4760 
4761   template <typename Numeric, unsigned Len, unsigned Dim>
4762   template <typename Accumulator>
4763   ArrayND<Numeric, Len, Dim> ArrayND<Numeric, Len, Dim>::derivative(const double inscale) const {
4764     if (!shapeIsKnown_)
4765       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"derivative\"");
4766     if (!dim_)
4767       throw npstat::NpstatInvalidArgument(
4768           "npstat::ArrayND::derivative method "
4769           "can not be used with array of 0 rank");
4770 
4771     const typename ProperDblFromCmpl<Accumulator>::type scale = inscale;
4772     const unsigned maxdim = CHAR_BIT * sizeof(unsigned long);
4773     if (dim_ >= maxdim)
4774       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::derivative: array rank is too large");
4775     const unsigned long maxcycle = 1UL << dim_;
4776 
4777     ArrayShape sh;
4778     sh.reserve(dim_);
4779     for (unsigned i = 0; i < dim_; ++i) {
4780       if (shape_[i] <= 1U)
4781         throw npstat::NpstatInvalidArgument(
4782             "In npstat::ArrayND::derivative: in some dimendions "
4783             "array size is too small");
4784       sh.push_back(shape_[i] - 1U);
4785     }
4786 
4787     ArrayND result(sh);
4788     const unsigned long rLen = result.length();
4789     for (unsigned long ilin = 0; ilin < rLen; ++ilin) {
4790       result.convertLinearIndex(ilin, &sh[0], dim_);
4791 
4792       Accumulator deriv = Accumulator();
4793       for (unsigned long icycle = 0UL; icycle < maxcycle; ++icycle) {
4794         unsigned long icell = 0UL;
4795         unsigned n1 = 0U;
4796         for (unsigned i = 0; i < dim_; ++i) {
4797           if (icycle & (1UL << i)) {
4798             ++n1;
4799             icell += strides_[i] * (sh[i] + 1);
4800           } else
4801             icell += strides_[i] * sh[i];
4802         }
4803         if ((dim_ - n1) % 2U)
4804           deriv -= data_[icell];
4805         else
4806           deriv += data_[icell];
4807       }
4808       result.data_[ilin] = static_cast<Numeric>(deriv * scale);
4809     }
4810 
4811     return result;
4812   }
4813 
4814   template <typename Numeric, unsigned Len, unsigned Dim>
4815   template <typename Accumulator>
4816   Accumulator ArrayND<Numeric, Len, Dim>::sumBelowLoop(const unsigned level,
4817                                                        unsigned long idx0,
4818                                                        const unsigned* limit) const {
4819     Accumulator cdf = Accumulator();
4820     const unsigned imax = limit[level] + 1U;
4821     if (level == dim_ - 1) {
4822       Numeric* base = data_ + idx0;
4823       for (unsigned i = 0; i < imax; ++i)
4824         cdf += base[i];
4825     } else {
4826       const unsigned long stride = strides_[level];
4827       for (unsigned i = 0; i < imax; ++i, idx0 += stride)
4828         cdf += sumBelowLoop<Accumulator>(level + 1, idx0, limit);
4829     }
4830     return cdf;
4831   }
4832 
4833   template <typename Numeric, unsigned Len, unsigned Dim>
4834   template <typename Accumulator>
4835   Accumulator ArrayND<Numeric, Len, Dim>::cdfValue(const unsigned* index, const unsigned indexLen) const {
4836     if (!shapeIsKnown_)
4837       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"cdfValue\"");
4838     if (!dim_)
4839       throw npstat::NpstatInvalidArgument(
4840           "npstat::ArrayND::cdfValue method "
4841           "can not be used with array of 0 rank");
4842     if (indexLen != dim_)
4843       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::cdfValue: incompatible index length");
4844     for (unsigned i = 0; i < indexLen; ++i)
4845       if (index[i] >= shape_[i])
4846         throw npstat::NpstatOutOfRange("In npstat::ArrayND::cdfValue: index out of range");
4847     return sumBelowLoop<Accumulator>(0, 0U, index);
4848   }
4849 
4850   template <typename Numeric, unsigned Len, unsigned Dim>
4851   template <typename Accumulator>
4852   ArrayND<Numeric, Len, Dim> ArrayND<Numeric, Len, Dim>::cdfArray(const double inscale) const {
4853     if (!shapeIsKnown_)
4854       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"cdfArray\"");
4855     if (!dim_)
4856       throw npstat::NpstatInvalidArgument(
4857           "npstat::ArrayND::cdfArray method "
4858           "can not be used with array of 0 rank");
4859 
4860     const proper_double scale = inscale;
4861     const unsigned maxdim = CHAR_BIT * sizeof(unsigned long);
4862     if (dim_ >= maxdim)
4863       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::cdfArray: array rank is too large");
4864     const unsigned long maxcycle = 1UL << dim_;
4865 
4866     ArrayShape sh;
4867     sh.reserve(dim_);
4868     for (unsigned i = 0; i < dim_; ++i)
4869       sh.push_back(shape_[i] + 1U);
4870 
4871     ArrayND<Accumulator> result(sh);
4872 
4873     unsigned* psh = &sh[0];
4874     const unsigned long len = result.length();
4875     for (unsigned long ipre = 0; ipre < len; ++ipre) {
4876       result.convertLinearIndex(ipre, psh, dim_);
4877       Accumulator deriv = Accumulator();
4878       bool has0 = false;
4879       for (unsigned i = 0; i < dim_; ++i)
4880         if (psh[i]-- == 0U) {
4881           has0 = true;
4882           break;
4883         }
4884       if (!has0) {
4885         for (unsigned long icycle = 0UL; icycle < maxcycle; ++icycle) {
4886           unsigned long icell = 0UL;
4887           unsigned n1 = 0U;
4888           for (unsigned i = 0; i < dim_; ++i) {
4889             if (icycle & (1UL << i)) {
4890               ++n1;
4891               icell += result.strides_[i] * (psh[i] + 1);
4892             } else
4893               icell += result.strides_[i] * psh[i];
4894           }
4895           if (n1 < dim_) {
4896             if ((dim_ - n1) % 2U)
4897               deriv += result.data_[icell];
4898             else
4899               deriv -= result.data_[icell];
4900           }
4901         }
4902         deriv += static_cast<Accumulator>(value(psh, dim_) * scale);
4903       }
4904       result.data_[ipre] = deriv;
4905     }
4906 
4907     // The "return" will convert Accumulator type into Numeric
4908     return result;
4909   }
4910 
4911   template <typename Numeric, unsigned Len, unsigned Dim>
4912   ArrayND<Numeric, Len, Dim> ArrayND<Numeric, Len, Dim>::transpose(const unsigned pos1, const unsigned pos2) const {
4913     if (!shapeIsKnown_)
4914       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"transpose\"");
4915     if (!(pos1 < dim_ && pos2 < dim_ && pos1 != pos2))
4916       throw npstat::NpstatInvalidArgument(
4917           "In npstat::ArrayND::transpose: "
4918           "incompatible transposition indices");
4919     if (dim_ == 2)
4920       return transpose();
4921     else {
4922       // Construct the new shape
4923       unsigned newshapeBuf[Dim];
4924       unsigned* newshape = makeBuffer(dim_, newshapeBuf, Dim);
4925       copyBuffer(newshape, shape_, dim_);
4926       std::swap(newshape[pos1], newshape[pos2]);
4927 
4928       // Form the result array
4929       ArrayND<Numeric, Len, Dim> result(newshape, dim_);
4930 
4931       // Fill the result array
4932       transposeLoop(0, pos1, pos2, 0UL, 0UL, result);
4933 
4934       destroyBuffer(newshape, newshapeBuf);
4935       return result;
4936     }
4937   }
4938 
4939   template <typename Numeric, unsigned Len, unsigned Dim>
4940   template <typename Num2, unsigned Len2, unsigned Dim2>
4941   void ArrayND<Numeric, Len, Dim>::multiMirror(ArrayND<Num2, Len2, Dim2>* out) const {
4942     assert(out);
4943     if (!shapeIsKnown_)
4944       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"multiMirror\"");
4945     if (!out->shapeIsKnown_)
4946       *out = ArrayND<Num2, Len2, Dim2>(doubleShape(shape()));
4947     if (dim_ != out->dim_)
4948       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::multiMirror: incompatible argument array rank");
4949 
4950     if (dim_) {
4951       const unsigned* dshape = out->shape_;
4952       for (unsigned i = 0; i < dim_; ++i)
4953         if (dshape[i] != shape_[i] * 2U)
4954           throw npstat::NpstatInvalidArgument(
4955               "In npstat::ArrayND::multiMirror: "
4956               "incompatible argument array shape");
4957 
4958       if (dim_ >= CHAR_BIT * sizeof(unsigned long))
4959         throw npstat::NpstatInvalidArgument(
4960             "In npstat::ArrayND::multiMirror: "
4961             "array rank is too large");
4962       const unsigned long maxcycle = 1UL << dim_;
4963       std::vector<unsigned> indexbuf(dim_ * 2U);
4964       unsigned* idx = &indexbuf[0];
4965       unsigned* mirror = idx + dim_;
4966 
4967       for (unsigned long ipt = 0; ipt < len_; ++ipt) {
4968         unsigned long l = ipt;
4969         for (unsigned i = 0; i < dim_; ++i) {
4970           idx[i] = l / strides_[i];
4971           l -= (idx[i] * strides_[i]);
4972         }
4973         for (unsigned long icycle = 0UL; icycle < maxcycle; ++icycle) {
4974           for (unsigned i = 0; i < dim_; ++i) {
4975             if (icycle & (1UL << i))
4976               mirror[i] = dshape[i] - idx[i] - 1U;
4977             else
4978               mirror[i] = idx[i];
4979           }
4980           out->value(mirror, dim_) = data_[ipt];
4981         }
4982       }
4983     } else
4984       out->localData_[0] = static_cast<Num2>(localData_[0]);
4985   }
4986 
4987   template <typename Numeric, unsigned Len, unsigned Dim>
4988   template <typename Num2, unsigned Len2, unsigned Dim2>
4989   void ArrayND<Numeric, Len, Dim>::rotate(const unsigned* shifts,
4990                                           const unsigned lenShifts,
4991                                           ArrayND<Num2, Len2, Dim2>* rotated) const {
4992     assert(rotated);
4993     if (!shapeIsKnown_)
4994       throw npstat::NpstatInvalidArgument("Initialize npstat::ArrayND before calling method \"rotate\"");
4995     // Can't rotate into itself -- it will be a mess
4996     if ((void*)rotated == (void*)this)
4997       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::rotate: can not rotate array into itself");
4998     if (!rotated->shapeIsKnown_)
4999       *rotated = *this;
5000     if (dim_ != rotated->dim_)
5001       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::rotate: incompatible argument array rank");
5002     if (lenShifts != dim_)
5003       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::rotate: incompatible dimensionality of shifts");
5004 
5005     if (dim_) {
5006       assert(shifts);
5007       if (dim_ > CHAR_BIT * sizeof(unsigned long))
5008         throw npstat::NpstatInvalidArgument("In npstat::ArrayND::rotate: array rank is too large");
5009       unsigned buf[CHAR_BIT * sizeof(unsigned long)];
5010       clearBuffer(buf, dim_);
5011       (const_cast<ArrayND*>(this))
5012           ->flatCircularLoop(0U, 0UL, 0UL, buf, shape_, shifts, *rotated, scast_assign_right<Numeric, Num2>());
5013     } else
5014       rotated->localData_[0] = static_cast<Num2>(localData_[0]);
5015   }
5016 
5017   template <typename Numeric, unsigned Len, unsigned Dim>
5018   template <typename Num2, unsigned Len2, unsigned Dim2>
5019   void ArrayND<Numeric, Len, Dim>::dotProductLoop(const unsigned level,
5020                                                   unsigned long idx0,
5021                                                   unsigned long idx1,
5022                                                   unsigned long idx2,
5023                                                   const ArrayND<Num2, Len2, Dim2>& r,
5024                                                   ArrayND& result) const {
5025     // idx0 -- this object
5026     // idx1 -- dot product argument
5027     // idx2 -- result
5028     if (level == result.dim_) {
5029       Numeric sum = Numeric();
5030       const unsigned imax = r.shape_[0];
5031       const unsigned rstride = r.strides_[0];
5032       const Numeric* l = data_ + idx0;
5033       const Num2* ri = r.data_ + idx1;
5034       for (unsigned i = 0; i < imax; ++i)
5035         sum += l[i] * ri[i * rstride];
5036       result.data_[idx2] = sum;
5037     } else {
5038       const unsigned imax = result.shape_[level];
5039       for (unsigned i = 0; i < imax; ++i) {
5040         dotProductLoop(level + 1, idx0, idx1, idx2, r, result);
5041         idx2 += result.strides_[level];
5042         if (level < dim_ - 1)
5043           idx0 += strides_[level];
5044         else
5045           idx1 += r.strides_[level + 2 - dim_];
5046       }
5047     }
5048   }
5049 
5050   template <typename Numeric, unsigned Len, unsigned Dim>
5051   template <typename Num2, unsigned Len2, unsigned Dim2>
5052   ArrayND<Numeric, Len, Dim> ArrayND<Numeric, Len, Dim>::dot(const ArrayND<Num2, Len2, Dim2>& r) const {
5053     if (!dim_)
5054       throw npstat::NpstatInvalidArgument(
5055           "npstat::ArrayND::dot method "
5056           "can not be used with array of 0 rank");
5057     if (!r.dim_)
5058       throw npstat::NpstatInvalidArgument(
5059           "npstat::ArrayND::dot method "
5060           "can not be used with argument array of 0 rank");
5061     if (shape_[dim_ - 1] != r.shape_[0])
5062       throw npstat::NpstatInvalidArgument("In npstat::ArrayND::dot: incompatible argument array shape");
5063 
5064     if (dim_ == 1 && r.dim_ == 1) {
5065       // Special case: the result is of 0 rank
5066       ArrayND<Numeric, Len, Dim> result(static_cast<unsigned*>(nullptr), 0U);
5067       Numeric sum = Numeric();
5068       const unsigned imax = shape_[0];
5069       for (unsigned i = 0; i < imax; ++i)
5070         sum += data_[i] * r.data_[i];
5071       result() = sum;
5072       return result;
5073     } else {
5074       unsigned newshapeBuf[2 * Dim];
5075       unsigned* newshape = makeBuffer(dim_ + r.dim_ - 2, newshapeBuf, 2 * Dim);
5076       copyBuffer(newshape, shape_, dim_ - 1);
5077       copyBuffer(newshape + (dim_ - 1), r.shape_ + 1, r.dim_ - 1);
5078       ArrayND<Numeric, Len, Dim> result(newshape, dim_ + r.dim_ - 2);
5079 
5080       dotProductLoop(0U, 0UL, 0UL, 0UL, r, result);
5081 
5082       destroyBuffer(newshape, newshapeBuf);
5083       return result;
5084     }
5085   }
5086 
5087   template <typename Numeric, unsigned Len, unsigned Dim>
5088   inline unsigned ArrayND<Numeric, Len, Dim>::span(const unsigned dim) const {
5089     if (dim >= dim_)
5090       throw npstat::NpstatOutOfRange("In npstat::ArrayND::span: dimension number is out of range");
5091     return shape_[dim];
5092   }
5093 
5094   template <typename Numeric, unsigned Len, unsigned Dim>
5095   unsigned ArrayND<Numeric, Len, Dim>::maximumSpan() const {
5096     unsigned maxspan = 0;
5097     for (unsigned i = 0; i < dim_; ++i)
5098       if (shape_[i] > maxspan)
5099         maxspan = shape_[i];
5100     return maxspan;
5101   }
5102 
5103   template <typename Numeric, unsigned Len, unsigned Dim>
5104   unsigned ArrayND<Numeric, Len, Dim>::minimumSpan() const {
5105     if (dim_ == 0)
5106       return 0U;
5107     unsigned minspan = shape_[0];
5108     for (unsigned i = 1;