ArrayND.icc
No OneTemporary
Actions

Size

222 KB

Referenced Files

None

Subscribers

None

ArrayND.icc
View Options

	#include <cmath>
	#include <climits>
	#include <algorithm>
	#include <sstream>
	#include <stdexcept>
	#include <set>

	#include "geners/GenericIO.hh"
	#include "geners/IOIsUnsigned.hh"

	#include "npstat/nm/allocators.hh"

	#include "npstat/nm/interpolate.hh"
	#include "npstat/nm/absDifference.hh"
	#include "npstat/nm/ComplexComparesFalse.hh"
	#include "npstat/nm/ComplexComparesAbs.hh"
	#include "npstat/nm/LinearMapper1d.hh"

	#define me_macro_check_loop_prerequisites(METHOD, INNERLOOP) /**/ \
	template<typename Numeric, unsigned Len, unsigned Dim> \
	template <typename Num2, unsigned Len2, unsigned Dim2, class Functor> \
	void ArrayND<Numeric,Len,Dim>:: METHOD ( \
	ArrayND<Num2, Len2, Dim2>& other, \
	const unsigned* thisCorner, \
	const unsigned* range, \
	const unsigned* otherCorner, \
	const unsigned arrLen, \
	Functor binaryFunct) \
	{ \
	if (!shapeIsKnown_) throw std::invalid_argument( \
	"Initialize npstat::ArrayND before calling method \"" \
	#METHOD "\""); \
	if (!other.shapeIsKnown_) throw std::invalid_argument( \
	"In npstat::ArrayND::" #METHOD ": uninitialized argument array");\
	if (dim_ != other.dim_) throw std::invalid_argument( \
	"In npstat::ArrayND::" #METHOD ": incompatible argument array rank");\
	if (arrLen != dim_) throw std::invalid_argument( \
	"In npstat::ArrayND::" #METHOD ": incompatible index length"); \
	if (dim_) \
	{ \
	assert(thisCorner); \
	assert(range); \
	assert(otherCorner); \
	INNERLOOP (0U, 0UL, 0UL, thisCorner, range, \
	otherCorner, other, binaryFunct); \
	} \
	else \
	binaryFunct(localData_[0], other.localData_[0]); \
	}

	namespace npstat {
	template<typename Numeric, unsigned Len, unsigned Dim>
	template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
	void ArrayND<Numeric,Len,Dim>::commonSubrangeLoop(
	unsigned level, unsigned long idx0,
	unsigned long idx1,
	const unsigned* thisCorner,
	const unsigned* range,
	const unsigned* otherCorner,
	ArrayND<Num2, Len2, Dim2>& r,
	Functor binaryFunct)
	{
	const unsigned imax = range[level];

	if (level == dim_ - 1)
	{
	Numeric* left = data_ + (idx0 + thisCorner[level]);
	Numeric* const lMax = data_ + (idx0 + shape_[level]);
	Num2* right = r.data_ + (idx1 + otherCorner[level]);
	Num2* const rMax = r.data_ + (idx1 + r.shape_[level]);

	for (unsigned i=0; i<imax && left<lMax && right<rMax; ++i)
	binaryFunct(left++, right++);
	}
	else
	{
	const unsigned long leftStride = strides_[level];
	const unsigned long leftMax = idx0 + shape_[level]*leftStride;
	idx0 += thisCorner[level]*leftStride;
	const unsigned long rightStride = r.strides_[level];
	const unsigned long rightMax = idx1 + r.shape_[level]*rightStride;
	idx1 += otherCorner[level]*rightStride;

	for (unsigned i=0; i<imax && idx0 < leftMax && idx1 < rightMax;
	++i, idx0 += leftStride, idx1 += rightStride)
	commonSubrangeLoop(level+1, idx0, idx1, thisCorner, range,
	otherCorner, r, binaryFunct);
	}
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
	void ArrayND<Numeric,Len,Dim>::dualCircularLoop(
	unsigned level, unsigned long idx0,
	unsigned long idx1,
	const unsigned* thisCorner,
	const unsigned* range,
	const unsigned* otherCorner,
	ArrayND<Num2, Len2, Dim2>& r,
	Functor binaryFunct)
	{
	const unsigned imax = range[level];
	const unsigned leftShift = thisCorner[level];
	const unsigned leftPeriod = shape_[level];
	const unsigned rightShift = otherCorner[level];
	const unsigned rightPeriod = r.shape_[level];

	if (level == dim_ - 1)
	{
	Numeric* left = data_ + idx0;
	Num2* right = r.data_ + idx1;
	for (unsigned i=0; i<imax; ++i)
	binaryFunct(left[(i+leftShift)%leftPeriod],
	right[(i+rightShift)%rightPeriod]);
	}
	else
	{
	const unsigned long leftStride = strides_[level];
	const unsigned long rightStride = r.strides_[level];
	for (unsigned i=0; i<imax; ++i)
	dualCircularLoop(
	level+1, idx0+((i+leftShift)%leftPeriod)*leftStride,
	idx1+((i+rightShift)%rightPeriod)*rightStride,
	thisCorner, range, otherCorner, r, binaryFunct);
	}
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
	void ArrayND<Numeric,Len,Dim>::flatCircularLoop(
	unsigned level, unsigned long idx0,
	unsigned long idx1,
	const unsigned* thisCorner,
	const unsigned* range,
	const unsigned* otherCorner,
	ArrayND<Num2, Len2, Dim2>& r,
	Functor binaryFunct)
	{
	const unsigned imax = range[level];
	const unsigned rightShift = otherCorner[level];
	const unsigned rightPeriod = r.shape_[level];

	if (level == dim_ - 1)
	{
	Numeric* left = data_ + (idx0 + thisCorner[level]);
	Numeric* const lMax = data_ + (idx0 + shape_[level]);
	Num2* right = r.data_ + idx1;

	for (unsigned i=0; i<imax && left<lMax; ++i)
	binaryFunct(*left++, right[(i+rightShift)%rightPeriod]);
	}
	else
	{
	const unsigned long leftStride = strides_[level];
	const unsigned long leftMax = idx0 + shape_[level]*leftStride;
	idx0 += thisCorner[level]*leftStride;
	const unsigned long rightStride = r.strides_[level];

	for (unsigned i=0; i<imax && idx0 < leftMax; ++i, idx0+=leftStride)
	flatCircularLoop(
	level+1, idx0,
	idx1+((i+rightShift)%rightPeriod)*rightStride,
	thisCorner, range, otherCorner, r, binaryFunct);
	}
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
	void ArrayND<Numeric,Len,Dim>::circularFlatLoop(
	unsigned level, unsigned long idx0,
	unsigned long idx1,
	const unsigned* thisCorner,
	const unsigned* range,
	const unsigned* otherCorner,
	ArrayND<Num2, Len2, Dim2>& r,
	Functor binaryFunct)
	{
	const unsigned imax = range[level];
	const unsigned leftShift = thisCorner[level];
	const unsigned leftPeriod = shape_[level];

	if (level == dim_ - 1)
	{
	Numeric* left = data_ + idx0;
	Num2* right = r.data_ + (idx1 + otherCorner[level]);
	Num2* const rMax = r.data_ + (idx1 + r.shape_[level]);

	for (unsigned i=0; i<imax && right<rMax; ++i)
	binaryFunct(left[(i+leftShift)%leftPeriod], *right++);
	}
	else
	{
	const unsigned long leftStride = strides_[level];
	const unsigned long rightStride = r.strides_[level];
	const unsigned long rightMax = idx1 + r.shape_[level]*rightStride;
	idx1 += otherCorner[level]*rightStride;

	for (unsigned i=0; i<imax && idx1<rightMax; ++i, idx1+=rightStride)
	circularFlatLoop(
	level+1, idx0+((i+leftShift)%leftPeriod)*leftStride,
	idx1, thisCorner, range, otherCorner, r, binaryFunct);
	}
	}

	me_macro_check_loop_prerequisites(jointSubrangeScan, commonSubrangeLoop)
	me_macro_check_loop_prerequisites(dualCircularScan, dualCircularLoop)
	me_macro_check_loop_prerequisites(flatCircularScan, flatCircularLoop)
	me_macro_check_loop_prerequisites(circularFlatScan, circularFlatLoop)

	template <typename Numeric, unsigned StackLen, unsigned StackDim>
	template <typename Num2, unsigned Len2, unsigned Dim2>
	Numeric ArrayND<Numeric,StackLen,StackDim>::marginalizeInnerLoop(
	unsigned long idx, const unsigned levelPr, unsigned long idxPr,
	const ArrayND<Num2,Len2,Dim2>& prior,
	const unsigned* indexMap) const
	{
	Numeric sum = Numeric();
	const unsigned long myStride = strides_[indexMap[levelPr]];
	const unsigned imax = prior.shape_[levelPr];
	assert(imax == shape_[indexMap[levelPr]]);
	if (levelPr == prior.dim_ - 1)
	{
	for (unsigned i=0; i<imax; ++i)
	sum += data_[idx+imyStride]prior.data_[idxPr++];
	}
	else
	{
	const unsigned long priorStride = prior.strides_[levelPr];
	for (unsigned i=0; i<imax; ++i)
	{
	sum += marginalizeInnerLoop(idx, levelPr+1U, idxPr,
	prior, indexMap);
	idx += myStride;
	idxPr += priorStride;
	}
	}
	return sum;
	}

	template <typename Numeric, unsigned StackLen, unsigned StackDim>
	template <typename Num2, unsigned Len2, unsigned Dim2>
	void ArrayND<Numeric,StackLen,StackDim>::marginalizeLoop(
	const unsigned level, unsigned long idx,
	const unsigned levelRes, unsigned long idxRes,
	const ArrayND<Num2,Len2,Dim2>& prior,
	const unsigned* indexMap, ArrayND& result) const
	{
	if (level == dim_)
	{
	const Numeric res = marginalizeInnerLoop(
	idx, 0U, 0UL, prior, indexMap);
	if (result.dim_)
	result.data_[idxRes] = res;
	else
	result.localData_[0] = res;
	}
	else
	{
	// Check if this level is mapped or not
	bool mapped = false;
	for (unsigned i=0; i<prior.dim_; ++i)
	if (level == indexMap[i])
	{
	mapped = true;
	break;
	}
	if (mapped)
	marginalizeLoop(level+1U, idx, levelRes, idxRes,
	prior, indexMap, result);
	else
	{
	const unsigned imax = shape_[level];
	const unsigned long myStride = strides_[level];
	const unsigned long resStride = result.strides_[levelRes];
	for (unsigned i=0; i<imax; ++i)
	{
	marginalizeLoop(level+1U, idx, levelRes+1U, idxRes,
	prior, indexMap, result);
	idx += myStride;
	idxRes += resStride;
	}
	}
	}
	}

	template<typename Numeric, unsigned StackLen, unsigned StackDim>
	template <typename Num2, unsigned Len2, unsigned Dim2>
	ArrayND<Numeric,StackLen,StackDim>
	ArrayND<Numeric,StackLen,StackDim>::marginalize(
	const ArrayND<Num2,Len2,Dim2>& prior,
	const unsigned* indexMap, const unsigned mapLen) const
	{
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"marginalize\"");
	if (!(prior.dim_ && prior.dim_ <= dim_)) throw std::invalid_argument(
	"In npstat::ArrayND::marginalize: incompatible argument array rank");
	const unsigned resultDim = dim_ - prior.dim_;

	// Check that the index map is reasonable
	if (mapLen != prior.dim_) throw std::invalid_argument(
	"In npstat::ArrayND::marginalize: incompatible index map length");
	assert(indexMap);
	for (unsigned i=0; i<mapLen; ++i)
	{
	const unsigned thisInd = indexMap[i];
	if (shape_[thisInd] != prior.shape_[i]) throw std::invalid_argument(
	"In npstat::ArrayND::marginalize: "
	"incompatible argument array dimensions");
	if (thisInd >= dim_) throw std::out_of_range(
	"In npstat::ArrayND::marginalize: index map entry out of range");
	for (unsigned j=0; j<i; ++j)
	if (indexMap[j] == thisInd) throw std::invalid_argument(
	"In npstat::ArrayND::marginalize: "
	"duplicate entry in the index map");
	}

	// Build the shape for the array of results
	ArrayShape newShape;
	newShape.reserve(resultDim);
	for (unsigned i=0; i<dim_; ++i)
	{
	bool mapped = false;
	for (unsigned j=0; j<mapLen; ++j)
	if (indexMap[j] == i)
	{
	mapped = true;
	break;
	}
	if (!mapped)
	newShape.push_back(shape_[i]);
	}

	ArrayND result(newShape);
	assert(result.dim_ == resultDim);
	bool calculated = false;
	if (resultDim == 0)
	{
	calculated = true;
	for (unsigned i=0; i<dim_; ++i)
	if (indexMap[i] != i)
	{
	calculated = false;
	break;
	}
	if (calculated)
	{
	Numeric sum = Numeric();
	for (unsigned long i=0; i<len_; ++i)
	sum += data_[i]*prior.data_[i];
	result.localData_[0] = sum;
	}
	}

	if (!calculated)
	marginalizeLoop(0U, 0UL, 0U, 0UL, prior, indexMap, result);

	return result;
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	void ArrayND<Numeric,Len,Dim>::buildFromShapePtr(
	const unsigned* sizes, const unsigned dim)
	{
	dim_ = dim;
	if (dim_)
	{
	assert(sizes);
	for (unsigned i=0; i<dim_; ++i)
	if (sizes[i] == 0)
	throw std::invalid_argument(
	"In npstat::ArrayND::buildFromShapePtr: "
	"detected span of zero");

	// Copy the array shape and figure out the array length
	shape_ = makeBuffer(dim_, localShape_, Dim);
	for (unsigned i=0; i<dim_; ++i)
	{
	shape_[i] = sizes[i];
	len_ *= shape_[i];
	}

	// Figure out the array strides
	buildStrides();

	// Allocate the data array
	data_ = makeBuffer(len_, localData_, Len);
	}
	else
	{
	localData_[0] = Numeric();
	data_ = localData_;
	}
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	void ArrayND<Numeric,Len,Dim>::buildExtFromShapePtr(
	const unsigned* sizes, const unsigned dim, Numeric* extdata)
	{
	if (dim)
	{
	assert(sizes);
	len_ = 1UL;
	shapeIsKnown_ = true;
	dataIsExternal_ = true;
	dim_ = dim;
	for (unsigned i=0; i<dim_; ++i)
	if (sizes[i] == 0)
	throw std::invalid_argument(
	"In npstat::ArrayND::buildExtFromShapePtr: "
	"detected span of zero");

	// Copy the array shape and figure out the array length
	shape_ = makeBuffer(dim_, localShape_, Dim);
	for (unsigned i=0; i<dim_; ++i)
	{
	shape_[i] = sizes[i];
	len_ *= shape_[i];
	}

	// Figure out the array strides
	buildStrides();

	// Assign the data
	data_ = extdata;
	}
	else
	throw std::invalid_argument(
	"In npstat::ArrayND::buildExtFromShapePtr: can not "
	"manage external memory using rank 0 arrays");
	}

	template<typename Numeric, unsigned StackLen, unsigned StackDim>
	template <typename Num2, unsigned Len2, unsigned Dim2>
	ArrayND<Numeric,StackLen,StackDim>::ArrayND(
	const ArrayND<Num2, Len2, Dim2>& slicedArray,
	const unsigned *fixedIndices, const unsigned nFixedIndices)
	: data_(0), strides_(0), shape_(0),
	len_(1UL), dim_(slicedArray.dim_ - nFixedIndices),
	shapeIsKnown_(true), dataIsExternal_(false)
	{
	if (nFixedIndices)
	{
	assert(fixedIndices);
	if (nFixedIndices > slicedArray.dim_) throw std::invalid_argument(
	"In npstat::ArrayND slicing constructor: too many fixed indices");
	if (!slicedArray.shapeIsKnown_) throw std::invalid_argument(
	"In npstat::ArrayND slicing constructor: "
	"uninitialized argument array");

	// Check that the fixed indices are within range
	for (unsigned j=0; j<nFixedIndices; ++j)
	if (fixedIndices[j] >= slicedArray.dim_)
	throw std::out_of_range("In npstat::ArrayND slicing "
	"constructor: fixed index out of range");

	// Build the shape for the slice
	shape_ = makeBuffer(dim_, localShape_, StackDim);
	unsigned idim = 0;
	for (unsigned i=0; i<slicedArray.dim_; ++i)
	{
	bool fixed = false;
	for (unsigned j=0; j<nFixedIndices; ++j)
	if (fixedIndices[j] == i)
	{
	fixed = true;
	break;
	}
	if (!fixed)
	{
	assert(idim < dim_);
	shape_[idim++] = slicedArray.shape_[i];
	}
	}
	assert(idim == dim_);

	if (dim_)
	{
	// Copy the array shape and figure out the array length
	for (unsigned i=0; i<dim_; ++i)
	len_ *= shape_[i];

	// Figure out the array strides
	buildStrides();

	// Allocate the data array
	data_ = makeBuffer(len_, localData_, StackLen);
	}
	else
	{
	localData_[0] = Numeric();
	data_ = localData_;
	}
	}
	else
	{
	new (this) ArrayND(slicedArray);
	}
	}

	template<typename Numeric, unsigned StackLen, unsigned StackDim>
	ArrayShape ArrayND<Numeric,StackLen,StackDim>::sliceShape(
	const unsigned *fixedIndices, const unsigned nFixedIndices) const
	{
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"sliceShape\"");
	if (nFixedIndices)
	{
	assert(fixedIndices);
	if (nFixedIndices > dim_)
	throw std::invalid_argument(
	"In npstat::ArrayND::sliceShape: too many fixed indices");
	for (unsigned j=0; j<nFixedIndices; ++j)
	if (fixedIndices[j] >= dim_)
	throw std::out_of_range("In npstat::ArrayND::sliceShape: "
	"fixed index out of range");
	ArrayShape sh;
	sh.reserve(dim_ - nFixedIndices);
	for (unsigned i=0; i<dim_; ++i)
	{
	bool fixed = false;
	for (unsigned j=0; j<nFixedIndices; ++j)
	if (fixedIndices[j] == i)
	{
	fixed = true;
	break;
	}
	if (!fixed)
	sh.push_back(shape_[i]);
	}
	return sh;
	}
	else
	return ArrayShape(shape_, shape_+dim_);
	}

	template<typename Numeric, unsigned StackLen, unsigned StackDim>
	template <typename Num2, unsigned Len2, unsigned Dim2>
	unsigned long ArrayND<Numeric,StackLen,StackDim>::verifySliceCompatibility(
	const ArrayND<Num2,Len2,Dim2>& slice,
	const unsigned *fixedIndices,
	const unsigned *fixedIndexValues,
	const unsigned nFixedIndices) const
	{
	if (!(nFixedIndices && nFixedIndices <= dim_))
	throw std::invalid_argument(
	"In npstat::ArrayND::verifySliceCompatibility: "
	"invalid number of fixed indices");
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling "
	"method \"verifySliceCompatibility\"");
	if (!slice.shapeIsKnown_) throw std::invalid_argument(
	"In npstat::ArrayND::verifySliceCompatibility: "
	"uninitialized argument array");
	if (slice.dim_ != dim_ - nFixedIndices) throw std::invalid_argument(
	"In npstat::ArrayND::verifySliceCompatibility: "
	"incompatible argument array rank");
	assert(fixedIndices);
	assert(fixedIndexValues);

	for (unsigned j=0; j<nFixedIndices; ++j)
	if (fixedIndices[j] >= dim_) throw std::out_of_range(
	"In npstat::ArrayND::verifySliceCompatibility: "
	"fixed index out of range");

	// Check slice shape compatibility
	unsigned long idx = 0UL;
	unsigned sliceDim = 0U;
	for (unsigned i=0; i<dim_; ++i)
	{
	bool fixed = false;
	for (unsigned j=0; j<nFixedIndices; ++j)
	if (fixedIndices[j] == i)
	{
	fixed = true;
	if (fixedIndexValues[j] >= shape_[i])
	throw std::out_of_range(
	"In npstat::ArrayND::verifySliceCompatibility: "
	"fixed index value out of range");
	idx += fixedIndexValues[j]*strides_[i];
	break;
	}
	if (!fixed)
	{
	if (shape_[i] != slice.shape_[sliceDim])
	throw std::invalid_argument(
	"In npstat::ArrayND::verifySliceCompatibility: "
	"incompatible argument array dimensions");
	++sliceDim;
	}
	}
	assert(sliceDim == slice.dim_);
	return idx;
	}

	template<typename Numeric, unsigned StackLen, unsigned StackDim>
	unsigned long
	ArrayND<Numeric,StackLen,StackDim>::verifyBufferSliceCompatibility(
	const unsigned long bufLen,
	const unsigned *fixedIndices,
	const unsigned *fixedIndexValues,
	const unsigned nFixedIndices,
	unsigned long* sliceStrides) const
	{
	if (!(nFixedIndices && nFixedIndices <= dim_))
	throw std::invalid_argument(
	"In npstat::ArrayND::verifyBufferSliceCompatibility: "
	"invalid number of fixed indices");
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling "
	"method \"verifyBufferSliceCompatibility\"");
	assert(fixedIndices);
	assert(fixedIndexValues);

	for (unsigned j=0; j<nFixedIndices; ++j)
	if (fixedIndices[j] >= dim_) throw std::out_of_range(
	"In npstat::ArrayND::verifyBufferSliceCompatibility: "
	"fixed index out of range");

	// Figure out slice shape. Store it, temporarily, in sliceStrides.
	unsigned long idx = 0UL;
	unsigned sliceDim = 0U;
	for (unsigned i=0; i<dim_; ++i)
	{
	bool fixed = false;
	for (unsigned j=0; j<nFixedIndices; ++j)
	if (fixedIndices[j] == i)
	{
	fixed = true;
	if (fixedIndexValues[j] >= shape_[i])
	throw std::out_of_range(
	"In npstat::ArrayND::verifyBufferSliceCompatibility:"
	" fixed index value out of range");
	idx += fixedIndexValues[j]*strides_[i];
	break;
	}
	if (!fixed)
	sliceStrides[sliceDim++] = shape_[i];
	}
	assert(sliceDim + nFixedIndices == dim_);

	// Convert the slice shape into slice strides
	unsigned long expectedBufLen = 1UL;
	if (sliceDim)
	{
	unsigned long shapeJ = sliceStrides[sliceDim - 1];
	sliceStrides[sliceDim - 1] = 1UL;
	for (unsigned j=sliceDim - 1; j>0; --j)
	{
	const unsigned long nextStride = sliceStrides[j]*shapeJ;
	shapeJ = sliceStrides[j - 1];
	sliceStrides[j - 1] = nextStride;
	}
	expectedBufLen = sliceStrides[0]*shapeJ;
	}
	if (expectedBufLen != bufLen)
	throw std::invalid_argument(
	"In npstat::ArrayND::verifyBufferSliceCompatibility: "
	"invalid memory buffer length");

	return idx;
	}

	template<typename Numeric, unsigned StackLen, unsigned StackDim>
	template <typename Num2, class Op>
	void ArrayND<Numeric,StackLen,StackDim>::jointSliceLoop(
	const unsigned level, const unsigned long idx0,
	const unsigned level1, const unsigned long idx1,
	Num2* sliceData, const unsigned long* sliceStrides,
	const unsigned *fixedIndices,
	const unsigned *fixedIndexValues,
	const unsigned nFixedIndices,
	Op fcn)
	{
	bool fixed = false;
	for (unsigned j=0; j<nFixedIndices; ++j)
	if (fixedIndices[j] == level)
	{
	fixed = true;
	break;
	}
	if (fixed)
	jointSliceLoop(level+1, idx0, level1, idx1,
	sliceData, sliceStrides, fixedIndices,
	fixedIndexValues, nFixedIndices, fcn);
	else
	{
	const unsigned imax = shape_[level];
	const unsigned long stride = strides_[level];

	if (level1 == dim_ - nFixedIndices - 1)
	{
	sliceData += idx1;
	Numeric* localData = data_ + idx0;
	for (unsigned i = 0; i<imax; ++i)
	fcn(localData[i*stride], sliceData[i]);
	}
	else
	{
	const unsigned long stride2 = sliceStrides[level1];
	for (unsigned i = 0; i<imax; ++i)
	jointSliceLoop(level+1, idx0+i*stride,
	level1+1, idx1+i*stride2,
	sliceData, sliceStrides, fixedIndices,
	fixedIndexValues, nFixedIndices, fcn);
	}
	}
	}

	template<typename Numeric, unsigned StackLen, unsigned StackDim>
	template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
	void ArrayND<Numeric,StackLen,StackDim>::jointSliceScan(
	ArrayND<Num2,Len2,Dim2>& slice,
	const unsigned *fixedIndices,
	const unsigned *fixedIndexValues,
	const unsigned nFixedIndices,
	Functor binaryFunct)
	{
	const unsigned long idx = verifySliceCompatibility(
	slice, fixedIndices, fixedIndexValues, nFixedIndices);
	if (slice.dim_)
	jointSliceLoop(0U, idx, 0U, 0UL, slice.data_, slice.strides_,
	fixedIndices, fixedIndexValues,
	nFixedIndices, binaryFunct);
	else
	binaryFunct(data_[idx], slice.localData_[0]);
	}

	template<typename Numeric, unsigned StackLen, unsigned StackDim>
	template <typename Num2, class Functor>
	void ArrayND<Numeric,StackLen,StackDim>::jointMemSliceScan(
	Num2* slice, const unsigned long len,
	const unsigned fixedIndices, const unsigned fixedIndexValues,
	unsigned nFixedIndices, Functor binaryFunct)
	{
	assert(slice);
	if (dim_ > CHAR_BIT*sizeof(unsigned long))
	throw std::invalid_argument(
	"In npstat::ArrayND::jointMemSliceScan: "
	"rank of this array is too large");
	unsigned long sliceStrides[CHAR_BIT*sizeof(unsigned long)];
	const unsigned long idx = verifyBufferSliceCompatibility(
	len, fixedIndices, fixedIndexValues, nFixedIndices, sliceStrides);
	if (dim_ > nFixedIndices)
	jointSliceLoop(0U, idx, 0U, 0UL, slice, sliceStrides,
	fixedIndices, fixedIndexValues,
	nFixedIndices, binaryFunct);
	else
	binaryFunct(data_[idx], *slice);
	}

	template<typename Numeric, unsigned StackLen, unsigned StackDim>
	template <typename Num2>
	void ArrayND<Numeric,StackLen,StackDim>::projectInnerLoop(
	const unsigned level, unsigned long idx0,
	unsigned* currentIndex,
	AbsArrayProjector<Numeric,Num2>& projector,
	const unsigned *projectedIndices,
	const unsigned nProjectedIndices) const
	{
	// level : dimension number among indices which are being iterated
	const unsigned idx = projectedIndices[level];
	const unsigned imax = shape_[idx];
	const unsigned long stride = strides_[idx];
	const bool last = (level == nProjectedIndices - 1);

	for (unsigned i = 0; i<imax; ++i)
	{
	currentIndex[idx] = i;
	if (last)
	projector.process(currentIndex, dim_, idx0, data_[idx0]);
	else
	projectInnerLoop(level+1, idx0, currentIndex, projector,
	projectedIndices, nProjectedIndices);
	idx0 += stride;
	}
	}

	template<typename Numeric, unsigned StackLen, unsigned StackDim>
	template <typename Num2, unsigned Len2, unsigned Dim2,
	typename Num3, class Op>
	void ArrayND<Numeric,StackLen,StackDim>::projectLoop(
	const unsigned level, const unsigned long idx0,
	const unsigned level1, const unsigned long idx1,
	unsigned* currentIndex,
	ArrayND<Num2,Len2,Dim2>* projection,
	AbsArrayProjector<Numeric,Num3>& projector,
	const unsigned *projectedIndices,
	const unsigned nProjectedIndices, Op fcn) const
	{
	// level : dimension number in this array
	// level1 : dimension number in the projection
	// idx0 : linear index in this array
	// idx1 : linear index in the projection
	// currentIndex : cycled over in this loop with the exception of the
	// dimensions which are iterated over to build the
	// projection
	if (level == dim_)
	{
	assert(level1 == projection->dim_);
	projector.clear();
	projectInnerLoop(0U, idx0, currentIndex, projector,
	projectedIndices, nProjectedIndices);
	if (projection->dim_)
	fcn(projection->data_[idx1], projector.result());
	else
	fcn(projection->localData_[0], projector.result());
	}
	else
	{
	bool iterated = false;
	for (unsigned j=0; j<nProjectedIndices; ++j)
	if (projectedIndices[j] == level)
	{
	iterated = true;
	break;
	}
	if (iterated)
	{
	// This index will be iterated over inside "projectInnerLoop"
	projectLoop(level+1, idx0, level1, idx1,
	currentIndex, projection, projector,
	projectedIndices, nProjectedIndices, fcn);
	}
	else
	{
	const unsigned imax = shape_[level];
	const unsigned long stride = strides_[level];
	// We will not be able to get here if projection->dim_ is 0.
	// Therefore, it is safe to access projection->strides_.
	const unsigned long stride2 = projection->strides_[level1];
	for (unsigned i = 0; i<imax; ++i)
	{
	currentIndex[level] = i;
	projectLoop(level+1, idx0+i*stride,
	level1+1, idx1+i*stride2,
	currentIndex, projection, projector,
	projectedIndices, nProjectedIndices, fcn);
	}
	}
	}
	}

	template<typename Numeric, unsigned StackLen, unsigned StackDim>
	template <typename Num2, unsigned Len2, unsigned Dim2>
	void ArrayND<Numeric,StackLen,StackDim>::verifyProjectionCompatibility(
	const ArrayND<Num2,Len2,Dim2>& projection,
	const unsigned *projectedIndices,
	const unsigned nProjectedIndices) const
	{
	if (!(nProjectedIndices && nProjectedIndices <= dim_))
	throw std::invalid_argument(
	"In npstat::ArrayND::verifyProjectionCompatibility: "
	"invalid number of projected indices");
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling "
	"method \"verifyProjectionCompatibility\"");
	if (!projection.shapeIsKnown_) throw std::invalid_argument(
	"In npstat::ArrayND::verifyProjectionCompatibility: "
	"uninitialized argument array");
	if (projection.dim_ != dim_ - nProjectedIndices)
	throw std::invalid_argument(
	"In npstat::ArrayND::verifyProjectionCompatibility: "
	"incompatible argument array rank");
	assert(projectedIndices);

	for (unsigned j=0; j<nProjectedIndices; ++j)
	if (projectedIndices[j] >= dim_) throw std::out_of_range(
	"In npstat::ArrayND::verifyProjectionCompatibility: "
	"projected index out of range");

	// Check projection shape compatibility
	unsigned sliceDim = 0U;
	for (unsigned i=0; i<dim_; ++i)
	{
	bool fixed = false;
	for (unsigned j=0; j<nProjectedIndices; ++j)
	if (projectedIndices[j] == i)
	{
	fixed = true;
	break;
	}
	if (!fixed)
	{
	if (shape_[i] != projection.shape_[sliceDim])
	throw std::invalid_argument(
	"In npstat::ArrayND::verifyProjectionCompatibility: "
	"incompatible argument array dimensions");
	++sliceDim;
	}
	}
	assert(sliceDim == projection.dim_);
	}

	template<typename Numeric, unsigned StackLen, unsigned StackDim>
	template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
	void ArrayND<Numeric,StackLen,StackDim>::project(
	ArrayND<Num2,Len2,Dim2>* projection,
	AbsArrayProjector<Numeric,Num3>& projector,
	const unsigned *projectedIndices,
	const unsigned nProjectedIndices) const
	{
	assert(projection);
	verifyProjectionCompatibility(*projection, projectedIndices,
	nProjectedIndices);
	unsigned ibuf[StackDim];
	unsigned* buf = makeBuffer(dim_, ibuf, StackDim);
	for (unsigned i=0; i<dim_; ++i)
	buf[i] = 0U;
	projectLoop(0U, 0UL, 0U, 0UL, buf, projection,
	projector, projectedIndices, nProjectedIndices,
	scast_assign_left<Num2,Num3>());
	destroyBuffer(buf, ibuf);
	}

	template<typename Numeric, unsigned StackLen, unsigned StackDim>
	template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
	void ArrayND<Numeric,StackLen,StackDim>::addToProjection(
	ArrayND<Num2,Len2,Dim2>* projection,
	AbsArrayProjector<Numeric,Num3>& projector,
	const unsigned *projectedIndices,
	const unsigned nProjectedIndices) const
	{
	assert(projection);
	verifyProjectionCompatibility(*projection, projectedIndices,
	nProjectedIndices);
	unsigned ibuf[StackDim];
	unsigned* buf = makeBuffer(dim_, ibuf, StackDim);
	for (unsigned i=0; i<dim_; ++i)
	buf[i] = 0U;
	projectLoop(0U, 0UL, 0U, 0UL, buf, projection,
	projector, projectedIndices, nProjectedIndices,
	scast_pluseq_left<Num2,Num3>());
	destroyBuffer(buf, ibuf);
	}

	template<typename Numeric, unsigned StackLen, unsigned StackDim>
	template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
	void ArrayND<Numeric,StackLen,StackDim>::subtractFromProjection(
	ArrayND<Num2,Len2,Dim2>* projection,
	AbsArrayProjector<Numeric,Num3>& projector,
	const unsigned *projectedIndices,
	const unsigned nProjectedIndices) const
	{
	assert(projection);
	verifyProjectionCompatibility(*projection, projectedIndices,
	nProjectedIndices);
	unsigned ibuf[StackDim];
	unsigned* buf = makeBuffer(dim_, ibuf, StackDim);
	for (unsigned i=0; i<dim_; ++i)
	buf[i] = 0U;
	projectLoop(0U, 0UL, 0U, 0UL, buf, projection,
	projector, projectedIndices, nProjectedIndices,
	scast_minuseq_left<Num2,Num3>());
	destroyBuffer(buf, ibuf);
	}

	template<typename Numeric, unsigned StackLen, unsigned StackDim>
	template <typename Num2>
	void ArrayND<Numeric,StackLen,StackDim>::projectInnerLoop2(
	const unsigned level, const unsigned long idx0,
	AbsVisitor<Numeric,Num2>& projector,
	const unsigned *projectedIndices,
	const unsigned nProjectedIndices) const
	{
	const unsigned idx = projectedIndices[level];
	const unsigned imax = shape_[idx];
	const unsigned long stride = strides_[idx];
	const bool last = (level == nProjectedIndices - 1);

	for (unsigned i = 0; i<imax; ++i)
	{
	if (last)
	projector.process(data_[idx0+i*stride]);
	else
	projectInnerLoop2(level+1, idx0+i*stride, projector,
	projectedIndices, nProjectedIndices);
	}
	}

	template<typename Numeric, unsigned StackLen, unsigned StackDim>
	template <typename Num2, unsigned Len2, unsigned Dim2,
	typename Num3, class Op>
	void ArrayND<Numeric,StackLen,StackDim>::projectLoop2(
	const unsigned level, const unsigned long idx0,
	const unsigned level1, const unsigned long idx1,
	ArrayND<Num2,Len2,Dim2>* projection,
	AbsVisitor<Numeric,Num3>& projector,
	const unsigned *projectedIndices,
	const unsigned nProjectedIndices, Op fcn) const
	{
	if (level == dim_)
	{
	assert(level1 == projection->dim_);
	projector.clear();
	projectInnerLoop2(0U, idx0, projector,
	projectedIndices, nProjectedIndices);
	if (projection->dim_)
	fcn(projection->data_[idx1], projector.result());
	else
	fcn(projection->localData_[0], projector.result());
	}
	else
	{
	bool fixed = false;
	for (unsigned j=0; j<nProjectedIndices; ++j)
	if (projectedIndices[j] == level)
	{
	fixed = true;
	break;
	}
	if (fixed)
	{
	projectLoop2(level+1, idx0, level1, idx1,
	projection, projector,
	projectedIndices, nProjectedIndices, fcn);
	}
	else
	{
	const unsigned imax = shape_[level];
	const unsigned long stride = strides_[level];
	const unsigned long stride2 = projection->strides_[level1];
	for (unsigned i = 0; i<imax; ++i)
	projectLoop2(level+1, idx0+i*stride,
	level1+1, idx1+i*stride2,
	projection, projector,
	projectedIndices, nProjectedIndices, fcn);
	}
	}
	}

	template<typename Numeric, unsigned StackLen, unsigned StackDim>
	template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
	void ArrayND<Numeric,StackLen,StackDim>::project(
	ArrayND<Num2,Len2,Dim2>* projection,
	AbsVisitor<Numeric,Num3>& projector,
	const unsigned *projectedIndices,
	const unsigned nProjectedIndices) const
	{
	assert(projection);
	verifyProjectionCompatibility(*projection, projectedIndices,
	nProjectedIndices);
	projectLoop2(0U, 0UL, 0U, 0UL, projection,
	projector, projectedIndices, nProjectedIndices,
	scast_assign_left<Num2,Num3>());
	}

	template<typename Numeric, unsigned StackLen, unsigned StackDim>
	template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
	void ArrayND<Numeric,StackLen,StackDim>::addToProjection(
	ArrayND<Num2,Len2,Dim2>* projection,
	AbsVisitor<Numeric,Num3>& projector,
	const unsigned *projectedIndices,
	const unsigned nProjectedIndices) const
	{
	assert(projection);
	verifyProjectionCompatibility(*projection, projectedIndices,
	nProjectedIndices);
	projectLoop2(0U, 0UL, 0U, 0UL, projection,
	projector, projectedIndices, nProjectedIndices,
	scast_pluseq_left<Num2,Num3>());
	}

	template<typename Numeric, unsigned StackLen, unsigned StackDim>
	template <typename Num2, unsigned Len2, unsigned Dim2, typename Num3>
	void ArrayND<Numeric,StackLen,StackDim>::subtractFromProjection(
	ArrayND<Num2,Len2,Dim2>* projection,
	AbsVisitor<Numeric,Num3>& projector,
	const unsigned *projectedIndices,
	const unsigned nProjectedIndices) const
	{
	assert(projection);
	verifyProjectionCompatibility(*projection, projectedIndices,
	nProjectedIndices);
	projectLoop2(0U, 0UL, 0U, 0UL, projection,
	projector, projectedIndices, nProjectedIndices,
	scast_minuseq_left<Num2,Num3>());
	}

	template<typename Numeric, unsigned StackLen, unsigned StackDim>
	template <typename Num2, class Functor>
	void ArrayND<Numeric,StackLen,StackDim>::scaleBySliceInnerLoop(
	const unsigned level, const unsigned long idx0,
	Num2& scale, const unsigned *projectedIndices,
	const unsigned nProjectedIndices, Functor binaryFunct)
	{
	const unsigned idx = projectedIndices[level];
	const unsigned imax = shape_[idx];
	const unsigned long stride = strides_[idx];

	if (level == nProjectedIndices - 1)
	{
	Numeric* data = data_ + idx0;
	for (unsigned i = 0; i<imax; ++i)
	binaryFunct(data[i*stride], scale);
	}
	else
	for (unsigned i = 0; i<imax; ++i)
	scaleBySliceInnerLoop(level+1, idx0+i*stride, scale,
	projectedIndices, nProjectedIndices,
	binaryFunct);
	}

	template<typename Numeric, unsigned StackLen, unsigned StackDim>
	template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
	void ArrayND<Numeric,StackLen,StackDim>::scaleBySliceLoop(
	unsigned level, unsigned long idx0,
	unsigned level1, unsigned long idx1,
	ArrayND<Num2,Len2,Dim2>& slice,
	const unsigned *projectedIndices,
	const unsigned nProjectedIndices,
	Functor binaryFunct)
	{
	if (level == dim_)
	{
	assert(level1 == slice.dim_);
	Num2& scaleFactor = slice.dim_ ? slice.data_[idx1] :
	slice.localData_[0];
	scaleBySliceInnerLoop(0U, idx0, scaleFactor, projectedIndices,
	nProjectedIndices, binaryFunct);
	}
	else
	{
	bool fixed = false;
	for (unsigned j=0; j<nProjectedIndices; ++j)
	if (projectedIndices[j] == level)
	{
	fixed = true;
	break;
	}
	if (fixed)
	{
	scaleBySliceLoop(level+1, idx0, level1, idx1, slice,
	projectedIndices, nProjectedIndices,
	binaryFunct);
	}
	else
	{
	const unsigned imax = shape_[level];
	const unsigned long stride = strides_[level];
	const unsigned long stride2 = slice.strides_[level1];
	for (unsigned i = 0; i<imax; ++i)
	scaleBySliceLoop(level+1, idx0+i*stride, level1+1,
	idx1+i*stride2, slice, projectedIndices,
	nProjectedIndices, binaryFunct);
	}
	}
	}

	template<typename Numeric, unsigned StackLen, unsigned StackDim>
	template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
	void ArrayND<Numeric,StackLen,StackDim>::jointScan(
	ArrayND<Num2, Len2, Dim2>& r, Functor binaryFunct)
	{
	if (!isShapeCompatible(r)) throw std::invalid_argument(
	"In npstat::ArrayND::jointScan: incompatible argument array shape");
	if (dim_)
	for (unsigned long i=0; i<len_; ++i)
	binaryFunct(data_[i], r.data_[i]);
	else
	binaryFunct(localData_[0], r.localData_[0]);
	}

	template<typename Numeric, unsigned StackLen, unsigned StackDim>
	template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
	void ArrayND<Numeric,StackLen,StackDim>::applySlice(
	ArrayND<Num2,Len2,Dim2>& slice,
	const unsigned *fixedIndices, const unsigned nFixedIndices,
	Functor binaryFunct)
	{
	if (nFixedIndices)
	{
	verifyProjectionCompatibility(slice, fixedIndices, nFixedIndices);
	if (slice.dim_ == 0U)
	for (unsigned long i=0; i<len_; ++i)
	binaryFunct(data_[i], slice.localData_[0]);
	else
	scaleBySliceLoop(0U, 0UL, 0U, 0UL, slice,
	fixedIndices, nFixedIndices, binaryFunct);
	}
	else
	jointScan(slice, binaryFunct);
	}

	template<typename Numeric, unsigned StackLen, unsigned StackDim>
	inline bool ArrayND<Numeric,StackLen,StackDim>::isCompatible(
	const unsigned* shape, unsigned dim) const
	{
	if (dim)
	assert(shape);
	if (!shapeIsKnown_)
	return false;
	if (dim_ != dim)
	return false;
	if (dim_)
	{
	for (unsigned i=0; i<dim_; ++i)
	if (shape_[i] != shape[i])
	return false;
	}
	else
	assert(len_ == 1UL);
	return true;
	}

	template<typename Numeric, unsigned StackLen, unsigned StackDim>
	inline bool ArrayND<Numeric,StackLen,StackDim>::isCompatible(
	const ArrayShape& shape) const
	{
	const unsigned sz = shape.size();
	return isCompatible(sz ? &shape[0] : 0, sz);
	}

	template<typename Numeric, unsigned StackLen, unsigned StackDim>
	template <typename Num2, unsigned Len2, unsigned Dim2>
	inline bool ArrayND<Numeric,StackLen,StackDim>::isShapeCompatible(
	const ArrayND<Num2,Len2,Dim2>& r) const
	{
	if (!shapeIsKnown_)
	return false;
	if (!r.shapeIsKnown_)
	return false;
	if (dim_ != r.dim_)
	return false;
	if (len_ != r.len_)
	return false;
	if (dim_)
	{
	assert(shape_);
	assert(r.shape_);
	for (unsigned i=0; i<dim_; ++i)
	if (shape_[i] != r.shape_[i])
	return false;
	}
	else
	assert(len_ == 1UL);
	return true;
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template<typename Num2, typename Integer>
	void ArrayND<Numeric,Len,Dim>::processSubrangeLoop(
	const unsigned level, unsigned long idx0,
	unsigned* currentIndex,
	AbsArrayProjector<Numeric,Num2>& f,
	const BoxND<Integer>& subrange) const
	{
	// Deal with possible negative limits first
	const Interval<Integer>& levelRange(subrange[level]);
	long long int iminl = static_cast<long long int>(levelRange.min());
	if (iminl < 0LL) iminl = 0LL;
	long long int imaxl = static_cast<long long int>(levelRange.max());
	if (imaxl < 0LL) imaxl = 0LL;

	// Now deal with possible out-of-range limits
	const unsigned imin = static_cast<unsigned>(iminl);
	unsigned imax = static_cast<unsigned>(imaxl);
	if (imax > shape_[level])
	imax = shape_[level];

	if (level == dim_ - 1)
	{
	idx0 += imin;
	for (unsigned i=imin; i<imax; ++i, ++idx0)
	{
	currentIndex[level] = i;
	f.process(currentIndex, dim_, idx0, data_[idx0]);
	}
	}
	else
	{
	const unsigned long stride = strides_[level];
	idx0 += imin*stride;
	for (unsigned i=imin; i<imax; ++i)
	{
	currentIndex[level] = i;
	processSubrangeLoop(level+1U, idx0, currentIndex, f, subrange);
	idx0 += stride;
	}
	}
	}

	template<typename Numeric, unsigned Len, unsigned StackDim>
	template <typename Num2, typename Integer>
	void ArrayND<Numeric,Len,StackDim>::processSubrange(
	AbsArrayProjector<Numeric,Num2>& f,
	const BoxND<Integer>& subrange) const
	{
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"processSubrange\"");
	if (!dim_) throw std::invalid_argument(
	"npstat::ArrayND::processSubrange method "
	"can not be used with array of 0 rank");
	if (dim_ != subrange.dim()) throw std::invalid_argument(
	"In npstat::ArrayND::processSubrange: incompatible subrange rank");
	unsigned ibuf[StackDim];
	unsigned* buf = makeBuffer(dim_, ibuf, StackDim);
	for (unsigned i=0; i<dim_; ++i)
	buf[i] = 0U;
	processSubrangeLoop(0U, 0UL, buf, f, subrange);
	destroyBuffer(buf, ibuf);
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template<typename Num2>
	inline ArrayND<Numeric,Len,Dim>& ArrayND<Numeric,Len,Dim>::setData(
	const Num2* data, const unsigned long dataLength)
	{
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"setData\"");
	if (dataLength != len_) throw std::invalid_argument(
	"In npstat::ArrayND::setData: incompatible input data length");
	copyBuffer(data_, data, dataLength);
	return *this;
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	void ArrayND<Numeric,Len,Dim>::buildStrides()
	{
	assert(dim_);
	if (strides_ == 0)
	strides_ = makeBuffer(dim_, localStrides_, Dim);
	strides_[dim_ - 1] = 1UL;
	for (unsigned j=dim_ - 1; j>0; --j)
	strides_[j - 1] = strides_[j]*shape_[j];
	}

	template<typename Numeric, unsigned StackLen, unsigned StackDim>
	inline ArrayND<Numeric,StackLen,StackDim>::ArrayND()
	: data_(0), strides_(0), shape_(0), len_(0UL), dim_(0U),
	shapeIsKnown_(false), dataIsExternal_(false)
	{
	localData_[0] = Numeric();
	data_ = localData_;
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	ArrayND<Numeric,Len,Dim>::ArrayND(const ArrayND& r)
	: data_(0), strides_(0), shape_(0),
	len_(r.len_), dim_(r.dim_), shapeIsKnown_(r.shapeIsKnown_),
	dataIsExternal_(false)
	{
	if (dim_)
	{
	// Copy the shape
	shape_ = makeBuffer(dim_, localShape_, Dim);
	copyBuffer(shape_, r.shape_, dim_);

	// Copy the strides
	strides_ = makeBuffer(dim_, localStrides_, Dim);
	copyBuffer(strides_, r.strides_, dim_);

	// Copy the data
	data_ = makeBuffer(len_, localData_, Len);
	copyBuffer(data_, r.data_, len_);
	}
	else
	{
	assert(len_ == 1UL);
	localData_[0] = r.localData_[0];
	data_ = localData_;
	}
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	ArrayND<Numeric,Len,Dim>::ArrayND(ArrayND&& r)
	: data_(0), strides_(0), shape_(0),
	len_(r.len_), dim_(r.dim_), shapeIsKnown_(r.shapeIsKnown_),
	dataIsExternal_(r.dataIsExternal_)
	{
	if (dim_)
	{
	// Copy the shape
	if (r.shape_ == r.localShape_)
	{
	shape_ = makeBuffer(dim_, localShape_, Dim);
	copyBuffer(shape_, r.shape_, dim_);
	}
	else
	{
	shape_ = r.shape_;
	r.shape_ = 0;
	}

	// Copy the strides
	if (r.strides_ == r.localStrides_)
	{
	strides_ = makeBuffer(dim_, localStrides_, Dim);
	copyBuffer(strides_, r.strides_, dim_);
	}
	else
	{
	strides_ = r.strides_;
	r.strides_ = 0;
	}

	// Copy the data
	if (dataIsExternal_)
	data_ = r.data_;
	else
	{
	if (r.data_ == r.localData_)
	{
	data_ = makeBuffer(len_, localData_, Len);
	copyBuffer(data_, r.data_, len_);
	}
	else
	{
	data_ = r.data_;
	r.data_ = 0;
	}
	}
	}
	else
	{
	assert(len_ == 1UL);
	assert(!dataIsExternal_);
	localData_[0] = r.localData_[0];
	data_ = localData_;
	}
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template <typename Num2, unsigned Len2, unsigned Dim2>
	ArrayND<Numeric,Len,Dim>::ArrayND(const ArrayND<Num2, Len2, Dim2>& r)
	: data_(0), strides_(0), shape_(0),
	len_(r.len_), dim_(r.dim_), shapeIsKnown_(r.shapeIsKnown_),
	dataIsExternal_(false)
	{
	if (dim_)
	{
	// Copy the shape
	shape_ = makeBuffer(dim_, localShape_, Dim);
	copyBuffer(shape_, r.shape_, dim_);

	// Copy the strides
	strides_ = makeBuffer(dim_, localStrides_, Dim);
	copyBuffer(strides_, r.strides_, dim_);

	// Copy the data
	data_ = makeBuffer(len_, localData_, Len);
	copyBuffer(data_, r.data_, len_);
	}
	else
	{
	assert(len_ == 1UL);
	localData_[0] = static_cast<Numeric>(r.localData_[0]);
	data_ = localData_;
	}
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template<typename Num2, unsigned Len2, unsigned Dim2, class Functor>
	ArrayND<Numeric,Len,Dim>::ArrayND(const ArrayND<Num2, Len2, Dim2>& r,
	Functor f)
	: data_(0), strides_(0), shape_(0),
	len_(r.len_), dim_(r.dim_), shapeIsKnown_(r.shapeIsKnown_),
	dataIsExternal_(false)
	{
	if (dim_)
	{
	// Copy the shape
	shape_ = makeBuffer(dim_, localShape_, Dim);
	copyBuffer(shape_, r.shape_, dim_);

	// Copy the strides
	strides_ = makeBuffer(dim_, localStrides_, Dim);
	copyBuffer(strides_, r.strides_, dim_);

	// Copy the data
	data_ = makeBuffer(len_, localData_, Len);
	for (unsigned long i=0; i<len_; ++i)
	data_[i] = static_cast<Numeric>(f(r.data_[i]));
	}
	else
	{
	assert(len_ == 1UL);
	localData_[0] = static_cast<Numeric>(f(r.localData_[0]));
	data_ = localData_;
	}
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
	void ArrayND<Numeric,Len,Dim>::copyRangeLoopFunct(
	const unsigned level, unsigned long idx0,
	unsigned long idx1,
	const ArrayND<Num2, Len2, Dim2>& r,
	const ArrayRange& range, Functor f)
	{
	const unsigned imax = shape_[level];
	if (level == dim_ - 1)
	{
	Numeric* to = data_ + idx0;
	const Num2* from = r.data_ + (idx1 + range[level].min());
	for (unsigned i=0; i<imax; ++i)
	to++ = static_cast<Numeric>(f(from++));
	}
	else
	{
	const unsigned long fromstride = r.strides_[level];
	const unsigned long tostride = strides_[level];
	idx1 += range[level].min()*fromstride;
	for (unsigned i=0; i<imax; ++i)
	{
	copyRangeLoopFunct(level+1, idx0, idx1, r, range, f);
	idx0 += tostride;
	idx1 += fromstride;
	}
	}
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template <typename Num2, unsigned Len2, unsigned Dim2>
	ArrayND<Numeric,Len,Dim>::ArrayND(
	const ArrayND<Num2, Len2, Dim2>& r, const ArrayRange& range)
	: data_(0), strides_(0), shape_(0),
	len_(r.len_), dim_(r.dim_), shapeIsKnown_(r.shapeIsKnown_),
	dataIsExternal_(false)
	{
	if (!range.isCompatible(r.shape_, r.dim_))
	throw std::invalid_argument(
	"In npstat::ArrayND subrange constructor: invalid subrange");
	if (dim_)
	{
	len_ = range.rangeSize();
	if (!len_)
	throw std::invalid_argument(
	"In npstat::ArrayND subrange constructor: empty subrange");

	// Figure out the shape
	shape_ = makeBuffer(dim_, localShape_, Dim);
	range.rangeLength(shape_, dim_);

	// Figure out the strides
	buildStrides();

	// Allocate the data array
	data_ = makeBuffer(len_, localData_, Len);

	// Copy the data
	if (dim_ > CHAR_BIT*sizeof(unsigned long))
	throw std::invalid_argument(
	"In npstat::ArrayND subrange constructor: "
	"input array rank is too large");
	unsigned lolim[CHAR_BIT*sizeof(unsigned long)];
	range.lowerLimits(lolim, dim_);
	unsigned toBuf[CHAR_BIT*sizeof(unsigned long)];
	clearBuffer(toBuf, dim_);
	(const_cast<ArrayND<Num2, Len2, Dim2>&>(r)).commonSubrangeLoop(
	0U, 0UL, 0UL, lolim, shape_, toBuf, *this,
	scast_assign_right<Num2,Numeric>());
	}
	else
	{
	assert(len_ == 1UL);
	localData_[0] = static_cast<Numeric>(r.localData_[0]);
	data_ = localData_;
	}
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
	ArrayND<Numeric,Len,Dim>::ArrayND(
	const ArrayND<Num2, Len2, Dim2>& r, const ArrayRange& range,
	Functor f)
	: data_(0), strides_(0), shape_(0),
	len_(r.len_), dim_(r.dim_), shapeIsKnown_(r.shapeIsKnown_),
	dataIsExternal_(false)
	{
	if (!range.isCompatible(r.shape_, r.dim_))
	throw std::invalid_argument(
	"In npstat::ArrayND transforming subrange constructor: "
	"incompatible subrange");
	if (dim_)
	{
	len_ = range.rangeSize();
	if (!len_)
	throw std::invalid_argument(
	"In npstat::ArrayND transforming subrange constructor: "
	"empty subrange");

	// Figure out the shape
	shape_ = makeBuffer(dim_, localShape_, Dim);
	for (unsigned i=0; i<dim_; ++i)
	shape_[i] = range[i].length();

	// Figure out the strides
	buildStrides();

	// Allocate the data array
	data_ = makeBuffer(len_, localData_, Len);

	// Transform the data
	copyRangeLoopFunct(0U, 0UL, 0UL, r, range, f);
	}
	else
	{
	assert(len_ == 1UL);
	localData_[0] = static_cast<Numeric>(f(r.localData_[0]));
	data_ = localData_;
	}
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	ArrayND<Numeric,Len,Dim>::ArrayND(const ArrayShape& sh)
	: data_(0), strides_(0), shape_(0), len_(1UL),
	shapeIsKnown_(true), dataIsExternal_(false)
	{
	const unsigned sz = sh.size();
	buildFromShapePtr(sz ? &sh[0] : 0, sz);
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	ArrayND<Numeric,Len,Dim>::ArrayND(const unsigned* sizes,
	const unsigned dim)
	: data_(0), strides_(0), shape_(0), len_(1UL),
	shapeIsKnown_(true), dataIsExternal_(false)
	{
	buildFromShapePtr(sizes, dim);
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	ArrayND<Numeric,Len,Dim>& ArrayND<Numeric,Len,Dim>::reshape(
	const unsigned* newShape, const unsigned dim)
	{
	if (!isCompatible(newShape, dim) \|\| dataIsExternal_)
	{
	uninitialize();
	shapeIsKnown_ = true;
	len_ = 1UL;
	buildFromShapePtr(newShape, dim);
	}
	return *this;
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template<typename Num2, unsigned Len2, unsigned Dim2>
	ArrayND<Numeric,Len,Dim>& ArrayND<Numeric,Len,Dim>::reshape(
	const ArrayND<Num2,Len2,Dim2>& r)
	{
	if (!r.shapeIsKnown_) throw std::invalid_argument(
	"In npstat::ArrayND::reshape: uninitialized argument array");
	return reshape(r.shape_, r.dim_);
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	ArrayND<Numeric,Len,Dim>& ArrayND<Numeric,Len,Dim>::reshape(
	const ArrayShape& newShape)
	{
	const unsigned sz = newShape.size();
	return reshape(sz ? &newShape[0] : 0, sz);
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	ArrayND<Numeric,Len,Dim>& ArrayND<Numeric,Len,Dim>::reshape(
	unsigned n0)
	{
	unsigned buf[1];
	buf[0] = n0;
	return reshape(buf, 1U);
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	ArrayND<Numeric,Len,Dim>& ArrayND<Numeric,Len,Dim>::reshape(
	unsigned n0, unsigned n1)
	{
	unsigned buf[2];
	buf[0] = n0;
	buf[1] = n1;
	return reshape(buf, 2U);
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	ArrayND<Numeric,Len,Dim>& ArrayND<Numeric,Len,Dim>::reshape(
	unsigned n0, unsigned n1, unsigned n2)
	{
	unsigned buf[3];
	buf[0] = n0;
	buf[1] = n1;
	buf[2] = n2;
	return reshape(buf, 3U);
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	ArrayND<Numeric,Len,Dim>& ArrayND<Numeric,Len,Dim>::reshape(
	unsigned n0, unsigned n1, unsigned n2, unsigned n3)
	{
	unsigned buf[4];
	buf[0] = n0;
	buf[1] = n1;
	buf[2] = n2;
	buf[3] = n3;
	return reshape(buf, 4U);
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	ArrayND<Numeric,Len,Dim>& ArrayND<Numeric,Len,Dim>::reshape(
	unsigned n0, unsigned n1, unsigned n2, unsigned n3, unsigned n4)
	{
	unsigned buf[5];
	buf[0] = n0;
	buf[1] = n1;
	buf[2] = n2;
	buf[3] = n3;
	buf[4] = n4;
	return reshape(buf, 5U);
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	ArrayND<Numeric,Len,Dim>& ArrayND<Numeric,Len,Dim>::reshape(
	unsigned n0, unsigned n1, unsigned n2, unsigned n3, unsigned n4,
	unsigned n5)
	{
	unsigned buf[6];
	buf[0] = n0;
	buf[1] = n1;
	buf[2] = n2;
	buf[3] = n3;
	buf[4] = n4;
	buf[5] = n5;
	return reshape(buf, 6U);
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	ArrayND<Numeric,Len,Dim>& ArrayND<Numeric,Len,Dim>::reshape(
	unsigned n0, unsigned n1, unsigned n2, unsigned n3, unsigned n4,
	unsigned n5, unsigned n6)
	{
	unsigned buf[7];
	buf[0] = n0;
	buf[1] = n1;
	buf[2] = n2;
	buf[3] = n3;
	buf[4] = n4;
	buf[5] = n5;
	buf[6] = n6;
	return reshape(buf, 7U);
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	ArrayND<Numeric,Len,Dim>& ArrayND<Numeric,Len,Dim>::reshape(
	unsigned n0, unsigned n1, unsigned n2, unsigned n3, unsigned n4,
	unsigned n5, unsigned n6, unsigned n7)
	{
	unsigned buf[8];
	buf[0] = n0;
	buf[1] = n1;
	buf[2] = n2;
	buf[3] = n3;
	buf[4] = n4;
	buf[5] = n5;
	buf[6] = n6;
	buf[7] = n7;
	return reshape(buf, 8U);
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	ArrayND<Numeric,Len,Dim>& ArrayND<Numeric,Len,Dim>::reshape(
	unsigned n0, unsigned n1, unsigned n2, unsigned n3, unsigned n4,
	unsigned n5, unsigned n6, unsigned n7, unsigned n8)
	{
	unsigned buf[9];
	buf[0] = n0;
	buf[1] = n1;
	buf[2] = n2;
	buf[3] = n3;
	buf[4] = n4;
	buf[5] = n5;
	buf[6] = n6;
	buf[7] = n7;
	buf[8] = n8;
	return reshape(buf, 9U);
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	ArrayND<Numeric,Len,Dim>& ArrayND<Numeric,Len,Dim>::reshape(
	unsigned n0, unsigned n1, unsigned n2, unsigned n3, unsigned n4,
	unsigned n5, unsigned n6, unsigned n7, unsigned n8, unsigned n9)
	{
	unsigned buf[10];
	buf[0] = n0;
	buf[1] = n1;
	buf[2] = n2;
	buf[3] = n3;
	buf[4] = n4;
	buf[5] = n5;
	buf[6] = n6;
	buf[7] = n7;
	buf[8] = n8;
	buf[9] = n9;
	return reshape(buf, 10U);
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	ArrayND<Numeric,Len,Dim>::ArrayND(const unsigned n0)
	: data_(0), strides_(0), shape_(0), len_(1UL),
	shapeIsKnown_(true), dataIsExternal_(false)
	{
	const unsigned dim = 1U;
	unsigned sizes[dim];
	sizes[0] = n0;
	buildFromShapePtr(sizes, dim);
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	ArrayND<Numeric,Len,Dim>::ArrayND(const unsigned n0,
	const unsigned n1)
	: data_(0), strides_(0), shape_(0), len_(1UL),
	shapeIsKnown_(true), dataIsExternal_(false)
	{
	const unsigned dim = 2U;
	unsigned sizes[dim];
	sizes[0] = n0;
	sizes[1] = n1;
	buildFromShapePtr(sizes, dim);
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	ArrayND<Numeric,Len,Dim>::ArrayND(const unsigned n0,
	const unsigned n1,
	const unsigned n2)
	: data_(0), strides_(0), shape_(0), len_(1UL),
	shapeIsKnown_(true), dataIsExternal_(false)
	{
	const unsigned dim = 3U;
	unsigned sizes[dim];
	sizes[0] = n0;
	sizes[1] = n1;
	sizes[2] = n2;
	buildFromShapePtr(sizes, dim);
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	ArrayND<Numeric,Len,Dim>::ArrayND(const unsigned n0,
	const unsigned n1,
	const unsigned n2,
	const unsigned n3)
	: data_(0), strides_(0), shape_(0), len_(1UL),
	shapeIsKnown_(true), dataIsExternal_(false)
	{
	const unsigned dim = 4U;
	unsigned sizes[dim];
	sizes[0] = n0;
	sizes[1] = n1;
	sizes[2] = n2;
	sizes[3] = n3;
	buildFromShapePtr(sizes, dim);
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	ArrayND<Numeric,Len,Dim>::ArrayND(const unsigned n0,
	const unsigned n1,
	const unsigned n2,
	const unsigned n3,
	const unsigned n4)
	: data_(0), strides_(0), shape_(0), len_(1UL),
	shapeIsKnown_(true), dataIsExternal_(false)
	{
	const unsigned dim = 5U;
	unsigned sizes[dim];
	sizes[0] = n0;
	sizes[1] = n1;
	sizes[2] = n2;
	sizes[3] = n3;
	sizes[4] = n4;
	buildFromShapePtr(sizes, dim);
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	ArrayND<Numeric,Len,Dim>::ArrayND(const unsigned n0,
	const unsigned n1,
	const unsigned n2,
	const unsigned n3,
	const unsigned n4,
	const unsigned n5)
	: data_(0), strides_(0), shape_(0), len_(1UL),
	shapeIsKnown_(true), dataIsExternal_(false)
	{
	const unsigned dim = 6U;
	unsigned sizes[dim];
	sizes[0] = n0;
	sizes[1] = n1;
	sizes[2] = n2;
	sizes[3] = n3;
	sizes[4] = n4;
	sizes[5] = n5;
	buildFromShapePtr(sizes, dim);
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	ArrayND<Numeric,Len,Dim>::ArrayND(const unsigned n0,
	const unsigned n1,
	const unsigned n2,
	const unsigned n3,
	const unsigned n4,
	const unsigned n5,
	const unsigned n6)
	: data_(0), strides_(0), shape_(0), len_(1UL),
	shapeIsKnown_(true), dataIsExternal_(false)
	{
	const unsigned dim = 7U;
	unsigned sizes[dim];
	sizes[0] = n0;
	sizes[1] = n1;
	sizes[2] = n2;
	sizes[3] = n3;
	sizes[4] = n4;
	sizes[5] = n5;
	sizes[6] = n6;
	buildFromShapePtr(sizes, dim);
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	ArrayND<Numeric,Len,Dim>::ArrayND(const unsigned n0,
	const unsigned n1,
	const unsigned n2,
	const unsigned n3,
	const unsigned n4,
	const unsigned n5,
	const unsigned n6,
	const unsigned n7)
	: data_(0), strides_(0), shape_(0), len_(1UL),
	shapeIsKnown_(true), dataIsExternal_(false)
	{
	const unsigned dim = 8U;
	unsigned sizes[dim];
	sizes[0] = n0;
	sizes[1] = n1;
	sizes[2] = n2;
	sizes[3] = n3;
	sizes[4] = n4;
	sizes[5] = n5;
	sizes[6] = n6;
	sizes[7] = n7;
	buildFromShapePtr(sizes, dim);
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	ArrayND<Numeric,Len,Dim>::ArrayND(const unsigned n0,
	const unsigned n1,
	const unsigned n2,
	const unsigned n3,
	const unsigned n4,
	const unsigned n5,
	const unsigned n6,
	const unsigned n7,
	const unsigned n8)
	: data_(0), strides_(0), shape_(0), len_(1UL),
	shapeIsKnown_(true), dataIsExternal_(false)
	{
	const unsigned dim = 9U;
	unsigned sizes[dim];
	sizes[0] = n0;
	sizes[1] = n1;
	sizes[2] = n2;
	sizes[3] = n3;
	sizes[4] = n4;
	sizes[5] = n5;
	sizes[6] = n6;
	sizes[7] = n7;
	sizes[8] = n8;
	buildFromShapePtr(sizes, dim);
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	ArrayND<Numeric,Len,Dim>::ArrayND(const unsigned n0,
	const unsigned n1,
	const unsigned n2,
	const unsigned n3,
	const unsigned n4,
	const unsigned n5,
	const unsigned n6,
	const unsigned n7,
	const unsigned n8,
	const unsigned n9)
	: data_(0), strides_(0), shape_(0), len_(1UL),
	shapeIsKnown_(true), dataIsExternal_(false)
	{
	const unsigned dim = 10U;
	unsigned sizes[dim];
	sizes[0] = n0;
	sizes[1] = n1;
	sizes[2] = n2;
	sizes[3] = n3;
	sizes[4] = n4;
	sizes[5] = n5;
	sizes[6] = n6;
	sizes[7] = n7;
	sizes[8] = n8;
	sizes[9] = n9;
	buildFromShapePtr(sizes, dim);
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template<typename Num1, unsigned Len1, unsigned Dim1,
	typename Num2, unsigned Len2, unsigned Dim2>
	void ArrayND<Numeric,Len,Dim>::outerProductLoop(
	const unsigned level, unsigned long idx0,
	unsigned long idx1, unsigned long idx2,
	const ArrayND<Num1, Len1, Dim1>& a1,
	const ArrayND<Num2, Len2, Dim2>& a2)
	{
	const unsigned imax = shape_[level];
	if (level == dim_ - 1)
	{
	for (unsigned i=0; i<imax; ++i)
	data_[idx0 + i] = a1.data_[idx1]*a2.data_[idx2 + i];
	}
	else
	{
	for (unsigned i=0; i<imax; ++i)
	{
	outerProductLoop(level+1, idx0, idx1, idx2, a1, a2);
	idx0 += strides_[level];
	if (level < a1.dim_)
	idx1 += a1.strides_[level];
	else
	idx2 += a2.strides_[level - a1.dim_];
	}
	}
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template<typename Num1, unsigned Len1, unsigned Dim1,
	typename Num2, unsigned Len2, unsigned Dim2>
	ArrayND<Numeric,Len,Dim>::ArrayND(const ArrayND<Num1, Len1, Dim1>& a1,
	const ArrayND<Num2, Len2, Dim2>& a2)
	: data_(0), strides_(0), shape_(0),
	len_(1UL), dim_(a1.dim_ + a2.dim_),
	shapeIsKnown_(true), dataIsExternal_(false)
	{
	if (!(a1.shapeIsKnown_ && a2.shapeIsKnown_))
	throw std::invalid_argument(
	"In npstat::ArrayND outer product constructor: "
	"uninitialized argument array");
	if (dim_)
	{
	shape_ = makeBuffer(dim_, localShape_, Dim);
	copyBuffer(shape_, a1.shape_, a1.dim_);
	copyBuffer(shape_+a1.dim_, a2.shape_, a2.dim_);

	for (unsigned i=0; i<dim_; ++i)
	{
	assert(shape_[i]);
	len_ *= shape_[i];
	}

	// Figure out the array strides
	buildStrides();

	// Allocate the data array
	data_ = makeBuffer(len_, localData_, Len);

	// Fill the data array
	if (a1.dim_ == 0)
	{
	for (unsigned long i=0; i<len_; ++i)
	data_[i] = a1.localData_[0] * a2.data_[i];
	}
	else if (a2.dim_ == 0)
	{
	for (unsigned long i=0; i<len_; ++i)
	data_[i] = a1.data_[i] * a2.localData_[0];
	}
	else
	outerProductLoop(0U, 0UL, 0UL, 0UL, a1, a2);
	}
	else
	{
	localData_[0] = a1.localData_[0] * a2.localData_[0];
	data_ = localData_;
	}
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline ArrayND<Numeric,Len,Dim>::~ArrayND()
	{
	if (!dataIsExternal_)
	destroyBuffer(data_, localData_);
	destroyBuffer(strides_, localStrides_);
	destroyBuffer(shape_, localShape_);
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	ArrayND<Numeric,Len,Dim>&
	ArrayND<Numeric,Len,Dim>::operator=(const ArrayND& r)
	{
	if (this == &r)
	return *this;
	if (shapeIsKnown_)
	{
	if (!r.shapeIsKnown_) throw std::invalid_argument(
	"In npstat::ArrayND assignment operator: "
	"uninitialized argument array");
	if (!isShapeCompatible(r)) throw std::invalid_argument(
	"In npstat::ArrayND assignment operator: "
	"incompatible argument array shape");
	if (dim_)
	{
	// It is possible that both arrays manage the same
	// external memory. This is the reason for the "if" below.
	if (data_ != r.data_)
	copyBuffer(data_, r.data_, len_);
	}
	else
	localData_[0] = r.localData_[0];
	}
	else
	{
	// This object is uninitialized. If the object on the
	// right is itself initialized, make an in-place copy.
	if (r.shapeIsKnown_)
	new (this) ArrayND(r);
	}
	return *this;
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	ArrayND<Numeric,Len,Dim>&
	ArrayND<Numeric,Len,Dim>::operator=(ArrayND&& r)
	{
	if (this == &r)
	return *this;
	if (shapeIsKnown_)
	{
	if (!r.shapeIsKnown_) throw std::invalid_argument(
	"In npstat::ArrayND assignment operator: "
	"uninitialized argument array");
	if (!isShapeCompatible(r)) throw std::invalid_argument(
	"In npstat::ArrayND assignment operator: "
	"incompatible argument array shape");
	if (dim_)
	{
	if (!dataIsExternal_)
	{
	destroyBuffer(data_, localData_);
	data_ = 0;
	}
	dataIsExternal_ = r.dataIsExternal_;
	if (dataIsExternal_)
	data_ = r.data_;
	else
	{
	if (r.data_ == r.localData_)
	{
	data_ = makeBuffer(len_, localData_, Len);
	copyBuffer(data_, r.data_, len_);
	}
	else
	{
	data_ = r.data_;
	r.data_ = 0;
	}
	}
	}
	else
	localData_[0] = r.localData_[0];
	}
	else
	{
	// This is where the move assignment is really useful
	if (r.shapeIsKnown_)
	{
	len_ = r.len_;
	dim_ = r.dim_;
	shapeIsKnown_ = true;
	dataIsExternal_ = r.dataIsExternal_;

	// Copy the shape
	if (r.shape_ == r.localShape_)
	{
	shape_ = makeBuffer(dim_, localShape_, Dim);
	copyBuffer(shape_, r.shape_, dim_);
	}
	else
	{
	shape_ = r.shape_;
	r.shape_ = 0;
	}

	// Copy the strides
	if (r.strides_ == r.localStrides_)
	{
	strides_ = makeBuffer(dim_, localStrides_, Dim);
	copyBuffer(strides_, r.strides_, dim_);
	}
	else
	{
	strides_ = r.strides_;
	r.strides_ = 0;
	}

	// Copy the data
	if (dataIsExternal_)
	data_ = r.data_;
	else
	{
	if (r.data_ == r.localData_)
	{
	data_ = makeBuffer(len_, localData_, Len);
	copyBuffer(data_, r.data_, len_);
	}
	else
	{
	data_ = r.data_;
	r.data_ = 0;
	}
	}
	}
	}
	return *this;
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template <typename Num2, unsigned Len2, unsigned Dim2>
	ArrayND<Numeric,Len,Dim>&
	ArrayND<Numeric,Len,Dim>::operator=(const ArrayND<Num2,Len2,Dim2>& r)
	{
	if ((void)this == (void)(&r))
	return *this;
	if (shapeIsKnown_)
	{
	if (!r.shapeIsKnown_) throw std::invalid_argument(
	"In npstat::ArrayND assignment operator: "
	"uninitialized argument array");
	if (!isShapeCompatible(r)) throw std::invalid_argument(
	"In npstat::ArrayND assignment operator: "
	"incompatible argument array shape");
	if (dim_)
	copyBuffer(data_, r.data_, len_);
	else
	localData_[0] = static_cast<Numeric>(r.localData_[0]);
	}
	else
	{
	// This object is uninitialized. If the object on the
	// right is itself initialized, make an in-place copy.
	if (r.shapeIsKnown_)
	new (this) ArrayND(r);
	}
	return *this;
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template <typename Num2, unsigned Len2, unsigned Dim2, class Functor>
	ArrayND<Numeric,Len,Dim>&
	ArrayND<Numeric,Len,Dim>::assign(const ArrayND<Num2,Len2,Dim2>& r,
	Functor f)
	{
	if (shapeIsKnown_)
	{
	if (!r.shapeIsKnown_) throw std::invalid_argument(
	"In npstat::ArrayND::assign: uninitialized argument array");
	if (!isShapeCompatible(r)) throw std::invalid_argument(
	"In npstat::ArrayND::assign: incompatible argument array shape");
	if (dim_)
	for (unsigned long i=0; i<len_; ++i)
	data_[i] = static_cast<Numeric>(f(r.data_[i]));
	else
	localData_[0] = static_cast<Numeric>(f(r.localData_[0]));
	}
	else
	{
	// This object is uninitialized. If the object on the
	// right is itself initialized, build new array in place.
	if (r.shapeIsKnown_)
	new (this) ArrayND(r, f);
	}
	return *this;
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline ArrayShape ArrayND<Numeric,Len,Dim>::shape() const
	{
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"shape\"");
	return ArrayShape(shape_, shape_+dim_);
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline ArrayRange ArrayND<Numeric,Len,Dim>::fullRange() const
	{
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"fullRange\"");
	ArrayRange range;
	if (dim_)
	{
	range.reserve(dim_);
	for (unsigned i=0; i<dim_; ++i)
	range.push_back(Interval<unsigned>(0U, shape_[i]));
	}
	return range;
	}

	template<class Arr>
	bool arrayIsDensity(const Arr& arr)
	{
	typedef typename Arr::value_type Numeric;

	if (!arr.isShapeKnown()) throw std::invalid_argument(
	"In npstat::arrayIsDensity: uninitialized array argument");
	const unsigned dim = arr.rank();
	const Numeric zero = Numeric();
	bool hasPositive = false;
	if (dim)
	{
	const unsigned long len = arr.length();
	for (unsigned long i=0; i<len; ++i)
	{
	// Don't make comparisons whose result can be
	// determined in advance by assuming that Numeric
	// is an unsigned type. Some compilers will
	// complain about it when this template is
	// instantiated with such a type.
	const Numeric& v = arr.linearValue(i);
	if (v == zero)
	continue;
	if (ComplexComparesFalse<Numeric>::less(zero, v))
	hasPositive = true;
	else
	return false;
	}
	}
	else
	{
	const Numeric& v = arr();
	hasPositive = ComplexComparesFalse<Numeric>::less(zero, v);
	}
	return hasPositive;
	}

	template<class Arr>
	bool arrayIsNonNegative(const Arr& arr)
	{
	typedef typename Arr::value_type Numeric;

	if (!arr.isShapeKnown()) throw std::invalid_argument(
	"In npstat::arrayIsNonNegative: uninitialized array argument");
	const unsigned dim = arr.rank();
	const Numeric zero = Numeric();
	if (dim)
	{
	const unsigned long len = arr.length();
	for (unsigned long i=0; i<len; ++i)
	{
	const Numeric& v = arr.linearValue(i);
	if (v == zero)
	continue;
	if (!ComplexComparesFalse<Numeric>::less(zero, v))
	return false;
	}
	}
	else
	{
	const Numeric& v = arr();
	if (v == zero)
	return true;
	if (!ComplexComparesFalse<Numeric>::less(zero, v))
	return false;
	}
	return true;
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	bool ArrayND<Numeric,Len,Dim>::isZero() const
	{
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"isZero\"");
	const Numeric zero = Numeric();
	if (dim_)
	{
	for (unsigned long i=0; i<len_; ++i)
	if (!(data_[i] == zero))
	return false;
	}
	else
	if (!(localData_[0] == zero))
	return false;
	return true;
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	bool ArrayND<Numeric,Len,Dim>::isNonZero() const
	{
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"isNonZero\"");
	const Numeric zero = Numeric();
	if (dim_)
	{
	for (unsigned long i=0; i<len_; ++i)
	if (data_[i] == zero)
	return false;
	}
	else
	if (localData_[0] == zero)
	return false;
	return true;
	}

	template<typename Numeric, unsigned Len, unsigned StackDim>
	ArrayRange ArrayND<Numeric,Len,StackDim>::nonZeroRange() const
	{
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling "
	"method \"nonZeroRange\"");
	ArrayRange range(dim_);
	if (dim_)
	{
	unsigned lolim[CHAR_BIT*sizeof(unsigned)];
	unsigned uplim[CHAR_BIT*sizeof(unsigned)];
	unsigned idx[CHAR_BIT*sizeof(unsigned)];

	for (unsigned i=0; i<dim_; ++i)
	{
	lolim[i] = shape_[i];
	uplim[i] = 0U;
	}
	const Numeric zero = Numeric();
	bool haveNon0 = false;
	for (unsigned long elem=0; elem<len_; ++elem)
	{
	if (!(data_[elem] == zero))
	{
	haveNon0 = true;
	convertLinearIndex(elem, &idx[0], dim_);
	for (unsigned i=0; i<dim_; ++i)
	{
	if (idx[i] < lolim[i])
	lolim[i] = idx[i];
	if (idx[i] >= uplim[i])
	uplim[i] = idx[i] + 1U;
	}
	}
	}

	if (!haveNon0)
	for (unsigned i=0; i<dim_; ++i)
	lolim[i] = 0U;

	for (unsigned i=0; i<dim_; ++i)
	range[i] = Interval<unsigned>(lolim[i], uplim[i]);
	}
	return range;
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	void ArrayND<Numeric,Len,Dim>::convertLinearIndex(
	unsigned long l, unsigned* idx, const unsigned idxLen) const
	{
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling "
	"method \"convertLinearIndex\"");
	if (!dim_) throw std::invalid_argument(
	"npstat::ArrayND::convertLinearIndex method "
	"can not be used with array of 0 rank");
	if (idxLen != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::convertLinearIndex: incompatible index length");
	if (l >= len_) throw std::out_of_range(
	"In npstat::ArrayND::convertLinearIndex: linear index out of range");
	assert(idx);

	for (unsigned i=0; i<dim_; ++i)
	{
	idx[i] = l / strides_[i];
	l -= (idx[i] * strides_[i]);
	}
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	unsigned long ArrayND<Numeric,Len,Dim>::linearIndex(
	const unsigned* index, unsigned idxLen) const
	{
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"linearIndex\"");
	if (!dim_) throw std::invalid_argument(
	"npstat::ArrayND::linearIndex method "
	"can not be used with array of 0 rank");
	if (idxLen != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::linearIndex: incompatible index length");
	assert(index);

	unsigned long idx = 0UL;
	for (unsigned i=0; i<dim_; ++i)
	{
	if (index[i] >= shape_[i])
	throw std::out_of_range(
	"In npstat::ArrayND::linearIndex: index out of range");
	idx += index[i]*strides_[i];
	}
	return idx;
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline Numeric& ArrayND<Numeric,Len,Dim>::value(
	const unsigned *index, const unsigned dim)
	{
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"value\"");
	if (dim != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::value: incompatible index length");
	if (dim)
	{
	assert(index);
	unsigned long idx = 0UL;
	for (unsigned i=0; i<dim_; ++i)
	idx += index[i]*strides_[i];
	return data_[idx];
	}
	else
	return localData_[0];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline const Numeric& ArrayND<Numeric,Len,Dim>::value(
	const unsigned *index, const unsigned dim) const
	{
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"value\"");
	if (dim != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::value: incompatible index length");
	if (dim)
	{
	assert(index);
	unsigned long idx = 0UL;
	for (unsigned i=0; i<dim_; ++i)
	idx += index[i]*strides_[i];
	return data_[idx];
	}
	else
	return localData_[0];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline Numeric& ArrayND<Numeric,Len,Dim>::linearValue(
	const unsigned long index)
	{
	return data_[index];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline const Numeric& ArrayND<Numeric,Len,Dim>::linearValue(
	const unsigned long index) const
	{
	return data_[index];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline Numeric& ArrayND<Numeric,Len,Dim>::linearValueAt(
	const unsigned long index)
	{
	if (index >= len_)
	throw std::out_of_range(
	"In npstat::ArrayND::linearValueAt: linear index out of range");
	return data_[index];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline const Numeric& ArrayND<Numeric,Len,Dim>::linearValueAt(
	const unsigned long index) const
	{
	if (index >= len_)
	throw std::out_of_range(
	"In npstat::ArrayND::linearValueAt: linear index out of range");
	return data_[index];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline unsigned ArrayND<Numeric,Len,Dim>::coordToIndex(
	const double x, const unsigned idim) const
	{
	if (x <= 0.0)
	return 0;
	else if (x >= static_cast<double>(shape_[idim] - 1))
	return shape_[idim] - 1;
	else
	return static_cast<unsigned>(std::floor(x + 0.5));
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline const Numeric& ArrayND<Numeric,Len,Dim>::closest(
	const double *x, const unsigned dim) const
	{
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"closest\"");
	if (dim != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::closest: incompatible data length");
	if (dim)
	{
	assert(x);
	unsigned long idx = 0UL;
	for (unsigned i=0; i<dim_; ++i)
	idx += coordToIndex(x[i], i)*strides_[i];
	return data_[idx];
	}
	else
	return localData_[0];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline Numeric& ArrayND<Numeric,Len,Dim>::closest(
	const double *x, const unsigned dim)
	{
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"closest\"");
	if (dim != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::closest: incompatible data length");
	if (dim)
	{
	assert(x);
	unsigned long idx = 0UL;
	for (unsigned i=0; i<dim_; ++i)
	idx += coordToIndex(x[i], i)*strides_[i];
	return data_[idx];
	}
	else
	return localData_[0];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline const Numeric& ArrayND<Numeric,Len,Dim>::valueAt(
	const unsigned *index, const unsigned dim) const
	{
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"valueAt\"");
	if (dim != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::valueAt: incompatible index length");
	if (dim)
	{
	assert(index);
	unsigned long idx = 0UL;
	for (unsigned i=0; i<dim_; ++i)
	{
	if (index[i] >= shape_[i]) throw std::out_of_range(
	"In npstat::ArrayND::valueAt: index out of range");
	idx += index[i]*strides_[i];
	}
	return data_[idx];
	}
	else
	return localData_[0];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline Numeric& ArrayND<Numeric,Len,Dim>::valueAt(
	const unsigned *index, const unsigned dim)
	{
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"valueAt\"");
	if (dim != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::valueAt: incompatible index length");
	if (dim)
	{
	assert(index);
	unsigned long idx = 0UL;
	for (unsigned i=0; i<dim_; ++i)
	{
	if (index[i] >= shape_[i]) throw std::out_of_range(
	"In npstat::ArrayND::valueAt: index out of range");
	idx += index[i]*strides_[i];
	}
	return data_[idx];
	}
	else
	return localData_[0];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline Numeric& ArrayND<Numeric,Len,Dim>::operator()()
	{
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"operator()\"");
	if (dim_) throw std::invalid_argument(
	"In npstat::ArrayND::operator(): wrong # of args (not rank 0 array)");
	return localData_[0];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline const Numeric& ArrayND<Numeric,Len,Dim>::operator()() const
	{
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"operator()\"");
	if (dim_) throw std::invalid_argument(
	"In npstat::ArrayND::operator(): wrong # of args (not rank 0 array)");
	return localData_[0];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline Numeric& ArrayND<Numeric,Len,Dim>::operator()(
	const unsigned i)
	{
	if (1U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::operator(): wrong # of args (not rank 1 array)");
	return data_[i];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline const Numeric& ArrayND<Numeric,Len,Dim>::operator()(
	const unsigned i) const
	{
	if (1U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::operator(): wrong # of args (not rank 1 array)");
	return data_[i];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	const Numeric& ArrayND<Numeric,Len,Dim>::at() const
	{
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"at\"");
	if (dim_) throw std::invalid_argument(
	"In npstat::ArrayND::at: wrong # of args (not rank 0 array)");
	return localData_[0];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	Numeric& ArrayND<Numeric,Len,Dim>::at()
	{
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"at\"");
	if (dim_) throw std::invalid_argument(
	"In npstat::ArrayND::at: wrong # of args (not rank 0 array)");
	return localData_[0];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	const Numeric& ArrayND<Numeric,Len,Dim>::at(
	const unsigned i0) const
	{
	if (1U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::at: wrong # of args (not rank 1 array)");
	if (i0 >= shape_[0]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 0 out of range (rank 1)");
	return data_[i0];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	Numeric& ArrayND<Numeric,Len,Dim>::at(
	const unsigned i0)
	{
	if (1U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::at: wrong # of args (not rank 1 array)");
	if (i0 >= shape_[0]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 0 out of range (rank 1)");
	return data_[i0];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline Numeric& ArrayND<Numeric,Len,Dim>::operator()(
	const unsigned i0,
	const unsigned i1)
	{
	if (2U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::operator(): wrong # of args (not rank 2 array)");
	return data_[i0*strides_[0] + i1];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline const Numeric& ArrayND<Numeric,Len,Dim>::operator()(
	const unsigned i0,
	const unsigned i1) const
	{
	if (2U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::operator(): wrong # of args (not rank 2 array)");
	return data_[i0*strides_[0] + i1];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	const Numeric& ArrayND<Numeric,Len,Dim>::at(
	const unsigned i0,
	const unsigned i1) const
	{
	if (2U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::at: wrong # of args (not rank 2 array)");
	if (i0 >= shape_[0]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 0 out of range (rank 2)");
	if (i1 >= shape_[1]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 1 out of range (rank 2)");
	return data_[i0*strides_[0] + i1];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	Numeric& ArrayND<Numeric,Len,Dim>::at(
	const unsigned i0,
	const unsigned i1)
	{
	if (2U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::at: wrong # of args (not rank 2 array)");
	if (i0 >= shape_[0]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 0 out of range (rank 2)");
	if (i1 >= shape_[1]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 1 out of range (rank 2)");
	return data_[i0*strides_[0] + i1];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline const Numeric& ArrayND<Numeric,Len,Dim>::operator()(
	const unsigned i0,
	const unsigned i1,
	const unsigned i2) const
	{
	if (3U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::operator(): wrong # of args (not rank 3 array)");
	return data_[i0strides_[0] + i1strides_[1] + i2];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline const Numeric& ArrayND<Numeric,Len,Dim>::operator()(
	const unsigned i0,
	const unsigned i1,
	const unsigned i2,
	const unsigned i3) const
	{
	if (4U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::operator(): wrong # of args (not rank 4 array)");
	return data_[i0strides_[0] + i1strides_[1] + i2*strides_[2] + i3];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline const Numeric& ArrayND<Numeric,Len,Dim>::operator()(
	const unsigned i0,
	const unsigned i1,
	const unsigned i2,
	const unsigned i3,
	const unsigned i4) const
	{
	if (5U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::operator(): wrong # of args (not rank 5 array)");
	return data_[i0strides_[0] + i1strides_[1] + i2*strides_[2] +
	i3*strides_[3] + i4];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline const Numeric& ArrayND<Numeric,Len,Dim>::operator()(
	const unsigned i0,
	const unsigned i1,
	const unsigned i2,
	const unsigned i3,
	const unsigned i4,
	const unsigned i5) const
	{
	if (6U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::operator(): wrong # of args (not rank 6 array)");
	return data_[i0strides_[0] + i1strides_[1] + i2*strides_[2] +
	i3strides_[3] + i4strides_[4] + i5];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline const Numeric& ArrayND<Numeric,Len,Dim>::operator()(
	const unsigned i0,
	const unsigned i1,
	const unsigned i2,
	const unsigned i3,
	const unsigned i4,
	const unsigned i5,
	const unsigned i6) const
	{
	if (7U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::operator(): wrong # of args (not rank 7 array)");
	return data_[i0strides_[0] + i1strides_[1] + i2*strides_[2] +
	i3strides_[3] + i4strides_[4] + i5*strides_[5] + i6];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline const Numeric& ArrayND<Numeric,Len,Dim>::operator()(
	const unsigned i0,
	const unsigned i1,
	const unsigned i2,
	const unsigned i3,
	const unsigned i4,
	const unsigned i5,
	const unsigned i6,
	const unsigned i7) const
	{
	if (8U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::operator(): wrong # of args (not rank 8 array)");
	return data_[i0strides_[0] + i1strides_[1] + i2*strides_[2] +
	i3strides_[3] + i4strides_[4] + i5*strides_[5] +
	i6*strides_[6] + i7];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline const Numeric& ArrayND<Numeric,Len,Dim>::operator()(
	const unsigned i0,
	const unsigned i1,
	const unsigned i2,
	const unsigned i3,
	const unsigned i4,
	const unsigned i5,
	const unsigned i6,
	const unsigned i7,
	const unsigned i8) const
	{
	if (9U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::operator(): wrong # of args (not rank 9 array)");
	return data_[i0strides_[0] + i1strides_[1] + i2*strides_[2] +
	i3strides_[3] + i4strides_[4] + i5*strides_[5] +
	i6strides_[6] + i7strides_[7] + i8];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline const Numeric& ArrayND<Numeric,Len,Dim>::operator()(
	const unsigned i0,
	const unsigned i1,
	const unsigned i2,
	const unsigned i3,
	const unsigned i4,
	const unsigned i5,
	const unsigned i6,
	const unsigned i7,
	const unsigned i8,
	const unsigned i9) const
	{
	if (10U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::operator(): wrong # of args (not rank 10 array)");
	return data_[i0strides_[0] + i1strides_[1] + i2*strides_[2] +
	i3strides_[3] + i4strides_[4] + i5*strides_[5] +
	i6strides_[6] + i7strides_[7] + i8*strides_[8] + i9];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline Numeric& ArrayND<Numeric,Len,Dim>::operator()(
	const unsigned i0,
	const unsigned i1,
	const unsigned i2)
	{
	if (3U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::operator(): wrong # of args (not rank 3 array)");
	return data_[i0strides_[0] + i1strides_[1] + i2];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline Numeric& ArrayND<Numeric,Len,Dim>::operator()(
	const unsigned i0,
	const unsigned i1,
	const unsigned i2,
	const unsigned i3)
	{
	if (4U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::operator(): wrong # of args (not rank 4 array)");
	return data_[i0strides_[0] + i1strides_[1] + i2*strides_[2] + i3];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline Numeric& ArrayND<Numeric,Len,Dim>::operator()(
	const unsigned i0,
	const unsigned i1,
	const unsigned i2,
	const unsigned i3,
	const unsigned i4)
	{
	if (5U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::operator(): wrong # of args (not rank 5 array)");
	return data_[i0strides_[0] + i1strides_[1] + i2*strides_[2] +
	i3*strides_[3] + i4];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline Numeric& ArrayND<Numeric,Len,Dim>::operator()(
	const unsigned i0,
	const unsigned i1,
	const unsigned i2,
	const unsigned i3,
	const unsigned i4,
	const unsigned i5)
	{
	if (6U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::operator(): wrong # of args (not rank 6 array)");
	return data_[i0strides_[0] + i1strides_[1] + i2*strides_[2] +
	i3strides_[3] + i4strides_[4] + i5];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline Numeric& ArrayND<Numeric,Len,Dim>::operator()(
	const unsigned i0,
	const unsigned i1,
	const unsigned i2,
	const unsigned i3,
	const unsigned i4,
	const unsigned i5,
	const unsigned i6)
	{
	if (7U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::operator(): wrong # of args (not rank 7 array)");
	return data_[i0strides_[0] + i1strides_[1] + i2*strides_[2] +
	i3strides_[3] + i4strides_[4] + i5*strides_[5] + i6];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline Numeric& ArrayND<Numeric,Len,Dim>::operator()(
	const unsigned i0,
	const unsigned i1,
	const unsigned i2,
	const unsigned i3,
	const unsigned i4,
	const unsigned i5,
	const unsigned i6,
	const unsigned i7)
	{
	if (8U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::operator(): wrong # of args (not rank 8 array)");
	return data_[i0strides_[0] + i1strides_[1] + i2*strides_[2] +
	i3strides_[3] + i4strides_[4] + i5*strides_[5] +
	i6*strides_[6] + i7];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline Numeric& ArrayND<Numeric,Len,Dim>::operator()(
	const unsigned i0,
	const unsigned i1,
	const unsigned i2,
	const unsigned i3,
	const unsigned i4,
	const unsigned i5,
	const unsigned i6,
	const unsigned i7,
	const unsigned i8)
	{
	if (9U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::operator(): wrong # of args (not rank 9 array)");
	return data_[i0strides_[0] + i1strides_[1] + i2*strides_[2] +
	i3strides_[3] + i4strides_[4] + i5*strides_[5] +
	i6strides_[6] + i7strides_[7] + i8];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline Numeric& ArrayND<Numeric,Len,Dim>::operator()(
	const unsigned i0,
	const unsigned i1,
	const unsigned i2,
	const unsigned i3,
	const unsigned i4,
	const unsigned i5,
	const unsigned i6,
	const unsigned i7,
	const unsigned i8,
	const unsigned i9)
	{
	if (10U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::operator(): wrong # of args (not rank 10 array)");
	return data_[i0strides_[0] + i1strides_[1] + i2*strides_[2] +
	i3strides_[3] + i4strides_[4] + i5*strides_[5] +
	i6strides_[6] + i7strides_[7] + i8*strides_[8] + i9];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	const Numeric& ArrayND<Numeric,Len,Dim>::at(
	const unsigned i0,
	const unsigned i1,
	const unsigned i2) const
	{
	if (3U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::at: wrong # of args (not rank 3 array)");
	if (i0 >= shape_[0]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 0 out of range (rank 3)");
	if (i1 >= shape_[1]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 1 out of range (rank 3)");
	if (i2 >= shape_[2]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 2 out of range (rank 3)");
	return data_[i0strides_[0] + i1strides_[1] + i2];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	const Numeric& ArrayND<Numeric,Len,Dim>::at(
	const unsigned i0,
	const unsigned i1,
	const unsigned i2,
	const unsigned i3) const
	{
	if (4U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::at: wrong # of args (not rank 4 array)");
	if (i0 >= shape_[0]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 0 out of range (rank 4)");
	if (i1 >= shape_[1]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 1 out of range (rank 4)");
	if (i2 >= shape_[2]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 2 out of range (rank 4)");
	if (i3 >= shape_[3]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 3 out of range (rank 4)");
	return data_[i0strides_[0] + i1strides_[1] + i2*strides_[2] + i3];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	const Numeric& ArrayND<Numeric,Len,Dim>::at(
	const unsigned i0,
	const unsigned i1,
	const unsigned i2,
	const unsigned i3,
	const unsigned i4) const
	{
	if (5U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::at: wrong # of args (not rank 5 array)");
	if (i0 >= shape_[0]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 0 out of range (rank 5)");
	if (i1 >= shape_[1]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 1 out of range (rank 5)");
	if (i2 >= shape_[2]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 2 out of range (rank 5)");
	if (i3 >= shape_[3]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 3 out of range (rank 5)");
	if (i4 >= shape_[4]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 4 out of range (rank 5)");
	return data_[i0strides_[0] + i1strides_[1] + i2*strides_[2] +
	i3*strides_[3] + i4];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	const Numeric& ArrayND<Numeric,Len,Dim>::at(
	const unsigned i0,
	const unsigned i1,
	const unsigned i2,
	const unsigned i3,
	const unsigned i4,
	const unsigned i5) const
	{
	if (6U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::at: wrong # of args (not rank 6 array)");
	if (i0 >= shape_[0]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 0 out of range (rank 6)");
	if (i1 >= shape_[1]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 1 out of range (rank 6)");
	if (i2 >= shape_[2]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 2 out of range (rank 6)");
	if (i3 >= shape_[3]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 3 out of range (rank 6)");
	if (i4 >= shape_[4]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 4 out of range (rank 6)");
	if (i5 >= shape_[5]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 5 out of range (rank 6)");
	return data_[i0strides_[0] + i1strides_[1] + i2*strides_[2] +
	i3strides_[3] + i4strides_[4] + i5];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	const Numeric& ArrayND<Numeric,Len,Dim>::at(
	const unsigned i0,
	const unsigned i1,
	const unsigned i2,
	const unsigned i3,
	const unsigned i4,
	const unsigned i5,
	const unsigned i6) const
	{
	if (7U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::at: wrong # of args (not rank 7 array)");
	if (i0 >= shape_[0]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 0 out of range (rank 7)");
	if (i1 >= shape_[1]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 1 out of range (rank 7)");
	if (i2 >= shape_[2]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 2 out of range (rank 7)");
	if (i3 >= shape_[3]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 3 out of range (rank 7)");
	if (i4 >= shape_[4]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 4 out of range (rank 7)");
	if (i5 >= shape_[5]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 5 out of range (rank 7)");
	if (i6 >= shape_[6]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 6 out of range (rank 7)");
	return data_[i0strides_[0] + i1strides_[1] + i2*strides_[2] +
	i3strides_[3] + i4strides_[4] + i5*strides_[5] + i6];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	const Numeric& ArrayND<Numeric,Len,Dim>::at(
	const unsigned i0,
	const unsigned i1,
	const unsigned i2,
	const unsigned i3,
	const unsigned i4,
	const unsigned i5,
	const unsigned i6,
	const unsigned i7) const
	{
	if (8U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::at: wrong # of args (not rank 8 array)");
	if (i0 >= shape_[0]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 0 out of range (rank 8)");
	if (i1 >= shape_[1]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 1 out of range (rank 8)");
	if (i2 >= shape_[2]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 2 out of range (rank 8)");
	if (i3 >= shape_[3]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 3 out of range (rank 8)");
	if (i4 >= shape_[4]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 4 out of range (rank 8)");
	if (i5 >= shape_[5]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 5 out of range (rank 8)");
	if (i6 >= shape_[6]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 6 out of range (rank 8)");
	if (i7 >= shape_[7]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 7 out of range (rank 8)");
	return data_[i0strides_[0] + i1strides_[1] + i2*strides_[2] +
	i3strides_[3] + i4strides_[4] + i5*strides_[5] +
	i6*strides_[6] + i7];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	const Numeric& ArrayND<Numeric,Len,Dim>::at(
	const unsigned i0,
	const unsigned i1,
	const unsigned i2,
	const unsigned i3,
	const unsigned i4,
	const unsigned i5,
	const unsigned i6,
	const unsigned i7,
	const unsigned i8) const
	{
	if (9U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::at: wrong # of args (not rank 9 array)");
	if (i0 >= shape_[0]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 0 out of range (rank 9)");
	if (i1 >= shape_[1]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 1 out of range (rank 9)");
	if (i2 >= shape_[2]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 2 out of range (rank 9)");
	if (i3 >= shape_[3]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 3 out of range (rank 9)");
	if (i4 >= shape_[4]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 4 out of range (rank 9)");
	if (i5 >= shape_[5]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 5 out of range (rank 9)");
	if (i6 >= shape_[6]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 6 out of range (rank 9)");
	if (i7 >= shape_[7]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 7 out of range (rank 9)");
	if (i8 >= shape_[8]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 8 out of range (rank 9)");
	return data_[i0strides_[0] + i1strides_[1] + i2*strides_[2] +
	i3strides_[3] + i4strides_[4] + i5*strides_[5] +
	i6strides_[6] + i7strides_[7] + i8];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	const Numeric& ArrayND<Numeric,Len,Dim>::at(
	const unsigned i0,
	const unsigned i1,
	const unsigned i2,
	const unsigned i3,
	const unsigned i4,
	const unsigned i5,
	const unsigned i6,
	const unsigned i7,
	const unsigned i8,
	const unsigned i9) const
	{
	if (10U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::at: wrong # of args (not rank 10 array)");
	if (i0 >= shape_[0]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 0 out of range (rank 10)");
	if (i1 >= shape_[1]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 1 out of range (rank 10)");
	if (i2 >= shape_[2]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 2 out of range (rank 10)");
	if (i3 >= shape_[3]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 3 out of range (rank 10)");
	if (i4 >= shape_[4]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 4 out of range (rank 10)");
	if (i5 >= shape_[5]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 5 out of range (rank 10)");
	if (i6 >= shape_[6]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 6 out of range (rank 10)");
	if (i7 >= shape_[7]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 7 out of range (rank 10)");
	if (i8 >= shape_[8]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 8 out of range (rank 10)");
	if (i9 >= shape_[9]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 9 out of range (rank 10)");
	return data_[i0strides_[0] + i1strides_[1] + i2*strides_[2] +
	i3strides_[3] + i4strides_[4] + i5*strides_[5] +
	i6strides_[6] + i7strides_[7] + i8*strides_[8] + i9];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	Numeric& ArrayND<Numeric,Len,Dim>::at(
	const unsigned i0,
	const unsigned i1,
	const unsigned i2)
	{
	if (3U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::at: wrong # of args (not rank 3 array)");
	if (i0 >= shape_[0]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 0 out of range (rank 3)");
	if (i1 >= shape_[1]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 1 out of range (rank 3)");
	if (i2 >= shape_[2]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 2 out of range (rank 3)");
	return data_[i0strides_[0] + i1strides_[1] + i2];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	Numeric& ArrayND<Numeric,Len,Dim>::at(
	const unsigned i0,
	const unsigned i1,
	const unsigned i2,
	const unsigned i3)
	{
	if (4U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::at: wrong # of args (not rank 4 array)");
	if (i0 >= shape_[0]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 0 out of range (rank 4)");
	if (i1 >= shape_[1]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 1 out of range (rank 4)");
	if (i2 >= shape_[2]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 2 out of range (rank 4)");
	if (i3 >= shape_[3]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 3 out of range (rank 4)");
	return data_[i0strides_[0] + i1strides_[1] + i2*strides_[2] + i3];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	Numeric& ArrayND<Numeric,Len,Dim>::at(
	const unsigned i0,
	const unsigned i1,
	const unsigned i2,
	const unsigned i3,
	const unsigned i4)
	{
	if (5U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::at: wrong # of args (not rank 5 array)");
	if (i0 >= shape_[0]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 0 out of range (rank 5)");
	if (i1 >= shape_[1]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 1 out of range (rank 5)");
	if (i2 >= shape_[2]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 2 out of range (rank 5)");
	if (i3 >= shape_[3]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 3 out of range (rank 5)");
	if (i4 >= shape_[4]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 4 out of range (rank 5)");
	return data_[i0strides_[0] + i1strides_[1] + i2*strides_[2] +
	i3*strides_[3] + i4];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	Numeric& ArrayND<Numeric,Len,Dim>::at(
	const unsigned i0,
	const unsigned i1,
	const unsigned i2,
	const unsigned i3,
	const unsigned i4,
	const unsigned i5)
	{
	if (6U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::at: wrong # of args (not rank 6 array)");
	if (i0 >= shape_[0]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 0 out of range (rank 6)");
	if (i1 >= shape_[1]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 1 out of range (rank 6)");
	if (i2 >= shape_[2]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 2 out of range (rank 6)");
	if (i3 >= shape_[3]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 3 out of range (rank 6)");
	if (i4 >= shape_[4]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 4 out of range (rank 6)");
	if (i5 >= shape_[5]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 5 out of range (rank 6)");
	return data_[i0strides_[0] + i1strides_[1] + i2*strides_[2] +
	i3strides_[3] + i4strides_[4] + i5];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	Numeric& ArrayND<Numeric,Len,Dim>::at(
	const unsigned i0,
	const unsigned i1,
	const unsigned i2,
	const unsigned i3,
	const unsigned i4,
	const unsigned i5,
	const unsigned i6)
	{
	if (7U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::at: wrong # of args (not rank 7 array)");
	if (i0 >= shape_[0]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 0 out of range (rank 7)");
	if (i1 >= shape_[1]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 1 out of range (rank 7)");
	if (i2 >= shape_[2]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 2 out of range (rank 7)");
	if (i3 >= shape_[3]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 3 out of range (rank 7)");
	if (i4 >= shape_[4]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 4 out of range (rank 7)");
	if (i5 >= shape_[5]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 5 out of range (rank 7)");
	if (i6 >= shape_[6]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 6 out of range (rank 7)");
	return data_[i0strides_[0] + i1strides_[1] + i2*strides_[2] +
	i3strides_[3] + i4strides_[4] + i5*strides_[5] + i6];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	Numeric& ArrayND<Numeric,Len,Dim>::at(
	const unsigned i0,
	const unsigned i1,
	const unsigned i2,
	const unsigned i3,
	const unsigned i4,
	const unsigned i5,
	const unsigned i6,
	const unsigned i7)
	{
	if (8U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::at: wrong # of args (not rank 8 array)");
	if (i0 >= shape_[0]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 0 out of range (rank 8)");
	if (i1 >= shape_[1]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 1 out of range (rank 8)");
	if (i2 >= shape_[2]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 2 out of range (rank 8)");
	if (i3 >= shape_[3]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 3 out of range (rank 8)");
	if (i4 >= shape_[4]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 4 out of range (rank 8)");
	if (i5 >= shape_[5]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 5 out of range (rank 8)");
	if (i6 >= shape_[6]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 6 out of range (rank 8)");
	if (i7 >= shape_[7]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 7 out of range (rank 8)");
	return data_[i0strides_[0] + i1strides_[1] + i2*strides_[2] +
	i3strides_[3] + i4strides_[4] + i5*strides_[5] +
	i6*strides_[6] + i7];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	Numeric& ArrayND<Numeric,Len,Dim>::at(
	const unsigned i0,
	const unsigned i1,
	const unsigned i2,
	const unsigned i3,
	const unsigned i4,
	const unsigned i5,
	const unsigned i6,
	const unsigned i7,
	const unsigned i8)
	{
	if (9U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::at: wrong # of args (not rank 9 array)");
	if (i0 >= shape_[0]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 0 out of range (rank 9)");
	if (i1 >= shape_[1]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 1 out of range (rank 9)");
	if (i2 >= shape_[2]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 2 out of range (rank 9)");
	if (i3 >= shape_[3]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 3 out of range (rank 9)");
	if (i4 >= shape_[4]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 4 out of range (rank 9)");
	if (i5 >= shape_[5]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 5 out of range (rank 9)");
	if (i6 >= shape_[6]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 6 out of range (rank 9)");
	if (i7 >= shape_[7]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 7 out of range (rank 9)");
	if (i8 >= shape_[8]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 8 out of range (rank 9)");
	return data_[i0strides_[0] + i1strides_[1] + i2*strides_[2] +
	i3strides_[3] + i4strides_[4] + i5*strides_[5] +
	i6strides_[6] + i7strides_[7] + i8];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	Numeric& ArrayND<Numeric,Len,Dim>::at(
	const unsigned i0,
	const unsigned i1,
	const unsigned i2,
	const unsigned i3,
	const unsigned i4,
	const unsigned i5,
	const unsigned i6,
	const unsigned i7,
	const unsigned i8,
	const unsigned i9)
	{
	if (10U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::at: wrong # of args (not rank 10 array)");
	if (i0 >= shape_[0]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 0 out of range (rank 10)");
	if (i1 >= shape_[1]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 1 out of range (rank 10)");
	if (i2 >= shape_[2]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 2 out of range (rank 10)");
	if (i3 >= shape_[3]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 3 out of range (rank 10)");
	if (i4 >= shape_[4]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 4 out of range (rank 10)");
	if (i5 >= shape_[5]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 5 out of range (rank 10)");
	if (i6 >= shape_[6]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 6 out of range (rank 10)");
	if (i7 >= shape_[7]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 7 out of range (rank 10)");
	if (i8 >= shape_[8]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 8 out of range (rank 10)");
	if (i9 >= shape_[9]) throw std::out_of_range(
	"In npstat::ArrayND::at: index 9 out of range (rank 10)");
	return data_[i0strides_[0] + i1strides_[1] + i2*strides_[2] +
	i3strides_[3] + i4strides_[4] + i5*strides_[5] +
	i6strides_[6] + i7strides_[7] + i8*strides_[8] + i9];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template<unsigned Len2, unsigned Dim2>
	double ArrayND<Numeric,Len,Dim>::maxAbsDifference(
	const ArrayND<Numeric,Len2,Dim2>& r) const
	{
	if (!isShapeCompatible(r)) throw std::invalid_argument(
	"In npstat::ArrayND::maxAbsDifference: "
	"incompatible argument array shape");
	if (dim_)
	{
	double maxd = 0.0;
	for (unsigned long i=0; i<len_; ++i)
	{
	const Numeric rval = r.data_[i];
	const double d = absDifference(data_[i], rval);
	if (d > maxd)
	maxd = d;
	}
	return maxd;
	}
	else
	{
	const Numeric rval = r.localData_[0];
	return absDifference(localData_[0], rval);
	}
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template<unsigned Len2, unsigned Dim2>
	bool ArrayND<Numeric,Len,Dim>::operator==(
	const ArrayND<Numeric,Len2,Dim2>& r) const
	{
	if (shapeIsKnown_ != r.shapeIsKnown_)
	return false;
	if (r.dim_ != dim_)
	return false;
	if (r.len_ != len_)
	return false;
	for (unsigned i=0; i<dim_; ++i)
	if (shape_[i] != r.shape_[i])
	return false;
	for (unsigned i=0; i<dim_; ++i)
	assert(strides_[i] == r.strides_[i]);
	for (unsigned long j=0; j<len_; ++j)
	if (data_[j] != r.data_[j])
	return false;
	return true;
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template<unsigned Len2, unsigned Dim2>
	inline bool ArrayND<Numeric,Len,Dim>::operator!=(
	const ArrayND<Numeric,Len2,Dim2>& r) const
	{
	return !(*this == r);
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template<typename Num2>
	ArrayND<Numeric,Len,Dim>
	ArrayND<Numeric,Len,Dim>::operator*(const Num2& r) const
	{
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"operator*\"");
	ArrayND<Numeric,Len,Dim> result(shape_, dim_);
	for (unsigned long i=0; i<len_; ++i)
	result.data_[i] = data_[i]*r;
	return result;
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template<typename Num2>
	ArrayND<Numeric,Len,Dim>
	ArrayND<Numeric,Len,Dim>::operator/(const Num2& r) const
	{
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"operator/\"");
	if (r == Num2()) throw std::runtime_error(
	"In npstat::ArrayND::operator/: division by zero");
	ArrayND<Numeric,Len,Dim> result(shape_, dim_);
	for (unsigned long i=0; i<len_; ++i)
	result.data_[i] = data_[i]/r;
	return result;
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline ArrayND<Numeric,Len,Dim> ArrayND<Numeric,Len,Dim>::operator+() const
	{
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"operator+\"");
	return *this;
	}

	// template<typename Numeric, unsigned Len, unsigned Dim>
	// ArrayND<Numeric,Len,Dim> ArrayND<Numeric,Len,Dim>::operator-() const
	// {
	// if (!shapeIsKnown_) throw std::invalid_argument(
	// "Initialize npstat::ArrayND before calling method \"operator-\"");
	// ArrayND<Numeric,Len,Dim> result(shape_, dim_);
	// for (unsigned long i=0; i<len_; ++i)
	// result.data_[i] = -data_[i];
	// return result;
	// }

	template<typename Numeric, unsigned Len, unsigned Dim>
	template <typename Num2>
	ArrayND<Numeric,Len,Dim>&
	ArrayND<Numeric,Len,Dim>::operator*=(const Num2& r)
	{
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"operator*=\"");
	for (unsigned long i=0; i<len_; ++i)
	data_[i] *= r;
	return *this;
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	ArrayND<Numeric,Len,Dim>&
	ArrayND<Numeric,Len,Dim>::makeNonNegative()
	{
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"makeNonNegative\"");
	const Numeric zero = Numeric();
	if (dim_)
	{
	for (unsigned long i=0; i<len_; ++i)
	if (!(ComplexComparesAbs<Numeric>::more(data_[i], zero)))
	data_[i] = zero;
	}
	else
	if (!(ComplexComparesAbs<Numeric>::more(localData_[0], zero)))
	localData_[0] = zero;
	return *this;
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	unsigned ArrayND<Numeric,Len,Dim>::makeCopulaSteps(
	const double tolerance, const unsigned nCycles)
	{
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"makeCopulaSteps\"");
	if (nCycles == 0U)
	return 0U;
	if (!dim_) throw std::invalid_argument(
	"npstat::ArrayND::makeCopulaSteps method "
	"can not be used with array of 0 rank");

	const Numeric zero = Numeric();
	for (unsigned long i=0; i<len_; ++i)
	if (!(ComplexComparesAbs<Numeric>::more(data_[i], zero)))
	data_[i] = zero;

	std::vector<Numeric*> axesPtrBuf(dim_);
	Numeric** axes = &axesPtrBuf[0];
	const Numeric one = static_cast<Numeric>(1);

	// Memory for the axis accumulators
	unsigned idxSum = 0;
	for (unsigned i=0; i<dim_; ++i)
	idxSum += shape_[i];
	std::vector<Numeric> axesBuf(idxSum);
	axes[0] = &axesBuf[0];
	for (unsigned i=1; i<dim_; ++i)
	axes[i] = axes[i-1] + shape_[i-1];

	// Accumulate axis projections
	unsigned icycle = 0;
	for (; icycle<nCycles; ++icycle)
	{
	for (unsigned i=0; i<idxSum; ++i)
	axesBuf[i] = zero;

	// Accumulate sums for each axis
	for (unsigned long idat=0; idat<len_; ++idat)
	{
	unsigned long l = idat;
	for (unsigned i=0; i<dim_; ++i)
	{
	const unsigned idx = l / strides_[i];
	l -= (idx * strides_[i]);
	axes[i][idx] += data_[idat];
	}
	}

	// Make averages out of sums
	bool withinTolerance = true;
	Numeric totalSum = zero;
	for (unsigned i=0; i<dim_; ++i)
	{
	Numeric axisSum = zero;
	const unsigned amax = shape_[i];
	for (unsigned a=0; a<amax; ++a)
	{
	if (axes[i][a] == zero)
	throw std::runtime_error(
	"In npstat::ArrayND::makeCopulaSteps: "
	"marginal density is zero");
	axisSum += axes[i][a];
	}
	totalSum += axisSum;
	const Numeric axisAverage = axisSum/static_cast<Numeric>(amax);
	for (unsigned a=0; a<amax; ++a)
	axes[i][a] /= axisAverage;
	for (unsigned a=0; a<amax && withinTolerance; ++a)
	{
	const double adelta = absDifference(axes[i][a], one);
	if (adelta > tolerance)
	withinTolerance = false;
	}
	}

	if (withinTolerance)
	break;

	const Numeric totalAverage = totalSum/
	static_cast<Numeric>(len_)/static_cast<Numeric>(dim_);

	// Run over all points again and divide by
	// the product of marginals
	for (unsigned long idat=0; idat<len_; ++idat)
	{
	unsigned long l = idat;
	for (unsigned i=0; i<dim_; ++i)
	{
	const unsigned idx = l / strides_[i];
	l -= (idx * strides_[i]);
	data_[idat] /= axes[i][idx];
	}
	data_[idat] /= totalAverage;
	}
	}

	return icycle;
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template <typename Num2>
	ArrayND<Numeric,Len,Dim>&
	ArrayND<Numeric,Len,Dim>::operator/=(const Num2& r)
	{
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"operator/=\"");
	if (r == Num2()) throw std::runtime_error(
	"In npstat::ArrayND::operator/=: division by zero");
	for (unsigned long i=0; i<len_; ++i)
	data_[i] /= r;
	return *this;
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template <typename Num2, unsigned Len2, unsigned Dim2>
	ArrayND<Numeric,Len,Dim>&
	ArrayND<Numeric,Len,Dim>::operator+=(const ArrayND<Num2,Len2,Dim2>& r)
	{
	if (!isShapeCompatible(r)) throw std::invalid_argument(
	"In npstat::ArrayND::operator+=: "
	"incompatible argument array shape");
	for (unsigned long i=0; i<len_; ++i)
	data_[i] += r.data_[i];
	return *this;
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template <typename Num2, unsigned Len2, unsigned Dim2>
	ArrayND<Numeric,Len,Dim>&
	ArrayND<Numeric,Len,Dim>::operator*=(const ArrayND<Num2,Len2,Dim2>& r)
	{
	if (!isShapeCompatible(r)) throw std::invalid_argument(
	"In npstat::ArrayND::operator*=: "
	"incompatible argument array shape");
	for (unsigned long i=0; i<len_; ++i)
	data_[i] *= r.data_[i];
	return *this;
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template <typename Num2, unsigned Len2, unsigned Dim2>
	ArrayND<Numeric,Len,Dim>&
	ArrayND<Numeric,Len,Dim>::operator/=(const ArrayND<Num2,Len2,Dim2>& r)
	{
	if (!isShapeCompatible(r)) throw std::invalid_argument(
	"In npstat::ArrayND::operator/=: "
	"incompatible argument array shape");
	if (!r.isNonZero()) throw std::invalid_argument(
	"In npstat::ArrayND::operator/=: "
	"division by zero encountered");
	for (unsigned long i=0; i<len_; ++i)
	data_[i] /= r.data_[i];
	return *this;
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template <typename Num3, typename Num2, unsigned Len2, unsigned Dim2>
	ArrayND<Numeric,Len,Dim>&
	ArrayND<Numeric,Len,Dim>::addmul(const ArrayND<Num2,Len2,Dim2>& r,
	const Num3& c)
	{
	if (!isShapeCompatible(r)) throw std::invalid_argument(
	"In npstat::ArrayND::addmul: "
	"incompatible argument array shape");
	for (unsigned long i=0; i<len_; ++i)
	data_[i] += r.data_[i]*c;
	return *this;
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template <typename Num2, unsigned Len2, unsigned Dim2>
	ArrayND<Numeric,Len,Dim>&
	ArrayND<Numeric,Len,Dim>::operator-=(const ArrayND<Num2,Len2,Dim2>& r)
	{
	if (!isShapeCompatible(r)) throw std::invalid_argument(
	"In npstat::ArrayND::operator-=: "
	"incompatible argument array shape");
	for (unsigned long i=0; i<len_; ++i)
	data_[i] -= r.data_[i];
	return *this;
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	Numeric ArrayND<Numeric,Len,Dim>::interpolate1(
	const double *coords, const unsigned dim) const
	{
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"interpolate1\"");
	if (dim != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::interpolate1: incompatible coordinate length");
	if (dim)
	{
	const unsigned maxdim = CHAR_BIT*sizeof(unsigned long);
	if (dim_ >= maxdim)
	throw std::invalid_argument(
	"In npstat::ArrayND::interpolate1: array rank is too large");

	double dx[maxdim];
	unsigned ix[maxdim];
	for (unsigned i=0; i<dim; ++i)
	{
	const double x = coords[i];
	if (x <= 0.0)
	{
	ix[i] = 0;
	dx[i] = 0.0;
	}
	else if (x >= static_cast<double>(shape_[i] - 1))
	{
	ix[i] = shape_[i] - 1;
	dx[i] = 0.0;
	}
	else
	{
	ix[i] = static_cast<unsigned>(std::floor(x));
	dx[i] = x - ix[i];
	}
	}

	Numeric sum = Numeric();
	const unsigned long maxcycle = 1UL << dim;
	for (unsigned long icycle=0UL; icycle<maxcycle; ++icycle)
	{
	double w = 1.0;
	unsigned long icell = 0UL;
	for (unsigned i=0; i<dim; ++i)
	{
	if (icycle & (1UL << i))
	{
	w *= dx[i];
	icell += strides_[i]*(ix[i] + 1U);
	}
	else
	{
	w *= (1.0 - dx[i]);
	icell += strides_[i]*ix[i];
	}
	}
	if (w > 0.0)
	sum += data_[icell]*static_cast<proper_double>(w);
	}
	return sum;
	}
	else
	return localData_[0];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	Numeric ArrayND<Numeric,Len,Dim>::interpolateLoop(
	const unsigned level, const double* coords, const Numeric* base) const
	{
	const unsigned npoints = shape_[level];
	const double x = coords[level];

	unsigned ix, npt = 1;
	double dx = 0.0;
	if (x < 0.0)
	ix = 0;
	else if (x > static_cast<double>(npoints - 1))
	ix = npoints - 1;
	else
	{
	ix = static_cast<unsigned>(std::floor(x));
	if (ix) --ix;
	unsigned imax = ix + 3;
	while (imax >= npoints)
	{
	if (ix) --ix;
	--imax;
	}
	dx = x - ix;
	npt = imax + 1 - ix;
	}
	assert(npt >= 1 && npt <= 4);

	Numeric fit[4];
	if (level < dim_ - 1)
	for (unsigned ipt=0; ipt<npt; ++ipt)
	fit[ipt] = interpolateLoop(level + 1, coords,
	base + (ix + ipt)*strides_[level]);

	const Numeric* const v = (level == dim_ - 1 ? base + ix : fit);
	switch (npt)
	{
	case 1:
	return v[0];
	case 2:
	return interpolate_linear(dx, v[0], v[1]);
	case 3:
	return interpolate_quadratic(dx, v[0], v[1], v[2]);
	case 4:
	return interpolate_cubic(dx, v[0], v[1], v[2], v[3]);
	default:
	assert(0);
	return Numeric();
	}
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline Numeric ArrayND<Numeric,Len,Dim>::interpolate3(
	const double* coords, const unsigned dim) const
	{
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"interpolate3\"");
	if (dim != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::interpolate3: incompatible coordinate length");
	if (dim)
	{
	assert(coords);
	return interpolateLoop(0, coords, data_);
	}
	else
	return localData_[0];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template<class Functor>
	ArrayND<Numeric,Len,Dim>& ArrayND<Numeric,Len,Dim>::apply(Functor f)
	{
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"apply\"");
	for (unsigned long i=0; i<len_; ++i)
	data_[i] = static_cast<Numeric>(f(data_[i]));
	return *this;
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template<class Functor>
	ArrayND<Numeric,Len,Dim>& ArrayND<Numeric,Len,Dim>::scanInPlace(
	Functor f)
	{
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"scanInPlace\"");
	for (unsigned long i=0; i<len_; ++i)
	f(data_[i]);
	return *this;
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	ArrayND<Numeric,Len,Dim>& ArrayND<Numeric,Len,Dim>::constFill(
	const Numeric c)
	{
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"constFill\"");
	for (unsigned long i=0; i<len_; ++i)
	data_[i] = c;
	return *this;
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline ArrayND<Numeric,Len,Dim>& ArrayND<Numeric,Len,Dim>::clear()
	{
	return constFill(Numeric());
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	ArrayND<Numeric,Len,Dim>& ArrayND<Numeric,Len,Dim>::uninitialize()
	{
	if (!dataIsExternal_)
	destroyBuffer(data_, localData_);
	destroyBuffer(strides_, localStrides_);
	destroyBuffer(shape_, localShape_);
	localData_[0] = Numeric();
	data_ = localData_;
	strides_ = 0;
	shape_ = 0;
	len_ = 0;
	dim_ = 0;
	shapeIsKnown_ = false;
	dataIsExternal_ = false;
	return *this;
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	ArrayND<Numeric,Len,Dim>& ArrayND<Numeric,Len,Dim>::makeUnit()
	{
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"makeUnit\"");
	if (dim_ < 2) throw std::invalid_argument(
	"npstat::ArrayND::makeUnit method "
	"can not be used with arrays of rank less than 2");
	constFill(Numeric());
	unsigned long stride = 0UL;
	const unsigned dimlen = shape_[0];
	for (unsigned i=0; i<dim_; ++i)
	{
	if (shape_[i] != dimlen) throw std::invalid_argument(
	"npstat::ArrayND::makeUnit method needs "
	"the array span to be the same in ech dimension");
	stride += strides_[i];
	}
	const Numeric one(static_cast<Numeric>(1));
	for (unsigned i=0; i<dimlen; ++i)
	data_[i*stride] = one;
	return *this;
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	void ArrayND<Numeric,Len,Dim>::linearFillLoop(
	const unsigned level, const double s0, const unsigned long idx,
	const double shift, const double* coeffs)
	{
	const unsigned imax = shape_[level];
	const double c = coeffs[level];
	if (level == dim_ - 1)
	{
	Numeric* d = &data_[idx];
	for (unsigned i=0; i<imax; ++i)
	{
	// Note that we want to add "shift" only at the
	// very end. This might improve the numerical
	// precision of the result.
	const double sum = s0 + c*i + shift;
	d[i] = static_cast<Numeric>(sum);
	}
	}
	else
	{
	const unsigned long stride = strides_[level];
	for (unsigned i=0; i<imax; ++i)
	linearFillLoop(level+1, s0 + ci, idx + istride,
	shift, coeffs);
	}
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	ArrayND<Numeric,Len,Dim>& ArrayND<Numeric,Len,Dim>::linearFill(
	const double* coeffs, const unsigned dimCoeffs, const double shift)
	{
	// Make sure the object has been initialized
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"linearFill\"");
	if (dim_ != dimCoeffs) throw std::invalid_argument(
	"In npstat::ArrayND::linearFill: incompatible number of coefficients");
	if (dim_)
	{
	assert(coeffs);
	linearFillLoop(0U, 0.0, 0UL, shift, coeffs);
	}
	else
	localData_[0] = static_cast<Numeric>(shift);
	return *this;
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template<class Functor>
	void ArrayND<Numeric,Len,Dim>::functorFillLoop(
	const unsigned level, const unsigned long idx,
	Functor f, unsigned* farg)
	{
	const unsigned imax = shape_[level];
	if (level == dim_ - 1)
	{
	Numeric* d = &data_[idx];
	const unsigned* myarg = farg;
	for (unsigned i = 0; i<imax; ++i)
	{
	farg[level] = i;
	d[i] = static_cast<Numeric>(f(myarg, dim_));
	}
	}
	else
	{
	const unsigned long stride = strides_[level];
	for (unsigned i = 0; i<imax; ++i)
	{
	farg[level] = i;
	functorFillLoop(level+1, idx + i*stride, f, farg);
	}
	}
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template <typename Accumulator>
	void ArrayND<Numeric,Len,Dim>::convertToLastDimCdfLoop(
	ArrayND* sumSlice, const unsigned level, unsigned long idx0,
	unsigned long idxSlice, const bool useTrapezoids)
	{
	static const proper_double half = 0.5;
	const unsigned imax = shape_[level];
	if (level == dim_ - 1)
	{
	Accumulator acc = Accumulator();
	Numeric* data = data_ + idx0;
	if (useTrapezoids)
	{
	Numeric oldval = Numeric();
	for (unsigned i = 0; i<imax; ++i)
	{
	acc += (data[i] + oldval)*half;
	oldval = data[i];
	data[i] = static_cast<Numeric>(acc);
	}
	acc += oldval*half;
	}
	else
	for (unsigned i = 0; i<imax; ++i)
	{
	acc += data[i];
	data[i] = static_cast<Numeric>(acc);
	}
	if (sumSlice->dim_)
	sumSlice->data_[idxSlice] = static_cast<Numeric>(acc);
	else
	sumSlice->localData_[0] = static_cast<Numeric>(acc);
	}
	else
	{
	const unsigned long stride = strides_[level];
	unsigned long sumStride = 0UL;
	if (sumSlice->dim_)
	sumStride = sumSlice->strides_[level];
	for (unsigned i = 0; i<imax; ++i)
	{
	convertToLastDimCdfLoop<Accumulator>(
	sumSlice, level+1, idx0, idxSlice, useTrapezoids);
	idx0 += stride;
	idxSlice += sumStride;
	}
	}
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template <typename Accumulator>
	inline void ArrayND<Numeric,Len,Dim>::convertToLastDimCdf(
	ArrayND* sumSlice, const bool useTrapezoids)
	{
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling "
	"method \"convertToLastDimCdf\"");
	if (!dim_) throw std::invalid_argument(
	"npstat::ArrayND::convertToLastDimCdf method "
	"can not be used with array of 0 rank");
	assert(sumSlice);
	if (!sumSlice->shapeIsKnown_) throw std::invalid_argument(
	"In npstat::ArrayND::convertToLastDimCdf: "
	"uninitialized argument array");
	convertToLastDimCdfLoop<Accumulator>(sumSlice, 0U, 0UL, 0UL,
	useTrapezoids);
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template <typename Num2, unsigned Len2, unsigned Dim2>
	void ArrayND<Numeric,Len,Dim>::coarseSum(
	const unsigned idim, const unsigned n, ArrayND<Num2,Len2,Dim2>* r) const
	{
	// Check that the input arguments make sense
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"coarseSum\"");
	if (idim >= dim_) throw std::invalid_argument(
	"In npstat::ArrayND::coarseSum: dimension number is out of range");
	if (!n) throw std::invalid_argument(
	"In npstat::ArrayND::coarseSum: number of divisions must be positive");

	assert(shape_);
	if (shape_[idim] % n) throw std::invalid_argument(
	"In npstat::ArrayND::coarseSum: dimension is not exactly divisible");
	assert(r);

	bool reinitResult = false;
	if (!r->shapeIsKnown_)
	reinitResult = true;
	else if (r->dim_ != dim_)
	reinitResult = true;
	else
	{
	assert(r->shape_);
	for (unsigned i=0; i<dim_ && !reinitResult; ++i)
	if (i == idim)
	reinitResult = shape_[i]/n != r->shape_[i];
	else
	reinitResult = shape_[i] != r->shape_[i];
	}
	if (reinitResult)
	{
	r->uninitialize();
	r->len_ = len_ / n;
	r->dim_ = dim_;
	r->shapeIsKnown_ = true;
	r->shape_ = makeBuffer(dim_, r->localShape_, Dim2);
	copyBuffer(r->shape_, shape_, dim_);
	r->shape_[idim] /= n;
	r->buildStrides();
	r->data_ = makeBuffer(len_, r->localData_, Len2);
	}

	switch (n)
	{
	case 1U:
	copyBuffer(r->data_, data_, len_);
	break;
	case 2U:
	coarseSum2(idim, r);
	break;
	case 3U:
	coarseSum3(idim, r);
	break;
	default:
	coarseSumN(idim, n, r);
	break;
	}
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template <typename Num2, unsigned Len2, unsigned Dim2>
	void ArrayND<Numeric,Len,Dim>::coarseAverage(
	const unsigned idim, const unsigned n, ArrayND<Num2,Len2,Dim2>* r) const
	{
	coarseSum(idim, n, r);
	if (n > 1U)
	*r /= static_cast<double>(n);
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template <typename Num2, unsigned Len2, unsigned Dim2>
	void ArrayND<Numeric,Len,Dim>::coarseSum2(
	const unsigned idim, ArrayND<Num2,Len2,Dim2>* r) const
	{
	Num2* result = r->data_;
	const unsigned long *rstrides = r->strides_;
	const unsigned long resultLen = r->length();

	// Cycle over result elements
	for (unsigned long ir=0; ir<resultLen; ++ir)
	{
	// Convert the result index into multidimensional index.
	// Simultaneously come up with the indices for this array.
	unsigned long l = ir, idx0 = 0UL, idx1 = 0UL;
	for (unsigned i=0; i<dim_; ++i)
	{
	const unsigned long idx = l / rstrides[i];
	l -= (idx * rstrides[i]);
	if (i == idim)
	{
	idx0 += 2ULidxstrides_[i];
	idx1 += (2ULidx + 1UL)strides_[i];
	}
	else
	{
	idx0 += idx*strides_[i];
	idx1 += idx*strides_[i];
	}
	}
	result[ir] = data_[idx0] + data_[idx1];
	}
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template <typename Num2, unsigned Len2, unsigned Dim2>
	void ArrayND<Numeric,Len,Dim>::coarseSum3(
	const unsigned idim, ArrayND<Num2,Len2,Dim2>* r) const
	{
	Num2* result = r->data_;
	const unsigned long *rstrides = r->strides_;
	const unsigned long resultLen = r->length();

	// Cycle over result elements
	for (unsigned long ir=0; ir<resultLen; ++ir)
	{
	// Convert the result index into multidimensional index.
	// Simultaneously come up with the indices for this array.
	unsigned long l = ir, idx0 = 0UL, idx1 = 0UL, idx2 = 0UL;
	for (unsigned i=0; i<dim_; ++i)
	{
	const unsigned long idx = l / rstrides[i];
	l -= (idx * rstrides[i]);
	if (i == idim)
	{
	idx0 += 3ULidxstrides_[i];
	idx1 += (3ULidx + 1UL)strides_[i];
	idx2 += (3ULidx + 2UL)strides_[i];
	}
	else
	{
	idx0 += idx*strides_[i];
	idx1 += idx*strides_[i];
	idx2 += idx*strides_[i];
	}
	}
	result[ir] = data_[idx0] + data_[idx1] + data_[idx2];
	}
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template <typename Num2, unsigned Len2, unsigned Dim2>
	void ArrayND<Numeric,Len,Dim>::coarseSumN(
	const unsigned idim, const unsigned n, ArrayND<Num2,Len2,Dim2>* r) const
	{
	const Num2 zero = Num2();
	Num2* result = r->data_;
	const unsigned long *rstrides = r->strides_;
	const unsigned long resultLen = r->length();

	unsigned long localidx[64];
	unsigned long* idxb = makeBuffer(n, localidx, 64);

	// Cycle over result elements
	for (unsigned long ir=0; ir<resultLen; ++ir)
	{
	// Convert the result index into multidimensional index.
	// Simultaneously come up with the indices for this array.
	unsigned long l = ir;
	for (unsigned iclean=0; iclean<n; ++iclean)
	idxb[iclean] = 0UL;
	for (unsigned i=0; i<dim_; ++i)
	{
	const unsigned long idx = l / rstrides[i];
	l -= (idx * rstrides[i]);
	if (i == idim)
	for (unsigned iclean=0; iclean<n; ++iclean)
	idxb[iclean] += (nidx + iclean)strides_[i];
	else
	for (unsigned iclean=0; iclean<n; ++iclean)
	idxb[iclean] += idx*strides_[i];
	}

	result[ir] = zero;
	for (unsigned iclean=0; iclean<n; ++iclean)
	result[ir] += data_[idxb[iclean]];
	}

	destroyBuffer(idxb, localidx);
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template<class Functor>
	ArrayND<Numeric,Len,Dim>& ArrayND<Numeric,Len,Dim>::functorFill(Functor f)
	{
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"functorFill\"");
	if (dim_)
	{
	unsigned localIndex[Dim];
	unsigned* index = makeBuffer(dim_, localIndex, Dim);
	functorFillLoop(0U, 0UL, f, index);
	destroyBuffer(index, localIndex);
	}
	else
	localData_[0] = static_cast<Numeric>(
	f(static_cast<unsigned*>(0), 0U));
	return *this;
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template <typename Num2, unsigned Len2, unsigned Dim2>
	bool ArrayND<Numeric,Len,Dim>::isClose(
	const ArrayND<Num2,Len2,Dim2>& r, const double eps) const
	{
	if (eps < 0.0) throw std::domain_error(
	"In npstat::ArrayND::isClose: tolerance must not be negative");
	if (!isShapeCompatible(r)) throw std::invalid_argument(
	"In npstat::ArrayND::isClose: incompatible argument array shape");
	if (dim_)
	{
	for (unsigned long i=0; i<len_; ++i)
	{
	const Numeric rval = r.data_[i];
	if (static_cast<double>(absDifference(data_[i], rval)) > eps)
	return false;
	}
	}
	else
	{
	const Numeric rval = r.localData_[0];
	if (static_cast<double>(absDifference(localData_[0], rval)) > eps)
	return false;
	}
	return true;
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template<typename Num2, unsigned Len2, unsigned Dim2>
	ArrayND<Numeric,Len,Dim> ArrayND<Numeric,Len,Dim>::outer(
	const ArrayND<Num2,Len2,Dim2>& r) const
	{
	return ArrayND<Numeric,Len,Dim>(*this, r);
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	void ArrayND<Numeric,Len,Dim>::contractLoop(
	unsigned thisLevel, const unsigned resLevel,
	const unsigned pos1, const unsigned pos2,
	unsigned long idxThis, unsigned long idxRes, ArrayND& result) const
	{
	while (thisLevel == pos1 \|\| thisLevel == pos2)
	++thisLevel;
	assert(thisLevel < dim_);

	if (resLevel == result.dim_ - 1)
	{
	const unsigned ncontract = shape_[pos1];
	const unsigned imax = result.shape_[resLevel];
	const unsigned long stride = strides_[pos1] + strides_[pos2];
	for (unsigned i=0; i<imax; ++i)
	{
	const Numeric* tmp = data_ + (idxThis + i*strides_[thisLevel]);
	Numeric sum = Numeric();
	for (unsigned j=0; j<ncontract; ++j)
	sum += tmp[j*stride];
	result.data_[idxRes + i] = sum;
	}
	}
	else
	{
	const unsigned imax = result.shape_[resLevel];
	assert(imax == shape_[thisLevel]);
	for (unsigned i=0; i<imax; ++i)
	{
	contractLoop(thisLevel+1, resLevel+1, pos1, pos2,
	idxThis, idxRes, result);
	idxThis += strides_[thisLevel];
	idxRes += result.strides_[resLevel];
	}
	}
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	ArrayND<Numeric,Len,Dim> ArrayND<Numeric,Len,Dim>::contract(
	const unsigned pos1, const unsigned pos2) const
	{
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"contract\"");
	if (!(pos1 < dim_ && pos2 < dim_ && pos1 != pos2))
	throw std::invalid_argument("In npstat::ArrayND::contract: "
	"incompatible contraction indices");
	if (shape_[pos1] != shape_[pos2])
	throw std::invalid_argument(
	"In npstat::ArrayND::contract: incompatible "
	"length of contracted dimensions");

	// Construct the new shape
	unsigned newshapeBuf[Dim];
	unsigned* newshape = makeBuffer(dim_ - 2, newshapeBuf, Dim);
	unsigned ishap = 0;
	for (unsigned i=0; i<dim_; ++i)
	if (i != pos1 && i != pos2)
	newshape[ishap++] = shape_[i];

	// Form the result array
	ArrayND<Numeric,Len,Dim> result(newshape, ishap);
	if (ishap)
	contractLoop(0, 0, pos1, pos2, 0UL, 0UL, result);
	else
	{
	// We are just calculating the trace
	Numeric sum = Numeric();
	const unsigned imax = shape_[0];
	const unsigned long stride = strides_[0] + strides_[1];
	for (unsigned i=0; i<imax; ++i)
	sum += data_[i*stride];
	result() = sum;
	}

	destroyBuffer(newshape, newshapeBuf);
	return result;
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	void ArrayND<Numeric,Len,Dim>::transposeLoop(
	const unsigned level, const unsigned pos1, const unsigned pos2,
	unsigned long idxThis, unsigned long idxRes, ArrayND& result) const
	{
	const unsigned imax = shape_[level];
	const unsigned long mystride = strides_[level];
	const unsigned relevel = level == pos1 ? pos2 :
	(level == pos2 ? pos1 : level);
	const unsigned long restride = result.strides_[relevel];
	const bool ready = (level == dim_ - 1);
	for (unsigned i=0; i<imax; ++i)
	{
	if (ready)
	result.data_[idxRes] = data_[idxThis];
	else
	transposeLoop(level+1, pos1, pos2, idxThis, idxRes, result);
	idxThis += mystride;
	idxRes += restride;
	}
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template<typename Num2>
	Num2 ArrayND<Numeric,Len,Dim>::sum() const
	{
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"sum\"");
	Num2 sum = Num2();
	for (unsigned long i=0; i<len_; ++i)
	sum += data_[i];
	return sum;
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template<typename Num2>
	Num2 ArrayND<Numeric,Len,Dim>::sumsq() const
	{
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"sumsq\"");
	Num2 sum = Num2();
	for (unsigned long i=0; i<len_; ++i)
	{
	const Num2 absval = absValue(data_[i]);
	sum += absval*absval;
	}
	return sum;
	}

	template<class Arr>
	typename Arr::value_type arrayMin(const Arr& arr)
	{
	typedef typename Arr::value_type Numeric;

	if (!arr.isShapeKnown()) throw std::invalid_argument(
	"In npstat::arrayMin: uninitialized array argument");
	if (arr.rank())
	{
	const unsigned long len = arr.length();
	const Numeric* minValue = &arr.linearValue(0);
	for (unsigned long i=1; i<len; ++i)
	{
	const Numeric& v = arr.linearValue(i);
	if (ComplexComparesAbs<Numeric>::less(v, *minValue))
	minValue = &v;
	}
	return *minValue;
	}
	else
	return arr();
	}

	template<class Arr>
	ArrayShape arrayArgmin(const Arr& arr)
	{
	typedef typename Arr::value_type Numeric;

	if (!arr.isShapeKnown()) throw std::invalid_argument(
	"In npstat::arrayMin: uninitialized array argument");
	ArrayShape sh(arr.rank());
	if (arr.rank())
	{
	const unsigned long len = arr.length();
	const Numeric* minValue = &arr.linearValue(0);
	unsigned long imin = 0;
	for (unsigned long i=1; i<len; ++i)
	{
	const Numeric& v = arr.linearValue(i);
	if (ComplexComparesAbs<Numeric>::less(v, *minValue))
	{
	minValue = &v;
	imin = i;
	}
	}
	arr.convertLinearIndex(imin, &sh[0], arr.rank());
	}
	return sh;
	}

	template<class Arr>
	typename Arr::value_type arrayMax(const Arr& arr)
	{
	typedef typename Arr::value_type Numeric;

	if (!arr.isShapeKnown()) throw std::invalid_argument(
	"In npstat::arrayMax: uninitialized array argument");
	if (arr.rank())
	{
	const unsigned long len = arr.length();
	const Numeric* maxValue = &arr.linearValue(0);
	for (unsigned long i=1; i<len; ++i)
	{
	const Numeric& v = arr.linearValue(i);
	if (ComplexComparesAbs<Numeric>::less(*maxValue, v))
	maxValue = &v;
	}
	return *maxValue;
	}
	else
	return arr();
	}

	template<class Arr>
	ArrayShape arrayArgmax(const Arr& arr)
	{
	typedef typename Arr::value_type Numeric;

	if (!arr.isShapeKnown()) throw std::invalid_argument(
	"In npstat::arrayMax: uninitialized array argument");
	ArrayShape sh(arr.rank());
	if (arr.rank())
	{
	const unsigned long len = arr.length();
	const Numeric* maxValue = &arr.linearValue(0);
	unsigned long imax = 0;
	for (unsigned long i=1; i<len; ++i)
	{
	const Numeric& v = arr.linearValue(i);
	if (ComplexComparesAbs<Numeric>::less(*maxValue, v))
	{
	maxValue = &v;
	imax = i;
	}
	}
	arr.convertLinearIndex(imax, &sh[0], arr.rank());
	}
	return sh;
	}

	template<class Arr>
	std::pair<typename Arr::value_type, typename Arr::value_type>
	arrayMinMax(const Arr& arr)
	{
	typedef typename Arr::value_type Numeric;

	if (!arr.isShapeKnown()) throw std::invalid_argument(
	"In npstat::arrayMinMax: uninitialized array argument");
	if (arr.rank())
	{
	const unsigned long len = arr.length();
	const Numeric* maxValue = &arr.linearValue(0);
	const Numeric* minValue = maxValue;
	for (unsigned long i=1; i<len; ++i)
	{
	const Numeric& v = arr.linearValue(i);
	if (ComplexComparesAbs<Numeric>::less(*maxValue, v))
	maxValue = &v;
	if (ComplexComparesAbs<Numeric>::less(v, *minValue))
	minValue = &v;
	}
	return std::pair<Numeric, Numeric>(minValue, maxValue);
	}
	else
	return std::pair<Numeric, Numeric>(arr(), arr());
	}

	// Faster function for 2d transpose
	template<typename Numeric, unsigned Len, unsigned Dim>
	ArrayND<Numeric,Len,Dim> ArrayND<Numeric,Len,Dim>::transpose() const
	{
	if (dim_ != 2) throw std::invalid_argument(
	"npstat::ArrayND::transpose method "
	"can not be used with arrays of rank other than 2");
	unsigned newshape[2];
	newshape[0] = shape_[1];
	newshape[1] = shape_[0];
	ArrayND<Numeric,Len,Dim> result(newshape, dim_);
	const unsigned imax = shape_[0];
	const unsigned jmax = shape_[1];
	for (unsigned i=0; i<imax; ++i)
	for (unsigned j=0; j<jmax; ++j)
	result.data_[jimax + i] = data_[ijmax + j];
	return result;
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template <typename Accumulator>
	ArrayND<Numeric,Len,Dim> ArrayND<Numeric,Len,Dim>::derivative(
	const double inscale) const
	{
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"derivative\"");
	if (!dim_) throw std::invalid_argument(
	"npstat::ArrayND::derivative method "
	"can not be used with array of 0 rank");

	const typename ProperDblFromCmpl<Accumulator>::type scale = inscale;
	const unsigned maxdim = CHAR_BIT*sizeof(unsigned long);
	if (dim_ >= maxdim) throw std::invalid_argument(
	"In npstat::ArrayND::derivative: array rank is too large");
	const unsigned long maxcycle = 1UL << dim_;

	ArrayShape sh;
	sh.reserve(dim_);
	for (unsigned i=0; i<dim_; ++i)
	{
	if (shape_[i] <= 1U)
	throw std::invalid_argument(
	"In npstat::ArrayND::derivative: in some dimendions "
	"array size is too small");
	sh.push_back(shape_[i] - 1U);
	}

	ArrayND result(sh);
	const unsigned long rLen = result.length();
	for (unsigned long ilin=0; ilin<rLen; ++ilin)
	{
	result.convertLinearIndex(ilin, &sh[0], dim_);

	Accumulator deriv = Accumulator();
	for (unsigned long icycle=0UL; icycle<maxcycle; ++icycle)
	{
	unsigned long icell = 0UL;
	unsigned n1 = 0U;
	for (unsigned i=0; i<dim_; ++i)
	{
	if (icycle & (1UL << i))
	{
	++n1;
	icell += strides_[i]*(sh[i] + 1);
	}
	else
	icell += strides_[i]*sh[i];
	}
	if ((dim_ - n1) % 2U)
	deriv -= data_[icell];
	else
	deriv += data_[icell];
	}
	result.data_[ilin] = static_cast<Numeric>(deriv*scale);
	}

	return result;
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template <typename Accumulator>
	Accumulator ArrayND<Numeric,Len,Dim>::sumBelowLoop(
	const unsigned level, unsigned long idx0,
	const unsigned* limit) const
	{
	Accumulator cdf = Accumulator();
	const unsigned imax = limit[level] + 1U;
	if (level == dim_ - 1)
	{
	Numeric* base = data_ + idx0;
	for (unsigned i=0; i<imax; ++i)
	cdf += base[i];
	}
	else
	{
	const unsigned long stride = strides_[level];
	for (unsigned i=0; i<imax; ++i, idx0+=stride)
	cdf += sumBelowLoop<Accumulator>(level+1, idx0, limit);
	}
	return cdf;
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template <typename Accumulator>
	Accumulator ArrayND<Numeric,Len,Dim>::cdfValue(
	const unsigned *index, const unsigned indexLen) const
	{
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"cdfValue\"");
	if (!dim_) throw std::invalid_argument(
	"npstat::ArrayND::cdfValue method "
	"can not be used with array of 0 rank");
	if (indexLen != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::cdfValue: incompatible index length");
	for (unsigned i=0; i<indexLen; ++i)
	if (index[i] >= shape_[i])
	throw std::out_of_range(
	"In npstat::ArrayND::cdfValue: index out of range");
	return sumBelowLoop<Accumulator>(0, 0U, index);
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template <typename Accumulator>
	ArrayND<Numeric,Len,Dim> ArrayND<Numeric,Len,Dim>::cdfArray(
	const double inscale) const
	{
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"cdfArray\"");
	if (!dim_) throw std::invalid_argument(
	"npstat::ArrayND::cdfArray method "
	"can not be used with array of 0 rank");

	const proper_double scale = inscale;
	const unsigned maxdim = CHAR_BIT*sizeof(unsigned long);
	if (dim_ >= maxdim)
	throw std::invalid_argument(
	"In npstat::ArrayND::cdfArray: array rank is too large");
	const unsigned long maxcycle = 1UL << dim_;

	ArrayShape sh;
	sh.reserve(dim_);
	for (unsigned i=0; i<dim_; ++i)
	sh.push_back(shape_[i] + 1U);

	ArrayND<Accumulator> result(sh);

	unsigned* psh = &sh[0];
	const unsigned long len = result.length();
	for (unsigned long ipre=0; ipre<len; ++ipre)
	{
	result.convertLinearIndex(ipre, psh, dim_);
	Accumulator deriv = Accumulator();
	bool has0 = false;
	for (unsigned i=0; i<dim_; ++i)
	if (psh[i]-- == 0U)
	{
	has0 = true;
	break;
	}
	if (!has0)
	{
	for (unsigned long icycle=0UL; icycle<maxcycle; ++icycle)
	{
	unsigned long icell = 0UL;
	unsigned n1 = 0U;
	for (unsigned i=0; i<dim_; ++i)
	{
	if (icycle & (1UL << i))
	{
	++n1;
	icell += result.strides_[i]*(psh[i] + 1);
	}
	else
	icell += result.strides_[i]*psh[i];
	}
	if (n1 < dim_)
	{
	if ((dim_ - n1) % 2U)
	deriv += result.data_[icell];
	else
	deriv -= result.data_[icell];
	}
	}
	deriv += static_cast<Accumulator>(value(psh, dim_)*scale);
	}
	result.data_[ipre] = deriv;
	}

	// The "return" will convert Accumulator type into Numeric
	return result;
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	ArrayND<Numeric,Len,Dim> ArrayND<Numeric,Len,Dim>::transpose(
	const unsigned pos1, const unsigned pos2) const
	{
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"transpose\"");
	if (!(pos1 < dim_ && pos2 < dim_ && pos1 != pos2))
	throw std::invalid_argument("In npstat::ArrayND::transpose: "
	"incompatible transposition indices");
	if (dim_ == 2)
	return transpose();
	else
	{
	// Construct the new shape
	unsigned newshapeBuf[Dim];
	unsigned *newshape = makeBuffer(dim_, newshapeBuf, Dim);
	copyBuffer(newshape, shape_, dim_);
	std::swap(newshape[pos1], newshape[pos2]);

	// Form the result array
	ArrayND<Numeric,Len,Dim> result(newshape, dim_);

	// Fill the result array
	transposeLoop(0, pos1, pos2, 0UL, 0UL, result);

	destroyBuffer(newshape, newshapeBuf);
	return result;
	}
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template<typename Num2, unsigned Len2, unsigned Dim2>
	void ArrayND<Numeric,Len,Dim>::mirror(
	const unsigned* mirrorDims, unsigned mirrorLen,
	ArrayND<Num2, Len2, Dim2>* out) const
	{
	assert(out);
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"mirror\"");
	if (mirrorLen > CHAR_BIT*sizeof(unsigned long)) throw std::invalid_argument(
	"In npstat::ArrayND::mirror: too many directions to mirror");
	if (!out->shapeIsKnown_)
	*out = ArrayND<Num2, Len2, Dim2>(shape_, dim_);
	if (!isShapeCompatible(*out)) throw std::invalid_argument(
	"In npstat::ArrayND::mirror: incompatible output array shape");

	unsigned nMirror = 0;
	unsigned toMirror[CHAR_BIT*sizeof(unsigned long)];
	if (mirrorLen > 0U)
	{
	assert(mirrorDims);
	std::set<unsigned> dimset(mirrorDims, mirrorDims+mirrorLen);
	for (std::set<unsigned>::const_iterator it = dimset.begin();
	it != dimset.end(); ++it)
	if (*it < dim_)
	toMirror[nMirror++] = *it;
	}
	if (!nMirror)
	{
	out = this;
	return;
	}

	unsigned idx[CHAR_BIT*sizeof(unsigned long)];
	for (unsigned long el=0; el<len_; ++el)
	{
	convertLinearIndex(el, &idx[0], dim_);
	for (unsigned im=0; im<nMirror; ++im)
	{
	const unsigned iswap = toMirror[im];
	idx[iswap] = shape_[iswap] - idx[iswap] - 1U;
	}
	out->value(&idx[0], dim_) = data_[el];
	}
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template<typename Num2, unsigned Len2, unsigned Dim2>
	void ArrayND<Numeric,Len,Dim>::multiMirror(
	ArrayND<Num2, Len2, Dim2>* out) const
	{
	assert(out);
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"multiMirror\"");
	if (!out->shapeIsKnown_)
	*out = ArrayND<Num2, Len2, Dim2>(doubleShape(shape()));
	if (dim_ != out->dim_) throw std::invalid_argument(
	"In npstat::ArrayND::multiMirror: incompatible argument array rank");

	if (dim_)
	{
	const unsigned *dshape = out->shape_;
	for (unsigned i=0; i<dim_; ++i)
	if (dshape[i] != shape_[i]*2U) throw std::invalid_argument(
	"In npstat::ArrayND::multiMirror: "
	"incompatible argument array shape");

	if (dim_ >= CHAR_BIT*sizeof(unsigned long))
	throw std::invalid_argument(
	"In npstat::ArrayND::multiMirror: "
	"array rank is too large");
	const unsigned long maxcycle = 1UL << dim_;
	std::vector<unsigned> indexbuf(dim_*2U);
	unsigned* idx = &indexbuf[0];
	unsigned* mirror = idx + dim_;

	for (unsigned long ipt=0; ipt<len_; ++ipt)
	{
	unsigned long l = ipt;
	for (unsigned i=0; i<dim_; ++i)
	{
	idx[i] = l / strides_[i];
	l -= (idx[i] * strides_[i]);
	}
	for (unsigned long icycle=0UL; icycle<maxcycle; ++icycle)
	{
	for (unsigned i=0; i<dim_; ++i)
	{
	if (icycle & (1UL << i))
	mirror[i] = dshape[i] - idx[i] - 1U;
	else
	mirror[i] = idx[i];
	}
	out->value(mirror, dim_) = data_[ipt];
	}
	}
	}
	else
	out->localData_[0] = static_cast<Num2>(localData_[0]);
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template<typename Num2, unsigned Len2, unsigned Dim2>
	void ArrayND<Numeric,Len,Dim>::rotate(
	const unsigned* shifts, const unsigned lenShifts,
	ArrayND<Num2, Len2, Dim2>* rotated) const
	{
	assert(rotated);
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"rotate\"");
	// Can't rotate into itself -- it will be a mess
	if ((void)rotated == (void)this) throw std::invalid_argument(
	"In npstat::ArrayND::rotate: can not rotate array into itself");
	if (!rotated->shapeIsKnown_)
	rotated = this;
	if (dim_ != rotated->dim_) throw std::invalid_argument(
	"In npstat::ArrayND::rotate: incompatible argument array rank");
	if (lenShifts != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::rotate: incompatible dimensionality of shifts");

	if (dim_)
	{
	assert(shifts);
	if (dim_ > CHAR_BIT*sizeof(unsigned long))
	throw std::invalid_argument(
	"In npstat::ArrayND::rotate: array rank is too large");
	unsigned buf[CHAR_BIT*sizeof(unsigned long)];
	clearBuffer(buf, dim_);
	(const_cast<ArrayND*>(this))->flatCircularLoop(
	0U, 0UL, 0UL, buf, shape_, shifts,
	*rotated, scast_assign_right<Numeric,Num2>());
	}
	else
	rotated->localData_[0] = static_cast<Num2>(localData_[0]);
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template<typename Num2, unsigned Len2, unsigned Dim2>
	void ArrayND<Numeric,Len,Dim>::dotProductLoop(
	const unsigned level, unsigned long idx0,
	unsigned long idx1, unsigned long idx2,
	const ArrayND<Num2, Len2, Dim2>& r,
	ArrayND& result) const
	{
	// idx0 -- this object
	// idx1 -- dot product argument
	// idx2 -- result
	if (level == result.dim_)
	{
	Numeric sum = Numeric();
	const unsigned imax = r.shape_[0];
	const unsigned rstride = r.strides_[0];
	const Numeric* l = data_ + idx0;
	const Num2* ri = r.data_ + idx1;
	for (unsigned i=0; i<imax; ++i)
	sum += l[i]ri[irstride];
	result.data_[idx2] = sum;
	}
	else
	{
	const unsigned imax = result.shape_[level];
	for (unsigned i=0; i<imax; ++i)
	{
	dotProductLoop(level+1, idx0, idx1, idx2, r, result);
	idx2 += result.strides_[level];
	if (level < dim_ - 1)
	idx0 += strides_[level];
	else
	idx1 += r.strides_[level + 2 - dim_];
	}
	}
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	template<typename Num2, unsigned Len2, unsigned Dim2>
	ArrayND<Numeric,Len,Dim> ArrayND<Numeric,Len,Dim>::dot(
	const ArrayND<Num2,Len2,Dim2>& r) const
	{
	if (!dim_) throw std::invalid_argument(
	"npstat::ArrayND::dot method "
	"can not be used with array of 0 rank");
	if (!r.dim_) throw std::invalid_argument(
	"npstat::ArrayND::dot method "
	"can not be used with argument array of 0 rank");
	if (shape_[dim_ - 1] != r.shape_[0]) throw std::invalid_argument(
	"In npstat::ArrayND::dot: incompatible argument array shape");

	if (dim_ == 1 && r.dim_ == 1)
	{
	// Special case: the result is of 0 rank
	ArrayND<Numeric,Len,Dim> result(static_cast<unsigned*>(0), 0U);
	Numeric sum = Numeric();
	const unsigned imax = shape_[0];
	for (unsigned i=0; i<imax; ++i)
	sum += data_[i]*r.data_[i];
	result() = sum;
	return result;
	}
	else
	{
	unsigned newshapeBuf[2*Dim];
	unsigned newshape = makeBuffer(dim_+r.dim_-2, newshapeBuf, 2Dim);
	copyBuffer(newshape, shape_, dim_-1);
	copyBuffer(newshape+(dim_-1), r.shape_+1, r.dim_-1);
	ArrayND<Numeric,Len,Dim> result(newshape, dim_+r.dim_-2);

	dotProductLoop(0U, 0UL, 0UL, 0UL, r, result);

	destroyBuffer(newshape, newshapeBuf);
	return result;
	}
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline unsigned ArrayND<Numeric,Len,Dim>::span(const unsigned dim) const
	{
	if (dim >= dim_)
	throw std::out_of_range(
	"In npstat::ArrayND::span: dimension number is out of range");
	return shape_[dim];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	unsigned ArrayND<Numeric,Len,Dim>::maximumSpan() const
	{
	unsigned maxspan = 0;
	for (unsigned i=0; i<dim_; ++i)
	if (shape_[i] > maxspan)
	maxspan = shape_[i];
	return maxspan;
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	unsigned ArrayND<Numeric,Len,Dim>::minimumSpan() const
	{
	if (dim_ == 0)
	return 0U;
	unsigned minspan = shape_[0];
	for (unsigned i=1; i<dim_; ++i)
	if (shape_[i] < minspan)
	minspan = shape_[i];
	return minspan;
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline const Numeric& ArrayND<Numeric,Len,Dim>::cl() const
	{
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"cl\"");
	if (dim_) throw std::invalid_argument(
	"In npstat::ArrayND::cl: wrong # of args (not rank 0 array)");
	return localData_[0];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline const Numeric& ArrayND<Numeric,Len,Dim>::cl(
	const double i0) const
	{
	if (1U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::cl: wrong # of args (not rank 1 array)");
	return data_[coordToIndex(i0, 0)];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline const Numeric& ArrayND<Numeric,Len,Dim>::cl(
	const double i0,
	const double i1) const
	{
	if (2U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::cl: wrong # of args (not rank 2 array)");
	return data_[coordToIndex(i0, 0)*strides_[0] +
	coordToIndex(i1, 1)];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline const Numeric& ArrayND<Numeric,Len,Dim>::cl(
	const double i0,
	const double i1,
	const double i2) const
	{
	if (3U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::cl: wrong # of args (not rank 3 array)");
	return data_[coordToIndex(i0, 0)*strides_[0] +
	coordToIndex(i1, 1)*strides_[1] +
	coordToIndex(i2, 2)];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline const Numeric& ArrayND<Numeric,Len,Dim>::cl(
	const double i0,
	const double i1,
	const double i2,
	const double i3) const
	{
	if (4U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::cl: wrong # of args (not rank 4 array)");
	return data_[coordToIndex(i0, 0)*strides_[0] +
	coordToIndex(i1, 1)*strides_[1] +
	coordToIndex(i2, 2)*strides_[2] +
	coordToIndex(i3, 3)];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline const Numeric& ArrayND<Numeric,Len,Dim>::cl(
	const double i0,
	const double i1,
	const double i2,
	const double i3,
	const double i4) const
	{
	if (5U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::cl: wrong # of args (not rank 5 array)");
	return data_[coordToIndex(i0, 0)*strides_[0] +
	coordToIndex(i1, 1)*strides_[1] +
	coordToIndex(i2, 2)*strides_[2] +
	coordToIndex(i3, 3)*strides_[3] +
	coordToIndex(i4, 4)];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline const Numeric& ArrayND<Numeric,Len,Dim>::cl(
	const double i0,
	const double i1,
	const double i2,
	const double i3,
	const double i4,
	const double i5) const
	{
	if (6U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::cl: wrong # of args (not rank 6 array)");
	return data_[coordToIndex(i0, 0)*strides_[0] +
	coordToIndex(i1, 1)*strides_[1] +
	coordToIndex(i2, 2)*strides_[2] +
	coordToIndex(i3, 3)*strides_[3] +
	coordToIndex(i4, 4)*strides_[4] +
	coordToIndex(i5, 5)];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline const Numeric& ArrayND<Numeric,Len,Dim>::cl(
	const double i0,
	const double i1,
	const double i2,
	const double i3,
	const double i4,
	const double i5,
	const double i6) const
	{
	if (7U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::cl: wrong # of args (not rank 7 array)");
	return data_[coordToIndex(i0, 0)*strides_[0] +
	coordToIndex(i1, 1)*strides_[1] +
	coordToIndex(i2, 2)*strides_[2] +
	coordToIndex(i3, 3)*strides_[3] +
	coordToIndex(i4, 4)*strides_[4] +
	coordToIndex(i5, 5)*strides_[5] +
	coordToIndex(i6, 6)];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline const Numeric& ArrayND<Numeric,Len,Dim>::cl(
	const double i0,
	const double i1,
	const double i2,
	const double i3,
	const double i4,
	const double i5,
	const double i6,
	const double i7) const
	{
	if (8U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::cl: wrong # of args (not rank 8 array)");
	return data_[coordToIndex(i0, 0)*strides_[0] +
	coordToIndex(i1, 1)*strides_[1] +
	coordToIndex(i2, 2)*strides_[2] +
	coordToIndex(i3, 3)*strides_[3] +
	coordToIndex(i4, 4)*strides_[4] +
	coordToIndex(i5, 5)*strides_[5] +
	coordToIndex(i6, 6)*strides_[6] +
	coordToIndex(i7, 7)];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline const Numeric& ArrayND<Numeric,Len,Dim>::cl(
	const double i0,
	const double i1,
	const double i2,
	const double i3,
	const double i4,
	const double i5,
	const double i6,
	const double i7,
	const double i8) const
	{
	if (9U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::cl: wrong # of args (not rank 9 array)");
	return data_[coordToIndex(i0, 0)*strides_[0] +
	coordToIndex(i1, 1)*strides_[1] +
	coordToIndex(i2, 2)*strides_[2] +
	coordToIndex(i3, 3)*strides_[3] +
	coordToIndex(i4, 4)*strides_[4] +
	coordToIndex(i5, 5)*strides_[5] +
	coordToIndex(i6, 6)*strides_[6] +
	coordToIndex(i7, 7)*strides_[7] +
	coordToIndex(i8, 8)];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline const Numeric& ArrayND<Numeric,Len,Dim>::cl(
	const double i0,
	const double i1,
	const double i2,
	const double i3,
	const double i4,
	const double i5,
	const double i6,
	const double i7,
	const double i8,
	const double i9) const
	{
	if (10U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::cl: wrong # of args (not rank 10 array)");
	return data_[coordToIndex(i0, 0)*strides_[0] +
	coordToIndex(i1, 1)*strides_[1] +
	coordToIndex(i2, 2)*strides_[2] +
	coordToIndex(i3, 3)*strides_[3] +
	coordToIndex(i4, 4)*strides_[4] +
	coordToIndex(i5, 5)*strides_[5] +
	coordToIndex(i6, 6)*strides_[6] +
	coordToIndex(i7, 7)*strides_[7] +
	coordToIndex(i8, 8)*strides_[8] +
	coordToIndex(i9, 9)];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline Numeric& ArrayND<Numeric,Len,Dim>::cl()
	{
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"cl\"");
	if (dim_) throw std::invalid_argument(
	"In npstat::ArrayND::cl: wrong # of args (not rank 0 array)");
	return localData_[0];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline Numeric& ArrayND<Numeric,Len,Dim>::cl(
	const double i0)
	{
	if (1U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::cl: wrong # of args (not rank 1 array)");
	return data_[coordToIndex(i0, 0)];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline Numeric& ArrayND<Numeric,Len,Dim>::cl(
	const double i0,
	const double i1)
	{
	if (2U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::cl: wrong # of args (not rank 2 array)");
	return data_[coordToIndex(i0, 0)*strides_[0] +
	coordToIndex(i1, 1)];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline Numeric& ArrayND<Numeric,Len,Dim>::cl(
	const double i0,
	const double i1,
	const double i2)
	{
	if (3U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::cl: wrong # of args (not rank 3 array)");
	return data_[coordToIndex(i0, 0)*strides_[0] +
	coordToIndex(i1, 1)*strides_[1] +
	coordToIndex(i2, 2)];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline Numeric& ArrayND<Numeric,Len,Dim>::cl(
	const double i0,
	const double i1,
	const double i2,
	const double i3)
	{
	if (4U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::cl: wrong # of args (not rank 4 array)");
	return data_[coordToIndex(i0, 0)*strides_[0] +
	coordToIndex(i1, 1)*strides_[1] +
	coordToIndex(i2, 2)*strides_[2] +
	coordToIndex(i3, 3)];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline Numeric& ArrayND<Numeric,Len,Dim>::cl(
	const double i0,
	const double i1,
	const double i2,
	const double i3,
	const double i4)
	{
	if (5U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::cl: wrong # of args (not rank 5 array)");
	return data_[coordToIndex(i0, 0)*strides_[0] +
	coordToIndex(i1, 1)*strides_[1] +
	coordToIndex(i2, 2)*strides_[2] +
	coordToIndex(i3, 3)*strides_[3] +
	coordToIndex(i4, 4)];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline Numeric& ArrayND<Numeric,Len,Dim>::cl(
	const double i0,
	const double i1,
	const double i2,
	const double i3,
	const double i4,
	const double i5)
	{
	if (6U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::cl: wrong # of args (not rank 6 array)");
	return data_[coordToIndex(i0, 0)*strides_[0] +
	coordToIndex(i1, 1)*strides_[1] +
	coordToIndex(i2, 2)*strides_[2] +
	coordToIndex(i3, 3)*strides_[3] +
	coordToIndex(i4, 4)*strides_[4] +
	coordToIndex(i5, 5)];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline Numeric& ArrayND<Numeric,Len,Dim>::cl(
	const double i0,
	const double i1,
	const double i2,
	const double i3,
	const double i4,
	const double i5,
	const double i6)
	{
	if (7U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::cl: wrong # of args (not rank 7 array)");
	return data_[coordToIndex(i0, 0)*strides_[0] +
	coordToIndex(i1, 1)*strides_[1] +
	coordToIndex(i2, 2)*strides_[2] +
	coordToIndex(i3, 3)*strides_[3] +
	coordToIndex(i4, 4)*strides_[4] +
	coordToIndex(i5, 5)*strides_[5] +
	coordToIndex(i6, 6)];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline Numeric& ArrayND<Numeric,Len,Dim>::cl(
	const double i0,
	const double i1,
	const double i2,
	const double i3,
	const double i4,
	const double i5,
	const double i6,
	const double i7)
	{
	if (8U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::cl: wrong # of args (not rank 8 array)");
	return data_[coordToIndex(i0, 0)*strides_[0] +
	coordToIndex(i1, 1)*strides_[1] +
	coordToIndex(i2, 2)*strides_[2] +
	coordToIndex(i3, 3)*strides_[3] +
	coordToIndex(i4, 4)*strides_[4] +
	coordToIndex(i5, 5)*strides_[5] +
	coordToIndex(i6, 6)*strides_[6] +
	coordToIndex(i7, 7)];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline Numeric& ArrayND<Numeric,Len,Dim>::cl(
	const double i0,
	const double i1,
	const double i2,
	const double i3,
	const double i4,
	const double i5,
	const double i6,
	const double i7,
	const double i8)
	{
	if (9U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::cl: wrong # of args (not rank 9 array)");
	return data_[coordToIndex(i0, 0)*strides_[0] +
	coordToIndex(i1, 1)*strides_[1] +
	coordToIndex(i2, 2)*strides_[2] +
	coordToIndex(i3, 3)*strides_[3] +
	coordToIndex(i4, 4)*strides_[4] +
	coordToIndex(i5, 5)*strides_[5] +
	coordToIndex(i6, 6)*strides_[6] +
	coordToIndex(i7, 7)*strides_[7] +
	coordToIndex(i8, 8)];
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	inline Numeric& ArrayND<Numeric,Len,Dim>::cl(
	const double i0,
	const double i1,
	const double i2,
	const double i3,
	const double i4,
	const double i5,
	const double i6,
	const double i7,
	const double i8,
	const double i9)
	{
	if (10U != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::cl: wrong # of args (not rank 10 array)");
	return data_[coordToIndex(i0, 0)*strides_[0] +
	coordToIndex(i1, 1)*strides_[1] +
	coordToIndex(i2, 2)*strides_[2] +
	coordToIndex(i3, 3)*strides_[3] +
	coordToIndex(i4, 4)*strides_[4] +
	coordToIndex(i5, 5)*strides_[5] +
	coordToIndex(i6, 6)*strides_[6] +
	coordToIndex(i7, 7)*strides_[7] +
	coordToIndex(i8, 8)*strides_[8] +
	coordToIndex(i9, 9)];
	}

	template<typename Numeric, unsigned StackLen, unsigned StackDim>
	const char* ArrayND<Numeric,StackLen,StackDim>::classname()
	{
	static const std::string name(
	gs::template_class_name<Numeric>("npstat::ArrayND"));
	return name.c_str();
	}

	template<typename Numeric, unsigned StackLen, unsigned StackDim>
	bool ArrayND<Numeric,StackLen,StackDim>::write(std::ostream& os) const
	{
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"write\"");
	gs::write_pod_vector(os, shape());
	return !os.fail() &&
	(dim_ ? gs::write_array(os, data_, len_) :
	gs::write_item(os, localData_[0], false));
	}

	template<typename Numeric, unsigned Len, unsigned Dim>
	void ArrayND<Numeric,Len,Dim>::restore(
	const gs::ClassId& id, std::istream& in, ArrayND* array)
	{
	static const gs::ClassId current(gs::ClassId::makeId<ArrayND<Numeric,Len,Dim> >());
	current.ensureSameId(id);

	ArrayShape rshape;
	gs::read_pod_vector(in, &rshape);
	if (in.fail()) throw gs::IOReadFailure(
	"In npstat::ArrayND::restore: input stream failure (checkpoint 0)");

	assert(array);
	array->uninitialize();
	array->dim_ = rshape.size();
	array->shapeIsKnown_ = true;
	array->len_ = 1UL;
	if (array->dim_)
	{
	array->shape_ = makeBuffer(array->dim_, array->localShape_, Dim);
	for (unsigned i=0; i<array->dim_; ++i)
	{
	array->shape_[i] = rshape[i];
	assert(array->shape_[i]);
	array->len_ *= array->shape_[i];
	}
	array->buildStrides();
	array->data_ = makeBuffer(array->len_, array->localData_, Len);
	gs::read_array(in, array->data_, array->len_);
	}
	else
	gs::restore_item(in, array->localData_, false);
	if (in.fail()) throw gs::IOReadFailure(
	"In npstat::ArrayND::restore: input stream failure (checkpoint 1)");
	}

	template<typename Numeric, unsigned Len, unsigned StackDim>
	template <typename Num2, unsigned Len2, unsigned Dim2>
	void ArrayND<Numeric,Len,StackDim>::exportSubrange(
	const unsigned* corner, const unsigned lenCorner,
	ArrayND<Num2, Len2, Dim2>* out) const
	{
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"exportSubrange\"");
	if (dim_ != lenCorner) throw std::invalid_argument(
	"In npstat::ArrayND::exportSubrange: incompatible corner index length");
	assert(out);
	if (!out->shapeIsKnown_) throw std::invalid_argument(
	"In npstat::ArrayND::exportSubrange: uninitialized argument array");
	if (out->dim_ != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::exportSubrange: incompatible argument array rank");

	if (dim_)
	{
	assert(corner);
	if (dim_ > CHAR_BIT*sizeof(unsigned long))
	throw std::invalid_argument(
	"In npstat::ArrayND::exportSubrange: "
	"array rank is too large");
	unsigned toBuf[CHAR_BIT*sizeof(unsigned long)];
	clearBuffer(toBuf, dim_);
	(const_cast<ArrayND*>(this))->commonSubrangeLoop(
	0U, 0UL, 0UL, corner, out->shape_, toBuf, *out,
	scast_assign_right<Numeric,Num2>());
	}
	else
	out->localData_[0] = static_cast<Num2>(localData_[0]);
	}

	template<typename Numeric, unsigned Len, unsigned StackDim>
	template <typename Num2, unsigned Len2, unsigned Dim2>
	void ArrayND<Numeric,Len,StackDim>::importSubrange(
	const unsigned* corner, const unsigned lenCorner,
	const ArrayND<Num2, Len2, Dim2>& from)
	{
	if (!shapeIsKnown_) throw std::invalid_argument(
	"Initialize npstat::ArrayND before calling method \"importSubrange\"");
	if (dim_ != lenCorner) throw std::invalid_argument(
	"In npstat::ArrayND::importSubrange: incompatible corner index length");
	if (!from.shapeIsKnown_) throw std::invalid_argument(
	"In npstat::ArrayND::importSubrange: uninitialized argument array");
	if (from.dim_ != dim_) throw std::invalid_argument(
	"In npstat::ArrayND::importSubrange: incompatible argument array rank");

	if (dim_)
	{
	assert(corner);
	if (dim_ > CHAR_BIT*sizeof(unsigned long))
	throw std::invalid_argument(
	"In npstat::ArrayND::importSubrange: "
	"array rank is too large");
	unsigned toBuf[CHAR_BIT*sizeof(unsigned long)];
	clearBuffer(toBuf, dim_);
	commonSubrangeLoop(0U, 0UL, 0UL, corner, from.shape_, toBuf,
	const_cast<ArrayND<Num2, Len2, Dim2>&>(from),
	scast_assign_left<Numeric,Num2>());
	}
	else
	localData_[0] = static_cast<Numeric>(from.localData_[0]);
	}

	template <typename Numeric>
	inline ArrayND<Numeric> externalMemArrayND(Numeric* data,
	const unsigned* shape,
	const unsigned dim)
	{
	ArrayND<Numeric> arr;
	arr.buildExtFromShapePtr(shape, dim, data);
	return arr;
	}

	template <typename Numeric>
	inline ArrayND<Numeric> externalMemArrayND(Numeric* data,
	const ArrayShape& shape)
	{
	const unsigned dim = shape.size();
	const unsigned* sh = dim ? &shape[0] : (unsigned*)0;
	ArrayND<Numeric> arr;
	arr.buildExtFromShapePtr(sh, dim, data);
	return arr;
	}

	template <typename Numeric>
	inline ArrayND<Numeric> externalMemArrayND(Numeric* data,
	const unsigned n0)
	{
	const unsigned dim = 1U;
	unsigned sizes[dim];
	sizes[0] = n0;
	ArrayND<Numeric> arr;
	arr.buildExtFromShapePtr(sizes, dim, data);
	return arr;
	}

	template <typename Numeric>
	inline ArrayND<Numeric> externalMemArrayND(Numeric* data,
	const unsigned n0,
	const unsigned n1)
	{
	const unsigned dim = 2U;
	unsigned sizes[dim];
	sizes[0] = n0;
	sizes[1] = n1;
	ArrayND<Numeric> arr;
	arr.buildExtFromShapePtr(sizes, dim, data);
	return arr;
	}

	template <typename Numeric>
	inline ArrayND<Numeric> externalMemArrayND(Numeric* data,
	const unsigned n0,
	const unsigned n1,
	const unsigned n2)
	{
	const unsigned dim = 3U;
	unsigned sizes[dim];
	sizes[0] = n0;
	sizes[1] = n1;
	sizes[2] = n2;
	ArrayND<Numeric> arr;
	arr.buildExtFromShapePtr(sizes, dim, data);
	return arr;
	}

	template <typename Numeric>
	inline ArrayND<Numeric> externalMemArrayND(Numeric* data,
	const unsigned n0,
	const unsigned n1,
	const unsigned n2,
	const unsigned n3)
	{
	const unsigned dim = 4U;
	unsigned sizes[dim];
	sizes[0] = n0;
	sizes[1] = n1;
	sizes[2] = n2;
	sizes[3] = n3;
	ArrayND<Numeric> arr;
	arr.buildExtFromShapePtr(sizes, dim, data);
	return arr;
	}

	template <typename Numeric>
	inline ArrayND<Numeric> externalMemArrayND(Numeric* data,
	const unsigned n0,
	const unsigned n1,
	const unsigned n2,
	const unsigned n3,
	const unsigned n4)
	{
	const unsigned dim = 5U;
	unsigned sizes[dim];
	sizes[0] = n0;
	sizes[1] = n1;
	sizes[2] = n2;
	sizes[3] = n3;
	sizes[4] = n4;
	ArrayND<Numeric> arr;
	arr.buildExtFromShapePtr(sizes, dim, data);
	return arr;
	}

	template <typename Numeric>
	inline ArrayND<Numeric> externalMemArrayND(Numeric* data,
	const unsigned n0,
	const unsigned n1,
	const unsigned n2,
	const unsigned n3,
	const unsigned n4,
	const unsigned n5)
	{
	const unsigned dim = 6U;
	unsigned sizes[dim];
	sizes[0] = n0;
	sizes[1] = n1;
	sizes[2] = n2;
	sizes[3] = n3;
	sizes[4] = n4;
	sizes[5] = n5;
	ArrayND<Numeric> arr;
	arr.buildExtFromShapePtr(sizes, dim, data);
	return arr;
	}

	template <typename Numeric>
	inline ArrayND<Numeric> externalMemArrayND(Numeric* data,
	const unsigned n0,
	const unsigned n1,
	const unsigned n2,
	const unsigned n3,
	const unsigned n4,
	const unsigned n5,
	const unsigned n6)
	{
	const unsigned dim = 7U;
	unsigned sizes[dim];
	sizes[0] = n0;
	sizes[1] = n1;
	sizes[2] = n2;
	sizes[3] = n3;
	sizes[4] = n4;
	sizes[5] = n5;
	sizes[6] = n6;
	ArrayND<Numeric> arr;
	arr.buildExtFromShapePtr(sizes, dim, data);
	return arr;
	}

	template <typename Numeric>
	inline ArrayND<Numeric> externalMemArrayND(Numeric* data,
	const unsigned n0,
	const unsigned n1,
	const unsigned n2,
	const unsigned n3,
	const unsigned n4,
	const unsigned n5,
	const unsigned n6,
	const unsigned n7)
	{
	const unsigned dim = 8U;
	unsigned sizes[dim];
	sizes[0] = n0;
	sizes[1] = n1;
	sizes[2] = n2;
	sizes[3] = n3;
	sizes[4] = n4;
	sizes[5] = n5;
	sizes[6] = n6;
	sizes[7] = n7;
	ArrayND<Numeric> arr;
	arr.buildExtFromShapePtr(sizes, dim, data);
	return arr;
	}

	template <typename Numeric>
	inline ArrayND<Numeric> externalMemArrayND(Numeric* data,
	const unsigned n0,
	const unsigned n1,
	const unsigned n2,
	const unsigned n3,
	const unsigned n4,
	const unsigned n5,
	const unsigned n6,
	const unsigned n7,
	const unsigned n8)
	{
	const unsigned dim = 9U;
	unsigned sizes[dim];
	sizes[0] = n0;
	sizes[1] = n1;
	sizes[2] = n2;
	sizes[3] = n3;
	sizes[4] = n4;
	sizes[5] = n5;
	sizes[6] = n6;
	sizes[7] = n7;
	sizes[8] = n8;
	ArrayND<Numeric> arr;
	arr.buildExtFromShapePtr(sizes, dim, data);
	return arr;
	}

	template <typename Numeric>
	inline ArrayND<Numeric> externalMemArrayND(Numeric* data,
	const unsigned n0,
	const unsigned n1,
	const unsigned n2,
	const unsigned n3,
	const unsigned n4,
	const unsigned n5,
	const unsigned n6,
	const unsigned n7,
	const unsigned n8,
	const unsigned n9)
	{
	const unsigned dim = 10U;
	unsigned sizes[dim];
	sizes[0] = n0;
	sizes[1] = n1;
	sizes[2] = n2;
	sizes[3] = n3;
	sizes[4] = n4;
	sizes[5] = n5;
	sizes[6] = n6;
	sizes[7] = n7;
	sizes[8] = n8;
	sizes[9] = n9;
	ArrayND<Numeric> arr;
	arr.buildExtFromShapePtr(sizes, dim, data);
	return arr;
	}

	template <typename Numeric>
	inline CPP11_auto_ptr<ArrayND<Numeric> > allocExternalMemArrayND(
	Numeric* data, const unsigned* shape, const unsigned dim)
	{
	CPP11_auto_ptr<ArrayND<Numeric> > arr(new ArrayND<Numeric>());
	arr->buildExtFromShapePtr(shape, dim, data);
	return arr;
	}

	template <typename Numeric>
	inline CPP11_auto_ptr<ArrayND<Numeric> > allocExternalMemArrayND(
	Numeric* data, const ArrayShape& shape)
	{
	const unsigned dim = shape.size();
	const unsigned* sh = dim ? &shape[0] : (unsigned*)0;
	CPP11_auto_ptr<ArrayND<Numeric> > arr(new ArrayND<Numeric>());
	arr->buildExtFromShapePtr(sh, dim, data);
	return arr;
	}
	}

	template <typename Numeric, unsigned StackLen, unsigned StackDim,
	typename Num2, unsigned Len2, unsigned Dim2>
	npstat::ArrayND<typename npstat::LongerType<Numeric,Num2>::type,
	npstat::BiggerUInt<StackLen,Len2>::value,
	npstat::BiggerUInt<StackDim,Dim2>::value>
	operator+(const npstat::ArrayND<Numeric,StackLen,StackDim>& l,
	const npstat::ArrayND<Num2,Len2,Dim2>& r)
	{
	if (!l.isShapeCompatible(r)) throw std::invalid_argument(
	"In ::operator+: incompatible shapes of ArrayND arguments");
	npstat::ArrayND<typename npstat::LongerType<Numeric,Num2>::type,
	npstat::BiggerUInt<StackLen,Len2>::value,
	npstat::BiggerUInt<StackDim,Dim2>::value> result(l.shape_, l.dim_);
	const unsigned long len = result.length();
	for (unsigned long i=0; i<len; ++i)
	result.data_[i] = l.data_[i] + r.data_[i];
	return result;
	}

	template <typename Numeric, unsigned StackLen, unsigned StackDim,
	typename Num2, unsigned Len2, unsigned Dim2>
	npstat::ArrayND<typename npstat::LongerType<Numeric,Num2>::type,
	npstat::BiggerUInt<StackLen,Len2>::value,
	npstat::BiggerUInt<StackDim,Dim2>::value>
	operator-(const npstat::ArrayND<Numeric,StackLen,StackDim>& l,
	const npstat::ArrayND<Num2,Len2,Dim2>& r)
	{
	if (!l.isShapeCompatible(r)) throw std::invalid_argument(
	"In ::operator-: incompatible shapes of ArrayND arguments");
	npstat::ArrayND<typename npstat::LongerType<Numeric,Num2>::type,
	npstat::BiggerUInt<StackLen,Len2>::value,
	npstat::BiggerUInt<StackDim,Dim2>::value> result(l.shape_, l.dim_);
	const unsigned long len = result.length();
	for (unsigned long i=0; i<len; ++i)
	result.data_[i] = l.data_[i] - r.data_[i];
	return result;
	}

	template <typename Numeric, unsigned StackLen, unsigned StackDim,
	typename Num2, unsigned Len2, unsigned Dim2>
	npstat::ArrayND<typename npstat::LongerType<Numeric,Num2>::type,
	npstat::BiggerUInt<StackLen,Len2>::value,
	npstat::BiggerUInt<StackDim,Dim2>::value>
	operator*(const npstat::ArrayND<Numeric,StackLen,StackDim>& l,
	const npstat::ArrayND<Num2,Len2,Dim2>& r)
	{
	if (!l.isShapeCompatible(r)) throw std::invalid_argument(
	"In ::operator*: incompatible shapes of ArrayND arguments");
	npstat::ArrayND<typename npstat::LongerType<Numeric,Num2>::type,
	npstat::BiggerUInt<StackLen,Len2>::value,
	npstat::BiggerUInt<StackDim,Dim2>::value> result(l.shape_, l.dim_);
	const unsigned long len = result.length();
	for (unsigned long i=0; i<len; ++i)
	result.data_[i] = l.data_[i]*r.data_[i];
	return result;
	}

	template <typename Numeric, unsigned StackLen, unsigned StackDim,
	typename Num2, unsigned Len2, unsigned Dim2>
	npstat::ArrayND<typename npstat::LongerType<Numeric,Num2>::type,
	npstat::BiggerUInt<StackLen,Len2>::value,
	npstat::BiggerUInt<StackDim,Dim2>::value>
	operator/(const npstat::ArrayND<Numeric,StackLen,StackDim>& l,
	const npstat::ArrayND<Num2,Len2,Dim2>& r)
	{
	if (!l.isShapeCompatible(r)) throw std::invalid_argument(
	"In ::operator/: incompatible shapes of ArrayND arguments");
	if (!r.isNonZero()) throw std::invalid_argument(
	"In ::operator/ for ArrayND: division by zero");
	npstat::ArrayND<typename npstat::LongerType<Numeric,Num2>::type,
	npstat::BiggerUInt<StackLen,Len2>::value,
	npstat::BiggerUInt<StackDim,Dim2>::value> result(l.shape_, l.dim_);
	const unsigned long len = result.length();
	for (unsigned long i=0; i<len; ++i)
	result.data_[i] = l.data_[i]/r.data_[i];
	return result;
	}

File Metadata

Mime Type: text/x-c
Expires: Tue, Sep 30, 4:48 AM (2 h, 59 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 6564826
Default Alt Text: ArrayND.icc (222 KB)

ArrayND.iccNo OneTemporaryActions

ArrayND.iccView Options

File Metadata

Event Timeline

ArrayND.icc
No OneTemporary
Actions

ArrayND.icc
View Options