Mercurial > hg > octave-lojdl
changeset 9513:9f870f73ab7d
implement built-in diff
author | Jaroslav Hajek <highegg@gmail.com> |
---|---|
date | Tue, 11 Aug 2009 09:08:12 +0200 |
parents | 56e850e3b06f |
children | af86991d8d47 |
files | liboctave/CNDArray.cc liboctave/CNDArray.h liboctave/ChangeLog liboctave/dNDArray.cc liboctave/dNDArray.h liboctave/fCNDArray.cc liboctave/fCNDArray.h liboctave/fNDArray.cc liboctave/fNDArray.h liboctave/intNDArray.cc liboctave/intNDArray.h liboctave/mx-inlines.cc scripts/ChangeLog scripts/general/Makefile.in scripts/general/diff.m src/ChangeLog src/data.cc |
diffstat | 17 files changed, 420 insertions(+), 151 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/CNDArray.cc +++ b/liboctave/CNDArray.cc @@ -672,6 +672,12 @@ } ComplexNDArray +ComplexNDArray::diff (octave_idx_type order, int dim) const +{ + return do_mx_diff_op<ComplexNDArray> (*this, dim, order, mx_inline_diff); +} + +ComplexNDArray ComplexNDArray::concat (const ComplexNDArray& rb, const Array<octave_idx_type>& ra_idx) { if (rb.numel () > 0)
--- a/liboctave/CNDArray.h +++ b/liboctave/CNDArray.h @@ -93,6 +93,8 @@ ComplexNDArray cummin (int dim = 0) const; ComplexNDArray cummin (ArrayN<octave_idx_type>& index, int dim = 0) const; + ComplexNDArray diff (octave_idx_type order = 1, int dim = 0) const; + ComplexNDArray& insert (const NDArray& a, octave_idx_type r, octave_idx_type c); ComplexNDArray& insert (const ComplexNDArray& a, octave_idx_type r, octave_idx_type c); ComplexNDArray& insert (const ComplexNDArray& a, const Array<octave_idx_type>& ra_idx);
--- a/liboctave/ChangeLog +++ b/liboctave/ChangeLog @@ -1,3 +1,20 @@ +2009-08-11 Jaroslav Hajek <highegg@gmail.com> + + * mx-inlines.cc (mx_inline_diff<T>): New overloaded template + function. + (get_extent_triplet): Use dim_vector::first_non_singleton. + (do_mx_diff_op): New template function. + * dNDArray.cc (NDArray::diff): New method. + * dNDArray.h: Declare it. + * fNDArray.cc (FloatNDArray::diff): New method. + * fNDArray.h: Declare it. + * CNDArray.cc (ComplexNDArray::diff): New method. + * CNDArray.h: Declare it. + * fCNDArray.cc (FloatComplexNDArray::diff): New method. + * fCNDArray.h: Declare it. + * intNDArray.cc (intNDArray<T>::diff): New method. + * intNDArray.h: Declare it. + 2009-08-10 Jaroslav Hajek <highegg@gmail.com> * dim-vector.h (dim_vector::first_non_singleton): New method.
--- a/liboctave/dNDArray.cc +++ b/liboctave/dNDArray.cc @@ -787,6 +787,12 @@ } NDArray +NDArray::diff (octave_idx_type order, int dim) const +{ + return do_mx_diff_op<NDArray> (*this, dim, order, mx_inline_diff); +} + +NDArray NDArray::concat (const NDArray& rb, const Array<octave_idx_type>& ra_idx) { if (rb.numel () > 0)
--- a/liboctave/dNDArray.h +++ b/liboctave/dNDArray.h @@ -105,6 +105,8 @@ NDArray cummin (int dim = 0) const; NDArray cummin (ArrayN<octave_idx_type>& index, int dim = 0) const; + NDArray diff (octave_idx_type order = 1, int dim = 0) const; + NDArray& insert (const NDArray& a, octave_idx_type r, octave_idx_type c); NDArray& insert (const NDArray& a, const Array<octave_idx_type>& ra_idx);
--- a/liboctave/fCNDArray.cc +++ b/liboctave/fCNDArray.cc @@ -667,6 +667,12 @@ } FloatComplexNDArray +FloatComplexNDArray::diff (octave_idx_type order, int dim) const +{ + return do_mx_diff_op<FloatComplexNDArray> (*this, dim, order, mx_inline_diff); +} + +FloatComplexNDArray FloatComplexNDArray::concat (const FloatComplexNDArray& rb, const Array<octave_idx_type>& ra_idx) { if (rb.numel () > 0)
--- a/liboctave/fCNDArray.h +++ b/liboctave/fCNDArray.h @@ -93,6 +93,8 @@ FloatComplexNDArray cummin (int dim = 0) const; FloatComplexNDArray cummin (ArrayN<octave_idx_type>& index, int dim = 0) const; + FloatComplexNDArray diff (octave_idx_type order = 1, int dim = 0) const; + FloatComplexNDArray& insert (const NDArray& a, octave_idx_type r, octave_idx_type c); FloatComplexNDArray& insert (const FloatComplexNDArray& a, octave_idx_type r, octave_idx_type c); FloatComplexNDArray& insert (const FloatComplexNDArray& a, const Array<octave_idx_type>& ra_idx);
--- a/liboctave/fNDArray.cc +++ b/liboctave/fNDArray.cc @@ -742,6 +742,12 @@ } FloatNDArray +FloatNDArray::diff (octave_idx_type order, int dim) const +{ + return do_mx_diff_op<FloatNDArray> (*this, dim, order, mx_inline_diff); +} + +FloatNDArray FloatNDArray::concat (const FloatNDArray& rb, const Array<octave_idx_type>& ra_idx) { if (rb.numel () > 0)
--- a/liboctave/fNDArray.h +++ b/liboctave/fNDArray.h @@ -102,6 +102,8 @@ FloatNDArray cummin (int dim = 0) const; FloatNDArray cummin (ArrayN<octave_idx_type>& index, int dim = 0) const; + FloatNDArray diff (octave_idx_type order = 1, int dim = 0) const; + FloatNDArray& insert (const FloatNDArray& a, octave_idx_type r, octave_idx_type c); FloatNDArray& insert (const FloatNDArray& a, const Array<octave_idx_type>& ra_idx);
--- a/liboctave/intNDArray.cc +++ b/liboctave/intNDArray.cc @@ -270,6 +270,13 @@ return do_mx_cumminmax_op<intNDArray<T> > (*this, idx_arg, dim, mx_inline_cummin); } +template <class T> +intNDArray<T> +intNDArray<T>::diff (octave_idx_type order, int dim) const +{ + return do_mx_diff_op<intNDArray<T> > (*this, dim, order, mx_inline_diff); +} + /* ;;; Local Variables: *** ;;; mode: C++ ***
--- a/liboctave/intNDArray.h +++ b/liboctave/intNDArray.h @@ -86,6 +86,8 @@ intNDArray sum (int dim) const; intNDArray cumsum (int dim) const; + intNDArray diff (octave_idx_type order = 1, int dim = 0) const; + intNDArray abs (void) const; intNDArray signum (void) const;
--- a/liboctave/mx-inlines.cc +++ b/liboctave/mx-inlines.cc @@ -844,6 +844,114 @@ OP_CUMMINMAX_FCNN (mx_inline_cummin) OP_CUMMINMAX_FCNN (mx_inline_cummax) +template <class T> +void mx_inline_diff (const T *v, T *r, octave_idx_type n, + octave_idx_type order) +{ + switch (order) + { + case 1: + for (octave_idx_type i = 0; i < n-1; i++) + r[i] = v[i+1] - v[i]; + break; + case 2: + { + T lst; + if (n > 1) + lst = v[1] - v[0]; + for (octave_idx_type i = 0; i < n-2; i++) + { + T dif = v[i+2] - v[i+1]; + r[i] = dif - lst; + lst = dif; + } + } + break; + default: + { + OCTAVE_LOCAL_BUFFER (T, buf, n-1); + + for (octave_idx_type i = 0; i < n-1; i++) + buf[i] = v[i+1] - v[i]; + + for (octave_idx_type o = 2; o <= order; o++) + { + for (octave_idx_type i = 0; i < n-o; i++) + buf[i] = buf[i+1] - buf[i]; + } + + for (octave_idx_type i = 0; i < n-order; i++) + r[i] = buf[i]; + } + } +} + +template <class T> +void mx_inline_diff (const T *v, T *r, + octave_idx_type m, octave_idx_type n, + octave_idx_type order) +{ + switch (order) + { + case 1: + for (octave_idx_type i = 0; i < m*(n-1); i++) + r[i] = v[i+m] - v[i]; + break; + case 2: + for (octave_idx_type i = 0; i < n-2; i++) + { + for (octave_idx_type j = i*m; j < i*m+m; j++) + r[j] = (v[j+m+m] - v[j+m]) + (v[j+m] - v[j]); + } + break; + default: + { + OCTAVE_LOCAL_BUFFER (T, buf, n-1); + + for (octave_idx_type j = 0; j < m; j++) + { + for (octave_idx_type i = 0; i < n-1; i++) + buf[i] = v[i*m+j+m] - v[i*m+j]; + + for (octave_idx_type o = 2; o <= order; o++) + { + for (octave_idx_type i = 0; i < n-o; i++) + buf[i] = buf[i+1] - buf[i]; + } + + for (octave_idx_type i = 0; i < n-order; i++) + r[i*m+j] = buf[i]; + } + } + } +} + +template <class T> +inline void +mx_inline_diff (const T *v, T *r, + octave_idx_type l, octave_idx_type n, octave_idx_type u, + octave_idx_type order) +{ + if (! n) return; + if (l == 1) + { + for (octave_idx_type i = 0; i < u; i++) + { + mx_inline_diff (v, r, n, order); + v += n; r += n-order; + } + } + else + { + for (octave_idx_type i = 0; i < u; i++) + { + mx_inline_diff (v, r, l, n, order); + v += l*n; + r += l*(n-order); + } + } +} + // Assistant function inline void @@ -861,10 +969,8 @@ else { if (dim < 0) - { - // find first non-singleton dim - for (dim = 0; dims(dim) == 1 && dim < ndims - 1; dim++) ; - } + dim = dims.first_non_singleton (); + // calculate extent triplet. l = 1, n = dims(dim), u = 1; for (octave_idx_type i = 0; i < dim; i++) @@ -1003,6 +1109,40 @@ return ret; } +template <class ArrayType> +inline ArrayType +do_mx_diff_op (const ArrayType& src, int dim, octave_idx_type order, + void (*mx_diff_op) (const typename ArrayType::element_type *, + typename ArrayType::element_type *, + octave_idx_type, octave_idx_type, octave_idx_type, + octave_idx_type)) +{ + octave_idx_type l, n, u; + if (order <= 0) + return src; + + dim_vector dims = src.dims (); + + get_extent_triplet (dims, dim, l, n, u); + if (dim >= dims.length ()) + dims.resize (dim+1, 1); + + if (dims(dim) <= order) + { + dims (dim) = 0; + return ArrayType (dims); + } + else + { + dims(dim) -= order; + } + + ArrayType ret (dims); + mx_diff_op (src.data (), ret.fortran_vec (), l, n, u, order); + + return ret; +} + #endif /*
--- a/scripts/ChangeLog +++ b/scripts/ChangeLog @@ -1,3 +1,8 @@ +2009-08-11 Jaroslav Hajek <highegg@gmail.com> + + * general/diff.m: Remove. + * general/Makefile.in: Update. + 2009-08-07 Jaroslav Hajek <highegg@gmail.com> * general/flipdim.m: Fix omitted check.
--- a/scripts/general/Makefile.in +++ b/scripts/general/Makefile.in @@ -36,7 +36,7 @@ SOURCES = __isequal__.m __splinen__.m accumarray.m arrayfun.m \ bicubic.m bitcmp.m bitget.m bitset.m blkdiag.m cart2pol.m \ cart2sph.m cellidx.m cell2mat.m celldisp.m circshift.m colon.m common_size.m \ - cplxpair.m cumtrapz.m dblquad.m deal.m del2.m diff.m display.m flipdim.m \ + cplxpair.m cumtrapz.m dblquad.m deal.m del2.m display.m flipdim.m \ fliplr.m flipud.m genvarname.m gradient.m idivide.m int2str.m \ interp1.m interp1q.m interp2.m interp3.m interpn.m interpft.m \ is_duplicate_entry.m isa.m isdefinite.m isdir.m isequal.m \
deleted file mode 100644 --- a/scripts/general/diff.m +++ /dev/null @@ -1,146 +0,0 @@ -## Copyright (C) 1995, 1996, 1999, 2000, 2002, 2004, 2005, 2006, 2007, -## 2008, 2009 Kurt Hornik -## -## This file is part of Octave. -## -## Octave is free software; you can redistribute it and/or modify it -## under the terms of the GNU General Public License as published by -## the Free Software Foundation; either version 3 of the License, or (at -## your option) any later version. -## -## Octave is distributed in the hope that it will be useful, but -## WITHOUT ANY WARRANTY; without even the implied warranty of -## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU -## General Public License for more details. -## -## You should have received a copy of the GNU General Public License -## along with Octave; see the file COPYING. If not, see -## <http://www.gnu.org/licenses/>. - -## -*- texinfo -*- -## @deftypefn {Function File} {} diff (@var{x}, @var{k}, @var{dim}) -## If @var{x} is a vector of length @var{n}, @code{diff (@var{x})} is the -## vector of first differences -## @tex -## $x_2 - x_1, \ldots{}, x_n - x_{n-1}$. -## @end tex -## @ifnottex -## @var{x}(2) - @var{x}(1), @dots{}, @var{x}(n) - @var{x}(n-1). -## @end ifnottex -## -## If @var{x} is a matrix, @code{diff (@var{x})} is the matrix of column -## differences along the first non-singleton dimension. -## -## The second argument is optional. If supplied, @code{diff (@var{x}, -## @var{k})}, where @var{k} is a non-negative integer, returns the -## @var{k}-th differences. It is possible that @var{k} is larger than -## then first non-singleton dimension of the matrix. In this case, -## @code{diff} continues to take the differences along the next -## non-singleton dimension. -## -## The dimension along which to take the difference can be explicitly -## stated with the optional variable @var{dim}. In this case the -## @var{k}-th order differences are calculated along this dimension. -## In the case where @var{k} exceeds @code{size (@var{x}, @var{dim})} -## then an empty matrix is returned. -## @end deftypefn - -## Author: KH <Kurt.Hornik@wu-wien.ac.at> -## Created: 2 February 1995 -## Adapted-By: jwe - -function x = diff (x, k, dim) - - if (nargin < 1 || nargin > 3) - print_usage (); - endif - - if (nargin < 2 || isempty(k)) - k = 1; - else - if (! (isscalar (k) && k == round (k) && k >= 0)) - error ("diff: k must be a nonnegative integer"); - elseif (k == 0) - return; - endif - endif - - nd = ndims (x); - sz = size (x); - if (nargin != 3) - %% Find the first non-singleton dimension - dim = 1; - while (dim < nd + 1 && sz (dim) == 1) - dim = dim + 1; - endwhile - if (dim > nd) - dim = 1; - endif - else - if (! (isscalar (dim) && dim == round (dim)) && dim > 0 && - dim < (nd + 1)) - error ("diff: dim must be an integer and valid dimension"); - endif - endif - - if (ischar (x)) - error ("diff: symbolic differentiation not (yet) supported"); - endif - - - if (nargin == 3) - if (sz (dim) <= k) - sz(dim) = 0; - x = zeros (sz); - else - n = sz (dim); - idx1 = cell (); - for i = 1:nd - idx1{i} = 1:sz(i); - endfor - idx2 = idx1; - for i = 1 : k; - idx1{dim} = 2 : (n - i + 1); - idx2{dim} = 1 : (n - i); - x = x(idx1{:}) - x(idx2{:}); - endfor - endif - else - if (sum (sz - 1) < k) - x = []; - else - idx1 = cell (); - for i = 1:nd - idx1{i} = 1:sz(i); - endfor - idx2 = idx1; - while (k) - n = sz (dim); - for i = 1 : min (k, n - 1) - idx1{dim} = 2 : (n - i + 1); - idx2{dim} = 1 : (n - i); - x = x(idx1{:}) - x(idx2{:}); - endfor - idx1{dim} = idx2{dim} = 1; - k = k - min (k, n - 1); - dim = dim + 1; - endwhile - endif - endif - -endfunction - -%!assert((diff ([1, 2, 3, 4]) == [1, 1, 1] -%! && diff ([1, 3, 7, 19], 2) == [2, 8] -%! && diff ([1, 2; 5, 4; 8, 7; 9, 6; 3, 1]) == [4, 2; 3, 3; 1, -1; -6, -5] -%! && diff ([1, 2; 5, 4; 8, 7; 9, 6; 3, 1], 3) == [-1, -5; -5, 0] -%! && isempty (diff (1)))); - -%!error diff ([1, 2; 3, 4], -1); - -%!error diff ("foo"); - -%!error diff (); - -%!error diff (1, 2, 3, 4); -
--- a/src/ChangeLog +++ b/src/ChangeLog @@ -1,3 +1,8 @@ +2009-08-11 Jaroslav Hajek <highegg@gmail.com> + + * data.cc (Fdiff): New built-in function. + (do_diff): New assistant function. + 2009-08-10 John W. Eaton <jwe@octave.org> * DLD-FUNCTIONS/dlmread.cc (Fdlmread): Perform tilde expansion on
--- a/src/data.cc +++ b/src/data.cc @@ -3,6 +3,7 @@ Copyright (C) 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009 John W. Eaton Copyright (C) 2009 Jaroslav Hajek +Copyright (C) 2009 VZLU Prague This file is part of Octave. @@ -6159,6 +6160,212 @@ #undef MAKE_INT_BRANCH +template <class SparseT> +static SparseT +do_sparse_diff (const SparseT& array, octave_idx_type order, + int dim) +{ + SparseT retval = array; + if (dim == 1) + { + octave_idx_type k = retval.columns (); + while (order > 0 && k > 0) + { + idx_vector col1 (':'), col2 (':'), sl1 (1, k), sl2 (0, k-1); + retval = retval.index (col1, sl1, 0) - retval.index (col2, sl2, 0); + assert (retval.columns () == k-1); + order--; + k--; + } + } + else + { + octave_idx_type k = retval.rows (); + while (order > 0 && k > 0) + { + idx_vector col1 (':'), col2 (':'), sl1 (1, k), sl2 (0, k-1); + retval = retval.index (sl1, col1, 0) - retval.index (sl2, col2, 0); + assert (retval.rows () == k-1); + order--; + k--; + } + } + + return retval; +} + +static octave_value +do_diff (const octave_value& array, octave_idx_type order, + int dim = -1) +{ + octave_value retval; + + const dim_vector& dv = array.dims (); + if (dim == -1) + { + dim = array.dims ().first_non_singleton (); + + // Bother Matlab. This behavior is really wicked. + if (dv(dim) <= order) + { + if (dv(dim) == 1) + retval = array.resize (dim_vector (0, 0)); + else + { + retval = array; + while (order > 0) + { + if (dim == dv.length ()) + { + retval = do_diff (array, order, dim - 1); + order = 0; + } + else if (dv(dim) == 1) + dim++; + else + { + retval = do_diff (array, dv(dim) - 1, dim); + order -= dv(dim) - 1; + dim++; + } + } + } + + return retval; + } + } + + if (array.is_integer_type ()) + { + if (array.is_int8_type ()) + retval = array.int8_array_value ().diff (order, dim); + else if (array.is_int16_type ()) + retval = array.int16_array_value ().diff (order, dim); + else if (array.is_int32_type ()) + retval = array.int32_array_value ().diff (order, dim); + else if (array.is_int64_type ()) + retval = array.int64_array_value ().diff (order, dim); + else if (array.is_uint8_type ()) + retval = array.uint8_array_value ().diff (order, dim); + else if (array.is_uint16_type ()) + retval = array.uint16_array_value ().diff (order, dim); + else if (array.is_uint32_type ()) + retval = array.uint32_array_value ().diff (order, dim); + else if (array.is_uint64_type ()) + retval = array.uint64_array_value ().diff (order, dim); + else + panic_impossible (); + } + else if (array.is_sparse_type ()) + { + if (array.is_complex_type ()) + retval = do_sparse_diff (array.sparse_complex_matrix_value (), order, dim); + else + retval = do_sparse_diff (array.sparse_matrix_value (), order, dim); + } + else if (array.is_single_type ()) + { + if (array.is_complex_type ()) + retval = array.float_complex_array_value ().diff (order, dim); + else + retval = array.float_array_value ().diff (order, dim); + } + else + { + if (array.is_complex_type ()) + retval = array.complex_array_value ().diff (order, dim); + else + retval = array.array_value ().diff (order, dim); + } + + return retval; +} + +DEFUN (diff, args, , + "-*- texinfo -*-\n\ +@deftypefn {Built-in Function} {} diff (@var{x}, @var{k}, @var{dim})\n\ +If @var{x} is a vector of length @var{n}, @code{diff (@var{x})} is the\n\ +vector of first differences\n\ +@tex\n\ + $x_2 - x_1, \\ldots{}, x_n - x_{n-1}$.\n\ +@end tex\n\ +@ifnottex\n\ + @var{x}(2) - @var{x}(1), @dots{}, @var{x}(n) - @var{x}(n-1).\n\ +@end ifnottex\n\ +\n\ +If @var{x} is a matrix, @code{diff (@var{x})} is the matrix of column\n\ +differences along the first non-singleton dimension.\n\ +\n\ +The second argument is optional. If supplied, @code{diff (@var{x},\n\ +@var{k})}, where @var{k} is a non-negative integer, returns the\n\ +@var{k}-th differences. It is possible that @var{k} is larger than\n\ +then first non-singleton dimension of the matrix. In this case,\n\ +@code{diff} continues to take the differences along the next\n\ +non-singleton dimension.\n\ +\n\ +The dimension along which to take the difference can be explicitly\n\ +stated with the optional variable @var{dim}. In this case the \n\ +@var{k}-th order differences are calculated along this dimension.\n\ +In the case where @var{k} exceeds @code{size (@var{x}, @var{dim})}\n\ +then an empty matrix is returned.\n\ +@end deftypefn") +{ + int nargin = args.length (); + octave_value retval; + + if (nargin < 1 || nargin > 3) + print_usage (); + else if (! args(0).is_numeric_type ()) + error ("diff: X must be numeric"); + + if (! error_state) + { + int dim = -1; + octave_idx_type order = 1; + if (nargin > 1) + { + if (args(1).is_scalar_type ()) + order = args(1).idx_type_value (true, false); + else if (! args(1).is_zero_by_zero ()) + error ("order must be a scalar or []"); + if (! error_state && order < 0) + error ("order must be non-negative"); + } + + if (nargin > 2) + { + dim = args(2).int_value (true, false); + if (! error_state && (dim < 1 || dim > args(0).ndims ())) + error ("needs a valid dimension"); + else + dim -= 1; + } + + if (! error_state) + retval = do_diff (args(0), order, dim); + } + + return retval; +} + +/* + +%!assert (diff ([1, 2, 3, 4]), [1, 1, 1]) +%!assert (diff ([1, 3, 7, 19], 2), [2, 8]) +%!assert (diff ([1, 2; 5, 4; 8, 7; 9, 6; 3, 1]), [4, 2; 3, 3; 1, -1; -6, -5]) +%!assert (diff ([1, 2; 5, 4; 8, 7; 9, 6; 3, 1], 3), [-1, -5; -5, 0]) +%!assert (isempty (diff (1))); + +%!error diff ([1, 2; 3, 4], -1); + +%!error diff ("foo"); + +%!error diff (); + +%!error diff (1, 2, 3, 4); + +*/ + /* ;;; Local Variables: *** ;;; mode: C++ ***