Array.h

Go to the documentation of this file.

/*

Copyright (C) 1993-2018 John W. Eaton
Copyright (C) 2008-2009 Jaroslav Hajek
Copyright (C) 2010 VZLU Prague

This file is part of Octave.

Octave is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

Octave is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with Octave; see the file COPYING.  If not, see
<https://www.gnu.org/licenses/>.

*/

#if ! defined (octave_Array_h)
#define octave_Array_h 1

#include "octave-config.h"

#include <cassert>
#include <cstddef>

#include <algorithm>
#include <iosfwd>
#include <string>

#include "dim-vector.h"
#include "idx-vector.h"
#include "lo-error.h"
#include "lo-traits.h"
#include "lo-utils.h"
#include "oct-refcount.h"
#include "oct-sort.h"
#include "quit.h"

//! N Dimensional Array with copy-on-write semantics.
//!
//! The Array class is at the root of Octave.  It provides a container
//! with an arbitrary number of dimensions.  The operator () provides
//! access to individual elements via subscript and linear indexing.
//! Indexing starts at 0.  Arrays are column-major order as in Fortran.
//!
//! @code{.cc}
//! // 3 D Array with 10 rows, 20 columns, and 5 pages, filled with 7.0
//! Array<double> A Array<double (dim_vector (10, 20, 5), 7.0);
//!
//! // set value for row 0, column 10, and page 3
//! A(0, 10, 3) = 2.5;
//!
//! // get value for row 1, column 2, and page 0
//! double v = A(1, 2, 0);
//!
//! // get value for 25th element (row 4, column 3, page 1)
//! double v = A(24);
//! @endcode
//!
//! ## Notes on STL compatibility
//!
//! ### size() and length()
//!
//! To access the total number of elements in an Array, use numel()
//! which is short for number of elements and is equivalent to the
//! Octave function with same name.
//!
//! @code{.cc}
//! Array<int> A (dim_vector (10, 20, 4), 1);
//!
//! octave_idx_type n = A.numel (); // returns 800 (10x20x4)
//!
//! octave_idx_type nr = A.size (0); // returns 10 (number of rows/dimension 0)
//! octave_idx_type nc = A.size (1); // returns 20 (number of columns)
//! octave_idx_type nc = A.size (2); // returns 4 (size of dimension 3)
//! octave_idx_type l6 = A.size (6); // returns 1 (implicit singleton dimension)
//!
//! // Alternatively, get a dim_vector which represents the dimensions.
//! dim_vector dims = A.dims ();
//! @endcode
//!
//! The methods size() and length() as they exist in the STL cause
//! confusion in the context of a N dimensional array.
//!
//! The size() of an array is the length of all dimensions.  In Octave,
//! the size() function returns a row vector with the length of each
//! dimension, or the size of a specific dimension.  Only the latter is
//! present in liboctave.
//!
//! Since there is more than 1 dimension, length() would not make sense
//! without expliciting which dimension.  If the function existed, which
//! length should it return?  Octave length() function returns the length
//! of the longest dimension which is an odd definition, only useful for
//! vectors and square matrices.  The alternatives numel(), rows(),
//! columns(), and size(d) are more explict and recommended.
//!
//! ### size_type
//!
//! Array::size_type is `octave_idx_type` which is a typedef for `int`
//! or `long int`, depending whether Octave was configured for 64-bit
//! indexing.
//!
//! This is a signed integer which may cause problems when mixed with
//! STL containers.  The reason is that Octave interacts with Fortran
//! routines, providing an interface many Fortran numeric libraries.
//!
//! ## Subclasses
//!
//! The following subclasses specializations, will be of most use:
//!   - Matrix: Array<double> with only 2 dimensions
//!   - ComplexMatrix: Array<std::complex<double>> with only 2 dimensions
//!   - boolNDArray: N dimensional Array<bool>
//!   - ColumnVector: Array<double> with 1 column
//!   - string_vector: Array<std::string> with 1 column
//!   - Cell: Array<octave_value>, equivalent to an Octave cell.

template <typename T>
class
Array
{
protected:

  //! The real representation of all arrays.
  class ArrayRep
  {
  public:

    T *data;
    octave_idx_type len;
    octave::refcount<int> count;

    ArrayRep (T *d, octave_idx_type l)
      : data (new T [l]), len (l), count (1)
    {
      std::copy_n (d, l, data);
    }

    template <typename U>
    ArrayRep (U *d, octave_idx_type l)
      : data (new T [l]), len (l), count (1)
    {
      std::copy_n (d, l, data);
    }

    // Use new instead of setting data to 0 so that fortran_vec and
    // data always return valid addresses, even for zero-size arrays.

    ArrayRep (void) : data (new T [0]), len (0), count (1) { }

    explicit ArrayRep (octave_idx_type n)
      : data (new T [n]), len (n), count (1) { }

    explicit ArrayRep (octave_idx_type n, const T& val)
      : data (new T [n]), len (n), count (1)
    {
      std::fill_n (data, n, val);
    }

    ArrayRep (const ArrayRep& a)
      : data (new T [a.len]), len (a.len), count (1)
    {
      std::copy_n (a.data, a.len, data);
    }

    ~ArrayRep (void) { delete [] data; }

    octave_idx_type numel (void) const { return len; }

  private:

    // No assignment!

    ArrayRep& operator = (const ArrayRep& a);
  };

  //--------------------------------------------------------------------

public:

  void make_unique (void)
  {
    if (rep->count > 1)
      {
        ArrayRep *r = new ArrayRep (slice_data, slice_len);

        if (--rep->count == 0)
          delete rep;

        rep = r;
        slice_data = rep->data;
      }
  }

  typedef T element_type;

  typedef T value_type;

  //! Used for operator(), and returned by numel() and size()
  //! (beware: signed integer)
  typedef octave_idx_type size_type;

  typedef typename ref_param<T>::type crefT;

  typedef bool (*compare_fcn_type) (typename ref_param<T>::type,
                                    typename ref_param<T>::type);

protected:

  dim_vector dimensions;

  typename Array<T>::ArrayRep *rep;

  // Rationale:
  // slice_data is a pointer to rep->data, denoting together with slice_len the
  // actual portion of the data referenced by this Array<T> object.  This
  // allows to make shallow copies not only of a whole array, but also of
  // contiguous subranges.  Every time rep is directly manipulated, slice_data
  // and slice_len need to be properly updated.

  T *slice_data;
  octave_idx_type slice_len;

  //! slice constructor
  Array (const Array<T>& a, const dim_vector& dv,
         octave_idx_type l, octave_idx_type u)
    : dimensions (dv), rep(a.rep), slice_data (a.slice_data+l), slice_len (u-l)
  {
    rep->count++;
    dimensions.chop_trailing_singletons ();
  }

private:

  static typename Array<T>::ArrayRep *nil_rep (void);

protected:

  //! For jit support
  Array (T *sdata, octave_idx_type slen, octave_idx_type *adims, void *arep)
    : dimensions (adims),
      rep (reinterpret_cast<typename Array<T>::ArrayRep *> (arep)),
      slice_data (sdata), slice_len (slen) { }

public:

  //! Empty ctor (0 by 0).
  Array (void)
    : dimensions (), rep (nil_rep ()), slice_data (rep->data),
      slice_len (rep->len)
  {
    rep->count++;
  }

  //! nD uninitialized ctor.
  explicit Array (const dim_vector& dv)
    : dimensions (dv),
      rep (new typename Array<T>::ArrayRep (dv.safe_numel ())),
      slice_data (rep->data), slice_len (rep->len)
  {
    dimensions.chop_trailing_singletons ();
  }

  //! nD initialized ctor.
  explicit Array (const dim_vector& dv, const T& val)
    : dimensions (dv),
      rep (new typename Array<T>::ArrayRep (dv.safe_numel ())),
      slice_data (rep->data), slice_len (rep->len)
  {
    fill (val);
    dimensions.chop_trailing_singletons ();
  }

  //! Reshape constructor.
  Array (const Array<T>& a, const dim_vector& dv);

  //! Constructor from standard library sequence containers.
  template<template <typename...> class Container>
  Array (const Container<T>& a, const dim_vector& dv);

  //! Type conversion case.
  template <typename U>
  Array (const Array<U>& a)
    : dimensions (a.dims ()),
      rep (new typename Array<T>::ArrayRep (a.data (), a.numel ())),
      slice_data (rep->data), slice_len (rep->len)
  { }

  //! No type conversion case.
  Array (const Array<T>& a)
    : dimensions (a.dimensions), rep (a.rep), slice_data (a.slice_data),
      slice_len (a.slice_len)
  {
    rep->count++;
  }

public:

  virtual ~Array (void)
  {
    if (--rep->count == 0)
      delete rep;
  }

  Array<T>& operator = (const Array<T>& a)
  {
    if (this != &a)
      {
        if (--rep->count == 0)
          delete rep;

        rep = a.rep;
        rep->count++;

        dimensions = a.dimensions;
        slice_data = a.slice_data;
        slice_len = a.slice_len;
      }

    return *this;
  }

  void fill (const T& val);

  void clear (void);
  void clear (const dim_vector& dv);

  void clear (octave_idx_type r, octave_idx_type c)
  { clear (dim_vector (r, c)); }

  // Number of elements in the array.  These are all synonyms.
  //@{
  //! Number of elements in the array.
  //! Synonymous with numel().
  //! @note This method is deprecated in favour of numel().
  OCTAVE_DEPRECATED (4.4, "use 'numel' instead")
  octave_idx_type capacity (void) const { return numel (); }

  //! Number of elements in the array.
  //!
  //! Synonymous with numel().
  //! @note This method is deprecated in favour of numel().
  //!
  //! @note
  //! This is @em not the same as @c %length() at the Octave interpreter.
  //! At the Octave interpreter, the function @c %length() returns the
  //! length of the greatest dimension.  This method returns the total
  //! number of elements.

  OCTAVE_DEPRECATED (4.4, "use 'numel' instead")
  octave_idx_type length (void) const { return numel (); }

  //! Number of elements in the array.
  //! Synonymous with numel().
  //! @note This method is deprecated in favour of numel().
  OCTAVE_DEPRECATED (4.4, "use 'numel' instead")
  octave_idx_type nelem (void) const { return numel (); }

  //! Number of elements in the array.
  octave_idx_type numel (void) const { return slice_len; }
  //@}

  //! Return the array as a column vector.
  Array<T> as_column (void) const
  {
    Array<T> retval (*this);
    if (dimensions.ndims () != 2 || dimensions(1) != 1)
      retval.dimensions = dim_vector (numel (), 1);

    return retval;
  }

  //! Return the array as a row vector.
  Array<T> as_row (void) const
  {
    Array<T> retval (*this);
    if (dimensions.ndims () != 2 || dimensions(0) != 1)
      retval.dimensions = dim_vector (1, numel ());

    return retval;
  }

  //! Return the array as a matrix.
  Array<T> as_matrix (void) const
  {
    Array<T> retval (*this);
    if (dimensions.ndims () != 2)
      retval.dimensions = dimensions.redim (2);

    return retval;
  }

  //! @name First dimension
  //!
  //! Get the first dimension of the array (number of rows)
  //@{
  octave_idx_type dim1 (void) const { return dimensions(0); }
  octave_idx_type rows (void) const { return dim1 (); }
  //@}

  //! @name Second dimension
  //!
  //! Get the second dimension of the array (number of columns)
  //@{
  octave_idx_type dim2 (void) const { return dimensions(1); }
  octave_idx_type cols (void) const { return dim2 (); }
  octave_idx_type columns (void) const { return dim2 (); }
  //@}

  //! @name Third dimension
  //!
  //! Get the third dimension of the array (number of pages)
  //@{
  octave_idx_type dim3 (void) const { return dimensions(2); }
  octave_idx_type pages (void) const { return dim3 (); }
  //@}

  //! Size of the specified dimension.
  //!
  //! Dimensions beyond the Array number of dimensions return 1 as
  //! those are implicit singleton dimensions.
  //!
  //! Equivalent to Octave's `size (A, DIM)`

  size_type size (const size_type d) const
  {
    // Should we throw for negative values?
    // Should >= ndims () be handled by dim_vector operator() instead ?
    return d >= ndims () ? 1 : dimensions(d);
  }

  size_t byte_size (void) const
  { return static_cast<size_t> (numel ()) * sizeof (T); }

  //! Return a const-reference so that dims ()(i) works efficiently.
  const dim_vector& dims (void) const { return dimensions; }

  //! Chop off leading singleton dimensions
  Array<T> squeeze (void) const;

  octave_idx_type compute_index (octave_idx_type i, octave_idx_type j) const;
  octave_idx_type compute_index (octave_idx_type i, octave_idx_type j,
                                 octave_idx_type k) const;
  octave_idx_type compute_index (const Array<octave_idx_type>& ra_idx) const;

  octave_idx_type compute_index_unchecked (const Array<octave_idx_type>& ra_idx)
  const
  { return dimensions.compute_index (ra_idx.data (), ra_idx.numel ()); }

  // No checking, even for multiple references, ever.

  T& xelem (octave_idx_type n) { return slice_data[n]; }
  crefT xelem (octave_idx_type n) const { return slice_data[n]; }

  T& xelem (octave_idx_type i, octave_idx_type j)
  { return xelem (dim1 ()*j+i); }
  crefT xelem (octave_idx_type i, octave_idx_type j) const
  { return xelem (dim1 ()*j+i); }

  T& xelem (octave_idx_type i, octave_idx_type j, octave_idx_type k)
  { return xelem (i, dim2 ()*k+j); }
  crefT xelem (octave_idx_type i, octave_idx_type j, octave_idx_type k) const
  { return xelem (i, dim2 ()*k+j); }

  T& xelem (const Array<octave_idx_type>& ra_idx)
  { return xelem (compute_index_unchecked (ra_idx)); }

  crefT xelem (const Array<octave_idx_type>& ra_idx) const
  { return xelem (compute_index_unchecked (ra_idx)); }

  // FIXME: would be nice to fix this so that we don't unnecessarily force
  //        a copy, but that is not so easy, and I see no clean way to do it.

  T& checkelem (octave_idx_type n);

  T& checkelem (octave_idx_type i, octave_idx_type j);

  T& checkelem (octave_idx_type i, octave_idx_type j, octave_idx_type k);

  T& checkelem (const Array<octave_idx_type>& ra_idx);

  T& elem (octave_idx_type n)
  {
    make_unique ();
    return xelem (n);
  }

  T& elem (octave_idx_type i, octave_idx_type j) { return elem (dim1 ()*j+i); }

  T& elem (octave_idx_type i, octave_idx_type j, octave_idx_type k)
  { return elem (i, dim2 ()*k+j); }

  T& elem (const Array<octave_idx_type>& ra_idx)
  { return Array<T>::elem (compute_index_unchecked (ra_idx)); }

  T& operator () (octave_idx_type n) { return elem (n); }
  T& operator () (octave_idx_type i, octave_idx_type j) { return elem (i, j); }
  T& operator () (octave_idx_type i, octave_idx_type j, octave_idx_type k)
  { return elem (i, j, k); }
  T& operator () (const Array<octave_idx_type>& ra_idx)
  { return elem (ra_idx); }

  crefT checkelem (octave_idx_type n) const;

  crefT checkelem (octave_idx_type i, octave_idx_type j) const;

  crefT checkelem (octave_idx_type i, octave_idx_type j,
                   octave_idx_type k) const;

  crefT checkelem (const Array<octave_idx_type>& ra_idx) const;

  crefT elem (octave_idx_type n) const { return xelem (n); }

  crefT elem (octave_idx_type i, octave_idx_type j) const
  { return xelem (i, j); }

  crefT elem (octave_idx_type i, octave_idx_type j, octave_idx_type k) const
  { return xelem (i, j, k); }

  crefT elem (const Array<octave_idx_type>& ra_idx) const
  { return Array<T>::xelem (compute_index_unchecked (ra_idx)); }

  crefT operator () (octave_idx_type n) const { return elem (n); }
  crefT operator () (octave_idx_type i, octave_idx_type j) const
  { return elem (i, j); }
  crefT operator () (octave_idx_type i, octave_idx_type j,
                     octave_idx_type k) const
  { return elem (i, j, k); }
  crefT operator () (const Array<octave_idx_type>& ra_idx) const
  { return elem (ra_idx); }

  // Fast extractors.  All of these produce shallow copies.

  //! Extract column: A(:,k+1).
  Array<T> column (octave_idx_type k) const;
  //! Extract page: A(:,:,k+1).
  Array<T> page (octave_idx_type k) const;

  //! Extract a slice from this array as a column vector: A(:)(lo+1:up).
  //! Must be 0 <= lo && up <= numel.  May be up < lo.
  Array<T> linear_slice (octave_idx_type lo, octave_idx_type up) const;

  Array<T> reshape (octave_idx_type nr, octave_idx_type nc) const
  { return Array<T> (*this, dim_vector (nr, nc)); }

  Array<T> reshape (const dim_vector& new_dims) const
  { return Array<T> (*this, new_dims); }

  Array<T> permute (const Array<octave_idx_type>& vec, bool inv = false) const;
  Array<T> ipermute (const Array<octave_idx_type>& vec) const
  { return permute (vec, true); }

  bool issquare (void) const { return (dim1 () == dim2 ()); }

  OCTAVE_DEPRECATED (4.4, "use 'issquare' instead")
  bool is_square (void) const
  { return issquare (); }

  bool isempty (void) const { return numel () == 0; }

  OCTAVE_DEPRECATED (4.4, "use 'isempty' instead")
  bool is_empty (void) const
  { return isempty (); }

  bool isvector (void) const { return dimensions.isvector (); }

  OCTAVE_DEPRECATED (4.4, "use 'isvector' instead")
  bool is_vector (void) const
  { return isvector (); }

  bool is_nd_vector (void) const { return dimensions.is_nd_vector (); }

  Array<T> transpose (void) const;
  Array<T> hermitian (T (*fcn) (const T&) = nullptr) const;

  const T * data (void) const { return slice_data; }

  const T * fortran_vec (void) const { return data (); }

  T * fortran_vec (void);

  bool is_shared (void) { return rep->count > 1; }

  int ndims (void) const { return dimensions.ndims (); }

  //@{
  //! Indexing without resizing.
  Array<T> index (const idx_vector& i) const;

  Array<T> index (const idx_vector& i, const idx_vector& j) const;

  Array<T> index (const Array<idx_vector>& ia) const;
  //@}

  virtual T resize_fill_value (void) const;

  //@{
  //! Resizing (with fill).
  void resize2 (octave_idx_type nr, octave_idx_type nc, const T& rfv);
  void resize2 (octave_idx_type nr, octave_idx_type nc)
  {
    resize2 (nr, nc, resize_fill_value ());
  }

  void resize1 (octave_idx_type n, const T& rfv);
  void resize1 (octave_idx_type n) { resize1 (n, resize_fill_value ()); }

  void resize (const dim_vector& dv, const T& rfv);
  void resize (const dim_vector& dv) { resize (dv, resize_fill_value ()); }
  //@}

  //@{
  //! Indexing with possible resizing and fill

  // FIXME: this is really a corner case, that should better be
  // handled directly in liboctinterp.

  Array<T> index (const idx_vector& i, bool resize_ok, const T& rfv) const;
  Array<T> index (const idx_vector& i, bool resize_ok) const
  {
    return index (i, resize_ok, resize_fill_value ());
  }

  Array<T> index (const idx_vector& i, const idx_vector& j, bool resize_ok,
                  const T& rfv) const;
  Array<T> index (const idx_vector& i, const idx_vector& j,
                  bool resize_ok) const
  {
    return index (i, j, resize_ok, resize_fill_value ());
  }

  Array<T> index (const Array<idx_vector>& ia, bool resize_ok,
                  const T& rfv) const;
  Array<T> index (const Array<idx_vector>& ia, bool resize_ok) const
  {
    return index (ia, resize_ok, resize_fill_value ());
  }
  //@}

  //@{
  //! Indexed assignment (always with resize & fill).
  void assign (const idx_vector& i, const Array<T>& rhs, const T& rfv);
  void assign (const idx_vector& i, const Array<T>& rhs)
  {
    assign (i, rhs, resize_fill_value ());
  }

  void assign (const idx_vector& i, const idx_vector& j, const Array<T>& rhs,
               const T& rfv);
  void assign (const idx_vector& i, const idx_vector& j, const Array<T>& rhs)
  {
    assign (i, j, rhs, resize_fill_value ());
  }

  void assign (const Array<idx_vector>& ia, const Array<T>& rhs, const T& rfv);
  void assign (const Array<idx_vector>& ia, const Array<T>& rhs)
  {
    assign (ia, rhs, resize_fill_value ());
  }
  //@}

  //@{
  //! Deleting elements.

  //! A(I) = [] (with a single subscript)
  void delete_elements (const idx_vector& i);

  //! A(:,...,I,...,:) = [] (>= 2 subscripts, one of them is non-colon)
  void delete_elements (int dim, const idx_vector& i);

  //! Dispatcher to the above two.
  void delete_elements (const Array<idx_vector>& ia);
  //@}

  //! Insert an array into another at a specified position.  If
  //! size (a) is [d1 d2 ... dN] and idx is [i1 i2 ... iN], this
  //! method is equivalent to x(i1:i1+d1-1, i2:i2+d2-1, ... ,
  //! iN:iN+dN-1) = a.
  Array<T>& insert (const Array<T>& a, const Array<octave_idx_type>& idx);

  //! This is just a special case for idx = [r c 0 ...]
  Array<T>& insert (const Array<T>& a, octave_idx_type r, octave_idx_type c);

  void maybe_economize (void)
  {
    if (rep->count == 1 && slice_len != rep->len)
      {
        ArrayRep *new_rep = new ArrayRep (slice_data, slice_len);
        delete rep;
        rep = new_rep;
        slice_data = rep->data;
      }
  }

  void print_info (std::ostream& os, const std::string& prefix) const;

  //! Give a pointer to the data in mex format.  Unsafe.  This function
  //! exists to support the MEX interface.  You should not use it
  //! anywhere else.
  void * mex_get_data (void) const { return const_cast<T *> (data ()); }

  Array<T> sort (int dim = 0, sortmode mode = ASCENDING) const;
  Array<T> sort (Array<octave_idx_type> &sidx, int dim = 0,
                 sortmode mode = ASCENDING) const;

  //! Ordering is auto-detected or can be specified.
  sortmode issorted (sortmode mode = UNSORTED) const;

  OCTAVE_DEPRECATED (4.4, "use 'issorted' instead")
  sortmode is_sorted (sortmode mode = UNSORTED) const
  { return issorted (mode); }

  //! Sort by rows returns only indices.
  Array<octave_idx_type> sort_rows_idx (sortmode mode = ASCENDING) const;

  //! Ordering is auto-detected or can be specified.
  sortmode is_sorted_rows (sortmode mode = UNSORTED) const;

  //! Do a binary lookup in a sorted array.  Must not contain NaNs.
  //! Mode can be specified or is auto-detected by comparing 1st and last element.
  octave_idx_type lookup (const T& value, sortmode mode = UNSORTED) const;

  //! Ditto, but for an array of values, specializing on the case when values
  //! are sorted.  NaNs get the value N.
  Array<octave_idx_type> lookup (const Array<T>& values,
                                 sortmode mode = UNSORTED) const;

  //! Count nonzero elements.
  octave_idx_type nnz (void) const;

  //! Find indices of (at most n) nonzero elements.  If n is specified,
  //! backward specifies search from backward.
  Array<octave_idx_type> find (octave_idx_type n = -1,
                               bool backward = false) const;

  //! Returns the n-th element in increasing order, using the same
  //! ordering as used for sort.  n can either be a scalar index or a
  //! contiguous range.
  Array<T> nth_element (const idx_vector& n, int dim = 0) const;

  //! Get the kth super or subdiagonal.  The zeroth diagonal is the
  //! ordinary diagonal.
  Array<T> diag (octave_idx_type k = 0) const;

  Array<T> diag (octave_idx_type m, octave_idx_type n) const;

  //! Concatenation along a specified (0-based) dimension, equivalent
  //! to cat().  dim = -1 corresponds to dim = 0 and dim = -2
  //! corresponds to dim = 1, but apply the looser matching rules of
  //! vertcat/horzcat.
  static Array<T>
  cat (int dim, octave_idx_type n, const Array<T> *array_list);

  //! Apply function fcn to each element of the Array<T>.  This function
  //! is optimized with a manually unrolled loop.
  template <typename U, typename F>
  Array<U>
  map (F fcn) const
  {
    octave_idx_type len = numel ();

    const T *m = data ();

    Array<U> result (dims ());
    U *p = result.fortran_vec ();

    octave_idx_type i;
    for (i = 0; i < len - 3; i += 4)
      {
        octave_quit ();

        p[i] = fcn (m[i]);
        p[i+1] = fcn (m[i+1]);
        p[i+2] = fcn (m[i+2]);
        p[i+3] = fcn (m[i+3]);
      }

    octave_quit ();

    for (; i < len; i++)
      p[i] = fcn (m[i]);

    return result;
  }

  //@{
  //! Overloads for function references.
  template <typename U>
  Array<U>
  map (U (&fcn) (T)) const
  { return map<U, U (&) (T)> (fcn); }

  template <typename U>
  Array<U>
  map (U (&fcn) (const T&)) const
  { return map<U, U (&) (const T&)> (fcn); }
  //@}

  //! Generic any/all test functionality with arbitrary predicate.
  template <typename F, bool zero>
  bool test (F fcn) const
  {
    return any_all_test<F, T, zero> (fcn, data (), numel ());
  }

  //@{
  //! Simpler calls.
  template <typename F>
  bool test_any (F fcn) const
  { return test<F, false> (fcn); }

  template <typename F>
  bool test_all (F fcn) const
  { return test<F, true> (fcn); }
  //@}

  //@{
  //! Overloads for function references.
  bool test_any (bool (&fcn) (T)) const
  { return test<bool (&) (T), false> (fcn); }

  bool test_any (bool (&fcn) (const T&)) const
  { return test<bool (&) (const T&), false> (fcn); }

  bool test_all (bool (&fcn) (T)) const
  { return test<bool (&) (T), true> (fcn); }

  bool test_all (bool (&fcn) (const T&)) const
  { return test<bool (&) (const T&), true> (fcn); }
  //@}

  template <typename U> friend class Array;

  //! Returns true if this->dims () == dv, and if so, replaces this->dimensions
  //! by a shallow copy of dv.  This is useful for maintaining several arrays
  //! with supposedly equal dimensions (e.g. structs in the interpreter).
  bool optimize_dimensions (const dim_vector& dv);

  //@{
  //! WARNING: Only call these functions from jit

  int * jit_ref_count (void) { return rep->count.get (); }

  T * jit_slice_data (void) const { return slice_data; }

  octave_idx_type * jit_dimensions (void) const { return dimensions.to_jit (); }

  void * jit_array_rep (void) const { return rep; }
  //@}

private:
  static void instantiation_guard ();
};

// We use a variadic template for template template parameter so that
// we don't have to specify all the template parameters and limit this
// to Container<T>. http://stackoverflow.com/a/20499809/1609556
template<typename T>
template<template <typename...> class Container>
Array<T>::Array (const Container<T>& a, const dim_vector& dv)
  : dimensions (dv), rep (new typename Array<T>::ArrayRep (dv.safe_numel ())),
    slice_data (rep->data), slice_len (rep->len)
{
  if (dimensions.safe_numel () != octave_idx_type (a.size ()))
    {
      std::string new_dims_str = dimensions.str ();

      (*current_liboctave_error_handler)
        ("reshape: can't reshape %i elements into %s array",
         a.size (), new_dims_str.c_str ());
    }

  octave_idx_type i = 0;
  for (const T& x : a)
    slice_data[i++] = x;

  dimensions.chop_trailing_singletons ();
}

//! This is a simple wrapper template that will subclass an Array<T>
//! type or any later type derived from it and override the default
//! non-const operator() to not check for the array's uniqueness.  It
//! is, however, the user's responsibility to ensure the array is
//! actually unaliased whenever elements are accessed.
template <typename ArrayClass>
class NoAlias : public ArrayClass
{
  typedef typename ArrayClass::element_type T;
public:
  NoAlias () : ArrayClass () { }

  // FIXME: this would be simpler once C++0x is available
  template <typename X>
    explicit NoAlias (X x) : ArrayClass (x) { }

  template <typename X, typename Y>
    explicit NoAlias (X x, Y y) : ArrayClass (x, y) { }

  template <typename X, typename Y, typename Z>
    explicit NoAlias (X x, Y y, Z z) : ArrayClass (x, y, z) { }

  T& operator () (octave_idx_type n)
  { return ArrayClass::xelem (n); }
  T& operator () (octave_idx_type i, octave_idx_type j)
  { return ArrayClass::xelem (i, j); }
  T& operator () (octave_idx_type i, octave_idx_type j, octave_idx_type k)
  { return ArrayClass::xelem (i, j, k); }
  T& operator () (const Array<octave_idx_type>& ra_idx)
  { return ArrayClass::xelem (ra_idx); }
};

template <typename T>
std::ostream&
operator << (std::ostream& os, const Array<T>& a);

#endif