// Template array classes /* Copyright (C) 1993-2015 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 . */ #if !defined (octave_Array_h) #define octave_Array_h 1 #include #include #include #include #include "dim-vector.h" #include "idx-vector.h" #include "lo-traits.h" #include "lo-utils.h" #include "oct-sort.h" #include "quit.h" #include "oct-refcount.h" //!Handles the reference counting for all the derived classes. template class Array { protected: //! The real representation of all arrays. class ArrayRep { public: T *data; octave_idx_type len; octave_refcount count; ArrayRep (T *d, octave_idx_type l) : data (new T [l]), len (l), count (1) { std::copy (d, d+l, data); } template ArrayRep (U *d, octave_idx_type l) : data (new T [l]), len (l), count (1) { std::copy (d, d+l, data); } ArrayRep (void) : data (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 (a.data, a.data + a.len, data); } ~ArrayRep (void) { delete [] data; } octave_idx_type length (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 typename ref_param::type crefT; typedef bool (*compare_fcn_type) (typename ref_param::type, typename ref_param::type); protected: dim_vector dimensions; typename Array::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 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& 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: typename Array::ArrayRep *nil_rep (void) const { // NR was originally allocated with new, but that does not seem // to be necessary since it will never be deleted. So just use // a static object instead. static typename Array::ArrayRep nr; return &nr; } protected: //! For jit support Array (T *sdata, octave_idx_type slen, octave_idx_type *adims, void *arep) : dimensions (adims), rep (reinterpret_cast::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::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::ArrayRep (dv.safe_numel ())), slice_data (rep->data), slice_len (rep->len) { fill (val); dimensions.chop_trailing_singletons (); } //! Reshape constructor. Array (const Array& a, const dim_vector& dv); //! Type conversion case. template Array (const Array& a) : dimensions (a.dims ()), rep (new typename Array::ArrayRep (a.data (), a.length ())), slice_data (rep->data), slice_len (rep->len) { } //! No type conversion case. Array (const Array& 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& operator = (const Array& 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 length(), nelem(), and numel(). octave_idx_type capacity (void) const { return slice_len; } //! Number of elements in the array. /*! Synonymous with capacity(), nelem(), and 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_idx_type length (void) const { return capacity (); } //! Number of elements in the array. //! Synonymous with capacity(), length(), and numel(). octave_idx_type nelem (void) const { return capacity (); } //! Number of elements in the array. //! Synonymous with capacity(), length(), and nelem(). octave_idx_type numel (void) const { return nelem (); } //@} //! Return the array as a column vector. Array as_column (void) const { Array retval (*this); if (dimensions.length () != 2 || dimensions(1) != 1) retval.dimensions = dim_vector (numel (), 1); return retval; } //! Return the array as a row vector. Array as_row (void) const { Array retval (*this); if (dimensions.length () != 2 || dimensions(0) != 1) retval.dimensions = dim_vector (1, numel ()); return retval; } //! Return the array as a matrix. Array as_matrix (void) const { Array retval (*this); if (dimensions.length () != 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_t byte_size (void) const { return static_cast (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 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& ra_idx) const; octave_idx_type compute_index_unchecked (const Array& ra_idx) const { return dimensions.compute_index (ra_idx.data (), ra_idx.length ()); } // 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& ra_idx) { return xelem (compute_index_unchecked (ra_idx)); } crefT xelem (const Array& 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& 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& ra_idx) { return Array::elem (compute_index_unchecked (ra_idx)); } #if defined (BOUNDS_CHECKING) T& operator () (octave_idx_type n) { return checkelem (n); } T& operator () (octave_idx_type i, octave_idx_type j) { return checkelem (i, j); } T& operator () (octave_idx_type i, octave_idx_type j, octave_idx_type k) { return checkelem (i, j, k); } T& operator () (const Array& ra_idx) { return checkelem (ra_idx); } #else 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& ra_idx) { return elem (ra_idx); } #endif 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& 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& ra_idx) const { return Array::xelem (compute_index_unchecked (ra_idx)); } #if defined (BOUNDS_CHECKING) crefT operator () (octave_idx_type n) const { return checkelem (n); } crefT operator () (octave_idx_type i, octave_idx_type j) const { return checkelem (i, j); } crefT operator () (octave_idx_type i, octave_idx_type j, octave_idx_type k) const { return checkelem (i, j, k); } crefT operator () (const Array& ra_idx) const { return checkelem (ra_idx); } #else 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& ra_idx) const { return elem (ra_idx); } #endif // Fast extractors. All of these produce shallow copies. // Warning: none of these do check bounds, unless BOUNDS_CHECKING is on! //! Extract column: A(:,k+1). Array column (octave_idx_type k) const; //! Extract page: A(:,:,k+1). Array 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 linear_slice (octave_idx_type lo, octave_idx_type up) const; Array reshape (octave_idx_type nr, octave_idx_type nc) const { return Array (*this, dim_vector (nr, nc)); } Array reshape (const dim_vector& new_dims) const { return Array (*this, new_dims); } Array permute (const Array& vec, bool inv = false) const; Array ipermute (const Array& vec) const { return permute (vec, true); } bool is_square (void) const { return (dim1 () == dim2 ()); } bool is_empty (void) const { return numel () == 0; } bool is_vector (void) const { return dimensions.is_vector (); } Array transpose (void) const; Array hermitian (T (*fcn) (const T&) = 0) 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.length (); } //@{ //! Indexing without resizing. Array index (const idx_vector& i) const; Array index (const idx_vector& i, const idx_vector& j) const; Array index (const Array& 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 index (const idx_vector& i, bool resize_ok, const T& rfv) const; Array index (const idx_vector& i, bool resize_ok) const { return index (i, resize_ok, resize_fill_value ()); } Array index (const idx_vector& i, const idx_vector& j, bool resize_ok, const T& rfv) const; Array index (const idx_vector& i, const idx_vector& j, bool resize_ok) const { return index (i, j, resize_ok, resize_fill_value ()); } Array index (const Array& ia, bool resize_ok, const T& rfv) const; Array index (const Array& 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& rhs, const T& rfv); void assign (const idx_vector& i, const Array& rhs) { assign (i, rhs, resize_fill_value ()); } void assign (const idx_vector& i, const idx_vector& j, const Array& rhs, const T& rfv); void assign (const idx_vector& i, const idx_vector& j, const Array& rhs) { assign (i, j, rhs, resize_fill_value ()); } void assign (const Array& ia, const Array& rhs, const T& rfv); void assign (const Array& ia, const Array& 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& 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& insert (const Array& a, const Array& idx); //! This is just a special case for idx = [r c 0 ...] Array& insert (const Array& 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 (data ()); } Array sort (int dim = 0, sortmode mode = ASCENDING) const; Array sort (Array &sidx, int dim = 0, sortmode mode = ASCENDING) const; //! Ordering is auto-detected or can be specified. sortmode is_sorted (sortmode mode = UNSORTED) const; //! Sort by rows returns only indices. Array sort_rows_idx (sortmode mode = ASCENDING) const; //! Ordering is auto-detected or can be specified. sortmode is_sorted_rows (sortmode mode = UNSORTED) const; //! @brief 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 lookup (const Array& 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 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 nth_element (const idx_vector& n, int dim = 0) const; //! Get the kth super or subdiagonal. The zeroth diagonal is the //! ordinary diagonal. Array diag (octave_idx_type k = 0) const; Array 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 cat (int dim, octave_idx_type n, const Array *array_list); //! Apply function fcn to each element of the Array. This function //! is optimised with a manually unrolled loop. template Array map (F fcn) const { octave_idx_type len = length (); const T *m = data (); Array 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 Array map (U (&fcn) (T)) const { return map (fcn); } template Array map (U (&fcn) (const T&)) const { return map (fcn); } //@} //! Generic any/all test functionality with arbitrary predicate. template bool test (F fcn) const { return any_all_test (fcn, data (), length ()); } //@{ //! Simpler calls. template bool test_any (F fcn) const { return test (fcn); } template bool test_all (F fcn) const { return test (fcn); } //@} //@{ //! Overloads for function references. bool test_any (bool (&fcn) (T)) const { return test (fcn); } bool test_any (bool (&fcn) (const T&)) const { return test (fcn); } bool test_all (bool (&fcn) (T)) const { return test (fcn); } bool test_all (bool (&fcn) (const T&)) const { return test (fcn); } //@} template 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 (); }; //! This is a simple wrapper template that will subclass an Array //! 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 class NoAlias : public ArrayClass { typedef typename ArrayClass::element_type T; public: NoAlias () : ArrayClass () { } // FIXME: this would be simpler once C++0x is available template explicit NoAlias (X x) : ArrayClass (x) { } template explicit NoAlias (X x, Y y) : ArrayClass (x, y) { } template 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& ra_idx) { return ArrayClass::xelem (ra_idx); } }; template std::ostream& operator << (std::ostream& os, const Array& a); #endif