Mercurial > hg > octave-lojdl
changeset 8290:7cbe01c21986
improve dense array indexing
author | Jaroslav Hajek <highegg@gmail.com> |
---|---|
date | Mon, 20 Oct 2008 16:54:28 +0200 |
parents | ac7f334d9652 |
children | 53f35799b235 |
files | liboctave/Array-C.cc liboctave/Array-d.cc liboctave/Array-f.cc liboctave/Array-fC.cc liboctave/Array-i.cc liboctave/Array-s.cc liboctave/Array-util.cc liboctave/Array-util.h liboctave/Array.cc liboctave/Array.h liboctave/Array2.h liboctave/Array3.h liboctave/ArrayN.h liboctave/ChangeLog liboctave/Sparse.cc liboctave/chMatrix.h liboctave/dim-vector.h liboctave/idx-vector.cc liboctave/idx-vector.h scripts/polynomial/polyval.m src/Cell.cc src/Cell.h src/ChangeLog src/DLD-FUNCTIONS/dlmread.cc src/ov-base-mat.cc src/ov-cell.cc src/ov-cx-mat.cc src/ov-flt-cx-mat.cc src/ov-list.cc src/pr-output.cc test/ChangeLog test/test_logical-wfi-f.m test/test_logical-wfi-t.m test/test_slice.m |
diffstat | 34 files changed, 2240 insertions(+), 3483 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/Array-C.cc +++ b/liboctave/Array-C.cc @@ -96,10 +96,10 @@ INSTANTIATE_ARRAY_AND_ASSIGN (Complex, OCTAVE_API); -INSTANTIATE_ARRAY_ASSIGN (Complex, double, OCTAVE_API); -INSTANTIATE_ARRAY_ASSIGN (Complex, int, OCTAVE_API); -INSTANTIATE_ARRAY_ASSIGN (Complex, short, OCTAVE_API); -INSTANTIATE_ARRAY_ASSIGN (Complex, char, OCTAVE_API); +INSTANTIATE_ARRAY_ASSIGN (Complex, double, OCTAVE_API) +INSTANTIATE_ARRAY_ASSIGN (Complex, int, OCTAVE_API) +INSTANTIATE_ARRAY_ASSIGN (Complex, short, OCTAVE_API) +INSTANTIATE_ARRAY_ASSIGN (Complex, char, OCTAVE_API) #include "Array2.h"
--- a/liboctave/Array-d.cc +++ b/liboctave/Array-d.cc @@ -388,9 +388,9 @@ INSTANTIATE_ARRAY_AND_ASSIGN (double, OCTAVE_API); -INSTANTIATE_ARRAY_ASSIGN (double, int, OCTAVE_API); -INSTANTIATE_ARRAY_ASSIGN (double, short, OCTAVE_API); -INSTANTIATE_ARRAY_ASSIGN (double, char, OCTAVE_API); +INSTANTIATE_ARRAY_ASSIGN (double, int, OCTAVE_API) +INSTANTIATE_ARRAY_ASSIGN (double, short, OCTAVE_API) +INSTANTIATE_ARRAY_ASSIGN (double, char, OCTAVE_API) #include "Array2.h"
--- a/liboctave/Array-f.cc +++ b/liboctave/Array-f.cc @@ -388,9 +388,9 @@ INSTANTIATE_ARRAY_AND_ASSIGN (float, OCTAVE_API); -INSTANTIATE_ARRAY_ASSIGN (float, int, OCTAVE_API); -INSTANTIATE_ARRAY_ASSIGN (float, short, OCTAVE_API); -INSTANTIATE_ARRAY_ASSIGN (float, char, OCTAVE_API); +INSTANTIATE_ARRAY_ASSIGN (float, int, OCTAVE_API) +INSTANTIATE_ARRAY_ASSIGN (float, short, OCTAVE_API) +INSTANTIATE_ARRAY_ASSIGN (float, char, OCTAVE_API) #include "Array2.h"
--- a/liboctave/Array-fC.cc +++ b/liboctave/Array-fC.cc @@ -96,10 +96,10 @@ INSTANTIATE_ARRAY_AND_ASSIGN (FloatComplex, OCTAVE_API); -INSTANTIATE_ARRAY_ASSIGN (FloatComplex, float, OCTAVE_API); -INSTANTIATE_ARRAY_ASSIGN (FloatComplex, int, OCTAVE_API); -INSTANTIATE_ARRAY_ASSIGN (FloatComplex, short, OCTAVE_API); -INSTANTIATE_ARRAY_ASSIGN (FloatComplex, char, OCTAVE_API); +INSTANTIATE_ARRAY_ASSIGN (FloatComplex, float, OCTAVE_API) +INSTANTIATE_ARRAY_ASSIGN (FloatComplex, int, OCTAVE_API) +INSTANTIATE_ARRAY_ASSIGN (FloatComplex, short, OCTAVE_API) +INSTANTIATE_ARRAY_ASSIGN (FloatComplex, char, OCTAVE_API) #include "Array2.h"
--- a/liboctave/Array-i.cc +++ b/liboctave/Array-i.cc @@ -45,8 +45,8 @@ INSTANTIATE_ARRAY_AND_ASSIGN (long long, OCTAVE_API); #endif -INSTANTIATE_ARRAY_ASSIGN (int, short, OCTAVE_API); -INSTANTIATE_ARRAY_ASSIGN (int, char, OCTAVE_API); +INSTANTIATE_ARRAY_ASSIGN (int, short, OCTAVE_API) +INSTANTIATE_ARRAY_ASSIGN (int, char, OCTAVE_API) INSTANTIATE_ARRAY_SORT (octave_int8); INSTANTIATE_ARRAY_SORT (octave_int16);
--- a/liboctave/Array-s.cc +++ b/liboctave/Array-s.cc @@ -35,7 +35,7 @@ INSTANTIATE_ARRAY_AND_ASSIGN (short, OCTAVE_API); -INSTANTIATE_ARRAY_ASSIGN (short, char, OCTAVE_API); +INSTANTIATE_ARRAY_ASSIGN (short, char, OCTAVE_API) #include "Array2.h"
--- a/liboctave/Array-util.cc +++ b/liboctave/Array-util.cc @@ -322,34 +322,6 @@ } bool -is_in (octave_idx_type num, const idx_vector& idx) -{ - octave_idx_type n = idx.capacity (); - - for (octave_idx_type i = 0; i < n; i++) - if (idx.elem (i) == num) - return true; - - return false; -} - -octave_idx_type -how_many_lgt (const octave_idx_type num, idx_vector& idxv) -{ - octave_idx_type retval = 0; - - octave_idx_type n = idxv.capacity (); - - for (octave_idx_type i = 0; i < n; i++) - { - if (num > idxv.elem (i)) - retval++; - } - - return retval; -} - -bool all_ones (const Array<octave_idx_type>& arr) { bool retval = true; @@ -413,62 +385,90 @@ return retval; } -dim_vector -short_freeze (Array<idx_vector>& ra_idx, const dim_vector& dimensions, - int resize_ok) +dim_vector zero_dims_inquire (const Array<idx_vector>& ia, + const dim_vector& rhdv) { - dim_vector retval; - - int n = ra_idx.length (); - - int n_dims = dimensions.length (); - - if (n == n_dims) - { - retval = freeze (ra_idx, dimensions, resize_ok); - } - else if (n < n_dims) + int ial = ia.length (), rhdvl = rhdv.length (); + dim_vector rdv; + rdv.resize (ial); + bool *scalar = new bool[ial], *colon = new bool[ial]; + // Mark scalars and colons, count non-scalar indices. + int nonsc = 0; + bool all_colons = true; + for (int i = 0; i < ial; i++) { - retval.resize (n); - - for (int i = 0; i < n - 1; i++) - retval(i) = ra_idx(i).freeze (dimensions(i), "dimension", resize_ok); - - int size_left = 1; - - for (int i = n - 1; i < n_dims; i++) - size_left *= dimensions(i); - - if (ra_idx(n-1).is_colon()) - { - retval(n-1) = size_left; - } - else - { - octave_idx_type last_ra_idx = ra_idx(n-1)(0); - for (octave_idx_type i = 1; i < ra_idx (n - 1).capacity (); i++) - last_ra_idx = (last_ra_idx > ra_idx(n-1)(i) ? last_ra_idx : - ra_idx(n-1)(i)); - - if (last_ra_idx < size_left) - { - retval(n-1) = ra_idx(n-1).freeze (size_left, - "dimension", resize_ok); - } - else - { - // Make it larger than it should be to get an error - // later. - - retval.resize (n_dims+1); - - (*current_liboctave_error_handler) - ("index exceeds N-d array dimensions"); - } - } + // FIXME: should we check for length() instead? + scalar[i] = ia(i).is_scalar (); + colon[i] = ia(i).is_colon (); + if (! scalar[i]) nonsc++; + if (! colon[i]) rdv(i) = ia(i).extent (0); + all_colons = all_colons && colon[i]; } - return retval; + bool match = false; + // If the number of nonscalar indices matches the dimensionality of RHS, + // we try an exact match, inquiring even singleton dimensions. + if (all_colons) + { + rdv = rhdv; + rdv.resize(ial, 1); + } + else if (nonsc == rhdvl) + { + for (int i = 0, j = 0; i < ial; i++) + { + if (scalar[i]) continue; + if (colon[i]) + rdv(i) = rhdv(j++); + } + } + else + { + dim_vector rhdv0 = rhdv; + rhdv0.chop_all_singletons (); + int rhdv0l = rhdv0.length (); + for (int i = 0, j = 0; i < ial; i++) + { + if (scalar[i]) continue; + if (colon[i]) + rdv(i) = (j < rhdv0l) ? rhdv0(j++) : 1; + } + } + + delete [] scalar; + delete [] colon; + + return rdv; +} + +dim_vector zero_dims_inquire (const idx_vector& i, const idx_vector& j, + const dim_vector& rhdv) +{ + bool icol = i.is_colon (), jcol = j.is_colon (); + dim_vector rdv; + if (icol && jcol && rhdv.length () == 2) + { + rdv(0) = rhdv(0); + rdv(1) = rhdv(1); + } + else + { + dim_vector rhdv0 = rhdv; + rhdv0.chop_all_singletons (); + int k = 0; + rdv(0) = i.extent (0); + if (icol) + rdv(0) = rhdv0(k++); + else if (! i.is_scalar ()) + k++; + rdv(1) = j.extent (0); + if (jcol) + rdv(1) = rhdv0(k++); + else if (! j.is_scalar ()) + k++; + } + + return rdv; } int
--- a/liboctave/Array-util.h +++ b/liboctave/Array-util.h @@ -65,10 +65,6 @@ extern OCTAVE_API bool all_colon_equiv (const Array<idx_vector>& ra_idx, const dim_vector& frozen_lengths); -extern OCTAVE_API bool is_in (octave_idx_type num, const idx_vector& idx); - -extern OCTAVE_API octave_idx_type how_many_lgt (const octave_idx_type num, idx_vector& idxv); - extern OCTAVE_API bool all_ones (const Array<octave_idx_type>& arr); extern OCTAVE_API Array<octave_idx_type> get_elt_idx (const Array<idx_vector>& ra_idx, @@ -76,9 +72,11 @@ extern OCTAVE_API Array<octave_idx_type> get_ra_idx (octave_idx_type idx, const dim_vector& dims); -extern OCTAVE_API dim_vector short_freeze (Array<idx_vector>& ra_idx, - const dim_vector& dimensions, - int resize_ok); +extern OCTAVE_API dim_vector zero_dims_inquire (const Array<idx_vector>& ia, + const dim_vector& rhdv); + +extern OCTAVE_API dim_vector zero_dims_inquire (const idx_vector& i, const idx_vector& j, + const dim_vector& rhdv); struct permute_vector
--- a/liboctave/Array.cc +++ b/liboctave/Array.cc @@ -3,6 +3,7 @@ Copyright (C) 1993, 1994, 1995, 1996, 1997, 2000, 2002, 2003, 2004, 2005, 2006, 2007 John W. Eaton +Copyright (C) 2008 Jaroslav Hajek <highegg@gmail.com> This file is part of Octave. @@ -32,6 +33,7 @@ #include <iostream> #include <sstream> #include <vector> +#include <algorithm> #include <new> #include "Array.h" @@ -44,7 +46,7 @@ template <class T> Array<T>::Array (const Array<T>& a, const dim_vector& dv) - : rep (a.rep), dimensions (dv), idx (0), idx_count (0) + : rep (a.rep), dimensions (dv) { rep->count++; @@ -58,8 +60,6 @@ { if (--rep->count <= 0) delete rep; - - delete [] idx; } template <class T> @@ -75,10 +75,6 @@ rep->count++; dimensions = a.dimensions; - - delete [] idx; - idx_count = 0; - idx = 0; } return *this; @@ -559,458 +555,896 @@ return retval; } -template <class T> -void -Array<T>::resize_no_fill (octave_idx_type n) +// Helper class for multi-d index reduction and recursive indexing/indexed assignment. +// Rationale: we could avoid recursion using a state machine instead. However, using +// recursion is much more amenable to possible parallelization in the future. +// Also, the recursion solution is cleaner and more understandable. +class rec_index_helper { - if (n < 0) + octave_idx_type *dim, *cdim; + idx_vector *idx; + int top; + +public: + rec_index_helper (const dim_vector& dv, const Array<idx_vector>& ia) { - (*current_liboctave_error_handler) - ("can't resize to negative dimension"); - return; + int n = ia.length (); + assert (n > 0 && (dv.length () == std::max (n, 2))); + + dim = new octave_idx_type [2*n]; + // A hack to avoid double allocation + cdim = dim + n; + idx = new idx_vector [n]; + top = 0; + + dim[0] = dv(0); + cdim[0] = 1; + idx[0] = ia(0); + + for (int i = 1; i < n; i++) + { + // Try reduction... + if (idx[top].maybe_reduce (dim[top], ia(i), dv(i))) + { + // Reduction successful, fold dimensions. + dim[top] *= dv(i); + } + else + { + // Unsuccessful, store index & cumulative dim. + top++; + idx[top] = ia(i); + dim[top] = dv(i); + cdim[top] = cdim[top-1] * dim[top-1]; + } + } + } + + ~rec_index_helper (void) { delete [] idx; delete [] dim; } + +private: + + // Recursive N-d indexing + template <class T> + T *do_index (const T *src, T *dest, int lev) const + { + if (lev == 0) + dest += idx[0].index (src, dim[0], dest); + else + { + octave_idx_type n = idx[lev].length (dim[lev]), d = cdim[lev]; + for (octave_idx_type i = 0; i < n; i++) + dest = do_index (src + d*idx[lev].xelem (i), dest, lev-1); + } + + return dest; + } + + // Recursive N-d indexed assignment + template <class T> + const T *do_assign (const T *src, T *dest, int lev) const + { + if (lev == 0) + src += idx[0].assign (src, dim[0], dest); + else + { + octave_idx_type n = idx[lev].length (dim[lev]), d = cdim[lev]; + for (octave_idx_type i = 0; i < n; i++) + src = do_assign (src, dest + d*idx[lev].xelem (i), lev-1); + } + + return src; } - if (n == length ()) - return; - - typename Array<T>::ArrayRep *old_rep = rep; - const T *old_data = data (); - octave_idx_type old_len = length (); - - rep = new typename Array<T>::ArrayRep (n); - - dimensions = dim_vector (n); - - if (n > 0 && old_data && old_len > 0) + // Recursive N-d indexed assignment + template <class T> + void do_fill (const T& val, T *dest, int lev) const { - octave_idx_type min_len = old_len < n ? old_len : n; - - for (octave_idx_type i = 0; i < min_len; i++) - xelem (i) = old_data[i]; + if (lev == 0) + idx[0].fill (val, dim[0], dest); + else + { + octave_idx_type n = idx[lev].length (dim[lev]), d = cdim[lev]; + for (octave_idx_type i = 0; i < n; i++) + do_fill (val, dest + d*idx[lev].xelem (i), lev-1); + } } - if (--old_rep->count <= 0) - delete old_rep; +public: + + template <class T> + void index (const T *src, T *dest) const { do_index (src, dest, top); } + + template <class T> + void assign (const T *src, T *dest) const { do_assign (src, dest, top); } + + template <class T> + void fill (const T& val, T *dest) const { do_fill (val, dest, top); } + +}; + +// Helper class for multi-d recursive resizing +// This handles resize () in an efficient manner, touching memory only +// once (apart from reinitialization) +class rec_resize_helper +{ + octave_idx_type *cext, *sext, *dext; + int n; + +public: + rec_resize_helper (const dim_vector& ndv, const dim_vector& odv) + { + int l = ndv.length (); + assert (odv.length () == l); + octave_idx_type ld = 1; + int i = 0; + for (; i < l-1 && ndv(i) == odv(i); i++) ld *= ndv(i); + n = l - i; + cext = new octave_idx_type[3*n]; + // Trick to avoid three allocations + sext = cext + n; + dext = sext + n; + + octave_idx_type sld = ld, dld = ld; + for (int j = 0; j < n; j++) + { + cext[j] = std::min (ndv(i+j), odv(i+j)); + sext[j] = sld *= odv(i+j); + dext[j] = dld *= ndv(i+j); + } + cext[0] *= ld; + } + + ~rec_resize_helper (void) { delete [] cext; } + +private: + // recursive resizing + template <class T> + void do_resize_fill (const T* src, T *dest, const T& rfv, int lev) const + { + if (lev == 0) + { + T* destc = std::copy (src, src + cext[0], dest); + std::fill (destc, dest + dext[0], rfv); + } + else + { + octave_idx_type sd = sext[lev-1], dd = dext[lev-1], k; + for (k = 0; k < cext[lev]; k++) + do_resize_fill (src + k * sd, dest + k * dd, rfv, lev - 1); + + std::fill (dest + k * dd, dest + dext[lev], rfv); + } + } +public: + template <class T> + void resize_fill (const T* src, T *dest, const T& rfv) const + { do_resize_fill (src, dest, rfv, n-1); } + +}; + +static void gripe_index_out_of_range (void) +{ + (*current_liboctave_error_handler) + ("A(I): Index exceeds matrix dimension."); } template <class T> -void -Array<T>::resize_no_fill (const dim_vector& dv) +Array<T> +Array<T>::index (const idx_vector& i) const { - octave_idx_type n = dv.length (); + octave_idx_type n = numel (); + Array<T> retval; - for (octave_idx_type i = 0; i < n; i++) + if (i.is_colon ()) { - if (dv(i) < 0) - { - (*current_liboctave_error_handler) - ("can't resize to negative dimension"); - return; - } + // A(:) produces a shallow copy as a column vector. + retval.dimensions = dim_vector (n, 1); + rep->count++; + retval.rep = rep; } - - bool same_size = true; - - if (dimensions.length () != n) + else if (i.extent (n) != n) { - same_size = false; + gripe_index_out_of_range (); } else { - for (octave_idx_type i = 0; i < n; i++) - { - if (dv(i) != dimensions(i)) - { - same_size = false; - break; - } - } + // FIXME: This is the only place where orig_dimensions are used. + dim_vector rd = i.orig_dimensions (); + octave_idx_type il = i.length (n); + + // FIXME: This is for Matlab compatibility. Matlab 2007 given `b = ones(3,1)' + // yields the following: + // b(zeros(0,0)) gives [] + // b(zeros(1,0)) gives zeros(0,1) + // b(zeros(0,1)) gives zeros(0,1) + // b(zeros(0,m)) gives zeros(0,m) + // b(zeros(m,0)) gives zeros(m,0) + // b(1:2) gives ones(2,1) + // b(ones(2)) gives ones(2) etc. + // As you can see, the behaviour is weird, but the tests end up pretty + // simple. Nah, I don't want to suggest that this is ad hoc :) Neither do + // I want to say that Matlab is a lousy piece of s...oftware. + if (ndims () == 2 && n != 1) + { + if (columns () == 1 && rd(0) == 1) + rd = dim_vector (il, 1); + else if (rows () == 1 && rd(1) == 1) + rd = dim_vector (1, il); + } + + // Don't use resize here to avoid useless initialization for POD types. + retval = Array<T> (rd); + + if (il != 0) + i.index (data (), n, retval.fortran_vec ()); + } + + return retval; +} + +template <class T> +Array<T> +Array<T>::index (const idx_vector& i, const idx_vector& j) const +{ + // Get dimensions, allowing Fortran indexing in the 2nd dim. + dim_vector dv = dimensions.redim (2); + octave_idx_type r = dv(0), c = dv(1); + Array<T> retval; + + if (i.is_colon () && j.is_colon ()) + { + // A(:,:) produces a shallow copy. + retval = Array<T> (*this, dv); + } + else if (i.extent (r) != r || j.extent (c) != c) + { + gripe_index_out_of_range (); + } + else + { + octave_idx_type n = numel (), il = i.length (r), jl = j.length (c); + + // Don't use resize here to avoid useless initialization for POD types. + retval = Array<T> (dim_vector (il, jl)); + + idx_vector ii (i); + + const T* src = data (); + T *dest = retval.fortran_vec (); + + if (ii.maybe_reduce (r, j, c)) + ii.index (src, n, dest); + else + { + for (octave_idx_type k = 0; k < jl; k++) + dest += i.index (src + r * j.xelem (k), r, dest); + } } - if (same_size) - return; - - typename Array<T>::ArrayRep *old_rep = rep; - const T *old_data = data (); - - octave_idx_type ts = get_size (dv); + return retval; +} - rep = new typename Array<T>::ArrayRep (ts); +template <class T> +Array<T> +Array<T>::index (const Array<idx_vector>& ia) const +{ + int ial = ia.length (); + Array<T> retval; - dim_vector dv_old = dimensions; - octave_idx_type dv_old_orig_len = dv_old.length (); - dimensions = dv; - octave_idx_type ts_old = get_size (dv_old); - - if (ts > 0 && ts_old > 0 && dv_old_orig_len > 0) + // FIXME: is this dispatching necessary? + if (ial == 1) + retval = index (ia(0)); + else if (ial == 2) + retval = index (ia(0), ia(1)); + else if (ial > 0) { - Array<octave_idx_type> ra_idx (dimensions.length (), 0); + // Get dimensions, allowing Fortran indexing in the last dim. + dim_vector dv = dimensions.redim (ial); - if (n > dv_old_orig_len) - { - dv_old.resize (n); + // Check for out of bounds conditions. + bool all_colons = true, mismatch = false; + for (int i = 0; i < ial; i++) + { + if (ia(i).extent (dv(i)) != dv(i)) + { + mismatch = true; + break; + } + else + all_colons = all_colons && ia(i).is_colon (); + } + - for (octave_idx_type i = dv_old_orig_len; i < n; i++) - dv_old.elem (i) = 1; - } + if (mismatch) + { + gripe_index_out_of_range (); + } + else if (all_colons) + { + // A(:,:,...,:) produces a shallow copy. + retval = Array<T> (*this, dv); + } + else + { + // Form result dimensions. + dim_vector rdv; + rdv.resize (ial); + for (int i = 0; i < ial; i++) rdv(i) = ia(i).length (dv(i)); + rdv.chop_trailing_singletons (); - for (octave_idx_type i = 0; i < ts; i++) - { - if (index_in_bounds (ra_idx, dv_old)) - rep->elem (i) = old_data[get_scalar_idx (ra_idx, dv_old)]; + // Don't use resize here to avoid useless initialization for POD types. + retval = Array<T> (rdv); - increment_index (ra_idx, dimensions); - } + // Prepare for recursive indexing + rec_index_helper rh (dv, ia); + + // Do it. + rh.index (data (), retval.fortran_vec ()); + } } - if (--old_rep->count <= 0) - delete old_rep; + return retval; } +// FIXME: the following is a common error message to resize, regardless of whether it's +// called from assign or elsewhere. It seems OK to me, but eventually the gripe can be +// specialized. Anyway, propagating various error messages into procedure is, IMHO, a +// nonsense. If anything, we should change error handling here (and throughout liboctave) +// to allow custom handling of errors +static void gripe_invalid_resize (void) +{ + (*current_liboctave_error_handler) + ("resize: Invalid resizing operation or ambiguous assignment to an out-of-bounds array element."); +} + +// The default fill value. Override if you want a different one. + +template <class T> +T Array<T>::resize_fill_value () +{ + return T (); +} + +// Yes, we could do resize using index & assign. However, that would possibly +// involve a lot more memory traffic than we actually need. + template <class T> void -Array<T>::resize_no_fill (octave_idx_type r, octave_idx_type c) +Array<T>::resize_fill (octave_idx_type n, const T& rfv) { - if (r < 0 || c < 0) + if (n >= 0 && ndims () == 2) { - (*current_liboctave_error_handler) - ("can't resize to negative dimension"); - return; - } - - int n = ndims (); - - if (n == 0) - dimensions = dim_vector (0, 0); - - assert (ndims () == 2); - - if (r == dim1 () && c == dim2 ()) - return; - - typename Array<T>::ArrayRep *old_rep = Array<T>::rep; - const T *old_data = data (); - - octave_idx_type old_d1 = dim1 (); - octave_idx_type old_d2 = dim2 (); - octave_idx_type old_len = length (); - - octave_idx_type ts = get_size (r, c); - - rep = new typename Array<T>::ArrayRep (ts); - - dimensions = dim_vector (r, c); - - if (ts > 0 && old_data && old_len > 0) - { - octave_idx_type min_r = old_d1 < r ? old_d1 : r; - octave_idx_type min_c = old_d2 < c ? old_d2 : c; - - for (octave_idx_type j = 0; j < min_c; j++) - for (octave_idx_type i = 0; i < min_r; i++) - xelem (i, j) = old_data[old_d1*j+i]; - } - - if (--old_rep->count <= 0) - delete old_rep; -} + dim_vector dv; + // This is driven by Matlab's behaviour of giving a *row* vector on + // some out-of-bounds assignments. Specifically, matlab allows a(i) with + // out-of-bouds i when a is either of 0x0, 1x0, 1x1, 0xN, and gives + // a row vector in all cases (yes, even the last one, search me why). + // Giving a column vector would make much more sense (given the way + // trailing singleton dims are treated), but hey, Matlab is not here to + // make sense, it's here to make money ;) + bool invalid = false; + if (rows () == 0 || rows () == 1) + dv = dim_vector (1, n); + else if (columns () == 1) + dv = dim_vector (n, 1); + else + invalid = true; + + if (invalid) + gripe_invalid_resize (); + else + { + octave_idx_type nx = numel (); + if (n != nx) + { + Array<T> tmp = Array<T> (dv); + T *dest = tmp.fortran_vec (); -template <class T> -void -Array<T>::resize_no_fill (octave_idx_type r, octave_idx_type c, octave_idx_type p) -{ - if (r < 0 || c < 0 || p < 0) - { - (*current_liboctave_error_handler) - ("can't resize to negative dimension"); - return; - } - - int n = ndims (); - - if (n == 0) - dimensions = dim_vector (0, 0, 0); - - assert (ndims () == 3); - - if (r == dim1 () && c == dim2 () && p == dim3 ()) - return; - - typename Array<T>::ArrayRep *old_rep = rep; - const T *old_data = data (); + octave_idx_type n0 = std::min (n, nx), n1 = n - n0; + dest = std::copy (data (), data () + n0, dest); + std::fill (dest, dest + n1, rfv); - octave_idx_type old_d1 = dim1 (); - octave_idx_type old_d2 = dim2 (); - octave_idx_type old_d3 = dim3 (); - octave_idx_type old_len = length (); - - octave_idx_type ts = get_size (get_size (r, c), p); - - rep = new typename Array<T>::ArrayRep (ts); - - dimensions = dim_vector (r, c, p); - - if (ts > 0 && old_data && old_len > 0) - { - octave_idx_type min_r = old_d1 < r ? old_d1 : r; - octave_idx_type min_c = old_d2 < c ? old_d2 : c; - octave_idx_type min_p = old_d3 < p ? old_d3 : p; - - for (octave_idx_type k = 0; k < min_p; k++) - for (octave_idx_type j = 0; j < min_c; j++) - for (octave_idx_type i = 0; i < min_r; i++) - xelem (i, j, k) = old_data[old_d1*(old_d2*k+j)+i]; + *this = tmp; + } + } } - - if (--old_rep->count <= 0) - delete old_rep; + else + gripe_invalid_resize (); } template <class T> void -Array<T>::resize_and_fill (octave_idx_type n, const T& val) +Array<T>::resize_fill (octave_idx_type r, octave_idx_type c, const T& rfv) { - if (n < 0) + if (r >= 0 && c >= 0 && ndims () == 2) { - (*current_liboctave_error_handler) - ("can't resize to negative dimension"); - return; + octave_idx_type rx = rows (), cx = columns (); + if (r != rx || c != cx) + { + Array<T> tmp = Array<T> (dim_vector (r, c)); + T *dest = tmp.fortran_vec (); + + octave_idx_type r0 = std::min (r, rx), r1 = r - r0; + octave_idx_type c0 = std::min (c, cx), c1 = c - c0; + const T *src = data (); + if (r == rx) + dest = std::copy (src, src + r * c0, dest); + else + { + for (octave_idx_type k = 0; k < c0; k++) + { + dest = std::copy (src, src + r0, dest); + src += rx; + std::fill (dest, dest + r1, rfv); + dest += r1; + } + } + + std::fill (dest, dest + r * c1, rfv); + + *this = tmp; + } + } + else + gripe_invalid_resize (); + +} + +template<class T> +void +Array<T>::resize_fill (const dim_vector& dv, const T& rfv) +{ + int dvl = dv.length (); + if (dvl == 1) + resize (dv(0), rfv); + else if (dvl == 2) + resize (dv(0), dv(1), rfv); + else if (dimensions != dv) + { + if (dimensions.length () <= dvl) + { + Array<T> tmp (dv); + // Prepare for recursive resizing. + rec_resize_helper rh (dv, dimensions.redim (dvl)); + + // Do it. + rh.resize_fill (data (), tmp.fortran_vec (), rfv); + *this = tmp; + } + else + gripe_invalid_resize (); + } +} + +template <class T> +Array<T> +Array<T>::index (const idx_vector& i, bool resize_ok, const T& rfv) const +{ + Array<T> tmp = *this; + if (resize_ok) + { + octave_idx_type n = numel (), nx = i.extent (n); + if (n != nx) + tmp.resize_fill (nx, rfv); + + if (tmp.numel () != nx) + return Array<T> (); } - if (n == length ()) - return; - - typename Array<T>::ArrayRep *old_rep = rep; - const T *old_data = data (); - octave_idx_type old_len = length (); - - rep = new typename Array<T>::ArrayRep (n); - - dimensions = dim_vector (n); + return tmp.index (i); +} - if (n > 0) +template <class T> +Array<T> +Array<T>::index (const idx_vector& i, const idx_vector& j, + bool resize_ok, const T& rfv) const +{ + Array<T> tmp = *this; + if (resize_ok) { - octave_idx_type min_len = old_len < n ? old_len : n; + dim_vector dv = dimensions.redim (2); + octave_idx_type r = dv(0), c = dv(1); + octave_idx_type rx = i.extent (r), cx = j.extent (c); + if (r != rx || c != cx) + tmp.resize_fill (rx, cx, rfv); - if (old_data && old_len > 0) - { - for (octave_idx_type i = 0; i < min_len; i++) - xelem (i) = old_data[i]; - } - - for (octave_idx_type i = old_len; i < n; i++) - xelem (i) = val; + if (tmp.rows () != rx || tmp.columns () != cx) + return Array<T> (); } - if (--old_rep->count <= 0) - delete old_rep; + return tmp.index (i, j); +} + +template <class T> +Array<T> +Array<T>::index (const Array<idx_vector>& ia, + bool resize_ok, const T& rfv) const +{ + Array<T> tmp = *this; + if (resize_ok) + { + int ial = ia.length (); + dim_vector dv = dimensions.redim (ial); + dim_vector dvx; dvx.resize (ial); + for (int i = 0; i < ial; i++) dvx(i) = ia(i).extent (dv (i)); + if (! (dvx == dv)) + tmp.resize_fill (dvx, rfv); + + if (tmp.dimensions != dvx) + return Array<T> (); + } + + return tmp.index (ia); +} + + +static void +gripe_invalid_assignment_size (void) +{ + (*current_liboctave_error_handler) + ("A(I) = X: X must have the same size as I"); +} + +static void +gripe_assignment_dimension_mismatch (void) +{ + (*current_liboctave_error_handler) + ("A(I,J,...) = X: dimensions mismatch"); } template <class T> void -Array<T>::resize_and_fill (octave_idx_type r, octave_idx_type c, const T& val) +Array<T>::assign (const idx_vector& i, const Array<T>& rhs, const T& rfv) { - if (r < 0 || c < 0) - { - (*current_liboctave_error_handler) - ("can't resize to negative dimension"); - return; - } - - if (ndims () == 0) - dimensions = dim_vector (0, 0); + octave_idx_type n = numel (), rhl = rhs.numel (); - assert (ndims () == 2); - - if (r == dim1 () && c == dim2 ()) - return; - - typename Array<T>::ArrayRep *old_rep = Array<T>::rep; - const T *old_data = data (); - - octave_idx_type old_d1 = dim1 (); - octave_idx_type old_d2 = dim2 (); - octave_idx_type old_len = length (); - - octave_idx_type ts = get_size (r, c); + if (rhl == 1 || i.length (n) == rhl) + { + octave_idx_type nx = i.extent (n); + // Try to resize first if necessary. + if (nx != n) + { + resize_fill (nx, rfv); + n = numel (); + } - rep = new typename Array<T>::ArrayRep (ts); - - dimensions = dim_vector (r, c); - - if (ts > 0) - { - octave_idx_type min_r = old_d1 < r ? old_d1 : r; - octave_idx_type min_c = old_d2 < c ? old_d2 : c; - - if (old_data && old_len > 0) - { - for (octave_idx_type j = 0; j < min_c; j++) - for (octave_idx_type i = 0; i < min_r; i++) - xelem (i, j) = old_data[old_d1*j+i]; - } - - for (octave_idx_type j = 0; j < min_c; j++) - for (octave_idx_type i = min_r; i < r; i++) - xelem (i, j) = val; - - for (octave_idx_type j = min_c; j < c; j++) - for (octave_idx_type i = 0; i < r; i++) - xelem (i, j) = val; + // If the resizing was ambiguous, resize has already griped. + if (nx == n) + { + if (i.is_colon ()) + { + // A(:) = X makes a full fill or a shallow copy. + if (rhl == 1) + fill (rhs(0)); + else + *this = rhs.reshape (dimensions); + } + else + { + if (rhl == 1) + i.fill (rhs(0), n, fortran_vec ()); + else + i.assign (rhs.data (), n, fortran_vec ()); + } + } } - - if (--old_rep->count <= 0) - delete old_rep; + else + gripe_invalid_assignment_size (); } template <class T> void -Array<T>::resize_and_fill (octave_idx_type r, octave_idx_type c, octave_idx_type p, const T& val) +Array<T>::assign (const idx_vector& i, const idx_vector& j, + const Array<T>& rhs, const T& rfv) { - if (r < 0 || c < 0 || p < 0) + // Get RHS extents, discarding singletons. + dim_vector rhdv = rhs.dims (); + // Get LHS extents, allowing Fortran indexing in the second dim. + dim_vector dv = dimensions.redim (2); + // Check for out-of-bounds and form resizing dimensions. + dim_vector rdv; + // In the special when all dimensions are zero, colons are allowed to inquire + // the shape of RHS. The rules are more obscure, so we solve that elsewhere. + if (dv.all_zero ()) + rdv = zero_dims_inquire (i, j, rhdv); + else { - (*current_liboctave_error_handler) - ("can't resize to negative dimension"); - return; + rdv(0) = i.extent (dv(0)); + rdv(1) = j.extent (dv(1)); } - if (ndims () == 0) - dimensions = dim_vector (0, 0, 0); - - assert (ndims () == 3); + bool isfill = rhs.numel () == 1; + octave_idx_type il = i.length (rdv(0)), jl = j.length (rdv(1)); + rhdv.chop_all_singletons (); + bool match = isfill || (rhdv.length () == 2 && il == rhdv(0) && jl == rhdv(1)); + match = match || (il == 1 && jl == rhdv(0) && rhdv(1) == 1); - if (r == dim1 () && c == dim2 () && p == dim3 ()) - return; - - typename Array<T>::ArrayRep *old_rep = rep; - const T *old_data = data (); + if (match) + { + // Resize if requested. + if (rdv != dv) + { + resize (rdv, rfv); + dv = dimensions; + } - octave_idx_type old_d1 = dim1 (); - octave_idx_type old_d2 = dim2 (); - octave_idx_type old_d3 = dim3 (); - - octave_idx_type old_len = length (); - - octave_idx_type ts = get_size (get_size (r, c), p); - - rep = new typename Array<T>::ArrayRep (ts); - - dimensions = dim_vector (r, c, p); - - if (ts > 0) - { - octave_idx_type min_r = old_d1 < r ? old_d1 : r; - octave_idx_type min_c = old_d2 < c ? old_d2 : c; - octave_idx_type min_p = old_d3 < p ? old_d3 : p; + // If the resizing was invalid, resize has already griped. + if (dv == rdv) + { + if (i.is_colon () && j.is_colon ()) + { + // A(:,:) = X makes a full fill or a shallow copy + if (isfill) + fill (rhs(0)); + else + *this = rhs.reshape (dimensions); + } + else + { + // The actual work. + octave_idx_type n = numel (), r = rows (), c = columns (); + idx_vector ii (i); - if (old_data && old_len > 0) - for (octave_idx_type k = 0; k < min_p; k++) - for (octave_idx_type j = 0; j < min_c; j++) - for (octave_idx_type i = 0; i < min_r; i++) - xelem (i, j, k) = old_data[old_d1*(old_d2*k+j)+i]; - - // FIXME -- if the copy constructor is expensive, this - // may win. Otherwise, it may make more sense to just copy the - // value everywhere when making the new ArrayRep. + const T* src = rhs.data (); + T *dest = fortran_vec (); - for (octave_idx_type k = 0; k < min_p; k++) - for (octave_idx_type j = min_c; j < c; j++) - for (octave_idx_type i = 0; i < min_r; i++) - xelem (i, j, k) = val; + // Try reduction first. + if (ii.maybe_reduce (r, j, c)) + { + if (isfill) + ii.fill (*src, n, dest); + else + ii.assign (src, n, dest); + } + else + { + if (isfill) + { + for (octave_idx_type k = 0; k < jl; k++) + i.fill (*src, r, dest + r * j.xelem (k)); + } + else + { + for (octave_idx_type k = 0; k < jl; k++) + src += i.assign (src, r, dest + r * j.xelem (k)); + } + } + } + + } - for (octave_idx_type k = 0; k < min_p; k++) - for (octave_idx_type j = 0; j < c; j++) - for (octave_idx_type i = min_r; i < r; i++) - xelem (i, j, k) = val; - - for (octave_idx_type k = min_p; k < p; k++) - for (octave_idx_type j = 0; j < c; j++) - for (octave_idx_type i = 0; i < r; i++) - xelem (i, j, k) = val; } - - if (--old_rep->count <= 0) - delete old_rep; + else + gripe_assignment_dimension_mismatch (); } template <class T> void -Array<T>::resize_and_fill (const dim_vector& dv, const T& val) +Array<T>::assign (const Array<idx_vector>& ia, + const Array<T>& rhs, const T& rfv) { - octave_idx_type n = dv.length (); + int ial = ia.length (); - for (octave_idx_type i = 0; i < n; i++) + // FIXME: is this dispatching necessary / desirable? + if (ial == 1) + assign (ia(0), rhs, rfv); + else if (ial == 2) + assign (ia(0), ia(1), rhs, rfv); + else if (ial > 0) { - if (dv(i) < 0) - { - (*current_liboctave_error_handler) - ("can't resize to negative dimension"); - return; - } - } + // Get RHS extents, discarding singletons. + dim_vector rhdv = rhs.dims (); + + // Get LHS extents, allowing Fortran indexing in the second dim. + dim_vector dv = dimensions.redim (ial); + + // Get the extents forced by indexing. + dim_vector rdv; + + // In the special when all dimensions are zero, colons are allowed to + // inquire the shape of RHS. The rules are more obscure, so we solve that + // elsewhere. + if (dv.all_zero ()) + rdv = zero_dims_inquire (ia, rhdv); + else + { + rdv.resize (ial); + for (int i = 0; i < ial; i++) + rdv(i) = ia(i).extent (dv(i)); + } + + // Check whether LHS and RHS match, up to singleton dims. + bool match = true, all_colons = true, isfill = rhs.numel () == 1; + + rhdv.chop_all_singletons (); + int j = 0, rhdvl = rhdv.length (); + for (int i = 0; i < ial; i++) + { + all_colons = all_colons && ia(i).is_colon (); + octave_idx_type l = ia(i).length (rdv(i)); + if (l == 1) continue; + match = match && j < rhdvl && l == rhdv(j++); + } + + match = match && (j == rhdvl || rhdv(j) == 1); + match = match || isfill; + + if (match) + { + // Resize first if necessary. + if (rdv != dv) + { + resize_fill (rdv, rfv); + dv = dimensions; + chop_trailing_singletons (); + } - bool same_size = true; + // If the resizing was invalid, resize has already griped. + if (dv == rdv) + { + if (all_colons) + { + // A(:,:,...,:) = X makes a full fill or a shallow copy. + if (isfill) + fill (rhs(0)); + else + *this = rhs.reshape (dimensions); + } + else + { + // Do the actual work. + + // Prepare for recursive indexing + rec_index_helper rh (dv, ia); - if (dimensions.length () != n) - { - same_size = false; + // Do it. + if (isfill) + rh.fill (rhs(0), fortran_vec ()); + else + rh.assign (rhs.data (), fortran_vec ()); + } + } + } + else + gripe_assignment_dimension_mismatch (); } - else +} + +template <class T> +void +Array<T>::delete_elements (const idx_vector& i) +{ + octave_idx_type n = numel (); + if (i.is_colon ()) + { + *this = Array<T> (); + } + else if (i.extent (n) != n) { - for (octave_idx_type i = 0; i < n; i++) - { - if (dv(i) != dimensions(i)) - { - same_size = false; - break; - } - } + gripe_index_out_of_range (); + } + else if (i.length (n) != 0) + { + octave_idx_type l, u; + + bool col_vec = ndims () == 2 && columns () == 1 && rows () != 1; + if (i.is_cont_range (n, l, u)) + { + // Special case deleting a contiguous range. + octave_idx_type m = n + l - u; + Array<T> tmp (dim_vector (col_vec ? m : 1, !col_vec ? m : 1)); + const T *src = data (); + T *dest = tmp.fortran_vec (); + dest = std::copy (src, src + l, dest); + dest = std::copy (src + u, src + n, dest); + *this = tmp; + } + else + { + // Use index. + *this = index (i.complement (n)); + } + } +} + +template <class T> +void +Array<T>::delete_elements (int dim, const idx_vector& i) +{ + if (dim > ndims ()) + { + (*current_liboctave_error_handler) + ("invalid dimension in delete_elements"); + return; } - if (same_size) - return; - - typename Array<T>::ArrayRep *old_rep = rep; - const T *old_data = data (); - - octave_idx_type len = get_size (dv); - - rep = new typename Array<T>::ArrayRep (len); + octave_idx_type n = dimensions (dim); + if (i.is_colon ()) + { + *this = Array<T> (); + } + else if (i.extent (n) != n) + { + gripe_index_out_of_range (); + } + else if (i.length (n) != 0) + { + octave_idx_type l, u; - dim_vector dv_old = dimensions; - octave_idx_type dv_old_orig_len = dv_old.length (); - dimensions = dv; + if (i.is_cont_range (n, l, u)) + { + // Special case deleting a contiguous range. + octave_idx_type nd = n + l - u, dl = 1, du = 1; + dim_vector rdv = dimensions; + rdv(dim) = nd; + for (int k = 0; k < dim; k++) dl *= dimensions(k); + for (int k = dim + 1; k < ndims (); k++) du *= dimensions(k); - if (len > 0 && dv_old_orig_len > 0) - { - Array<octave_idx_type> ra_idx (dimensions.length (), 0); - - if (n > dv_old_orig_len) - { - dv_old.resize (n); + // Special case deleting a contiguous range. + Array<T> tmp = Array<T> (rdv); + const T *src = data (); + T *dest = tmp.fortran_vec (); + l *= dl; u *= dl; n *= dl; + for (octave_idx_type k = 0; k < du; k++) + { + dest = std::copy (src, src + l, dest); + dest = std::copy (src + u, src + n, dest); + src += n; + } - for (octave_idx_type i = dv_old_orig_len; i < n; i++) - dv_old.elem (i) = 1; - } + *this = tmp; + } + else + { + // Use index. + Array<idx_vector> ia (ndims (), idx_vector::colon); + ia (dim) = i.complement (n); + *this = index (ia); + } + } +} - for (octave_idx_type i = 0; i < len; i++) - { - if (index_in_bounds (ra_idx, dv_old)) - rep->elem (i) = old_data[get_scalar_idx (ra_idx, dv_old)]; - else - rep->elem (i) = val; - - increment_index (ra_idx, dimensions); - } +template <class T> +void +Array<T>::delete_elements (const Array<idx_vector>& ia) +{ + if (ia.length () == 1) + delete_elements (ia(0)); + else + { + int len = ia.length (), k, dim = -1; + for (k = 0; k < len; k++) + { + if (! ia(k).is_colon ()) + { + if (dim < 0) + dim = k; + else + break; + } + } + if (dim < 0) + { + dim_vector dv = dimensions; + dv(0) = 0; + *this = Array<T> (dv); + } + else if (k == len) + { + delete_elements (dim, ia(dim)); + } + else + { + (*current_liboctave_error_handler) + ("A null assignment can only have one non-colon index."); + } } - else - for (octave_idx_type i = 0; i < len; i++) - rep->elem (i) = val; - if (--old_rep->count <= 0) - delete old_rep; } +// FIXME: Remove these methods or implement them using assign. + template <class T> Array<T>& Array<T>::insert (const Array<T>& a, octave_idx_type r, octave_idx_type c) @@ -1194,6 +1628,7 @@ return *this; } + template <class T> Array<T> Array<T>::transpose (void) const @@ -1406,1178 +1841,6 @@ } template <class T> -void -Array<T>::clear_index (void) const -{ - delete [] idx; - idx = 0; - idx_count = 0; -} - -template <class T> -void -Array<T>::set_index (const idx_vector& idx_arg) const -{ - int nd = ndims (); - - if (! idx && nd > 0) - idx = new idx_vector [nd]; - - if (idx_count < nd) - { - idx[idx_count++] = idx_arg; - } - else - { - idx_vector *new_idx = new idx_vector [idx_count+1]; - - for (int i = 0; i < idx_count; i++) - new_idx[i] = idx[i]; - - new_idx[idx_count++] = idx_arg; - - delete [] idx; - - idx = new_idx; - } -} - -template <class T> -void -Array<T>::maybe_delete_elements (idx_vector& idx_arg) -{ - switch (ndims ()) - { - case 1: - maybe_delete_elements_1 (idx_arg); - break; - - case 2: - maybe_delete_elements_2 (idx_arg); - break; - - default: - (*current_liboctave_error_handler) - ("Array<T>::maybe_delete_elements: invalid operation"); - break; - } -} - -template <class T> -void -Array<T>::maybe_delete_elements_1 (idx_vector& idx_arg) -{ - octave_idx_type len = length (); - - if (len == 0) - return; - - if (idx_arg.is_colon_equiv (len, 1)) - resize_no_fill (0); - else - { - int num_to_delete = idx_arg.length (len); - - if (num_to_delete != 0) - { - octave_idx_type new_len = len; - - octave_idx_type iidx = 0; - - for (octave_idx_type i = 0; i < len; i++) - if (i == idx_arg.elem (iidx)) - { - iidx++; - new_len--; - - if (iidx == num_to_delete) - break; - } - - if (new_len > 0) - { - T *new_data = new T [new_len]; - - octave_idx_type ii = 0; - iidx = 0; - for (octave_idx_type i = 0; i < len; i++) - { - if (iidx < num_to_delete && i == idx_arg.elem (iidx)) - iidx++; - else - { - new_data[ii] = xelem (i); - ii++; - } - } - - if (--rep->count <= 0) - delete rep; - - rep = new typename Array<T>::ArrayRep (new_data, new_len); - - dimensions.resize (1); - dimensions(0) = new_len; - } - else - (*current_liboctave_error_handler) - ("A(idx) = []: index out of range"); - } - } -} - -template <class T> -void -Array<T>::maybe_delete_elements_2 (idx_vector& idx_arg) -{ - assert (ndims () == 2); - - octave_idx_type nr = dim1 (); - octave_idx_type nc = dim2 (); - - if (idx_arg.is_colon ()) - { - // A(:) = [] always gives 0-by-0 matrix, even if A was empty. - resize_no_fill (0, 0); - return; - } - - octave_idx_type n; - if (nr == 1) - n = nc; - else if (nc == 1) - n = nr; - else if (! idx_arg.orig_empty ()) - { - // Reshape to row vector for Matlab compatibility. - - n = nr * nc; - nr = 1; - nc = n; - } - else - return; - - idx_arg.sort (true); - - if (idx_arg.is_colon_equiv (n, 1)) - { - if (nr == 1) - resize_no_fill (1, 0); - else if (nc == 1) - resize_no_fill (0, 1); - return; - } - - octave_idx_type num_to_delete = idx_arg.length (n); - - if (num_to_delete != 0) - { - octave_idx_type new_n = n; - - octave_idx_type iidx = 0; - - for (octave_idx_type i = 0; i < n; i++) - if (i == idx_arg.elem (iidx)) - { - iidx++; - new_n--; - - if (iidx == num_to_delete) - break; - } - - if (new_n > 0) - { - T *new_data = new T [new_n]; - - octave_idx_type ii = 0; - iidx = 0; - for (octave_idx_type i = 0; i < n; i++) - { - if (iidx < num_to_delete && i == idx_arg.elem (iidx)) - iidx++; - else - { - new_data[ii] = xelem (i); - - ii++; - } - } - - if (--(Array<T>::rep)->count <= 0) - delete Array<T>::rep; - - Array<T>::rep = new typename Array<T>::ArrayRep (new_data, new_n); - - dimensions.resize (2); - - if (nr == 1) - { - dimensions(0) = 1; - dimensions(1) = new_n; - } - else - { - dimensions(0) = new_n; - dimensions(1) = 1; - } - } - else - (*current_liboctave_error_handler) - ("A(idx) = []: index out of range"); - } -} - -template <class T> -void -Array<T>::maybe_delete_elements (idx_vector& idx_i, idx_vector& idx_j) -{ - assert (ndims () == 2); - - octave_idx_type nr = dim1 (); - octave_idx_type nc = dim2 (); - - if (idx_i.is_colon () && idx_j.is_colon ()) - { - // A special case: A(:,:). Matlab gives 0-by-nc here, but perhaps we - // should not? - resize_no_fill (0, nc); - } - else if (idx_i.is_colon ()) - { - idx_j.sort (true); // sort in advance to speed-up the following check - - if (idx_j.is_colon_equiv (nc, 1)) - resize_no_fill (nr, 0); - else - { - octave_idx_type num_to_delete = idx_j.length (nc); - - if (num_to_delete != 0) - { - octave_idx_type new_nc = nc; - - octave_idx_type iidx = 0; - - for (octave_idx_type j = 0; j < nc; j++) - if (j == idx_j.elem (iidx)) - { - iidx++; - new_nc--; - - if (iidx == num_to_delete) - break; - } - - if (new_nc > 0) - { - T *new_data = new T [nr * new_nc]; - - octave_idx_type jj = 0; - iidx = 0; - for (octave_idx_type j = 0; j < nc; j++) - { - if (iidx < num_to_delete && j == idx_j.elem (iidx)) - iidx++; - else - { - for (octave_idx_type i = 0; i < nr; i++) - new_data[nr*jj+i] = xelem (i, j); - jj++; - } - } - - if (--(Array<T>::rep)->count <= 0) - delete Array<T>::rep; - - Array<T>::rep = new typename Array<T>::ArrayRep (new_data, nr * new_nc); - - dimensions.resize (2); - dimensions(1) = new_nc; - } - else - (*current_liboctave_error_handler) - ("A(idx) = []: index out of range"); - } - } - } - else if (idx_j.is_colon ()) - { - idx_i.sort (true); // sort in advance to speed-up the following check - - if (idx_i.is_colon_equiv (nr, 1)) - resize_no_fill (0, nc); - else - { - octave_idx_type num_to_delete = idx_i.length (nr); - - if (num_to_delete != 0) - { - octave_idx_type new_nr = nr; - - octave_idx_type iidx = 0; - - for (octave_idx_type i = 0; i < nr; i++) - if (i == idx_i.elem (iidx)) - { - iidx++; - new_nr--; - - if (iidx == num_to_delete) - break; - } - - if (new_nr > 0) - { - T *new_data = new T [new_nr * nc]; - - octave_idx_type ii = 0; - iidx = 0; - for (octave_idx_type i = 0; i < nr; i++) - { - if (iidx < num_to_delete && i == idx_i.elem (iidx)) - iidx++; - else - { - for (octave_idx_type j = 0; j < nc; j++) - new_data[new_nr*j+ii] = xelem (i, j); - ii++; - } - } - - if (--(Array<T>::rep)->count <= 0) - delete Array<T>::rep; - - Array<T>::rep = new typename Array<T>::ArrayRep (new_data, new_nr * nc); - - dimensions.resize (2); - dimensions(0) = new_nr; - } - else - (*current_liboctave_error_handler) - ("A(idx) = []: index out of range"); - } - } - } - else if (! (idx_i.orig_empty () || idx_j.orig_empty ())) - { - (*current_liboctave_error_handler) - ("a null assignment can have only one non-colon index"); - } -} - -template <class T> -void -Array<T>::maybe_delete_elements (idx_vector&, idx_vector&, idx_vector&) -{ - assert (0); -} - -template <class T> -void -Array<T>::maybe_delete_elements (Array<idx_vector>& ra_idx) -{ - octave_idx_type n_idx = ra_idx.length (); - - // Special case matrices - if (ndims () == 2) - { - if (n_idx == 1) - { - maybe_delete_elements (ra_idx (0)); - return; - } - else if (n_idx == 2) - { - maybe_delete_elements (ra_idx (0), ra_idx (1)); - return; - } - } - - dim_vector lhs_dims = dims (); - - int n_lhs_dims = lhs_dims.length (); - - if (lhs_dims.all_zero ()) - return; - - if (n_idx == 1 && ra_idx(0).is_colon ()) - { - resize (dim_vector (0, 0)); - return; - } - - if (n_idx > n_lhs_dims) - { - for (int i = n_idx; i < n_lhs_dims; i++) - { - // Ensure that extra indices are either colon or 1. - - if (! ra_idx(i).is_colon_equiv (1, 1)) - { - (*current_liboctave_error_handler) - ("index exceeds array dimensions"); - return; - } - } - - ra_idx.resize (n_lhs_dims); - - n_idx = n_lhs_dims; - } - - Array<int> idx_is_colon (n_idx, 0); - - Array<int> idx_is_colon_equiv (n_idx, 0); - - // Initialization of colon arrays. - - for (octave_idx_type i = 0; i < n_idx; i++) - { - if (ra_idx(i).orig_empty ()) return; - idx_is_colon_equiv(i) = ra_idx(i).is_colon_equiv (lhs_dims(i), 1); - - idx_is_colon(i) = ra_idx(i).is_colon (); - } - - bool idx_ok = true; - - // Check for index out of bounds. - - for (octave_idx_type i = 0 ; i < n_idx - 1; i++) - { - if (! (idx_is_colon(i) || idx_is_colon_equiv(i))) - { - ra_idx(i).sort (true); - - if (ra_idx(i).max () > lhs_dims(i)) - { - (*current_liboctave_error_handler) - ("index exceeds array dimensions"); - - idx_ok = false; - break; - } - else if (ra_idx(i).min () < 0) // I believe this is checked elsewhere - { - (*current_liboctave_error_handler) - ("index must be one or larger"); - - idx_ok = false; - break; - } - } - } - - if (n_idx <= n_lhs_dims) - { - octave_idx_type last_idx = ra_idx(n_idx-1).max (); - - octave_idx_type sum_el = lhs_dims(n_idx-1); - - for (octave_idx_type i = n_idx; i < n_lhs_dims; i++) - sum_el *= lhs_dims(i); - - if (last_idx > sum_el - 1) - { - (*current_liboctave_error_handler) - ("index exceeds array dimensions"); - - idx_ok = false; - } - } - - if (idx_ok) - { - if (n_idx > 1 - && (all_ones (idx_is_colon))) - { - // A(:,:,:) -- we are deleting elements in all dimensions, so - // the result is [](0x0x0). - - dim_vector newdim = dims (); - newdim(0) = 0; - - resize (newdim); - } - - else if (n_idx > 1 - && num_ones (idx_is_colon) == n_idx - 1 - && num_ones (idx_is_colon_equiv) == n_idx) - { - // A(:,:,j) -- we are deleting elements in one dimension by - // enumerating them. - // - // If we enumerate all of the elements, we should have zero - // elements in that dimension with the same number of elements - // in the other dimensions that we started with. - - dim_vector temp_dims; - temp_dims.resize (n_idx); - - for (octave_idx_type i = 0; i < n_idx; i++) - { - if (idx_is_colon (i)) - temp_dims(i) = lhs_dims(i); - else - temp_dims(i) = 0; - } - - resize (temp_dims); - } - else if (n_idx > 1 && num_ones (idx_is_colon) == n_idx - 1) - { - // We have colons in all indices except for one. - // This index tells us which slice to delete - - if (n_idx < n_lhs_dims) - { - // Collapse dimensions beyond last index. - - if (! (ra_idx(n_idx-1).is_colon ())) - (*current_liboctave_warning_with_id_handler) - ("Octave:fortran-indexing", - "fewer indices than dimensions for N-d array"); - - for (octave_idx_type i = n_idx; i < n_lhs_dims; i++) - lhs_dims(n_idx-1) *= lhs_dims(i); - - lhs_dims.resize (n_idx); - - // Reshape *this. - dimensions = lhs_dims; - } - - int non_col = 0; - - // Find the non-colon column. - - for (octave_idx_type i = 0; i < n_idx; i++) - { - if (! idx_is_colon(i)) - non_col = i; - } - - // The length of the non-colon dimension. - - octave_idx_type non_col_dim = lhs_dims (non_col); - - octave_idx_type num_to_delete = ra_idx(non_col).length (lhs_dims (non_col)); - - if (num_to_delete > 0) - { - int temp = lhs_dims.num_ones (); - - if (non_col_dim == 1) - temp--; - - if (temp == n_idx - 1 && num_to_delete == non_col_dim) - { - // We have A with (1x1x4), where A(1,:,1:4) - // Delete all (0x0x0) - - dim_vector zero_dims (n_idx, 0); - - resize (zero_dims); - } - else - { - // New length of non-colon dimension - // (calculated in the next for loop) - - octave_idx_type new_dim = non_col_dim; - - octave_idx_type iidx = 0; - - for (octave_idx_type j = 0; j < non_col_dim; j++) - if (j == ra_idx(non_col).elem (iidx)) - { - iidx++; - - new_dim--; - - if (iidx == num_to_delete) - break; - } - - // Creating the new nd array after deletions. - - if (new_dim > 0) - { - // Calculate number of elements in new array. - - octave_idx_type num_new_elem=1; - - for (int i = 0; i < n_idx; i++) - { - if (i == non_col) - num_new_elem *= new_dim; - - else - num_new_elem *= lhs_dims(i); - } - - T *new_data = new T [num_new_elem]; - - Array<octave_idx_type> result_idx (n_lhs_dims, 0); - - dim_vector new_lhs_dim = lhs_dims; - - new_lhs_dim(non_col) = new_dim; - - octave_idx_type num_elem = 1; - - octave_idx_type numidx = 0; - - octave_idx_type n = length (); - - for (int i = 0; i < n_lhs_dims; i++) - if (i != non_col) - num_elem *= lhs_dims(i); - - num_elem *= ra_idx(non_col).capacity (); - - for (octave_idx_type i = 0; i < n; i++) - { - if (numidx < num_elem - && is_in (result_idx(non_col), ra_idx(non_col))) - numidx++; - - else - { - Array<octave_idx_type> temp_result_idx = result_idx; - - octave_idx_type num_lgt = how_many_lgt (result_idx(non_col), - ra_idx(non_col)); - - temp_result_idx(non_col) -= num_lgt; - - octave_idx_type kidx - = ::compute_index (temp_result_idx, new_lhs_dim); - - new_data[kidx] = xelem (result_idx); - } - - increment_index (result_idx, lhs_dims); - } - - if (--rep->count <= 0) - delete rep; - - rep = new typename Array<T>::ArrayRep (new_data, - num_new_elem); - - dimensions = new_lhs_dim; - } - } - } - } - else if (n_idx == 1) - { - // This handle cases where we only have one index (not - // colon). The index denotes which elements we should - // delete in the array which can be of any dimension. We - // return a column vector, except for the case where we are - // operating on a row vector. The elements are numerated - // column by column. - // - // A(3,3,3)=2; - // A(3:5) = []; A(6)=[] - - octave_idx_type lhs_numel = numel (); - - idx_vector idx_vec = ra_idx(0); - - idx_vec.freeze (lhs_numel, 0, true); - - idx_vec.sort (true); - - octave_idx_type num_to_delete = idx_vec.length (lhs_numel); - - if (num_to_delete > 0) - { - octave_idx_type new_numel = lhs_numel - num_to_delete; - - T *new_data = new T[new_numel]; - - Array<octave_idx_type> lhs_ra_idx (ndims (), 0); - - octave_idx_type ii = 0; - octave_idx_type iidx = 0; - - for (octave_idx_type i = 0; i < lhs_numel; i++) - { - if (iidx < num_to_delete && i == idx_vec.elem (iidx)) - { - iidx++; - } - else - { - new_data[ii++] = xelem (lhs_ra_idx); - } - - increment_index (lhs_ra_idx, lhs_dims); - } - - if (--(Array<T>::rep)->count <= 0) - delete Array<T>::rep; - - Array<T>::rep = new typename Array<T>::ArrayRep (new_data, new_numel); - - dimensions.resize (2); - - if (lhs_dims.length () == 2 && lhs_dims(1) == 1) - { - dimensions(0) = new_numel; - dimensions(1) = 1; - } - else - { - dimensions(0) = 1; - dimensions(1) = new_numel; - } - } - } - else if (num_ones (idx_is_colon) < n_idx) - { - (*current_liboctave_error_handler) - ("a null assignment can have only one non-colon index"); - } - } -} - -template <class T> -Array<T> -Array<T>::value (void) const -{ - Array<T> retval; - - int n_idx = index_count (); - - if (n_idx == 2) - { - idx_vector *tmp = get_idx (); - - idx_vector idx_i = tmp[0]; - idx_vector idx_j = tmp[1]; - - retval = index (idx_i, idx_j); - } - else if (n_idx == 1) - { - retval = index (idx[0]); - } - else - { - clear_index (); - - (*current_liboctave_error_handler) - ("Array<T>::value: invalid number of indices specified"); - } - - clear_index (); - - return retval; -} - -template <class T> -Array<T> -Array<T>::index (idx_vector& idx_arg, int resize_ok, const T& rfv) const -{ - Array<T> retval; - - dim_vector dv = idx_arg.orig_dimensions (); - - if (dv.length () > 2 || ndims () > 2) - retval = indexN (idx_arg, resize_ok, rfv); - else - { - switch (ndims ()) - { - case 1: - retval = index1 (idx_arg, resize_ok, rfv); - break; - - case 2: - retval = index2 (idx_arg, resize_ok, rfv); - break; - - default: - (*current_liboctave_error_handler) - ("invalid array (internal error)"); - break; - } - } - - return retval; -} - -template <class T> -Array<T> -Array<T>::index1 (idx_vector& idx_arg, int resize_ok, const T& rfv) const -{ - Array<T> retval; - - octave_idx_type len = length (); - - octave_idx_type n = idx_arg.freeze (len, "vector", resize_ok); - - if (idx_arg) - { - if (idx_arg.is_colon_equiv (len)) - { - retval = *this; - } - else if (n == 0) - { - retval.resize_no_fill (0); - } - else - { - retval.resize_no_fill (n); - - for (octave_idx_type i = 0; i < n; i++) - { - octave_idx_type ii = idx_arg.elem (i); - if (ii >= len) - retval.elem (i) = rfv; - else - retval.elem (i) = elem (ii); - } - } - } - - // idx_vector::freeze() printed an error message for us. - - return retval; -} - -template <class T> -Array<T> -Array<T>::index2 (idx_vector& idx_arg, int resize_ok, const T& rfv) const -{ - Array<T> retval; - - assert (ndims () == 2); - - octave_idx_type nr = dim1 (); - octave_idx_type nc = dim2 (); - - octave_idx_type orig_len = nr * nc; - - dim_vector idx_orig_dims = idx_arg.orig_dimensions (); - - octave_idx_type idx_orig_rows = idx_arg.orig_rows (); - octave_idx_type idx_orig_columns = idx_arg.orig_columns (); - - if (idx_arg.is_colon ()) - { - // Fast magic colon processing. - - octave_idx_type result_nr = nr * nc; - octave_idx_type result_nc = 1; - - retval = Array<T> (*this, dim_vector (result_nr, result_nc)); - } - else if (nr == 1 && nc == 1) - { - Array<T> tmp = Array<T>::index1 (idx_arg, resize_ok); - - octave_idx_type len = tmp.length (); - - if (len >= idx_orig_dims.numel ()) - retval = Array<T> (tmp, idx_orig_dims); - } - else if (nr == 1 || nc == 1) - { - // If indexing a vector with a matrix, return value has same - // shape as the index. Otherwise, it has same orientation as - // indexed object. - - Array<T> tmp = Array<T>::index1 (idx_arg, resize_ok); - - octave_idx_type len = tmp.length (); - - if (idx_orig_rows == 1 || idx_orig_columns == 1) - { - if (nr == 1) - retval = Array<T> (tmp, dim_vector (1, len)); - else - retval = Array<T> (tmp, dim_vector (len, 1)); - } - else if (len >= idx_orig_dims.numel ()) - retval = Array<T> (tmp, idx_orig_dims); - } - else - { - (*current_liboctave_warning_with_id_handler) - ("Octave:fortran-indexing", "single index used for matrix"); - - // This code is only for indexing matrices. The vector - // cases are handled above. - - idx_arg.freeze (nr * nc, "matrix", resize_ok); - - if (idx_arg) - { - octave_idx_type result_nr = idx_orig_rows; - octave_idx_type result_nc = idx_orig_columns; - - retval.resize_no_fill (result_nr, result_nc); - - octave_idx_type k = 0; - for (octave_idx_type j = 0; j < result_nc; j++) - { - for (octave_idx_type i = 0; i < result_nr; i++) - { - octave_idx_type ii = idx_arg.elem (k++); - if (ii >= orig_len) - retval.elem (i, j) = rfv; - else - { - octave_idx_type fr = ii % nr; - octave_idx_type fc = (ii - fr) / nr; - retval.elem (i, j) = elem (fr, fc); - } - } - } - } - // idx_vector::freeze() printed an error message for us. - } - - return retval; -} - -template <class T> -Array<T> -Array<T>::indexN (idx_vector& ra_idx, int resize_ok, const T& rfv) const -{ - Array<T> retval; - - dim_vector dv = dims (); - - int n_dims = dv.length (); - - octave_idx_type orig_len = dv.numel (); - - dim_vector idx_orig_dims = ra_idx.orig_dimensions (); - - if (ra_idx.is_colon ()) - { - // Fast magic colon processing. - - retval = Array<T> (*this, dim_vector (orig_len, 1)); - } - else - { - bool vec_equiv = vector_equivalent (dv); - - if (! (vec_equiv || ra_idx.is_colon ())) - (*current_liboctave_warning_with_id_handler) - ("Octave:fortran-indexing", "single index used for N-d array"); - - octave_idx_type frozen_len - = ra_idx.freeze (orig_len, "nd-array", resize_ok); - - if (ra_idx) - { - dim_vector result_dims; - - if (vec_equiv && ! orig_len == 1) - { - result_dims = dv; - - for (int i = 0; i < n_dims; i++) - { - if (result_dims(i) != 1) - { - // All but this dim should be one. - result_dims(i) = frozen_len; - break; - } - } - } - else - result_dims = idx_orig_dims; - - result_dims.chop_trailing_singletons (); - - retval.resize (result_dims); - - octave_idx_type n = result_dims.numel (); - - octave_idx_type k = 0; - - for (octave_idx_type i = 0; i < n; i++) - { - octave_idx_type ii = ra_idx.elem (k++); - - if (ii >= orig_len) - retval.elem (i) = rfv; - else - retval.elem (i) = elem (ii); - } - } - } - - return retval; -} - -template <class T> -Array<T> -Array<T>::index (idx_vector& idx_i, idx_vector& idx_j, int resize_ok, - const T& rfv) const -{ - Array<T> retval; - - if (ndims () != 2) - { - Array<idx_vector> ra_idx (2); - ra_idx(0) = idx_i; - ra_idx(1) = idx_j; - return index (ra_idx, resize_ok, rfv); - } - - octave_idx_type nr = dim1 (); - octave_idx_type nc = dim2 (); - - octave_idx_type n = idx_i.freeze (nr, "row", resize_ok); - octave_idx_type m = idx_j.freeze (nc, "column", resize_ok); - - if (idx_i && idx_j) - { - if (idx_i.orig_empty () || idx_j.orig_empty () || n == 0 || m == 0) - { - retval.resize_no_fill (n, m); - } - else if (idx_i.is_colon_equiv (nr) && idx_j.is_colon_equiv (nc)) - { - retval = *this; - } - else - { - retval.resize_no_fill (n, m); - - for (octave_idx_type j = 0; j < m; j++) - { - octave_idx_type jj = idx_j.elem (j); - for (octave_idx_type i = 0; i < n; i++) - { - octave_idx_type ii = idx_i.elem (i); - if (ii >= nr || jj >= nc) - retval.elem (i, j) = rfv; - else - retval.elem (i, j) = elem (ii, jj); - } - } - } - } - - // idx_vector::freeze() printed an error message for us. - - return retval; -} - -template <class T> -Array<T> -Array<T>::index (Array<idx_vector>& ra_idx, int resize_ok, const T& rfv) const -{ - // This function handles all calls with more than one idx. - // For (3x3x3), the call can be A(2,5), A(2,:,:), A(3,2,3) etc. - - Array<T> retval; - - int n_dims = dimensions.length (); - - // Remove trailing singletons in ra_idx, but leave at least ndims - // elements. - - octave_idx_type ra_idx_len = ra_idx.length (); - - bool trim_trailing_singletons = true; - for (octave_idx_type j = ra_idx_len; j > n_dims; j--) - { - idx_vector iidx = ra_idx (ra_idx_len-1); - if (iidx.capacity () == 1 && trim_trailing_singletons) - ra_idx_len--; - else - trim_trailing_singletons = false; - - if (! resize_ok) - { - for (octave_idx_type i = 0; i < iidx.capacity (); i++) - if (iidx (i) != 0) - { - (*current_liboctave_error_handler) - ("index exceeds N-d array dimensions"); - - return retval; - } - } - } - - ra_idx.resize (ra_idx_len); - - dim_vector new_dims = dims (); - dim_vector frozen_lengths; - - if (!ra_idx (ra_idx_len - 1).orig_empty () && ra_idx_len < n_dims) - frozen_lengths = short_freeze (ra_idx, dimensions, resize_ok); - else - { - new_dims.resize (ra_idx_len, 1); - frozen_lengths = freeze (ra_idx, new_dims, resize_ok); - } - - if (all_ok (ra_idx)) - { - if (any_orig_empty (ra_idx) || frozen_lengths.any_zero ()) - { - frozen_lengths.chop_trailing_singletons (); - - retval.resize (frozen_lengths); - } - else if (frozen_lengths.length () == n_dims - && all_colon_equiv (ra_idx, dimensions)) - { - retval = *this; - } - else - { - dim_vector frozen_lengths_for_resize = frozen_lengths; - - frozen_lengths_for_resize.chop_trailing_singletons (); - - retval.resize (frozen_lengths_for_resize); - - octave_idx_type n = retval.length (); - - Array<octave_idx_type> result_idx (ra_idx.length (), 0); - - Array<octave_idx_type> elt_idx; - - for (octave_idx_type i = 0; i < n; i++) - { - elt_idx = get_elt_idx (ra_idx, result_idx); - - octave_idx_type numelem_elt = get_scalar_idx (elt_idx, new_dims); - - if (numelem_elt >= length () || numelem_elt < 0) - retval.elem (i) = rfv; - else - retval.elem (i) = elem (numelem_elt); - - increment_index (result_idx, frozen_lengths); - - } - } - } - - return retval; -} - -template <class T> bool ascending_compare (T a, T b) { @@ -2851,7 +2114,7 @@ if (nnr == 1) { octave_idx_type n = nnc + std::abs (k); - d = Array<T> (dim_vector (n, n), resize_fill_value (T ())); + d = Array<T> (dim_vector (n, n), resize_fill_value ()); for (octave_idx_type i = 0; i < nnc; i++) d.xelem (i+roff, i+coff) = elem (0, i); @@ -2859,7 +2122,7 @@ else { octave_idx_type n = nnr + std::abs (k); - d = Array<T> (dim_vector (n, n), resize_fill_value (T ())); + d = Array<T> (dim_vector (n, n), resize_fill_value ()); for (octave_idx_type i = 0; i < nnr; i++) d.xelem (i+roff, i+coff) = elem (i, 0); @@ -2870,1011 +2133,6 @@ return d; } -// FIXME -- this is a mess. - -template <class LT, class RT> -int -assign (Array<LT>& lhs, const Array<RT>& rhs, const LT& rfv) -{ - int n_idx = lhs.index_count (); - - // kluge... - if (lhs.ndims () == 0) - lhs.resize_no_fill (0, 0); - - return (lhs.ndims () == 2 - && (n_idx == 1 - || (n_idx < 3 - && rhs.ndims () == 2 - && rhs.rows () == 0 && rhs.columns () == 0))) - ? assign2 (lhs, rhs, rfv) : assignN (lhs, rhs, rfv); -} - -template <class LT, class RT> -int -assign1 (Array<LT>& lhs, const Array<RT>& rhs, const LT& rfv) -{ - int retval = 1; - - idx_vector *tmp = lhs.get_idx (); - - idx_vector lhs_idx = tmp[0]; - - octave_idx_type lhs_len = lhs.length (); - octave_idx_type rhs_len = rhs.length (); - - octave_idx_type n = lhs_idx.freeze (lhs_len, "vector", true); - - if (n != 0) - { - dim_vector lhs_dims = lhs.dims (); - - if (lhs_len != 0 - || lhs_dims.all_zero () - || (lhs_dims.length () == 2 && lhs_dims(0) < 2)) - { - if (rhs_len == n || rhs_len == 1) - { - lhs.make_unique (); - - octave_idx_type max_idx = lhs_idx.max () + 1; - if (max_idx > lhs_len) - lhs.resize_and_fill (max_idx, rfv); - } - - if (rhs_len == n) - { - lhs.make_unique (); - - if (lhs_idx.is_colon ()) - { - for (octave_idx_type i = 0; i < n; i++) - lhs.xelem (i) = rhs.elem (i); - } - else - { - for (octave_idx_type i = 0; i < n; i++) - { - octave_idx_type ii = lhs_idx.elem (i); - lhs.xelem (ii) = rhs.elem (i); - } - } - } - else if (rhs_len == 1) - { - lhs.make_unique (); - - RT scalar = rhs.elem (0); - - if (lhs_idx.is_colon ()) - { - for (octave_idx_type i = 0; i < n; i++) - lhs.xelem (i) = scalar; - } - else - { - for (octave_idx_type i = 0; i < n; i++) - { - octave_idx_type ii = lhs_idx.elem (i); - lhs.xelem (ii) = scalar; - } - } - } - else - { - lhs.clear_index (); - - (*current_liboctave_error_handler) - ("A(I) = X: X must be a scalar or a vector with same length as I"); - - retval = 0; - } - } - else - { - lhs.clear_index (); - - (*current_liboctave_error_handler) - ("A(I) = X: unable to resize A"); - - retval = 0; - } - } - else if (lhs_idx.is_colon ()) - { - dim_vector lhs_dims = lhs.dims (); - - if (lhs_dims.all_zero ()) - { - lhs.make_unique (); - - lhs.resize_no_fill (rhs_len); - - for (octave_idx_type i = 0; i < rhs_len; i++) - lhs.xelem (i) = rhs.elem (i); - } - else if (rhs_len != lhs_len) - { - lhs.clear_index (); - - (*current_liboctave_error_handler) - ("A(:) = X: A must be the same size as X"); - } - } - else if (! (rhs_len == 1 || rhs_len == 0)) - { - lhs.clear_index (); - - (*current_liboctave_error_handler) - ("A([]) = X: X must also be an empty matrix or a scalar"); - - retval = 0; - } - - lhs.clear_index (); - - return retval; -} - -#define MAYBE_RESIZE_LHS \ - do \ - { \ - octave_idx_type max_row_idx = idx_i_is_colon ? rhs_nr : idx_i.max () + 1; \ - octave_idx_type max_col_idx = idx_j_is_colon ? rhs_nc : idx_j.max () + 1; \ - \ - octave_idx_type new_nr = max_row_idx > lhs_nr ? max_row_idx : lhs_nr; \ - octave_idx_type new_nc = max_col_idx > lhs_nc ? max_col_idx : lhs_nc; \ - \ - lhs.resize_and_fill (new_nr, new_nc, rfv); \ - } \ - while (0) - -template <class LT, class RT> -int -assign2 (Array<LT>& lhs, const Array<RT>& rhs, const LT& rfv) -{ - int retval = 1; - - int n_idx = lhs.index_count (); - - octave_idx_type lhs_nr = lhs.rows (); - octave_idx_type lhs_nc = lhs.cols (); - - Array<RT> xrhs = rhs; - - octave_idx_type rhs_nr = xrhs.rows (); - octave_idx_type rhs_nc = xrhs.cols (); - - if (xrhs.ndims () > 2) - { - xrhs = xrhs.squeeze (); - - dim_vector dv_tmp = xrhs.dims (); - - switch (dv_tmp.length ()) - { - case 1: - // FIXME -- this case should be unnecessary, because - // squeeze should always return an object with 2 dimensions. - if (rhs_nr == 1) - rhs_nc = dv_tmp.elem (0); - break; - - case 2: - rhs_nr = dv_tmp.elem (0); - rhs_nc = dv_tmp.elem (1); - break; - - default: - lhs.clear_index (); - (*current_liboctave_error_handler) - ("Array<T>::assign2: Dimension mismatch"); - return 0; - } - } - - bool rhs_is_scalar = rhs_nr == 1 && rhs_nc == 1; - - idx_vector *tmp = lhs.get_idx (); - - idx_vector idx_i; - idx_vector idx_j; - - if (n_idx > 1) - idx_j = tmp[1]; - - if (n_idx > 0) - idx_i = tmp[0]; - - if (n_idx == 2) - { - octave_idx_type n = idx_i.freeze (lhs_nr, "row", true); - octave_idx_type m = idx_j.freeze (lhs_nc, "column", true); - - int idx_i_is_colon = idx_i.is_colon (); - int idx_j_is_colon = idx_j.is_colon (); - - if (lhs_nr == 0 && lhs_nc == 0) - { - if (idx_i_is_colon) - n = rhs_nr; - - if (idx_j_is_colon) - m = rhs_nc; - } - - if (idx_i && idx_j) - { - if (rhs_is_scalar && n >= 0 && m >= 0) - { - // No need to do anything if either of the indices - // are empty. - - if (n > 0 && m > 0) - { - lhs.make_unique (); - - MAYBE_RESIZE_LHS; - - RT scalar = xrhs.elem (0, 0); - - for (octave_idx_type j = 0; j < m; j++) - { - octave_idx_type jj = idx_j.elem (j); - for (octave_idx_type i = 0; i < n; i++) - { - octave_idx_type ii = idx_i.elem (i); - lhs.xelem (ii, jj) = scalar; - } - } - } - } - else if ((n == 1 || m == 1) - && (rhs_nr == 1 || rhs_nc == 1) - && n * m == rhs_nr * rhs_nc) - { - lhs.make_unique (); - - MAYBE_RESIZE_LHS; - - if (n > 0 && m > 0) - { - octave_idx_type k = 0; - - for (octave_idx_type j = 0; j < m; j++) - { - octave_idx_type jj = idx_j.elem (j); - for (octave_idx_type i = 0; i < n; i++) - { - octave_idx_type ii = idx_i.elem (i); - lhs.xelem (ii, jj) = xrhs.elem (k++); - } - } - } - } - else if (n == rhs_nr && m == rhs_nc) - { - lhs.make_unique (); - - MAYBE_RESIZE_LHS; - - if (n > 0 && m > 0) - { - for (octave_idx_type j = 0; j < m; j++) - { - octave_idx_type jj = idx_j.elem (j); - for (octave_idx_type i = 0; i < n; i++) - { - octave_idx_type ii = idx_i.elem (i); - lhs.xelem (ii, jj) = xrhs.elem (i, j); - } - } - } - } - else if (n == 0 && m == 0) - { - if (! (rhs_is_scalar || (rhs_nr == 0 || rhs_nc == 0))) - { - lhs.clear_index (); - - (*current_liboctave_error_handler) - ("A([], []) = X: X must be an empty matrix or a scalar"); - - retval = 0; - } - } - else - { - lhs.clear_index (); - - (*current_liboctave_error_handler) - ("A(I, J) = X: X must be a scalar or the number of elements in I must match the number of rows in X and the number of elements in J must match the number of columns in X"); - - retval = 0; - } - } - // idx_vector::freeze() printed an error message for us. - } - else if (n_idx == 1) - { - int lhs_is_empty = lhs_nr == 0 || lhs_nc == 0; - - if (lhs_is_empty || (lhs_nr == 1 && lhs_nc == 1)) - { - octave_idx_type lhs_len = lhs.length (); - - idx_i.freeze (lhs_len, 0, true); - - if (idx_i) - { - if (lhs_is_empty - && idx_i.is_colon () - && ! (rhs_nr == 1 || rhs_nc == 1)) - { - (*current_liboctave_warning_with_id_handler) - ("Octave:fortran-indexing", - "A(:) = X: X is not a vector or scalar"); - } - else - { - octave_idx_type idx_nr = idx_i.orig_rows (); - octave_idx_type idx_nc = idx_i.orig_columns (); - - if (! (rhs_nr == idx_nr && rhs_nc == idx_nc)) - (*current_liboctave_warning_with_id_handler) - ("Octave:fortran-indexing", - "A(I) = X: X does not have same shape as I"); - } - - if (assign1 (lhs, xrhs, rfv)) - { - octave_idx_type len = lhs.length (); - - if (len > 0) - { - // The following behavior is much simplified - // over previous versions of Octave. It - // seems to be compatible with Matlab. - - lhs.dimensions = dim_vector (1, lhs.length ()); - } - } - else - retval = 0; - } - // idx_vector::freeze() printed an error message for us. - } - else if (lhs_nr == 1) - { - idx_i.freeze (lhs_nc, "vector", true); - - if (idx_i) - { - if (assign1 (lhs, xrhs, rfv)) - lhs.dimensions = dim_vector (1, lhs.length ()); - else - retval = 0; - } - // idx_vector::freeze() printed an error message for us. - } - else if (lhs_nc == 1) - { - idx_i.freeze (lhs_nr, "vector", true); - - if (idx_i) - { - if (assign1 (lhs, xrhs, rfv)) - lhs.dimensions = dim_vector (lhs.length (), 1); - else - retval = 0; - } - // idx_vector::freeze() printed an error message for us. - } - else - { - if (! idx_i.is_colon ()) - (*current_liboctave_warning_with_id_handler) - ("Octave:fortran-indexing", "single index used for matrix"); - - octave_idx_type len = idx_i.freeze (lhs_nr * lhs_nc, "matrix"); - - if (idx_i) - { - if (len == 0) - { - if (! (rhs_is_scalar || (rhs_nr == 0 || rhs_nc == 0))) - { - lhs.clear_index (); - - (*current_liboctave_error_handler) - ("A([]) = X: X must be an empty matrix or scalar"); - - retval = 0; - } - } - else if (len == rhs_nr * rhs_nc) - { - lhs.make_unique (); - - if (idx_i.is_colon ()) - { - for (octave_idx_type i = 0; i < len; i++) - lhs.xelem (i) = xrhs.elem (i); - } - else - { - for (octave_idx_type i = 0; i < len; i++) - { - octave_idx_type ii = idx_i.elem (i); - lhs.xelem (ii) = xrhs.elem (i); - } - } - } - else if (rhs_is_scalar) - { - lhs.make_unique (); - - RT scalar = rhs.elem (0, 0); - - if (idx_i.is_colon ()) - { - for (octave_idx_type i = 0; i < len; i++) - lhs.xelem (i) = scalar; - } - else - { - for (octave_idx_type i = 0; i < len; i++) - { - octave_idx_type ii = idx_i.elem (i); - lhs.xelem (ii) = scalar; - } - } - } - else - { - lhs.clear_index (); - - (*current_liboctave_error_handler) - ("A(I) = X: X must be a scalar or a matrix with the same size as I"); - - retval = 0; - } - } - // idx_vector::freeze() printed an error message for us. - } - } - else - { - (*current_liboctave_error_handler) - ("invalid number of indices for matrix expression"); - - retval = 0; - } - - lhs.clear_index (); - - return retval; -} - -template <class LT, class RT> -int -assignN (Array<LT>& lhs, const Array<RT>& rhs, const LT& rfv) -{ - int retval = 1; - - dim_vector rhs_dims = rhs.dims (); - - octave_idx_type rhs_dims_len = rhs_dims.length (); - - bool rhs_is_scalar = is_scalar (rhs_dims); - - int n_idx = lhs.index_count (); - - idx_vector *idx_vex = lhs.get_idx (); - - Array<idx_vector> idx = conv_to_array (idx_vex, n_idx); - - if (n_idx == 0) - { - lhs.clear_index (); - - (*current_liboctave_error_handler) - ("invalid number of indices for matrix expression"); - - retval = 0; - } - else if (n_idx == 1) - { - idx_vector iidx = idx(0); - int iidx_is_colon = iidx.is_colon (); - - if (! iidx_is_colon) - (*current_liboctave_warning_with_id_handler) - ("Octave:fortran-indexing", "single index used for N-d array"); - - octave_idx_type lhs_len = lhs.length (); - - octave_idx_type len = iidx.freeze (lhs_len, "N-d arrray"); - - if (iidx) - { - if (len == 0) - { - if (! (rhs_dims.all_ones () || rhs_dims.any_zero ())) - { - lhs.clear_index (); - - (*current_liboctave_error_handler) - ("A([]) = X: X must be an empty matrix or scalar"); - - retval = 0; - } - } - else if (len == rhs.length ()) - { - lhs.make_unique (); - - if (iidx_is_colon) - { - for (octave_idx_type i = 0; i < len; i++) - lhs.xelem (i) = rhs.elem (i); - } - else - { - for (octave_idx_type i = 0; i < len; i++) - { - octave_idx_type ii = iidx.elem (i); - - lhs.xelem (ii) = rhs.elem (i); - } - } - } - else if (rhs_is_scalar) - { - RT scalar = rhs.elem (0); - - lhs.make_unique (); - - if (iidx_is_colon) - { - for (octave_idx_type i = 0; i < len; i++) - lhs.xelem (i) = scalar; - } - else - { - for (octave_idx_type i = 0; i < len; i++) - { - octave_idx_type ii = iidx.elem (i); - - lhs.xelem (ii) = scalar; - } - } - } - else - { - lhs.clear_index (); - - (*current_liboctave_error_handler) - ("A(I) = X: X must be a scalar or a matrix with the same size as I"); - - retval = 0; - } - - // idx_vector::freeze() printed an error message for us. - } - } - else - { - // Maybe expand to more dimensions. - - dim_vector lhs_dims = lhs.dims (); - - octave_idx_type lhs_dims_len = lhs_dims.length (); - - dim_vector final_lhs_dims = lhs_dims; - - dim_vector frozen_len; - - octave_idx_type orig_lhs_dims_len = lhs_dims_len; - - bool orig_empty = lhs_dims.all_zero (); - - if (n_idx < lhs_dims_len) - { - // Collapse dimensions beyond last index. Note that we - // delay resizing LHS until we know that the assignment will - // succeed. - - if (! (idx(n_idx-1).is_colon ())) - (*current_liboctave_warning_with_id_handler) - ("Octave:fortran-indexing", - "fewer indices than dimensions for N-d array"); - - for (int i = n_idx; i < lhs_dims_len; i++) - lhs_dims(n_idx-1) *= lhs_dims(i); - - lhs_dims.resize (n_idx); - - lhs_dims_len = lhs_dims.length (); - } - - // Resize. - - dim_vector new_dims; - new_dims.resize (n_idx); - - if (orig_empty) - { - if (rhs_is_scalar) - { - for (int i = 0; i < n_idx; i++) - { - if (idx(i).is_colon ()) - new_dims(i) = 1; - else - new_dims(i) = idx(i).orig_empty () ? 0 : idx(i).max () + 1; - } - } - else if (is_vector (rhs_dims)) - { - int ncolon = 0; - int fcolon = 0; - octave_idx_type new_dims_numel = 1; - int new_dims_vec = 0; - for (int i = 0; i < n_idx; i++) - if (idx(i).is_colon ()) - { - ncolon ++; - if (ncolon == 1) - fcolon = i; - } - else - { - octave_idx_type cap = idx(i).capacity (); - new_dims_numel *= cap; - if (cap != 1) - new_dims_vec ++; - } - - if (ncolon == n_idx) - { - new_dims = rhs_dims; - new_dims.resize (n_idx); - for (int i = rhs_dims_len; i < n_idx; i++) - new_dims (i) = 1; - } - else - { - octave_idx_type rhs_dims_numel = rhs_dims.numel (); - - for (int i = 0; i < n_idx; i++) - new_dims(i) = idx(i).orig_empty () ? 0 : idx(i).max () + 1; - - if (new_dims_numel != rhs_dims_numel && - ncolon > 0 && new_dims_numel == 1) - { - if (ncolon == rhs_dims_len) - { - int k = 0; - for (int i = 0; i < n_idx; i++) - if (idx(i).is_colon ()) - new_dims (i) = rhs_dims (k++); - } - else - new_dims (fcolon) = rhs_dims_numel; - } - else if (new_dims_numel != rhs_dims_numel || new_dims_vec > 1) - { - lhs.clear_index (); - - (*current_liboctave_error_handler) - ("A(IDX-LIST) = RHS: mismatched index and RHS dimension"); - return retval; - } - } - } - else - { - int k = 0; - for (int i = 0; i < n_idx; i++) - { - if (idx(i).is_colon ()) - { - if (k < rhs_dims_len) - new_dims(i) = rhs_dims(k++); - else - new_dims(i) = 1; - } - else - { - octave_idx_type nelem = idx(i).capacity (); - - if (nelem >= 1 - && (k < rhs_dims_len && nelem == rhs_dims(k))) - k++; - else if (nelem != 1) - { - lhs.clear_index (); - - (*current_liboctave_error_handler) - ("A(IDX-LIST) = RHS: mismatched index and RHS dimension"); - return retval; - } - new_dims(i) = idx(i).orig_empty () ? 0 : - idx(i).max () + 1; - } - } - } - } - else - { - for (int i = 0; i < n_idx; i++) - { - // We didn't start out with all zero dimensions, so if - // index is a colon, it refers to the current LHS - // dimension. Otherwise, it is OK to enlarge to a - // dimension given by the largest index, but if that - // index is a colon the new dimension is singleton. - - if (i < lhs_dims_len - && (idx(i).is_colon () - || idx(i).orig_empty () - || idx(i).max () < lhs_dims(i))) - new_dims(i) = lhs_dims(i); - else if (! idx(i).is_colon ()) - new_dims(i) = idx(i).max () + 1; - else - new_dims(i) = 1; - } - } - - if (retval != 0) - { - if (! orig_empty - && n_idx < orig_lhs_dims_len - && new_dims(n_idx-1) != lhs_dims(n_idx-1)) - { - // We reshaped and the last dimension changed. This has to - // be an error, because we don't know how to undo that - // later... - - lhs.clear_index (); - - (*current_liboctave_error_handler) - ("array index %d (= %d) for assignment requires invalid resizing operation", - n_idx, new_dims(n_idx-1)); - - retval = 0; - } - else - { - // Determine final dimensions for LHS and reset the - // current size of the LHS. Note that we delay actually - // resizing LHS until we know that the assignment will - // succeed. - - if (n_idx < orig_lhs_dims_len) - { - for (int i = 0; i < n_idx-1; i++) - final_lhs_dims(i) = new_dims(i); - } - else - final_lhs_dims = new_dims; - - lhs_dims_len = new_dims.length (); - - frozen_len = freeze (idx, new_dims, true); - - for (int i = 0; i < idx.length (); i++) - { - if (! idx(i)) - { - retval = 0; - lhs.clear_index (); - return retval; - } - } - - if (rhs_is_scalar) - { - lhs.make_unique (); - - if (n_idx < orig_lhs_dims_len) - lhs = lhs.reshape (lhs_dims); - - lhs.resize_and_fill (new_dims, rfv); - - if (! final_lhs_dims.any_zero ()) - { - RT scalar = rhs.elem (0); - - if (n_idx == 1) - { - idx_vector iidx = idx(0); - - octave_idx_type len = frozen_len(0); - - if (iidx.is_colon ()) - { - for (octave_idx_type i = 0; i < len; i++) - lhs.xelem (i) = scalar; - } - else - { - for (octave_idx_type i = 0; i < len; i++) - { - octave_idx_type ii = iidx.elem (i); - - lhs.xelem (ii) = scalar; - } - } - } - else if (lhs_dims_len == 2 && n_idx == 2) - { - idx_vector idx_i = idx(0); - idx_vector idx_j = idx(1); - - octave_idx_type i_len = frozen_len(0); - octave_idx_type j_len = frozen_len(1); - - if (idx_i.is_colon()) - { - for (octave_idx_type j = 0; j < j_len; j++) - { - octave_idx_type off = new_dims (0) * - idx_j.elem (j); - for (octave_idx_type i = 0; i < i_len; i++) - lhs.xelem (i + off) = scalar; - } - } - else - { - for (octave_idx_type j = 0; j < j_len; j++) - { - octave_idx_type off = new_dims (0) * - idx_j.elem (j); - for (octave_idx_type i = 0; i < i_len; i++) - { - octave_idx_type ii = idx_i.elem (i); - lhs.xelem (ii + off) = scalar; - } - } - } - } - else - { - octave_idx_type n = Array<LT>::get_size (frozen_len); - - Array<octave_idx_type> result_idx (lhs_dims_len, 0); - - for (octave_idx_type i = 0; i < n; i++) - { - Array<octave_idx_type> elt_idx = get_elt_idx (idx, result_idx); - - lhs.xelem (elt_idx) = scalar; - - increment_index (result_idx, frozen_len); - } - } - } - } - else - { - // RHS is matrix or higher dimension. - - octave_idx_type n = Array<LT>::get_size (frozen_len); - - if (n != rhs.numel ()) - { - lhs.clear_index (); - - (*current_liboctave_error_handler) - ("A(IDX-LIST) = X: X must be a scalar or size of X must equal number of elements indexed by IDX-LIST"); - - retval = 0; - } - else - { - lhs.make_unique (); - - if (n_idx < orig_lhs_dims_len) - lhs = lhs.reshape (lhs_dims); - - lhs.resize_and_fill (new_dims, rfv); - - if (! final_lhs_dims.any_zero ()) - { - if (n_idx == 1) - { - idx_vector iidx = idx(0); - - octave_idx_type len = frozen_len(0); - - if (iidx.is_colon ()) - { - for (octave_idx_type i = 0; i < len; i++) - lhs.xelem (i) = rhs.elem (i); - } - else - { - for (octave_idx_type i = 0; i < len; i++) - { - octave_idx_type ii = iidx.elem (i); - - lhs.xelem (ii) = rhs.elem (i); - } - } - } - else if (lhs_dims_len == 2 && n_idx == 2) - { - idx_vector idx_i = idx(0); - idx_vector idx_j = idx(1); - - octave_idx_type i_len = frozen_len(0); - octave_idx_type j_len = frozen_len(1); - octave_idx_type k = 0; - - if (idx_i.is_colon()) - { - for (octave_idx_type j = 0; j < j_len; j++) - { - octave_idx_type off = new_dims (0) * - idx_j.elem (j); - for (octave_idx_type i = 0; - i < i_len; i++) - lhs.xelem (i + off) = rhs.elem (k++); - } - } - else - { - for (octave_idx_type j = 0; j < j_len; j++) - { - octave_idx_type off = new_dims (0) * - idx_j.elem (j); - for (octave_idx_type i = 0; i < i_len; i++) - { - octave_idx_type ii = idx_i.elem (i); - lhs.xelem (ii + off) = rhs.elem (k++); - } - } - } - - } - else - { - n = Array<LT>::get_size (frozen_len); - - Array<octave_idx_type> result_idx (lhs_dims_len, 0); - - for (octave_idx_type i = 0; i < n; i++) - { - Array<octave_idx_type> elt_idx = get_elt_idx (idx, result_idx); - - lhs.xelem (elt_idx) = rhs.elem (i); - - increment_index (result_idx, frozen_len); - } - } - } - } - } - } - } - - lhs.clear_index (); - - if (retval != 0) - lhs = lhs.reshape (final_lhs_dims); - } - - if (retval != 0) - lhs.chop_trailing_singletons (); - - lhs.clear_index (); - - return retval; -} - template <class T> void Array<T>::print_info (std::ostream& os, const std::string& prefix) const
--- a/liboctave/Array.h +++ b/liboctave/Array.h @@ -3,6 +3,7 @@ Copyright (C) 1993, 1994, 1995, 1996, 1997, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 John W. Eaton +Copyright (C) 2008 Jaroslav Hajek <highegg@gmail.com> This file is part of Octave. @@ -31,23 +32,15 @@ #include <iostream> #include "dim-vector.h" +#include "idx-vector.h" #include "lo-utils.h" #include "oct-sort.h" #include "quit.h" -class idx_vector; - // One dimensional array class. Handles the reference counting for // all the derived classes. template <class T> -T -resize_fill_value (const T& x) -{ - return x; -} - -template <class T> class Array { @@ -148,16 +141,12 @@ protected: - mutable idx_vector *idx; - mutable int idx_count; - Array (T *d, octave_idx_type n) - : rep (new typename Array<T>::ArrayRep (d, n)), dimensions (n), - idx (0), idx_count (0) { } + : rep (new typename Array<T>::ArrayRep (d, n)), dimensions (n) { } Array (T *d, const dim_vector& dv) : rep (new typename Array<T>::ArrayRep (d, get_size (dv))), - dimensions (dv), idx (0), idx_count (0) { } + dimensions (dv) { } private: @@ -184,16 +173,13 @@ public: Array (void) - : rep (nil_rep ()), dimensions (), - idx (0), idx_count (0) { rep->count++; } + : rep (nil_rep ()), dimensions () { rep->count++; } explicit Array (octave_idx_type n) - : rep (new typename Array<T>::ArrayRep (n)), dimensions (n), - idx (0), idx_count (0) { } + : rep (new typename Array<T>::ArrayRep (n)), dimensions (n) { } explicit Array (octave_idx_type n, const T& val) - : rep (new typename Array<T>::ArrayRep (n)), dimensions (n), - idx (0), idx_count (0) + : rep (new typename Array<T>::ArrayRep (n)), dimensions (n) { fill (val); } @@ -202,13 +188,13 @@ template <class U> Array (const Array<U>& a) : rep (new typename Array<T>::ArrayRep (coerce (a.data (), a.length ()), a.length ())), - dimensions (a.dimensions), idx (0), idx_count (0) + dimensions (a.dimensions) { } // No type conversion case. Array (const Array<T>& a) - : rep (a.rep), dimensions (a.dimensions), idx (0), idx_count (0) + : rep (a.rep), dimensions (a.dimensions) { rep->count++; } @@ -217,11 +203,11 @@ Array (const dim_vector& dv) : rep (new typename Array<T>::ArrayRep (get_size (dv))), - dimensions (dv), idx (0), idx_count (0) { } + dimensions (dv) { } Array (const dim_vector& dv, const T& val) : rep (new typename Array<T>::ArrayRep (get_size (dv))), - dimensions (dv), idx (0), idx_count (0) + dimensions (dv) { fill (val); } @@ -419,43 +405,6 @@ Array<T> ipermute (const Array<octave_idx_type>& vec) const { return permute (vec, true); } - void resize_no_fill (octave_idx_type n); - void resize_and_fill (octave_idx_type n, const T& val); - - // !!! WARNING !!! -- the following resize_no_fill and - // resize_and_fill functions are public because template friends - // don't work properly with versions of gcc earlier than 3.3. You - // should use these functions only in classes that are derived - // from Array<T>. - - // protected: - - void resize_no_fill (octave_idx_type r, octave_idx_type c); - void resize_and_fill (octave_idx_type r, octave_idx_type c, const T& val); - - void resize_no_fill (octave_idx_type r, octave_idx_type c, octave_idx_type p); - void resize_and_fill (octave_idx_type r, octave_idx_type c, octave_idx_type p, const T& val); - - void resize_no_fill (const dim_vector& dv); - void resize_and_fill (const dim_vector& dv, const T& val); - -public: - - void resize (octave_idx_type n) { resize_no_fill (n); } - - void resize (octave_idx_type n, const T& val) { resize_and_fill (n, val); } - - void resize (const dim_vector& dv) { resize_no_fill (dv); } - - void resize (const dim_vector& dv, const T& val) - { resize_and_fill (dv, val); } - - Array<T>& insert (const Array<T>& a, octave_idx_type r, octave_idx_type c); - Array<T>& insert2 (const Array<T>& a, octave_idx_type r, octave_idx_type c); - Array<T>& insertN (const Array<T>& a, octave_idx_type r, octave_idx_type c); - - Array<T>& insert (const Array<T>& a, const Array<octave_idx_type>& idx); - bool is_square (void) const { return (dim1 () == dim2 ()); } bool is_empty (void) const { return numel () == 0; } @@ -482,47 +431,98 @@ void maybe_delete_dims (void); - void clear_index (void) const; + // Indexing without resizing. + + Array<T> index (const idx_vector& i) const; - void set_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; - int index_count (void) const { return idx_count; } + static T resize_fill_value (); + + // Resizing (with fill). - idx_vector *get_idx (void) const { return idx; } + void resize_fill (octave_idx_type n, const T& rfv); - void maybe_delete_elements (idx_vector& i); + void resize_fill (octave_idx_type nr, octave_idx_type nc, const T& rfv); + + void resize_fill (const dim_vector& dv, const T& rfv); - void maybe_delete_elements_1 (idx_vector& i); + // Resizing with default fill. + // Rationale: + // These use the default fill value rather than leaving memory uninitialized. + // Resizing without fill leaves the resulting array in a rather weird state, + // where part of the data is initialized an part isn't. - void maybe_delete_elements_2 (idx_vector& i); + void resize (octave_idx_type n) + { resize_fill (n, resize_fill_value ()); } - void maybe_delete_elements (idx_vector& i, idx_vector& j); + // FIXME: This method cannot be defined here because it would clash with + // void resize (octave_idx_type, const T&) + // (these become undistinguishable when T = octave_idx_type). + // In the future, I think the resize (.., const T& rfv) overloads should go away + // in favor of using resize_fill. - void maybe_delete_elements (idx_vector& i, idx_vector& j, idx_vector& k); + // void resize (octave_idx_type nr, octave_idx_type nc) + // { resize_fill (nr, nc, resize_fill_value ()); } - void maybe_delete_elements (Array<idx_vector>& ra_idx); + void resize (dim_vector dv) + { resize_fill (dv, resize_fill_value ()); } + + // FIXME: these are here for backward compatibility. They should go away in favor + // of using resize_fill directly. + void resize (octave_idx_type n, const T& rfv) + { resize_fill (n, static_cast<T> (rfv)); } - Array<T> value (void) const; + void resize (octave_idx_type nr, octave_idx_type nc, const T& rfv) + { resize_fill (nr, nc, rfv); } - Array<T> index (idx_vector& i, int resize_ok = 0, - const T& rfv = resize_fill_value (T ())) const; + void resize (dim_vector dv, const T& rfv) + { resize_fill (dv, rfv); } + + // Indexing with possible resizing and fill + // FIXME: This is really a corner case, that should better be handled + // directly in liboctinterp. - Array<T> index1 (idx_vector& i, int resize_ok = 0, - const T& rfv = resize_fill_value (T ())) const; + Array<T> index (const idx_vector& i, bool resize_ok, + const T& rfv = resize_fill_value ()) const; + + Array<T> index (const idx_vector& i, const idx_vector& j, + bool resize_ok, const T& rfv = resize_fill_value ()) const; - Array<T> index2 (idx_vector& i, int resize_ok = 0, - const T& rfv = resize_fill_value (T ())) const; + Array<T> index (const Array<idx_vector>& ia, + bool resize_ok, const T& rfv = resize_fill_value ()) const; + + // Indexed assignment (always with resize & fill). + + void assign (const idx_vector& i, const Array<T>& rhs, + const T& rfv = resize_fill_value ()); - Array<T> indexN (idx_vector& i, int resize_ok = 0, - const T& rfv = resize_fill_value (T ())) const; + void assign (const idx_vector& i, const idx_vector& j, const Array<T>& rhs, + const T& rfv = resize_fill_value ()); + + void assign (const Array<idx_vector>& ia, const Array<T>& rhs, + const T& rfv = resize_fill_value ()); + + // Deleting elements. + + // A(I) = [] (with a single subscript) + void delete_elements (const idx_vector& i); - Array<T> index (idx_vector& i, idx_vector& j, int resize_ok = 0, - const T& rfv = resize_fill_value (T ())) const; + // 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); - Array<T> index (Array<idx_vector>& ra_idx, int resize_ok = 0, - const T& rfv = resize_fill_value (T ())) const; + // FIXME: are these required? What exactly are they supposed to do?. - // static T resize_fill_value (void) { return T (); } + Array<T>& insert (const Array<T>& a, octave_idx_type r, octave_idx_type c); + Array<T>& insert2 (const Array<T>& a, octave_idx_type r, octave_idx_type c); + Array<T>& insertN (const Array<T>& a, octave_idx_type r, octave_idx_type c); + + Array<T>& insert (const Array<T>& a, const Array<octave_idx_type>& idx); void print_info (std::ostream& os, const std::string& prefix) const; @@ -558,48 +558,18 @@ } }; -// NOTE: these functions should be friends of the Array<T> class and -// Array<T>::dimensions should be protected, not public, but we can't -// do that because of bugs in gcc prior to 3.3. - -template <class LT, class RT> -/* friend */ int -assign (Array<LT>& lhs, const Array<RT>& rhs, const LT& rfv); - -template <class LT, class RT> -/* friend */ int -assign1 (Array<LT>& lhs, const Array<RT>& rhs, const LT& rfv); - -template <class LT, class RT> -/* friend */ int -assign2 (Array<LT>& lhs, const Array<RT>& rhs, const LT& rfv); - -template <class LT, class RT> -/* friend */ int -assignN (Array<LT>& lhs, const Array<RT>& rhs, const LT& rfv); +#define INSTANTIATE_ARRAY(T, API) \ + template class API Array<T> -template <class LT, class RT> -int -assign (Array<LT>& lhs, const Array<RT>& rhs) -{ - return assign (lhs, rhs, resize_fill_value (LT ())); -} +// FIXME: These are here for compatibility. In current implementation, only +// homogeneous array assignments are actually instantiated. I think heterogeneous +// indexed assignments are rare enough to be implemented via conversion first. +// This decision may still be revised, that's why these macros stay here. +#define INSTANTIATE_ARRAY_AND_ASSIGN(T, API) \ + INSTANTIATE_ARRAY(T, API) -#define INSTANTIATE_ARRAY_ASSIGN(LT, RT, API) \ - template API int assign (Array<LT>&, const Array<RT>&, const LT&); \ - template API int assign1 (Array<LT>&, const Array<RT>&, const LT&); \ - template API int assign2 (Array<LT>&, const Array<RT>&, const LT&); \ - template API int assignN (Array<LT>&, const Array<RT>&, const LT&); \ - template API int assign (Array<LT>&, const Array<RT>&) - - -#define INSTANTIATE_ARRAY(T, API) \ - template class API Array<T>; \ - template API T resize_fill_value (const T&); \ - -#define INSTANTIATE_ARRAY_AND_ASSIGN(T, API) \ - INSTANTIATE_ARRAY (T, API); \ - INSTANTIATE_ARRAY_ASSIGN (T, T, API) +#define INSTANTIATE_ARRAY_ASSIGN(LT, RT, API) + // do nothing #define INSTANTIATE_ARRAY_SORT(T) \ template class octave_sort<T>; \
--- a/liboctave/Array2.h +++ b/liboctave/Array2.h @@ -92,10 +92,11 @@ return *this; } - void resize (octave_idx_type r, octave_idx_type c) { this->resize_no_fill (r, c); } + void resize (octave_idx_type r, octave_idx_type c) + { Array<T>::resize_fill (r, c, Array<T>::resize_fill_value ()); } void resize (octave_idx_type r, octave_idx_type c, const T& val) - { this->resize_and_fill (r, c, val); } + { Array<T>::resize_fill (r, c, val); } Array2<T>& insert (const Array2<T>& a, octave_idx_type r, octave_idx_type c) {
--- a/liboctave/Array3.h +++ b/liboctave/Array3.h @@ -71,10 +71,11 @@ return *this; } - void resize (octave_idx_type r, octave_idx_type c, octave_idx_type p) { this->resize_no_fill (r, c, p); } + void resize (octave_idx_type r, octave_idx_type c, octave_idx_type p) + { Array<T>::resize (dim_vector (r, c, p)); } void resize (octave_idx_type r, octave_idx_type c, octave_idx_type p, const T& val) - { this->resize_and_fill (r, c, p, val); } + { Array<T>::resize_fill (dim_vector (r, c, p), val); } Array3<T> sort (octave_idx_type dim = 0, sortmode mode = ASCENDING) const {
--- a/liboctave/ArrayN.h +++ b/liboctave/ArrayN.h @@ -93,12 +93,6 @@ ArrayN<T> ipermute (const Array<octave_idx_type>& vec) const { return Array<T>::ipermute (vec); } - void resize (const dim_vector& dv) - { this->resize_no_fill (dv); } - - void resize (const dim_vector& dv, const T& val) - { Array<T>::resize (dv, val); } - ArrayN<T> squeeze (void) const { return Array<T>::squeeze (); } ArrayN<T> transpose (void) const { return Array<T>::transpose (); } @@ -117,21 +111,21 @@ } ArrayN<T> index (idx_vector& i, int resize_ok = 0, - const T& rfv = resize_fill_value (T ())) const + const T& rfv = Array<T>::resize_fill_value ()) const { Array<T> tmp = Array<T>::index (i, resize_ok, rfv); return ArrayN<T> (tmp, tmp.dims ()); } ArrayN<T> index (idx_vector& i, idx_vector& j, int resize_ok = 0, - const T& rfv = resize_fill_value (T ())) const + const T& rfv = Array<T>::resize_fill_value ()) const { Array<T> tmp = Array<T>::index (i, j, resize_ok, rfv); return ArrayN<T> (tmp, tmp.dims ()); } ArrayN<T> index (Array<idx_vector>& ra_idx, int resize_ok = 0, - const T& rfv = resize_fill_value (T ())) const + const T& rfv = Array<T>::resize_fill_value ()) const { Array<T> tmp = Array<T>::index (ra_idx, resize_ok, rfv); return ArrayN<T> (tmp, tmp.dims ());
--- a/liboctave/ChangeLog +++ b/liboctave/ChangeLog @@ -1,3 +1,28 @@ +2008-10-28 Jaroslav Hajek <highegg@gmail.com> + + * Array-C.cc Array-d.cc Array-f.cc Array-fC.cc Array-i.cc Array-s.cc: + Don't use semicolon after INSTANTIATE_ARRAY_ASSIGN. + * Array-util.h (zero_dims_inquire): New declarations. + (is_in, how_many_lgt, short_freeze): Remove declarations. + * Array-util.cc (zero_dims_inquire): New functions. + (is_in, how_many_lgt, short_freeze): Remove functions. + * Array.cc (Array<T>::index, Array<T>::resize_fill, Array<T>::resize, + Array<T>::assign, Array<T>::delete_elements): + Rewrite. + * Array.h (Array<T>::index, Array<T>::resize_fill, Array<T>::resize, + Array<T>::assign, Array<T>::delete_elements): + Rewrite interfaces. + * Array2.h (Array2<T>::resize): Call Array<T>::resize_fill. + * Array3.h (Array3<T>::resize): Call Array<T>::resize_fill. + * ArrayN.h (ArrayN<T>::resize): Remove declarations. + (ArrayN<T>::index): Fix call to resize_fill_value. + * Sparse.cc (assign, assign1): Use zero-based indices. + * chMatrix.h: Include mx-op-defs.h + * dim-vector.h (dim_vector::any_neg, dim_vector::chop_all_singletons, + dim_vector::redim): New member functions. + * idx-vector.cc: Mostly rewrite. + * idx-vector.h: Mostly rewrite. + 2008-10-29 Michael Goffioul <michael.goffioul@gmail.com> * lo-specfun.cc (cbesj, cbesy, cbesi, cbesk, cbesh1, cbesh2): Do not
--- a/liboctave/Sparse.cc +++ b/liboctave/Sparse.cc @@ -2507,7 +2507,7 @@ for (octave_idx_type i = 0; i < n; i++) { - new_idx.xelem(i) = sidx[i]->i + 1; + new_idx.xelem(i) = sidx[i]->i; rhs_idx[i] = sidx[i]->idx; } @@ -3061,7 +3061,7 @@ for (octave_idx_type i = 0; i < n; i++) { - new_idx.xelem(i) = sidx[i]->i + 1; + new_idx.xelem(i) = sidx[i]->i; rhs_idx_i[i] = sidx[i]->idx; } @@ -3100,7 +3100,7 @@ for (octave_idx_type i = 0; i < m; i++) { - new_idx.xelem(i) = sidx[i]->i + 1; + new_idx.xelem(i) = sidx[i]->i; rhs_idx_j[i] = sidx[i]->idx; } @@ -3318,7 +3318,7 @@ for (octave_idx_type i = 0; i < len; i++) { - new_idx.xelem(i) = sidx[i]->i + 1; + new_idx.xelem(i) = sidx[i]->i; rhs_idx[i] = sidx[i]->idx; }
--- a/liboctave/chMatrix.h +++ b/liboctave/chMatrix.h @@ -29,6 +29,7 @@ #include "MArray2.h" #include "mx-defs.h" +#include "mx-op-defs.h" #include "str-vec.h" class
--- a/liboctave/dim-vector.h +++ b/liboctave/dim-vector.h @@ -130,6 +130,18 @@ } } + void chop_all_singletons (void) + { + int j = 0; + for (int i = 0; i < ndims; i++) + { + if (dims[i] != 1) + dims[j++] = dims[i]; + } + if (j == 1) dims[1] = 1; + ndims = j > 2 ? j : 2; + } + private: // No assignment! @@ -309,12 +321,26 @@ return retval; } + bool any_neg (void) const + { + int n_dims = length (), i; + for (i = 0; i < n_dims; i++) + if (elem (i) < 0) break; + return i < n_dims; + } + void chop_trailing_singletons (void) { make_unique (); rep->chop_trailing_singletons (); } + void chop_all_singletons (void) + { + make_unique (); + rep->chop_all_singletons (); + } + dim_vector squeeze (void) const { dim_vector new_dims = *this; @@ -439,6 +465,44 @@ return true; } + + // Forces certain dimensionality, preserving numel (). Missing dimensions are + // set to 1, redundant are folded into the trailing one. If n = 1, the result + // is 2d and the second dim is 1 (dim_vectors are always at least 2D). + // If the original dimensions were all zero, the padding value is zero. + dim_vector redim (int n) const + { + int n_dims = length (); + if (n_dims == n) + return *this; + else + { + dim_vector retval; + retval.resize (n == 1 ? 2 : n, 1); + + bool zeros = true; + for (int i = 0; i < n && i < n_dims; i++) + { + retval(i) = elem (i); + zeros = zeros && elem (i) == 0; + } + + if (n < n_dims) + { + octave_idx_type k = 1; + for (int i = n; i < n_dims; i++) + k *= elem (i); + retval(n - 1) *= k; + } + else if (zeros) + { + for (int i = n_dims; i < n; i++) + retval.elem (i) = 0; + } + + return retval; + } + } }; static inline bool
--- a/liboctave/idx-vector.cc +++ b/liboctave/idx-vector.cc @@ -2,6 +2,7 @@ Copyright (C) 1993, 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2002, 2003, 2004, 2005, 2006, 2007 John W. Eaton +Copyright (C) 2008 Jaroslav Hajek This file is part of Octave. @@ -29,587 +30,539 @@ #include <iostream> +#include "idx-vector.h" +#include "Array.h" #include "Range.h" -#include "boolNDArray.h" -#include "dColVector.h" -#include "dNDArray.h" -#include "idx-vector.h" #include "lo-error.h" #include "lo-mappers.h" -#define IDX_VEC_REP idx_vector::idx_vector_rep +static void gripe_invalid_index () +{ + (*current_liboctave_error_handler) + ("subscript indices must be either positive integers or logicals."); +} + +static void gripe_invalid_range () +{ + (*current_liboctave_error_handler) + ("invalid range used as index."); +} + +static void gripe_index_out_of_range () +{ + (*current_liboctave_error_handler) + ("internal error: idx_vector index out of range."); +} + +idx_vector::idx_colon_rep::idx_colon_rep (char c) +{ + if (c != ':') + { + (*current_liboctave_error_handler) + ("internal error: invalid character converted to idx_vector. Must be ':'."); + err = true; + } +} + +octave_idx_type +idx_vector::idx_colon_rep::checkelem (octave_idx_type i) const +{ + if (i < 0) + { + gripe_index_out_of_range (); + return 0; + } + else + return i; +} -IDX_VEC_REP::idx_vector_rep (const IDX_VEC_REP& a) - : data (0), len (a.len), num_zeros (a.num_zeros), num_ones (a.num_ones), - range_base (a.range_base), range_step (a.range_step), - max_val (a.max_val), min_val (a.min_val), - frozen_at_z_len (a.frozen_at_z_len), frozen_len (a.frozen_len), - colon (a.colon), range(a.range), initialized (a.initialized), - frozen (a.frozen), colon_equiv_checked (a.colon_equiv_checked), - colon_equiv (a.colon_equiv), orig_dims (a.orig_dims) +std::ostream& +idx_vector::idx_colon_rep::print (std::ostream& os) const +{ + return os << ":"; +} + +idx_vector::idx_range_rep::idx_range_rep (octave_idx_type _start, octave_idx_type _limit, + octave_idx_type _step) + : start(_start), len (_step ? (_limit - _start) / _step : -1), step (_step) +{ + if (len < 0) + { + gripe_invalid_range (); + err = true; + } + else if (start < 0) + { + gripe_invalid_index (); + err = true; + } +} + +static void gripe_non_int_range () { - if (len > 0) + (*current_liboctave_error_handler) + ("If a range is used as subscript, all elements are expected to be integers."); +} + +idx_vector::idx_range_rep::idx_range_rep (const Range& r) + : start (0), len (r.nelem ()), step (1) +{ + if (len < 0) { - if (! range) - { - data = new octave_idx_type [len]; - - for (octave_idx_type i = 0; i < len; i++) - data[i] = a.data[i]; - } + gripe_invalid_range (); + err = true; + } + else if (len > 0) + { + if (r.all_elements_are_ints ()) + { + start = static_cast<octave_idx_type> (r.base ()) - 1; + step = static_cast<octave_idx_type> (r.inc ()); + } + else + { + gripe_non_int_range (); + err = true; + } } } octave_idx_type -IDX_VEC_REP::tree_to_mat_idx (double x, bool& conversion_error) +idx_vector::idx_range_rep::checkelem (octave_idx_type i) const { - octave_idx_type retval = -1; - - conversion_error = false; - - if (D_NINT (x) != x) + if (i < 0 || i >= len) { - (*current_liboctave_error_handler) - ("expecting integer index, found %f", x); - - conversion_error = true; + gripe_index_out_of_range (); + return 0; } else - retval = static_cast<octave_idx_type> (x - 1); - - return retval; + return start + i*step; } -static inline bool -idx_is_inf_or_nan (double x) +idx_vector::idx_base_rep * +idx_vector::idx_range_rep::sort_uniq_clone (bool) { - bool retval = false; - - if (xisnan (x)) + if (step < 0) + return new idx_range_rep (start + (len - 1)*step, len, -step, DIRECT); + else { - (*current_liboctave_error_handler) ("NaN invalid as index"); - retval = true; + count++; + return this; } - else if (xisinf (x)) - { - (*current_liboctave_error_handler) ("Inf invalid as index"); - retval = true; - } - - return retval; } -IDX_VEC_REP::idx_vector_rep (const ColumnVector& v) - : data (0), len (v.length ()), num_zeros (0), num_ones (0), - range_base (0), range_step (0), max_val (0), min_val (0), count (1), - frozen_at_z_len (0), frozen_len (0), colon (0), range(0), - initialized (0), frozen (0), colon_equiv_checked (0), - colon_equiv (0), orig_dims (len, 1) +std::ostream& +idx_vector::idx_range_rep::print (std::ostream& os) const { - if (len == 0) - { - initialized = 1; - return; - } - else - { - data = new octave_idx_type [len]; - - bool conversion_error = false; - - for (octave_idx_type i = 0; i < len; i++) - { - double d = v.elem (i); - - if (idx_is_inf_or_nan (d)) - return; - else - data[i] = tree_to_mat_idx (d, conversion_error); - - if (conversion_error) - return; - } - } - - init_state (); + os << start << ':' << step << ':' << start + len*step; + return os; } -IDX_VEC_REP::idx_vector_rep (const NDArray& nda) - : data (0), len (nda.length ()), num_zeros (0), num_ones (0), - range_base (0), range_step (0), max_val (0), min_val (0), count (1), - frozen_at_z_len (0), frozen_len (0), colon (0), range(0), - initialized (0), frozen (0), colon_equiv_checked (0), - colon_equiv (0), orig_dims (nda.dims ()) +inline octave_idx_type +convert_index (octave_idx_type i, bool& conv_error, + octave_idx_type& ext) { - if (len == 0) + if (i <= 0) + conv_error = true; + if (ext < i) ext = i; + return i - 1; +} + +inline octave_idx_type +convert_index (double x, bool& conv_error, octave_idx_type& ext) +{ + octave_idx_type i; + if (xisnan (x) || xisinf (x)) { - initialized = 1; - return; + i = 0; + conv_error = true; } else { - octave_idx_type k = 0; - data = new octave_idx_type [len]; - - bool conversion_error = false; - - for (octave_idx_type i = 0; i < len; i++) - { - double d = nda.elem (i); - - if (idx_is_inf_or_nan (d)) - return; - else - data[k++] = tree_to_mat_idx (d, conversion_error); - - if (conversion_error) - return; - } + i = static_cast<octave_idx_type> (x); + if (static_cast<double> (i) != x) + conv_error = true; } - init_state (); + return convert_index (i, conv_error, ext); } -IDX_VEC_REP::idx_vector_rep (const Range& r) - : data (0), len (r.nelem ()), num_zeros (0), num_ones (0), - range_base (0), range_step (0), max_val (0), min_val (0), - count (1), frozen_at_z_len (0), frozen_len (0), colon (0), - range(1), initialized (0), frozen (0), - colon_equiv_checked (0), colon_equiv (0), orig_dims (1, len) +inline octave_idx_type +convert_index (float x, bool& conv_error, octave_idx_type& ext) { - if (len < 0) - { - (*current_liboctave_error_handler) ("invalid range used as index"); - return; - } - else if (len == 0) - { - initialized = 1; - return; - } - - if (r.all_elements_are_ints ()) - { - range_base = static_cast<octave_idx_type> (r.base () - 1); - range_step = static_cast<octave_idx_type> (r.inc ()); - - init_state (); - } - else - (*current_liboctave_error_handler) - ("expecting integer index, found non integer Range"); + return convert_index (static_cast<double> (x), conv_error, ext); } -IDX_VEC_REP::idx_vector_rep (double d) - : data (0), len (1), num_zeros (0), num_ones (0), - range_base (0), range_step (0), max_val (0), min_val (0), - count (1), frozen_at_z_len (0), frozen_len (0), colon (0), - range(1), initialized (0), frozen (0), colon_equiv_checked (0), - colon_equiv (0), orig_dims (1, 1) +template <class T> +inline octave_idx_type +convert_index (octave_int<T> x, bool& conv_error, + octave_idx_type& ext) { - if (idx_is_inf_or_nan (d)) - return; - else - { - bool conversion_error = false; - - range_base = tree_to_mat_idx (d, conversion_error); - range_step = 1; - - if (conversion_error) - return; - } - - init_state (); + octave_idx_type i = octave_int<octave_idx_type> (x).value (); + return convert_index (i, conv_error, ext); } -IDX_VEC_REP::idx_vector_rep (octave_idx_type i) - : data (0), len (1), num_zeros (0), num_ones (0), - range_base (tree_to_mat_idx (i)), range_step (1), - max_val (0), min_val (0), count (1), frozen_at_z_len (0), - frozen_len (0), colon (0), range(1), initialized (0), - frozen (0), colon_equiv_checked (0), colon_equiv (0), orig_dims (1, 1) +template <class T> +idx_vector::idx_scalar_rep::idx_scalar_rep (T x) { - init_state (); + octave_idx_type dummy; + data = convert_index (x, err, dummy); + if (err) gripe_invalid_index (); } -IDX_VEC_REP::idx_vector_rep (char c) - : data (0), len (0), num_zeros (0), num_ones (0), range_base (0), - range_step (0), max_val (0), min_val (0), count (1), - frozen_at_z_len (0), frozen_len (0), colon (1), range(0), - initialized (0), frozen (0), colon_equiv_checked (0), - colon_equiv (0), orig_dims (0, 0) +idx_vector::idx_scalar_rep::idx_scalar_rep (octave_idx_type i) + : data (i) { - assert (c == ':'); - - init_state (); -} - -IDX_VEC_REP::idx_vector_rep (bool b) - : data (0), len (b ? 1 : 0), num_zeros (0), num_ones (0), range_base (0), - range_step (0), max_val (0), min_val (0), count (1), - frozen_at_z_len (0), frozen_len (0), colon (0), range(0), - initialized (0), frozen (0), colon_equiv_checked (0), - colon_equiv (0), orig_dims (len, len) -{ - if (len == 0) - initialized = 1; - else + if (data < 0) { - data = new octave_idx_type [len]; - data[0] = 0; - init_state (); + gripe_invalid_index (); + err = true; } } -IDX_VEC_REP::idx_vector_rep (const boolNDArray& bnda) - : data (0), len (bnda.nnz ()), num_zeros (0), num_ones (0), - range_base (0), range_step (0), max_val (0), min_val (0), - count (1), frozen_at_z_len (0), frozen_len (0), colon (0), - range(0), initialized (0), frozen (0), - colon_equiv_checked (0), colon_equiv (0), orig_dims () +octave_idx_type +idx_vector::idx_scalar_rep::checkelem (octave_idx_type i) const +{ + if (i != 0) gripe_index_out_of_range (); + return data; +} + +std::ostream& idx_vector::idx_scalar_rep::print (std::ostream& os) const +{ + return os << data; +} + +template <class T> +idx_vector::idx_vector_rep::idx_vector_rep (const Array<T>& nda) + : data (0), len (nda.numel ()), ext (0), aowner (0), orig_dims (nda.dims ()) +{ + if (len != 0) + { + octave_idx_type *d = new octave_idx_type[len]; + for (octave_idx_type i = 0; i < len; i++) + d[i] = convert_index (nda.xelem (i), err, ext); + data = d; + + if (err) gripe_invalid_index (); + } +} + +// Note that this makes a shallow copy of the index array. + +idx_vector::idx_vector_rep::idx_vector_rep (const Array<octave_idx_type>& inda) + : data (inda.data ()), len (inda.numel ()), ext (0), + aowner (new Array<octave_idx_type> (inda)), orig_dims (inda.dims ()) { + if (len != 0) + { + octave_idx_type max = -1; + for (octave_idx_type i = 0; i < len; i++) + { + octave_idx_type k = inda.xelem (i); + if (k < 0) + err = true; + else if (k > max) + max = k; + } + + ext = max + 1; + + if (err) gripe_invalid_index (); + } +} + +idx_vector::idx_vector_rep::idx_vector_rep (bool b) + : data (0), len (b ? 1 : 0), ext (0), aowner (0), orig_dims (len, len) +{ + if (len != 0) + { + octave_idx_type *d = new octave_idx_type [1]; + d[0] = 0; + data = d; + ext = 1; + } +} + +idx_vector::idx_vector_rep::idx_vector_rep (const Array<bool>& bnda) + : data (0), len (0), ext (0), aowner (0), orig_dims () +{ + for (octave_idx_type i = 0, l = bnda.numel (); i < l; i++) + if (bnda.xelem (i)) len++; + dim_vector dv = bnda.dims (); orig_dims = ((dv.length () == 2 && dv(0) == 1) ? dim_vector (1, len) : orig_dims = dim_vector (len, 1)); - if (len == 0) - initialized = 1; - else + if (len != 0) { - data = new octave_idx_type [len]; + octave_idx_type *d = new octave_idx_type [len]; octave_idx_type ntot = bnda.length (); - for (octave_idx_type i = 0, k = 0; i < ntot && k < len; i++) - if (bnda.elem (i)) - data[k++] = i; + octave_idx_type k = 0; + for (octave_idx_type i = 0; i < ntot; i++) + if (bnda.xelem (i)) d[k++] = i; - init_state (); + data = d; + + ext = k; } } -IDX_VEC_REP& -IDX_VEC_REP::operator = (const IDX_VEC_REP& a) -{ - if (this != &a) - { - delete [] data; - len = a.len; - - if (a.data) - { - data = new octave_idx_type [len]; - - for (octave_idx_type i = 0; i < len; i++) - data[i] = a.data[i]; - } - else - data = 0; - - num_zeros = a.num_zeros; - num_ones = a.num_ones; - range_base = a.range_base; - range_step = a.range_step; - max_val = a.max_val; - min_val = a.min_val; - frozen_at_z_len = a.frozen_at_z_len; - frozen_len = a.frozen_len; - colon = a.colon; - range = a.range; - initialized = a.initialized; - frozen = a.frozen; - colon_equiv_checked = a.colon_equiv_checked; - colon_equiv = a.colon_equiv; - orig_dims = a.orig_dims; - } - - return *this; -} - -void -IDX_VEC_REP::init_state (void) -{ - num_zeros = 0; - num_ones = 0; - - if (colon) - { - min_val = 0; - max_val = 0; - } - else if (range) - { - if (range_step >= 0) - { - min_val = range_base; - max_val = (len > 0) ? range_base + (len-1)*range_step : range_base; - } - else - { - max_val = range_base; - min_val = (len > 0) ? range_base + (len-1)*range_step : range_base; - } - - if ((range_base <= 0 && range_step > 0) - || (range_base >= 0 && range_step < 0)) - num_zeros = 1; - - if ((range_base <= 1 && range_step > 0) - || (range_base >= 1 && range_step < 0)) - num_zeros = 0; - } +idx_vector::idx_vector_rep::~idx_vector_rep (void) +{ + if (aowner) + delete aowner; else - { - min_val = max_val = data[0]; - - octave_idx_type i = 0; - do - { - if (data[i] == -1) - num_zeros++; - else if (data[i] == 0) - num_ones++; - - if (data[i] > max_val) - max_val = data[i]; - - if (data[i] < min_val) - min_val = data[i]; - } - while (++i < len); - } - - initialized = 1; + delete [] data; } octave_idx_type -IDX_VEC_REP::checkelem (octave_idx_type n) const +idx_vector::idx_vector_rep::checkelem (octave_idx_type n) const { if (n < 0 || n >= len) { - (*current_liboctave_error_handler) ("idx-vector: index out of range"); + gripe_invalid_index (); return 0; } - return elem (n); -} - -static inline int -intcmp (const void *ii, const void *jj) -{ - return (*(static_cast<const octave_idx_type *> (ii)) - *(static_cast<const octave_idx_type *> (jj))); -} - -static inline void -sort_data (octave_idx_type *d, octave_idx_type l) -{ - qsort (d, l, sizeof (octave_idx_type), intcmp); -} - -static inline octave_idx_type -make_uniq (octave_idx_type *d, octave_idx_type l) -{ - if (l < 2) - return l; - - octave_idx_type k = 0; - for (octave_idx_type ii = 1; ii < l; ii++) - { - if (d[ii] != d[k]) - { - k++; - d[k] = d[ii]; - } - } - return k+1; -} - -static inline octave_idx_type * -copy_data (const octave_idx_type *d, octave_idx_type l) -{ - octave_idx_type *new_data = new octave_idx_type [l]; - - for (octave_idx_type ii = 0; ii < l; ii++) - new_data[ii] = d[ii]; - - return new_data; + return xelem (n); } -int -IDX_VEC_REP::is_colon_equiv (octave_idx_type n, int sort_uniq) +idx_vector::idx_base_rep * +idx_vector::idx_vector_rep::sort_uniq_clone (bool uniq) { - if (! colon_equiv_checked) - { - if (colon) - { - colon_equiv = 1; - } - else if (range) - { - colon_equiv = (range_base == 0 - && len == n - && (range_step == 1 - || (range_step == -1 && sort_uniq))); - } - else if (static_cast<octave_idx_type> (len) > 1) - { - if (sort_uniq) - { - octave_idx_type *tmp_data = copy_data (data, len); - - sort_data (tmp_data, len); - - octave_idx_type tmp_len = make_uniq (tmp_data, len); + octave_idx_type *new_data = new octave_idx_type[len]; + std::copy (data, data + len, new_data); + std::sort (new_data, new_data + len); + octave_idx_type new_len; + if (uniq) + new_len = std::unique (new_data, new_data + len) - new_data; + else + new_len = len; - colon_equiv = (tmp_len == n - && tmp_data[0] == 0 - && tmp_data[tmp_len-1] == tmp_len - 1); - - delete [] tmp_data; - } - else - { - if (len == n) - { - colon_equiv = 1; - - for (octave_idx_type ii = 0; ii < n; ii++) - if (data[ii] != ii) - { - colon_equiv = 0; - break; - } - } - } - } - else - colon_equiv = (len == n && (n == 0 || (n == 1 && data[0] == 0))); - - colon_equiv_checked = 1; - } - - return colon_equiv; + return new idx_vector_rep (new_data, new_len, ext, orig_dims, DIRECT); } -void -IDX_VEC_REP::sort (bool uniq) +std::ostream& +idx_vector::idx_vector_rep::print (std::ostream& os) const { - if (range && len) - { - if (range_step < 0) - { - range_base += (len-1)*(range_step); - range_step = -range_step; - } - } - else if (len > 1) - { - sort_data (data, len); - - if (uniq) - len = make_uniq (data, len); - } -} - -void -IDX_VEC_REP::shorten (octave_idx_type n) -{ - if (n > 0 && n <= len) - len = n; - else - (*current_liboctave_error_handler) - ("idx_vector::shorten: internal error!"); -} - -std::ostream& -IDX_VEC_REP::print (std::ostream& os) const -{ - if (colon) - os << "colon" << std::endl; - else if (range) - os << "range_base: " << range_base - << ", range_step: " << range_step << std::endl; - else - { - for (octave_idx_type ii = 0; ii < len; ii++) - os << data[ii] << "\n"; - } + os << '['; + for (octave_idx_type ii = 0; ii < len - 1; ii++) + os << data[ii] << ',' << ' '; + if (len > 0) os << data[len-1]; os << ']'; return os; } -octave_idx_type -IDX_VEC_REP::freeze (octave_idx_type z_len, const char *tag, bool resize_ok) -{ - if (frozen) - return frozen_len; +const idx_vector idx_vector::colon (new idx_vector::idx_colon_rep ()); - frozen_len = -1; +bool idx_vector::maybe_reduce (octave_idx_type n, const idx_vector& j, + octave_idx_type nj) +{ + bool reduced = false; - if (colon) - frozen_len = z_len; - else + // Empty index always reduces. + if (rep->length (n) == 0) { - if (len == 0) - frozen_len = 0; - else - { - max_val = max (); - min_val = min (); - - if (min_val < 0) - { - if (tag) - (*current_liboctave_error_handler) - ("invalid %s index = %d", tag, min_val+1); - else - (*current_liboctave_error_handler) - ("invalid index = %d", min_val+1); + *this = idx_vector (); + return true; + } - initialized = 0; - } - else if (! resize_ok && max_val >= z_len) - { - if (tag) - (*current_liboctave_error_handler) - ("invalid %s index = %d", tag, max_val+1); - else - (*current_liboctave_error_handler) - ("invalid index = %d", max_val+1); - - initialized = 0; - } - else - { - if (max_val >= z_len) - { - if (tag) - (*current_liboctave_warning_with_id_handler) - ("Octave:resize-on-range-error", - "resizing object with %s index = %d out of bounds", - tag, max_val+1); - else - (*current_liboctave_warning_with_id_handler) - ("Octave:resize-on-range-error", - "resizing object with index = %d out of bounds", - max_val+1); - } - - frozen_len = length (z_len); - } - } + switch (j.idx_class ()) + { + case class_colon: + if (rep->is_colon_equiv (n)) + { + // (:,:) reduces to (:) + *this = colon; + reduced = true; + } + else if (rep->idx_class () == class_scalar) + { + // (i,:) reduces to a range. + idx_scalar_rep * r = dynamic_cast<idx_scalar_rep *> (rep); + octave_idx_type k = r->get_data (); + *this = new idx_range_rep (k, nj, n, DIRECT); + reduced = true; + } + break; + case class_range: + if (rep->is_colon_equiv (n)) + { + // (:,i:j) reduces to a range (the step must be 1) + idx_range_rep * rj = dynamic_cast<idx_range_rep *> (j.rep); + if (rj->get_step () == 1) + { + octave_idx_type s = rj->get_start (), l = rj->length (nj); + *this = new idx_range_rep (s * n, l * n, 1, DIRECT); + reduced = true; + } + } + else if (rep->idx_class () == class_scalar) + { + // (k,i:d:j) reduces to a range. + idx_scalar_rep * r = dynamic_cast<idx_scalar_rep *> (rep); + idx_range_rep * rj = dynamic_cast<idx_range_rep *> (j.rep); + octave_idx_type k = r->get_data (); + octave_idx_type s = rj->get_start (), l = rj->length (nj); + octave_idx_type t = rj->get_step (); + *this = new idx_range_rep (n * s + k, l, n * t, DIRECT); + reduced = true; + } + break; + case class_scalar: + switch (rep->idx_class ()) + { + case class_scalar: + { + // (i,j) reduces to a single index. + idx_scalar_rep * r = dynamic_cast<idx_scalar_rep *> (rep); + idx_scalar_rep * rj = dynamic_cast<idx_scalar_rep *> (j.rep); + octave_idx_type k = r->get_data () + n * rj->get_data (); + *this = new idx_scalar_rep (k, DIRECT); + reduced = true; + } + break; + case class_range: + { + // (i:d:j,k) reduces to a range. + idx_range_rep * r = dynamic_cast<idx_range_rep *> (rep); + idx_scalar_rep * rj = dynamic_cast<idx_scalar_rep *> (j.rep); + octave_idx_type s = r->get_start (), l = r->length (nj); + octave_idx_type t = r->get_step (); + octave_idx_type k = rj->get_data (); + *this = new idx_range_rep (n * k + s, l, t, DIRECT); + reduced = true; + } + break; + case class_colon: + { + // (:,k) reduces to a range. + idx_scalar_rep * rj = dynamic_cast<idx_scalar_rep *> (j.rep); + octave_idx_type k = rj->get_data (); + *this = new idx_range_rep (n * k, n, 1, DIRECT); + reduced = true; + } + break; + default: + break; + } + break; + default: + break; } - frozen = 1; + return reduced; +} + +bool idx_vector::is_cont_range (octave_idx_type n, + octave_idx_type& l, octave_idx_type& u) const +{ + bool res = false; + switch (rep->idx_class ()) + { + case class_colon: + l = 0; u = n; + res = true; + break; + case class_range: + { + idx_range_rep * r = dynamic_cast<idx_range_rep *> (rep); + if (r->get_step () == 1) + { + l = r->get_start (); + u = l + r->length (n); + res = true; + } + } + break; + case class_scalar: + { + idx_scalar_rep * r = dynamic_cast<idx_scalar_rep *> (rep); + l = r->get_data (); + u = l + 1; + res = true; + } + break; + default: + break; + } + + return res; +} + +idx_vector +idx_vector::complement (octave_idx_type n) const +{ + + bool *left = new bool[n]; + + std::fill (left, left + n, true); + + octave_idx_type cnt = n; - frozen_at_z_len = z_len ? z_len : len; + for (octave_idx_type i = 0, len = length (); i < len; i++) + { + octave_idx_type k = xelem (i); + if (k < n && left[k]) + { + left[k] = false; + cnt--; + } + } + + octave_idx_type len = cnt, *data = new octave_idx_type[len]; + for (octave_idx_type i = 0, j = 0; i < n; i++) + if (left[i]) data[j++] = i; + + return new idx_vector_rep (data, len, data[len-1], dim_vector (1, len), DIRECT); +} + +octave_idx_type +idx_vector::freeze (octave_idx_type z_len, const char *tag, bool resize_ok) +{ + if (! resize_ok && extent (z_len) > z_len) + { + (*current_liboctave_error_handler) + ("invalid matrix index = %d", extent (z_len)); + rep->err = true; + chkerr (); + } - return frozen_len; + return length (z_len); +} + +octave_idx_type +idx_vector::ones_count () const +{ + octave_idx_type n = 0; + if (is_colon ()) + n = 1; + else + for (octave_idx_type i = 0; i < length (1); i++) + if (xelem (i) == 0) n++; + return n; } +// Instantiate the octave_int constructors we want. +#define INSTANTIATE_SCALAR_VECTOR_REP_CONST(T) \ + template idx_vector::idx_scalar_rep::idx_scalar_rep (T); \ + template idx_vector::idx_vector_rep::idx_vector_rep (const Array<T>&); + +INSTANTIATE_SCALAR_VECTOR_REP_CONST (float) +INSTANTIATE_SCALAR_VECTOR_REP_CONST (double) +INSTANTIATE_SCALAR_VECTOR_REP_CONST (octave_int8) +INSTANTIATE_SCALAR_VECTOR_REP_CONST (octave_int16) +INSTANTIATE_SCALAR_VECTOR_REP_CONST (octave_int32) +INSTANTIATE_SCALAR_VECTOR_REP_CONST (octave_int64) +INSTANTIATE_SCALAR_VECTOR_REP_CONST (octave_uint8) +INSTANTIATE_SCALAR_VECTOR_REP_CONST (octave_uint16) +INSTANTIATE_SCALAR_VECTOR_REP_CONST (octave_uint32) +INSTANTIATE_SCALAR_VECTOR_REP_CONST (octave_uint64) + /* ;;; Local Variables: *** ;;; mode: C++ ***
--- a/liboctave/idx-vector.h +++ b/liboctave/idx-vector.h @@ -2,6 +2,7 @@ Copyright (C) 1993, 1994, 1995, 1996, 1997, 1998, 2000, 2002, 2003, 2004, 2005, 2006, 2007 John W. Eaton +Copyright (C) 2008 Jaroslav Hajek This file is part of Octave. @@ -25,199 +26,382 @@ #define octave_idx_vector_h 1 #include <iostream> +#include <algorithm> +#include <cassert> #include "dim-vector.h" #include "oct-inttypes.h" -#include "intNDArray.h" + +template<class T> class Array; +class Range; -class ColumnVector; -class boolNDArray; -class NDArray; -class Range; +// Design rationale: +// idx_vector is a reference-counting, polymorphic pointer, that can contain +// 4 types of index objects: a magic colon, a range, a scalar, or an index vector. +// Polymorphic methods for single element access are provided, as well as +// templates implementing "early dispatch", i.e. hoisting the checks for index +// type out of loops. class OCTAVE_API idx_vector { +public: + + enum idx_class_type + { + class_invalid = -1, + class_colon = 0, + class_range, + class_scalar, + class_vector + }; + private: - class - OCTAVE_API - idx_vector_rep + class idx_base_rep + { + public: + idx_base_rep (void) : count (1), err (false) { } + + virtual ~idx_base_rep (void) { } + + // Non-range-checking element query. + virtual octave_idx_type xelem (octave_idx_type i) const = 0; + + // Range-checking element query. + virtual octave_idx_type checkelem (octave_idx_type i) const = 0; + + // Length of the index vector. + virtual octave_idx_type length (octave_idx_type n) const = 0; + + // The maximum index + 1. The actual dimension is passed in. + virtual octave_idx_type extent (octave_idx_type n) const = 0; + + // Index class. + virtual idx_class_type idx_class () const { return class_invalid; } + + // Sorts, maybe uniqifies, and returns a clone object pointer. + virtual idx_base_rep *sort_uniq_clone (bool uniq = false) = 0; + + // Checks whether the index is colon or a range equivalent to colon. + virtual bool is_colon_equiv (octave_idx_type) const + { return false; } + + // The original dimensions of this object (used when subscribing by matrices). + virtual dim_vector orig_dimensions () const + { return dim_vector (); } + + // i/o + virtual std::ostream& print (std::ostream& os) const = 0; + + int count; + + bool err; + + private: + + // No copying! + idx_base_rep (const idx_base_rep&); + }; + + // The magic colon index. + class idx_colon_rep : public idx_base_rep { public: + idx_colon_rep (void) { } - idx_vector_rep (void) - : data (0), len (0), num_zeros (0), num_ones (0), - range_base (0), range_step (0), max_val (0), - min_val (0), count (1), frozen_at_z_len (0), frozen_len (0), - colon (0), range(0), initialized (0), frozen (0), - colon_equiv_checked (0), colon_equiv (0), orig_dims () { } + idx_colon_rep (char c); + + octave_idx_type xelem (octave_idx_type i) const + { return i; } + + octave_idx_type checkelem (octave_idx_type i) const; - idx_vector_rep (const ColumnVector& v); + octave_idx_type length (octave_idx_type n) const + { return n; } - idx_vector_rep (const NDArray& nda); + octave_idx_type extent (octave_idx_type n) const + { return n; } + + idx_class_type idx_class () const { return class_colon; } - template <class U> - idx_vector_rep (const intNDArray<U>& inda) - : data (0), len (inda.length ()), num_zeros (0), num_ones (0), - range_base (0), range_step (0), max_val (0), min_val (0), - count (1), frozen_at_z_len (0), frozen_len (0), colon (0), - range(0), initialized (0), frozen (0), - colon_equiv_checked (0), colon_equiv (0), orig_dims (inda.dims ()) - { - if (len == 0) - { - initialized = 1; - return; - } - else - { - data = new octave_idx_type [len]; + idx_base_rep *sort_uniq_clone (bool = false) + { count++; return this; } + + bool is_colon_equiv (octave_idx_type) const + { return true; } + + std::ostream& print (std::ostream& os) const; - for (octave_idx_type i = 0; i < len; i++) - data[i] = tree_to_mat_idx (inda.elem (i)); - } - - init_state (); - } + private: - idx_vector_rep (const Range& r); - - idx_vector_rep (double d); + // No copying! + idx_colon_rep (const idx_colon_rep& idx); + }; - idx_vector_rep (octave_idx_type i); - - idx_vector_rep (char c); - - idx_vector_rep (bool b); + // To distinguish the "direct" constructors that blindly trust the data. + enum direct { DIRECT }; - template <class U> - idx_vector_rep (const octave_int<U>& i) - : data (0), len (1), num_zeros (0), num_ones (0), - range_base (0), range_step (0), max_val (0), min_val (0), - count (1), frozen_at_z_len (0), frozen_len (0), colon (0), - range(0), initialized (0), frozen (0), - colon_equiv_checked (0), colon_equiv (0), orig_dims (1, 1) - { - data = new octave_idx_type [len]; - - data[0] = tree_to_mat_idx (i); + // The integer range index. + class idx_range_rep : public idx_base_rep + { + public: + idx_range_rep (octave_idx_type _start, octave_idx_type _len, + octave_idx_type _step, direct) + : idx_base_rep (), start(_start), len(_len), step(_step) { } - init_state (); - } - - idx_vector_rep (const boolNDArray& bnda); - - idx_vector_rep (const idx_vector_rep& a); + idx_range_rep (void) + : start(0), len(1), step(1) { } - ~idx_vector_rep (void) { delete [] data; } - - idx_vector_rep& operator = (const idx_vector_rep& a); + // Zero-based constructor. + idx_range_rep (octave_idx_type _start, octave_idx_type _limit, + octave_idx_type _step); - int ok (void) { return initialized; } + idx_range_rep (const Range&); - octave_idx_type capacity (void) const { return len; } - octave_idx_type length (octave_idx_type colon_len) const { return colon ? colon_len : len; } + octave_idx_type xelem (octave_idx_type i) const + { return start + i * step; } - octave_idx_type elem (octave_idx_type n) const - { - return colon ? n : (range ? range_base + range_step*n : data[n]); - } + octave_idx_type checkelem (octave_idx_type i) const; + + octave_idx_type length (octave_idx_type) const + { return len; } - octave_idx_type checkelem (octave_idx_type n) const; - octave_idx_type operator () (octave_idx_type n) const { return checkelem (n); } + octave_idx_type extent (octave_idx_type n) const + { return std::max (n, (start + 1 + (step < 0 ? 0 : step * (len - 1)))); } - octave_idx_type max (void) const { return max_val; } - octave_idx_type min (void) const { return min_val; } - - octave_idx_type zeros_count (void) const { return num_zeros; } - octave_idx_type ones_count (void) const { return num_ones; } + idx_class_type idx_class () const { return class_range; } - int is_colon (void) const { return colon; } - int is_colon_equiv (octave_idx_type n, int sort_uniq); + idx_base_rep *sort_uniq_clone (bool uniq = false); - void sort (bool uniq); - - octave_idx_type orig_rows (void) const { return orig_dims(0); } - octave_idx_type orig_columns (void) const { return orig_dims(1); } + bool is_colon_equiv (octave_idx_type n) const + { return start == 0 && step == 1 && len == n; } - dim_vector orig_dimensions (void) const { return orig_dims; } - - // other stuff + dim_vector orig_dimensions () const + { return dim_vector (1, len); } - void shorten (octave_idx_type n); // Unsafe. Avoid at all cost. + octave_idx_type get_start () const { return start; } - octave_idx_type freeze (octave_idx_type z_len, const char *tag, bool resize_ok); - - // i/o + octave_idx_type get_step () const { return step; } std::ostream& print (std::ostream& os) const; - octave_idx_type *data; - octave_idx_type len; - octave_idx_type num_zeros; - octave_idx_type num_ones; - octave_idx_type range_base; - octave_idx_type range_step; - octave_idx_type max_val; - octave_idx_type min_val; + private: + + // No copying! + idx_range_rep (const idx_range_rep& idx); + + octave_idx_type start, len, step; + + }; + + // The integer scalar index. + class idx_scalar_rep : public idx_base_rep + { + public: + idx_scalar_rep (octave_idx_type i, direct) + : data (i) { } + + idx_scalar_rep (void) + : data (0) { } + + // Zero-based constructor. + idx_scalar_rep (octave_idx_type i); + + template <class T> + idx_scalar_rep (T x); + + octave_idx_type xelem (octave_idx_type) const + { return data; } + + octave_idx_type checkelem (octave_idx_type i) const; + + octave_idx_type length (octave_idx_type) const + { return 1; } + + octave_idx_type extent (octave_idx_type n) const + { return std::max (n, data + 1); } + + idx_class_type idx_class () const { return class_scalar; } + + idx_base_rep *sort_uniq_clone (bool = false) + { count++; return this; } + + bool is_colon_equiv (octave_idx_type n) const + { return n == 1 && data == 0; } + + dim_vector orig_dimensions () const + { return dim_vector (1, 1); } + + octave_idx_type get_data () const { return data; } + + std::ostream& print (std::ostream& os) const; + + private: + + // No copying! + idx_scalar_rep (const idx_scalar_rep& idx); + + octave_idx_type data; + + }; - int count; + // The integer vector index. + class idx_vector_rep : public idx_base_rep + { + public: + // Direct constructor. + idx_vector_rep (octave_idx_type *_data, octave_idx_type _len, + octave_idx_type _ext, const dim_vector& od, direct) + : data (_data), len (_len), ext (_ext), aowner (0), orig_dims (od) { } + + idx_vector_rep (void) + : data (0), len (0), aowner (0) + { } - octave_idx_type frozen_at_z_len; - octave_idx_type frozen_len; + // Zero-based constructor. + idx_vector_rep (const Array<octave_idx_type>& inda); + + template <class T> + idx_vector_rep (const Array<T>&); + + idx_vector_rep (bool); + + idx_vector_rep (const Array<bool>&); + + ~idx_vector_rep (void); + + octave_idx_type xelem (octave_idx_type i) const + { return data[i]; } + + octave_idx_type checkelem (octave_idx_type i) const; - unsigned int colon : 1; - unsigned int range : 1; - unsigned int initialized : 1; - unsigned int frozen : 1; - unsigned int colon_equiv_checked : 1; - unsigned int colon_equiv : 1; + octave_idx_type length (octave_idx_type) const + { return len; } + + octave_idx_type extent (octave_idx_type n) const + { return std::max (n, ext); } + + idx_class_type idx_class () const { return class_vector; } + + idx_base_rep *sort_uniq_clone (bool uniq = false); + + dim_vector orig_dimensions () const + { return orig_dims; } + + const octave_idx_type *get_data () const { return data; } + + std::ostream& print (std::ostream& os) const; + + private: + + // No copying! + idx_vector_rep (const idx_vector_rep& idx); + + const octave_idx_type *data; + octave_idx_type len, ext; + + // This is a trick to allow user-given zero-based arrays to be used as indices + // without copying. If the following pointer is nonzero, we do not own the data, + // but rather have an Array<octave_idx_type> object that provides us the data. + // Note that we need a pointer because we deferred the Array<T> declaration and + // we do not want it yet to be defined. + + Array<octave_idx_type> *aowner; dim_vector orig_dims; - - void init_state (void); + }; + + + idx_vector (idx_base_rep *r) : rep (r) { } - octave_idx_type tree_to_mat_idx (double x, bool& conversion_error); + // The shared empty vector representation (for fast default constructor) + static idx_vector_rep *nil_rep () + { + static idx_vector_rep ivr; + return &ivr; + } - octave_idx_type tree_to_mat_idx (octave_idx_type i) { return i - 1; } + // The shared empty vector representation with the error flag set. + static idx_vector_rep *err_rep () + { + static idx_vector_rep ivr; + ivr.err = true; + return &ivr; + } - template <class U> octave_idx_type tree_to_mat_idx (const octave_int<U>& i) - { return i.value () - 1; } - }; + // If there was an error in constructing the rep, replace it with empty vector + // for safety. + void chkerr (void) + { + if (rep->err) + { + if (--rep->count == 0) + delete rep; + rep = err_rep (); + rep->count++; + } + } public: - idx_vector (void) : rep (new idx_vector_rep ()) { } + // Fast empty constructor. + idx_vector (void) : rep (nil_rep ()) { rep->count++; } - idx_vector (const ColumnVector& v) : rep (new idx_vector_rep (v)) { } + // Zero-based constructors (for use from C++). + idx_vector (octave_idx_type i) : rep (new idx_scalar_rep (i)) + { chkerr (); } + + idx_vector (octave_idx_type start, octave_idx_type limit, + octave_idx_type step = 1) + : rep (new idx_range_rep (start, limit, step)) + { chkerr (); } - idx_vector (const NDArray& nda) : rep (new idx_vector_rep (nda)) { } + idx_vector (const Array<octave_idx_type>& inda) + : rep (new idx_vector_rep (inda)) + { chkerr (); } - template <class U> - idx_vector (const intNDArray<U>& inda) : rep (new idx_vector_rep (inda)) { } + // Colon is best constructed by simply copying (or referencing) this member. + static const idx_vector colon; - idx_vector (const Range& r) : rep (new idx_vector_rep (r)) { } + // or passing ':' here + idx_vector (char c) : rep (new idx_colon_rep (c)) { chkerr (); } + + // Conversion constructors (used by interpreter). - idx_vector (double d) : rep (new idx_vector_rep (d)) { } + template <class T> + idx_vector (octave_int<T> x) : rep (new idx_scalar_rep (x)) { chkerr (); } + + idx_vector (double x) : rep (new idx_scalar_rep (x)) { chkerr (); } - idx_vector (octave_idx_type i) : rep (new idx_vector_rep (i)) { } + idx_vector (float x) : rep (new idx_scalar_rep (x)) { chkerr (); } - idx_vector (char c) : rep (new idx_vector_rep (c)) { } + // A scalar bool does not necessarily map to scalar index. + idx_vector (bool x) : rep (new idx_vector_rep (x)) { chkerr (); } - idx_vector (bool b) : rep (new idx_vector_rep (b)) { } + template <class T> + idx_vector (const Array<octave_int<T> >& nda) : rep (new idx_vector_rep (nda)) + { chkerr (); } + + idx_vector (const Array<double>& nda) : rep (new idx_vector_rep (nda)) + { chkerr (); } - template <class U> - idx_vector (const octave_int<U>& i) : rep (new idx_vector_rep (i)) { } + idx_vector (const Array<float>& nda) : rep (new idx_vector_rep (nda)) + { chkerr (); } - idx_vector (const boolNDArray& bnda) : rep (new idx_vector_rep (bnda)) { } + idx_vector (const Array<bool>& nda) : rep (new idx_vector_rep (nda)) + { chkerr (); } + + idx_vector (const Range& r) + : rep (new idx_range_rep (r)) + { chkerr (); } idx_vector (const idx_vector& a) : rep (a.rep) { rep->count++; } ~idx_vector (void) { - if (--rep->count <= 0) + if (--rep->count == 0) delete rep; } @@ -225,7 +409,7 @@ { if (this != &a) { - if (--rep->count <= 0) + if (--rep->count == 0) delete rep; rep = a.rep; @@ -234,54 +418,329 @@ return *this; } - operator bool () const { return rep->ok (); } + idx_class_type idx_class (void) const { return rep->idx_class (); } - octave_idx_type capacity (void) const { return rep->capacity (); } - octave_idx_type length (octave_idx_type cl) const { return rep->length (cl); } + octave_idx_type length (octave_idx_type n = 0) const + { return rep->length (n); } - octave_idx_type elem (octave_idx_type n) const { return rep->elem (n); } - octave_idx_type checkelem (octave_idx_type n) const { return rep->checkelem (n); } - octave_idx_type operator () (octave_idx_type n) const { return rep->operator () (n); } + octave_idx_type extent (octave_idx_type n) const + { return rep->extent (n); } - octave_idx_type max (void) const { return rep->max (); } - octave_idx_type min (void) const { return rep->min (); } + octave_idx_type xelem (octave_idx_type n) const + { return rep->xelem (n); } + + octave_idx_type checkelem (octave_idx_type n) const + { return rep->checkelem (n); } - int one_zero_only (void) const { return 0; } - octave_idx_type zeros_count (void) const { return rep->zeros_count (); } - octave_idx_type ones_count (void) const { return rep->ones_count (); } + octave_idx_type operator () (octave_idx_type n) const + { +#if defined (BOUNDS_CHECKING) + return rep->checkelem (n); +#else + return rep->xelem (n); +#endif + } + + operator bool () const + { return ! rep->err; } - int is_colon (void) const { return rep->is_colon (); } - int is_colon_equiv (octave_idx_type n, int sort_uniq = 0) const - { return rep->is_colon_equiv (n, sort_uniq); } + bool is_colon (void) const + { return rep->idx_class () == class_colon; } + + bool is_scalar (void) const + { return rep->idx_class () == class_scalar; } - void sort (bool uniq = false) { rep->sort (uniq); } + bool is_colon_equiv (octave_idx_type n) const + { return rep->is_colon_equiv (n); } - octave_idx_type orig_rows (void) const { return rep->orig_rows (); } - octave_idx_type orig_columns (void) const { return rep->orig_columns (); } + idx_vector sorted (bool uniq = false) const + { return idx_vector (rep->sort_uniq_clone (uniq)); } dim_vector orig_dimensions (void) const { return rep->orig_dimensions (); } + octave_idx_type orig_rows () const + { return orig_dimensions () (0); } + + octave_idx_type orig_columns () const + { return orig_dimensions () (1); } + int orig_empty (void) const { return (! is_colon () && orig_dimensions().any_zero ()); } - // Unsafe. Avoid at all cost. - void shorten (octave_idx_type n) { rep->shorten (n); } - // i/o - octave_idx_type freeze (octave_idx_type z_len, const char *tag, bool resize_ok = false) - { return rep->freeze (z_len, tag, resize_ok); } - std::ostream& print (std::ostream& os) const { return rep->print (os); } friend std::ostream& operator << (std::ostream& os, const idx_vector& a) { return a.print (os); } + // Slice with specializations. No checking of bounds! + // + // This is equivalent to the following loop (but much faster): + // + // for (octave_idx_type i = 0; i < idx->length (n); i++) + // dest[i] = src[idx(i)]; + // return i; + // + template <class T> + octave_idx_type + index (const T *src, octave_idx_type n, T *dest) const + { + octave_idx_type len = rep->length (n); + switch (rep->idx_class ()) + { + case class_colon: + std::copy (src, src + len, dest); + break; + case class_range: + { + idx_range_rep * r = dynamic_cast<idx_range_rep *> (rep); + octave_idx_type start = r->get_start (), step = r->get_step (); + const T *ssrc = src + start; + if (step == 1) + std::copy (ssrc, ssrc + len, dest); + else if (step == -1) + std::reverse_copy (ssrc - len + 1, ssrc + 1, dest); + else + { + for (octave_idx_type i = 0, j = 0; i < len; i++, j += step) + dest[i] = ssrc[j]; + } + } + break; + case class_scalar: + { + idx_scalar_rep * r = dynamic_cast<idx_scalar_rep *> (rep); + dest[0] = src[r->get_data ()]; + } + break; + case class_vector: + { + idx_vector_rep * r = dynamic_cast<idx_vector_rep *> (rep); + const octave_idx_type *data = r->get_data (); + for (octave_idx_type i = 0; i < len; i++) + dest[i] = src[data[i]]; + } + break; + default: + assert (false); + break; + } + + return len; + } + + // Slice assignment with specializations. No checking of bounds! + // + // This is equivalent to the following loop (but much faster): + // + // for (octave_idx_type i = 0; i < idx->length (n); i++) + // dest[idx(i)] = src[i]; + // return i; + // + template <class T> + octave_idx_type + assign (const T *src, octave_idx_type n, T *dest) const + { + octave_idx_type len = rep->length (n); + switch (rep->idx_class ()) + { + case class_colon: + std::copy (src, src + len, dest); + break; + case class_range: + { + idx_range_rep * r = dynamic_cast<idx_range_rep *> (rep); + octave_idx_type start = r->get_start (), step = r->get_step (); + T *sdest = dest + start; + if (step == 1) + std::copy (src, src + len, sdest); + else if (step == -1) + std::reverse_copy (src, src + len, sdest - len + 1); + else + { + for (octave_idx_type i = 0, j = 0; i < len; i++, j += step) + sdest[j] = src[i]; + } + } + break; + case class_scalar: + { + idx_scalar_rep * r = dynamic_cast<idx_scalar_rep *> (rep); + dest[r->get_data ()] = src[0]; + } + break; + case class_vector: + { + idx_vector_rep * r = dynamic_cast<idx_vector_rep *> (rep); + const octave_idx_type *data = r->get_data (); + for (octave_idx_type i = 0; i < len; i++) + dest[data[i]] = src[i]; + } + break; + default: + assert (false); + break; + } + + return len; + } + + // Slice fill with specializations. No checking of bounds! + // + // This is equivalent to the following loop (but much faster): + // + // for (octave_idx_type i = 0; i < idx->length (n); i++) + // dest[idx(i)] = val; + // return i; + // + template <class T> + octave_idx_type + fill (const T& val, octave_idx_type n, T *dest) const + { + octave_idx_type len = rep->length (n); + switch (rep->idx_class ()) + { + case class_colon: + std::fill (dest, dest + len, val); + break; + case class_range: + { + idx_range_rep * r = dynamic_cast<idx_range_rep *> (rep); + octave_idx_type start = r->get_start (), step = r->get_step (); + T *sdest = dest + start; + if (step == 1) + std::fill (sdest, sdest + len, val); + else if (step == -1) + std::fill (sdest - len + 1, sdest + 1, val); + else + { + for (octave_idx_type i = 0, j = 0; i < len; i++, j += step) + sdest[j] = val; + } + } + break; + case class_scalar: + { + idx_scalar_rep * r = dynamic_cast<idx_scalar_rep *> (rep); + dest[r->get_data ()] = val; + } + break; + case class_vector: + { + idx_vector_rep * r = dynamic_cast<idx_vector_rep *> (rep); + const octave_idx_type *data = r->get_data (); + for (octave_idx_type i = 0; i < len; i++) + dest[data[i]] = val; + } + break; + default: + assert (false); + break; + } + + return len; + } + + // Generic breakable indexed loop. The loop body should be encapsulated in a + // single functor body. + // This is equivalent to the following loop (but faster, at least for simple + // inlined bodies): + // + // for (octave_idx_type i = 0; i < idx->length (n); i++) + // if (body (i, idx(i))) break; + // return i; + // + + template <class Functor> + octave_idx_type + loop (octave_idx_type n, const Functor& body) const + { + octave_idx_type len = rep->length (n), ret; + switch (rep->idx_class ()) + { + case class_colon: + { + octave_idx_type i; + for (i = 0; i < len && body (i); i++) ; + ret = i; + } + break; + case class_range: + { + idx_range_rep * r = dynamic_cast<idx_range_rep *> (rep); + octave_idx_type start = r->get_start (), step = r->get_step (); + octave_idx_type i, j; + if (step == 1) + for (i = start, j = start + len; i < j && body (i); i++) ; + else if (step == -1) + for (i = start, j = start - len; i > j && body (i); i--) ; + else + for (i = 0, j = start; i < len && body (j); i++, j += step) ; + ret = i; + } + break; + case class_scalar: + { + idx_scalar_rep * r = dynamic_cast<idx_scalar_rep *> (rep); + ret = body (r->get_data ()) ? 1 : 0; + } + break; + case class_vector: + { + idx_vector_rep * r = dynamic_cast<idx_vector_rep *> (rep); + octave_idx_type *data = r->get_data (); + octave_idx_type i; + for (i = 0; i < len && body (data[i]); i++) ; + ret = i; + } + break; + default: + assert (false); + break; + } + + return ret; + } + + // Rationale: + // This method is the key to "smart indexing". When indexing cartesian + // arrays, sometimes consecutive index vectors can be reduced into a single + // index. If rows (A) = k and i.maybe_reduce (j) gives k, then A(i,j)(:) is + // equal to A(k)(:). + + // If the next index can be reduced, returns true and updates this. + bool maybe_reduce (octave_idx_type n, const idx_vector& j, + octave_idx_type nj); + + bool is_cont_range (octave_idx_type n, + octave_idx_type& l, octave_idx_type& u) const; + + idx_vector + complement (octave_idx_type n) const; + + // FIXME: These are here for compatibility. They should be removed when no + // longer in use. + + octave_idx_type elem (octave_idx_type n) const + { return (*this) (n); } + + bool is_colon_equiv (octave_idx_type n, int) const + { return is_colon_equiv (n); } + + octave_idx_type freeze (octave_idx_type z_len, const char *tag, bool resize_ok = false); + + void sort (bool uniq = false) + { *this = sorted (uniq); } + + octave_idx_type ones_count () const; + + octave_idx_type max () const { return extent (1) - 1; } + private: - idx_vector_rep *rep; + idx_base_rep *rep; - void init_state (void) { rep->init_state (); } }; #endif
--- a/scripts/polynomial/polyval.m +++ b/scripts/polynomial/polyval.m @@ -71,7 +71,7 @@ k = numel (x); x = (x - mu(1)) / mu(2); A = (x(:) * ones (1, n+1)) .^ (ones (k, 1) * (n:-1:0)); - y(:) = A * p(:); + y = A * p(:); y = reshape (y, size (x)); if (nargout == 2)
--- a/src/Cell.cc +++ b/src/Cell.cc @@ -164,10 +164,14 @@ const octave_value& fill_val) { - for (octave_idx_type i = 0; i < idx_arg.length (); i++) - set_index (idx_arg(i).index_vector ()); + octave_idx_type len = idx_arg.length (); + + Array<idx_vector> ra_idx (len); - ::assign (*this, rhs, fill_val); + for (octave_idx_type i = 0; i < len; i++) + ra_idx(i) = idx_arg(i).index_vector (); + + Array<octave_value>::assign (ra_idx, rhs, fill_val); return *this; } @@ -176,11 +180,14 @@ Cell::delete_elements (const octave_value_list& idx_arg) { - Array<idx_vector> ra_idx (idx_arg.length ()); - for (octave_idx_type i = 0; i < idx_arg.length (); i++) + octave_idx_type len = idx_arg.length (); + + Array<idx_vector> ra_idx (len); + + for (octave_idx_type i = 0; i < len; i++) ra_idx.xelem (i) = idx_arg(i).index_vector (); - maybe_delete_elements (ra_idx); + Array<octave_value>::delete_elements (ra_idx); return *this; }
--- a/src/Cell.h +++ b/src/Cell.h @@ -89,11 +89,22 @@ const octave_value& rfv = resize_fill_value ()) const { return Cell (ArrayN<octave_value>::index (ra_idx, resize_ok, rfv)); } + // FIXME: This seems necessary for octave_base_mat<Cell>::delete_elements + // to work, but I don't understand why. + void delete_elements (const Array<idx_vector>& ia) + { ArrayN<octave_value>::delete_elements (ia); } + Cell& delete_elements (const octave_value_list& idx); Cell& assign (const octave_value_list& idx, const Cell& rhs, const octave_value& fill_val = octave_value ()); + // FIXME: This seems necessary for octave_base_mat<Cell>::assign + // to work, but I don't understand why. + void assign (const Array<idx_vector>& ia, const Array<octave_value>& rhs, + const octave_value& fill_val = octave_value ()) + { ArrayN<octave_value>::assign (ia, rhs, fill_val); } + Cell reshape (const dim_vector& new_dims) const { return ArrayN<octave_value>::reshape (new_dims); }
--- a/src/ChangeLog +++ b/src/ChangeLog @@ -1,3 +1,23 @@ +2008-10-28 Jaroslav Hajek <highegg@gmail.com> + + * Cell.h (Cell::assign (const Array<idx_vector>&, ...), + Cell::delete_elements (const Array<idx_vector>&, ...)): + New member functions. + * Cell.cc (Cell::assign (const octave_value_list&, ...), + Cell::delete_elements (const octave_value_list&, ...)): + Call Array<T>::assign. + * DLD-FUNCTIONS/dlmread.cc: Call Array<T>::resize_fill. + * ov-base-mat.cc (octave_base_matrix::assign): Call Array<T>::assign. + (octave_base_matrix::delete_elements):: Call Array<T>::delete_elements. + * ov-cell.cc (Fcell): Call Array<T>::chop_trailing_singletons, + simplify. + * ov-cx-mat.cc (octave_complex_matrix::assign): Call Array<T>::assign. + * ov-flt-cx-mat.cc (octave_float_complex_matrix::assign): Call + Array<T>::assign. + * ov-list.cc (octave_list::subsasgn): Call Array<T>::assign. + * pr-output.cc (PRINT_ND_ARRAY): Use zero-based indices. + + 2008-10-29 Thorsten Meyer <thorsten.meyier@gmx.de> * data.cc (Fcolumns): Remove "and" from @seealso string.
--- a/src/DLD-FUNCTIONS/dlmread.cc +++ b/src/DLD-FUNCTIONS/dlmread.cc @@ -301,9 +301,9 @@ // Use resize_and_fill for the case of not-equal // length rows. if (iscmplx) - cdata.resize_and_fill (r, c, 0); + cdata.resize_fill (r, c, 0); else - rdata.resize_and_fill (r, c, 0); + rdata.resize_fill (r, c, 0); rmax = r; cmax = c; }
--- a/src/ov-base-mat.cc +++ b/src/ov-base-mat.cc @@ -195,11 +195,12 @@ { octave_idx_type len = idx.length (); - for (octave_idx_type i = 0; i < len; i++) - matrix.set_index (idx(i).index_vector ()); + Array<idx_vector> ra_idx (len); - ::assign (matrix, rhs, MT::resize_fill_value ()); + for (octave_idx_type i = 0; i < len; i++) + ra_idx(i) = idx(i).index_vector (); + matrix.assign (ra_idx, rhs, MT::resize_fill_value ()); // Invalidate the matrix type typ.invalidate_type (); @@ -216,7 +217,7 @@ for (octave_idx_type i = 0; i < len; i++) ra_idx(i) = idx(i).index_vector (); - matrix.maybe_delete_elements (ra_idx); + matrix.delete_elements (ra_idx); // Invalidate the matrix type typ.invalidate_type ();
--- a/src/ov-cell.cc +++ b/src/ov-cell.cc @@ -1017,23 +1017,12 @@ if (! error_state) { - int ndim = dims.length (); + dims.chop_trailing_singletons (); check_dimensions (dims, "cell"); if (! error_state) - { - switch (ndim) - { - case 1: - retval = Cell (dims(0), dims(0), Matrix ()); - break; - - default: - retval = Cell (dims, Matrix ()); - break; - } - } + retval = Cell (dims, Matrix ()); } return retval;
--- a/src/ov-cx-mat.cc +++ b/src/ov-cx-mat.cc @@ -122,10 +122,12 @@ { octave_idx_type len = idx.length (); + Array<idx_vector> ra_idx (len); + for (octave_idx_type i = 0; i < len; i++) - matrix.set_index (idx(i).index_vector ()); + ra_idx(i) = idx(i).index_vector (); - ::assign (matrix, rhs); + matrix.assign (ra_idx, rhs); } bool
--- a/src/ov-flt-cx-mat.cc +++ b/src/ov-flt-cx-mat.cc @@ -108,14 +108,16 @@ void octave_float_complex_matrix::assign (const octave_value_list& idx, - const FloatNDArray& rhs) + const FloatNDArray& rhs) { octave_idx_type len = idx.length (); + Array<idx_vector> ra_idx (len); + for (octave_idx_type i = 0; i < len; i++) - matrix.set_index (idx(i).index_vector ()); + ra_idx(i) = idx(i).index_vector (); - ::assign (matrix, rhs); + matrix.assign (ra_idx, rhs); } bool
--- a/src/ov-list.cc +++ b/src/ov-list.cc @@ -206,12 +206,7 @@ { octave_value_list i = idx.front (); - octave_idx_type len = i.length (); - - for (octave_idx_type k = 0; k < len; k++) - data.set_index (i(k).index_vector ()); - - ::assign (data, Cell (t_rhs), Cell::resize_fill_value ()); + data.assign (i, Cell (t_rhs), Cell::resize_fill_value ()); count++; retval = octave_value (this);
--- a/src/pr-output.cc +++ b/src/pr-output.cc @@ -1747,7 +1747,7 @@ idx(1) = idx_vector (':'); \ \ for (int k = 2; k < ndims; k++) \ - idx(k) = idx_vector (ra_idx(k) + 1); \ + idx(k) = idx_vector (ra_idx(k)); \ \ octave_value page \ = MAT_T (Array2<ELT_T> (nda.index (idx), nr, nc)); \ @@ -2291,7 +2291,7 @@ idx(1) = idx_vector (':'); for (int k = 2; k < ndims; k++) - idx(k) = idx_vector (ra_idx(k) + 1); + idx(k) = idx_vector (ra_idx(k)); Array2<std::string> page (nda.index (idx), nr, nc); @@ -2558,7 +2558,7 @@ idx(1) = idx_vector (':'); for (int k = 2; k < ndims; k++) - idx(k) = idx_vector (ra_idx(k) + 1); + idx(k) = idx_vector (ra_idx(k)); Array2<T> page (nda.index (idx), nr, nc); @@ -2663,7 +2663,7 @@ idx(1) = idx_vector (':'); for (int k = 2; k < ndims; k++) - idx(k) = idx_vector (ra_idx(k) + 1); + idx(k) = idx_vector (ra_idx(k)); Array2<T> page (nda.index (idx), nr, nc);
--- a/test/ChangeLog +++ b/test/ChangeLog @@ -1,3 +1,9 @@ +2008-10-28 Jaroslav Hajek <highegg@gmail.com> + + * test_logical-wfi-f.m: Fix error messages. + * test_logical-wfi-t.m: Fix error messages. + * test_slice.m: Fix error messages. + 2008-09-26 Jaroslav Hajek <highegg@gmail.com> * test_null_assign.m: More test for null assignments.
--- a/test/test_logical-wfi-f.m +++ b/test/test_logical-wfi-f.m @@ -48,7 +48,7 @@ %! warning ("off", "Octave:fortran-indexing"); %!shared a %! a = 2; -%!error <invalid vector index> a(logical ([1,1])); +%!error <A\(I\): Index exceeds matrix dimension\.> a(logical ([1,1])); %! warning ("wfi.state", "Octave:fortran-indexing"); %% test/octave.test/logical-wfi-f/v-1.m
--- a/test/test_logical-wfi-t.m +++ b/test/test_logical-wfi-t.m @@ -48,7 +48,7 @@ %! warning ("on", "Octave:fortran-indexing"); %!shared a %! a = 2; -%!error <invalid vector index.*> a(logical ([1,1])); +%!error <A\(I\): Index exceeds matrix dimension\.> a(logical ([1,1])); %! warning ("wfi.state", "Octave:fortran-indexing"); %% test/octave.test/logical-wfi-t/v-1.m
--- a/test/test_slice.m +++ b/test/test_slice.m @@ -38,9 +38,9 @@ ## size = [2 0] %!assert(set_slice([2 0], 11, []), zeros([2 0])); -%!error <A\(I\) = X: unable to resize A> set_slice([2 0], 11, 1) -%!error <A\(I\) = X: unable to resize A> set_slice([2 0], 11, 2) -%!error <A\(I\) = X: unable to resize A> set_slice([2 0], 11, 3) +%!error <resize: Invalid.*> set_slice([2 0], 11, 1) +%!error <resize: Invalid.*> set_slice([2 0], 11, 2) +%!error <resize: Invalid.*> set_slice([2 0], 11, 3) %!assert(set_slice([2 0], 21, []), zeros([2 0])); %!assert(set_slice([2 0], 21, 1), zeros([2 0])); %!assert(set_slice([2 0], 21, 2), zeros([2 0])); @@ -152,8 +152,8 @@ %!assert(set_slice([2 2], 11, 2), [1 1;2 1]); %!assert(set_slice([2 2], 11, 3), [1 2;1 1]); %!assert(set_slice([2 2], 11, 4), [1 1;1 2]); -%!error <invalid matrix index = 5> set_slice([2 2], 11, 5) -%!error <invalid matrix index = 6> set_slice([2 2], 11, 6) +%!error <resize: Invalid.*> set_slice([2 2], 11, 5) +%!error <resize: Invalid.*> set_slice([2 2], 11, 6) %!assert(set_slice([2 2], 21, []), ones([2 2])); %!assert(set_slice([2 2], 21, 1), [2 2;1 1]); %!assert(set_slice([2 2], 21, 2), [1 1;2 2]);