Mercurial > hg > octave-max
changeset 4844:9f7ef92b50b0
[project @ 2004-04-02 17:26:53 by jwe]
line wrap: on
line diff
--- a/liboctave/CNDArray.cc +++ b/liboctave/CNDArray.cc @@ -683,6 +683,192 @@ return ::cat_ra(*this, ra_arg, dim, iidx, move); } +static const Complex Complex_NaN_result (octave_NaN, octave_NaN); + +ComplexNDArray +ComplexNDArray::max (int dim) const +{ + ArrayN<int> dummy_idx; + return max (dummy_idx, dim); +} + +ComplexNDArray +ComplexNDArray::max (ArrayN<int>& idx_arg, int dim) const +{ + dim_vector dv = dims (); + dim_vector dr = dims (); + + if (dv.numel () == 0 || dim > dv.length () || dim < 0) + return ComplexNDArray (); + + dr(dim) = 1; + + ComplexNDArray result (dr); + idx_arg.resize (dr); + + int x_stride = 1; + int x_len = dv(dim); + for (int i = 0; i < dim; i++) + x_stride *= dv(i); + + for (int i = 0; i < dr.numel (); i++) + { + int x_offset; + if (x_stride == 1) + x_offset = i * x_len; + else + { + int x_offset2 = 0; + x_offset = i; + while (x_offset >= x_stride) + { + x_offset -= x_stride; + x_offset2++; + } + x_offset += x_offset2 * x_stride * x_len; + } + + int idx_j; + + Complex tmp_max; + + double abs_max = octave_NaN; + + for (idx_j = 0; idx_j < x_len; idx_j++) + { + tmp_max = elem (idx_j * x_stride + x_offset); + + if (! octave_is_NaN_or_NA (tmp_max)) + { + abs_max = ::abs(tmp_max); + break; + } + } + + for (int j = idx_j+1; j < x_len; j++) + { + Complex tmp = elem (j * x_stride + x_offset); + + if (octave_is_NaN_or_NA (tmp)) + continue; + + double abs_tmp = ::abs (tmp); + + if (abs_tmp > abs_max) + { + idx_j = j; + tmp_max = tmp; + abs_max = abs_tmp; + } + } + + if (octave_is_NaN_or_NA (tmp_max)) + { + result.elem (i) = Complex_NaN_result; + idx_arg.elem (i) = 0; + } + else + { + result.elem (i) = tmp_max; + idx_arg.elem (i) = idx_j; + } + } + + return result; +} + +ComplexNDArray +ComplexNDArray::min (int dim) const +{ + ArrayN<int> dummy_idx; + return min (dummy_idx, dim); +} + +ComplexNDArray +ComplexNDArray::min (ArrayN<int>& idx_arg, int dim) const +{ + dim_vector dv = dims (); + dim_vector dr = dims (); + + if (dv.numel () == 0 || dim > dv.length () || dim < 0) + return ComplexNDArray (); + + dr(dim) = 1; + + ComplexNDArray result (dr); + idx_arg.resize (dr); + + int x_stride = 1; + int x_len = dv(dim); + for (int i = 0; i < dim; i++) + x_stride *= dv(i); + + for (int i = 0; i < dr.numel (); i++) + { + int x_offset; + if (x_stride == 1) + x_offset = i * x_len; + else + { + int x_offset2 = 0; + x_offset = i; + while (x_offset >= x_stride) + { + x_offset -= x_stride; + x_offset2++; + } + x_offset += x_offset2 * x_stride * x_len; + } + + int idx_j; + + Complex tmp_min; + + double abs_min = octave_NaN; + + for (idx_j = 0; idx_j < x_len; idx_j++) + { + tmp_min = elem (idx_j * x_stride + x_offset); + + if (! octave_is_NaN_or_NA (tmp_min)) + { + abs_min = ::abs(tmp_min); + break; + } + } + + for (int j = idx_j+1; j < x_len; j++) + { + Complex tmp = elem (j * x_stride + x_offset); + + if (octave_is_NaN_or_NA (tmp)) + continue; + + double abs_tmp = ::abs (tmp); + + if (abs_tmp < abs_min) + { + idx_j = j; + tmp_min = tmp; + abs_min = abs_tmp; + } + } + + if (octave_is_NaN_or_NA (tmp_min)) + { + result.elem (i) = Complex_NaN_result; + idx_arg.elem (i) = 0; + } + else + { + result.elem (i) = tmp_min; + idx_arg.elem (i) = idx_j; + } + } + + return result; +} + NDArray ComplexNDArray::abs (void) const { @@ -836,6 +1022,141 @@ return is; } +// XXX FIXME XXX -- it would be nice to share code among the min/max +// functions below. + +#define EMPTY_RETURN_CHECK(T) \ + if (nel == 0) \ + return T (dv); + +ComplexNDArray +min (const Complex& c, const ComplexNDArray& m) +{ + dim_vector dv = m.dims (); + int nel = dv.numel (); + + EMPTY_RETURN_CHECK (ComplexNDArray); + + ComplexNDArray result (dv); + + for (int i = 0; i < nel; i++) + { + OCTAVE_QUIT; + result (i) = xmin (c, m (i)); + } + + return result; +} + +ComplexNDArray +min (const ComplexNDArray& m, const Complex& c) +{ + dim_vector dv = m.dims (); + int nel = dv.numel (); + + EMPTY_RETURN_CHECK (ComplexNDArray); + + ComplexNDArray result (dv); + + for (int i = 0; i < nel; i++) + { + OCTAVE_QUIT; + result (i) = xmin (c, m (i)); + } + + return result; +} + +ComplexNDArray +min (const ComplexNDArray& a, const ComplexNDArray& b) +{ + dim_vector dv = a.dims (); + int nel = dv.numel (); + + if (dv != b.dims ()) + { + (*current_liboctave_error_handler) + ("two-arg min expecting args of same size"); + return ComplexNDArray (); + } + + EMPTY_RETURN_CHECK (ComplexNDArray); + + ComplexNDArray result (dv); + + for (int i = 0; i < nel; i++) + { + OCTAVE_QUIT; + result (i) = xmin (a (i), b (i)); + } + + return result; +} + +ComplexNDArray +max (const Complex& c, const ComplexNDArray& m) +{ + dim_vector dv = m.dims (); + int nel = dv.numel (); + + EMPTY_RETURN_CHECK (ComplexNDArray); + + ComplexNDArray result (dv); + + for (int i = 0; i < nel; i++) + { + OCTAVE_QUIT; + result (i) = xmax (c, m (i)); + } + + return result; +} + +ComplexNDArray +max (const ComplexNDArray& m, const Complex& c) +{ + dim_vector dv = m.dims (); + int nel = dv.numel (); + + EMPTY_RETURN_CHECK (ComplexNDArray); + + ComplexNDArray result (dv); + + for (int i = 0; i < nel; i++) + { + OCTAVE_QUIT; + result (i) = xmax (c, m (i)); + } + + return result; +} + +ComplexNDArray +max (const ComplexNDArray& a, const ComplexNDArray& b) +{ + dim_vector dv = a.dims (); + int nel = dv.numel (); + + if (dv != b.dims ()) + { + (*current_liboctave_error_handler) + ("two-arg max expecting args of same size"); + return ComplexNDArray (); + } + + EMPTY_RETURN_CHECK (ComplexNDArray); + + ComplexNDArray result (dv); + + for (int i = 0; i < nel; i++) + { + OCTAVE_QUIT; + result (i) = xmax (a (i), b (i)); + } + + return result; +} + NDS_CMP_OPS(ComplexNDArray, real, Complex, real) NDS_BOOL_OPS(ComplexNDArray, Complex, 0.0)
--- a/liboctave/CNDArray.h +++ b/liboctave/CNDArray.h @@ -89,6 +89,11 @@ ComplexNDArray sumsq (int dim = -1) const; int cat (const ComplexNDArray& ra_arg, int dim, int iidx, int move); + ComplexNDArray max (int dim = 0) const; + ComplexNDArray max (ArrayN<int>& index, int dim = 0) const; + ComplexNDArray min (int dim = 0) const; + ComplexNDArray min (ArrayN<int>& index, int dim = 0) const; + ComplexNDArray& insert (const NDArray& a, int r, int c); ComplexNDArray& insert (const ComplexNDArray& a, int r, int c); @@ -130,6 +135,14 @@ : MArrayN<Complex> (d, dv) { } }; +extern ComplexNDArray min (const Complex& c, const ComplexNDArray& m); +extern ComplexNDArray min (const ComplexNDArray& m, const Complex& c); +extern ComplexNDArray min (const ComplexNDArray& a, const ComplexNDArray& b); + +extern ComplexNDArray max (const Complex& c, const ComplexNDArray& m); +extern ComplexNDArray max (const ComplexNDArray& m, const Complex& c); +extern ComplexNDArray max (const ComplexNDArray& a, const ComplexNDArray& b); + NDS_CMP_OP_DECLS (ComplexNDArray, Complex) NDS_BOOL_OP_DECLS (ComplexNDArray, Complex)
--- a/liboctave/ChangeLog +++ b/liboctave/ChangeLog @@ -1,3 +1,21 @@ +2004-04-02 David Bateman <dbateman@free.fr> + + * lo-specfun.cc (besselj, bessely, besseli, besselk, besselh1, + besselh2, airy, biry, betainc, gammainc, do_bessel): + New NDArray versions. + (SN_BESSEL, NS_BESSEL, NN_BESSEL): New macros. + * lo-specfun.h (besselj, bessely, besseli, besselk, besselh1, + besselh2, airy, biry, betainc, gammainc): Provide decls. + + * dNDArray.cc (NDArray::min, NDArray::max, min, max): + New functions. + * dNDArray.h (NDArray::min, NDArray::max, min, max): Provide decls. + + * CNDArray.cc (ComplexNDArray::min, ComplexNDArray::max, min, max): + New functions. + * CNDArray.h (ComplexNDArray::min, ComplexNDArray::max, min, max): + Provide decls. + 2004-03-17 David Hoover <jazzdaq@yahoo.com> * DASPK.cc (DASPK::do_integrate): Always add n*n elements to the
--- a/liboctave/dNDArray.cc +++ b/liboctave/dNDArray.cc @@ -656,6 +656,156 @@ MX_ND_REAL_OP_REDUCTION (+= elem (iter_idx), 0); } +NDArray +NDArray::max (int dim) const +{ + ArrayN<int> dummy_idx; + return max (dummy_idx, dim); +} + +NDArray +NDArray::max (ArrayN<int>& idx_arg, int dim) const +{ + dim_vector dv = dims (); + dim_vector dr = dims (); + + if (dv.numel () == 0 || dim > dv.length () || dim < 0) + return NDArray (); + + dr(dim) = 1; + + NDArray result (dr); + idx_arg.resize (dr); + + int x_stride = 1; + int x_len = dv(dim); + for (int i = 0; i < dim; i++) + x_stride *= dv(i); + + for (int i = 0; i < dr.numel (); i++) + { + int x_offset; + if (x_stride == 1) + x_offset = i * x_len; + else + { + int x_offset2 = 0; + x_offset = i; + while (x_offset >= x_stride) + { + x_offset -= x_stride; + x_offset2++; + } + x_offset += x_offset2 * x_stride * x_len; + } + + int idx_j; + + double tmp_max = octave_NaN; + + for (idx_j = 0; idx_j < x_len; idx_j++) + { + tmp_max = elem (idx_j * x_stride + x_offset); + + if (! octave_is_NaN_or_NA (tmp_max)) + break; + } + + for (int j = idx_j+1; j < x_len; j++) + { + double tmp = elem (j * x_stride + x_offset); + + if (octave_is_NaN_or_NA (tmp)) + continue; + else if (tmp > tmp_max) + { + idx_j = j; + tmp_max = tmp; + } + } + + result.elem (i) = tmp_max; + idx_arg.elem (i) = octave_is_NaN_or_NA (tmp_max) ? 0 : idx_j; + } + + return result; +} + +NDArray +NDArray::min (int dim) const +{ + ArrayN<int> dummy_idx; + return min (dummy_idx, dim); +} + +NDArray +NDArray::min (ArrayN<int>& idx_arg, int dim) const +{ + dim_vector dv = dims (); + dim_vector dr = dims (); + + if (dv.numel () == 0 || dim > dv.length () || dim < 0) + return NDArray (); + + dr(dim) = 1; + + NDArray result (dr); + idx_arg.resize (dr); + + int x_stride = 1; + int x_len = dv(dim); + for (int i = 0; i < dim; i++) + x_stride *= dv(i); + + for (int i = 0; i < dr.numel (); i++) + { + int x_offset; + if (x_stride == 1) + x_offset = i * x_len; + else + { + int x_offset2 = 0; + x_offset = i; + while (x_offset >= x_stride) + { + x_offset -= x_stride; + x_offset2++; + } + x_offset += x_offset2 * x_stride * x_len; + } + + int idx_j; + + double tmp_min = octave_NaN; + + for (idx_j = 0; idx_j < x_len; idx_j++) + { + tmp_min = elem (idx_j * x_stride + x_offset); + + if (! octave_is_NaN_or_NA (tmp_min)) + break; + } + + for (int j = idx_j+1; j < x_len; j++) + { + double tmp = elem (j * x_stride + x_offset); + + if (octave_is_NaN_or_NA (tmp)) + continue; + else if (tmp < tmp_min) + { + idx_j = j; + tmp_min = tmp; + } + } + + result.elem (i) = tmp_min; + idx_arg.elem (i) = octave_is_NaN_or_NA (tmp_min) ? 0 : idx_j; + } + + return result; +} + int NDArray::cat (const NDArray& ra_arg, int dim, int iidx, int move) { @@ -776,6 +926,141 @@ return is; } +// XXX FIXME XXX -- it would be nice to share code among the min/max +// functions below. + +#define EMPTY_RETURN_CHECK(T) \ + if (nel == 0) \ + return T (dv); + +NDArray +min (double d, const NDArray& m) +{ + dim_vector dv = m.dims (); + int nel = dv.numel (); + + EMPTY_RETURN_CHECK (NDArray); + + NDArray result (dv); + + for (int i = 0; i < nel; i++) + { + OCTAVE_QUIT; + result (i) = xmin (d, m (i)); + } + + return result; +} + +NDArray +min (const NDArray& m, double d) +{ + dim_vector dv = m.dims (); + int nel = dv.numel (); + + EMPTY_RETURN_CHECK (NDArray); + + NDArray result (dv); + + for (int i = 0; i < nel; i++) + { + OCTAVE_QUIT; + result (i) = xmin (d, m (i)); + } + + return result; +} + +NDArray +min (const NDArray& a, const NDArray& b) +{ + dim_vector dv = a.dims (); + int nel = dv.numel (); + + if (dv != b.dims ()) + { + (*current_liboctave_error_handler) + ("two-arg min expecting args of same size"); + return NDArray (); + } + + EMPTY_RETURN_CHECK (NDArray); + + NDArray result (dv); + + for (int i = 0; i < nel; i++) + { + OCTAVE_QUIT; + result (i) = xmin (a (i), b (i)); + } + + return result; +} + +NDArray +max (double d, const NDArray& m) +{ + dim_vector dv = m.dims (); + int nel = dv.numel (); + + EMPTY_RETURN_CHECK (NDArray); + + NDArray result (dv); + + for (int i = 0; i < nel; i++) + { + OCTAVE_QUIT; + result (i) = xmax (d, m (i)); + } + + return result; +} + +NDArray +max (const NDArray& m, double d) +{ + dim_vector dv = m.dims (); + int nel = dv.numel (); + + EMPTY_RETURN_CHECK (NDArray); + + NDArray result (dv); + + for (int i = 0; i < nel; i++) + { + OCTAVE_QUIT; + result (i) = xmax (d, m (i)); + } + + return result; +} + +NDArray +max (const NDArray& a, const NDArray& b) +{ + dim_vector dv = a.dims (); + int nel = dv.numel (); + + if (dv != b.dims ()) + { + (*current_liboctave_error_handler) + ("two-arg max expecting args of same size"); + return NDArray (); + } + + EMPTY_RETURN_CHECK (NDArray); + + NDArray result (dv); + + for (int i = 0; i < nel; i++) + { + OCTAVE_QUIT; + result (i) = xmax (a (i), b (i)); + } + + return result; +} + NDS_CMP_OPS(NDArray, , double, ) NDS_BOOL_OPS(NDArray, double, 0.0)
--- a/liboctave/dNDArray.h +++ b/liboctave/dNDArray.h @@ -87,7 +87,12 @@ NDArray sum (int dim = -1) const; NDArray sumsq (int dim = -1) const; int cat (const NDArray& ra_arg, int dim, int iidx, int move); - + + NDArray max (int dim = 0) const; + NDArray max (ArrayN<int>& index, int dim = 0) const; + NDArray min (int dim = 0) const; + NDArray min (ArrayN<int>& index, int dim = 0) const; + NDArray abs (void) const; ComplexNDArray fourier (int dim = 1) const; @@ -125,6 +130,14 @@ NDArray (double *d, const dim_vector& dv) : MArrayN<double> (d, dv) { } }; +extern NDArray min (double d, const NDArray& m); +extern NDArray min (const NDArray& m, double d); +extern NDArray min (const NDArray& a, const NDArray& b); + +extern NDArray max (double d, const NDArray& m); +extern NDArray max (const NDArray& m, double d); +extern NDArray max (const NDArray& a, const NDArray& b); + NDS_CMP_OP_DECLS (NDArray, double) NDS_BOOL_OP_DECLS (NDArray, double)
--- a/liboctave/lo-specfun.cc +++ b/liboctave/lo-specfun.cc @@ -29,6 +29,8 @@ #include "CMatrix.h" #include "dRowVector.h" #include "dMatrix.h" +#include "dNDArray.h" +#include "CNDArray.h" #include "f77-fcn.h" #include "lo-error.h" #include "lo-ieee.h" @@ -607,6 +609,62 @@ return retval; } +static inline ComplexNDArray +do_bessel (fptr f, const char *, double alpha, const ComplexNDArray& x, + bool scaled, ArrayN<int>& ierr) +{ + dim_vector dv = x.dims (); + int nel = dv.numel (); + ComplexNDArray retval (dv); + + ierr.resize (dv); + + for (int i = 0; i < nel; i++) + retval(i) = f (x(i), alpha, (scaled ? 2 : 1), ierr(i)); + + return retval; +} + +static inline ComplexNDArray +do_bessel (fptr f, const char *, const NDArray& alpha, const Complex& x, + bool scaled, ArrayN<int>& ierr) +{ + dim_vector dv = alpha.dims (); + int nel = dv.numel (); + ComplexNDArray retval (dv); + + ierr.resize (dv); + + for (int i = 0; i < nel; i++) + retval(i) = f (x, alpha(i), (scaled ? 2 : 1), ierr(i)); + + return retval; +} + +static inline ComplexNDArray +do_bessel (fptr f, const char *fn, const NDArray& alpha, + const ComplexNDArray& x, bool scaled, ArrayN<int>& ierr) +{ + dim_vector dv = x.dims (); + ComplexNDArray retval; + + if (dv == alpha.dims ()) + { + int nel = dv.numel (); + + retval.resize (dv); + ierr.resize (dv); + + for (int i = 0; i < nel; i++) + retval(i) = f (x(i), alpha(i), (scaled ? 2 : 1), ierr(i)); + } + else + (*current_liboctave_error_handler) + ("%s: the sizes of alpha and x must conform", fn); + + return retval; +} + static inline ComplexMatrix do_bessel (fptr f, const char *, const RowVector& alpha, const ComplexColumnVector& x, bool scaled, Array2<int>& ierr) @@ -656,6 +714,30 @@ return do_bessel (fcn, #name, alpha, x, scaled, ierr); \ } +#define SN_BESSEL(name, fcn) \ + ComplexNDArray \ + name (double alpha, const ComplexNDArray& x, bool scaled, \ + ArrayN<int>& ierr) \ + { \ + return do_bessel (fcn, #name, alpha, x, scaled, ierr); \ + } + +#define NS_BESSEL(name, fcn) \ + ComplexNDArray \ + name (const NDArray& alpha, const Complex& x, bool scaled, \ + ArrayN<int>& ierr) \ + { \ + return do_bessel (fcn, #name, alpha, x, scaled, ierr); \ + } + +#define NN_BESSEL(name, fcn) \ + ComplexNDArray \ + name (const NDArray& alpha, const ComplexNDArray& x, bool scaled, \ + ArrayN<int>& ierr) \ + { \ + return do_bessel (fcn, #name, alpha, x, scaled, ierr); \ + } + #define RC_BESSEL(name, fcn) \ ComplexMatrix \ name (const RowVector& alpha, const ComplexColumnVector& x, bool scaled, \ @@ -669,6 +751,9 @@ SM_BESSEL (name, fcn) \ MS_BESSEL (name, fcn) \ MM_BESSEL (name, fcn) \ + SN_BESSEL (name, fcn) \ + NS_BESSEL (name, fcn) \ + NN_BESSEL (name, fcn) \ RC_BESSEL (name, fcn) ALL_BESSEL (besselj, zbesj) @@ -778,6 +863,36 @@ return retval; } +ComplexNDArray +airy (const ComplexNDArray& z, bool deriv, bool scaled, ArrayN<int>& ierr) +{ + dim_vector dv = z.dims (); + int nel = dv.numel (); + ComplexNDArray retval (dv); + + ierr.resize (dv); + + for (int i = 0; i < nel; i++) + retval (i) = airy (z(i), deriv, scaled, ierr(i)); + + return retval; +} + +ComplexNDArray +biry (const ComplexNDArray& z, bool deriv, bool scaled, ArrayN<int>& ierr) +{ + dim_vector dv = z.dims (); + int nel = dv.numel (); + ComplexNDArray retval (dv); + + ierr.resize (dv); + + for (int i = 0; i < nel; i++) + retval (i) = biry (z(i), deriv, scaled, ierr(i)); + + return retval; +} + static void gripe_betainc_nonconformant (int r1, int c1, int r2, int c2, int r3, int c3) @@ -787,6 +902,21 @@ r1, c1, r2, c2, r3, c3); } +static dim_vector null_dims (0); + +static void +gripe_betainc_nonconformant (const dim_vector& d1, const dim_vector& d2, + const dim_vector& d3) +{ + std::string d1_str = d1.str (); + std::string d2_str = d2.str (); + std::string d3_str = d3.str (); + + (*current_liboctave_error_handler) + ("betainc: nonconformant arguments (x is %s, a is %s, b is %s)", + d1_str.c_str (), d2_str.c_str (), d3_str.c_str ()); +} + double betainc (double x, double a, double b) { @@ -850,6 +980,56 @@ return retval; } +NDArray +betainc (double x, double a, const NDArray& b) +{ + dim_vector dv = b.dims (); + int nel = dv.numel (); + + NDArray retval (dv); + + for (int i = 0; i < nel; i++) + retval (i) = betainc (x, a, b(i)); + + return retval; +} + +NDArray +betainc (double x, const NDArray& a, double b) +{ + dim_vector dv = a.dims (); + int nel = dv.numel (); + + NDArray retval (dv); + + for (int i = 0; i < nel; i++) + retval (i) = betainc (x, a(i), b); + + return retval; +} + +NDArray +betainc (double x, const NDArray& a, const NDArray& b) +{ + NDArray retval; + dim_vector dv = a.dims (); + + if (dv == b.dims ()) + { + int nel = dv.numel (); + + retval.resize (dv); + + for (int i = 0; i < nel; i++) + retval (i) = betainc (x, a(i), b(i)); + } + else + gripe_betainc_nonconformant (dim_vector (0), dv, b.dims ()); + + return retval; +} + + Matrix betainc (const Matrix& x, double a, double b) { @@ -943,6 +1123,83 @@ return retval; } +NDArray +betainc (const NDArray& x, double a, double b) +{ + dim_vector dv = x.dims (); + int nel = dv.numel (); + + NDArray retval (dv); + + for (int i = 0; i < nel; i++) + retval (i) = betainc (x(i), a, b); + + return retval; +} + +NDArray +betainc (const NDArray& x, double a, const NDArray& b) +{ + NDArray retval; + dim_vector dv = x.dims (); + + if (dv == b.dims ()) + { + int nel = dv.numel (); + + retval.resize (dv); + + for (int i = 0; i < nel; i++) + retval (i) = betainc (x(i), a, b(i)); + } + else + gripe_betainc_nonconformant (dv, dim_vector (0), b.dims ()); + + return retval; +} + +NDArray +betainc (const NDArray& x, const NDArray& a, double b) +{ + NDArray retval; + dim_vector dv = x.dims (); + + if (dv == a.dims ()) + { + int nel = dv.numel (); + + retval.resize (dv); + + for (int i = 0; i < nel; i++) + retval (i) = betainc (x(i), a(i), b); + } + else + gripe_betainc_nonconformant (dv, a.dims (), dim_vector (0)); + + return retval; +} + +NDArray +betainc (const NDArray& x, const NDArray& a, const NDArray& b) +{ + NDArray retval; + dim_vector dv = x.dims (); + + if (dv == a.dims () && dv == b.dims ()) + { + int nel = dv.numel (); + + retval.resize (dv); + + for (int i = 0; i < nel; i++) + retval (i) = betainc (x(i), a(i), b(i)); + } + else + gripe_betainc_nonconformant (dv, a.dims (), b.dims ()); + + return retval; +} + // XXX FIXME XXX -- there is still room for improvement here... double @@ -1058,6 +1315,98 @@ return retval; } +NDArray +gammainc (double x, const NDArray& a) +{ + dim_vector dv = a.dims (); + int nel = dv.numel (); + + NDArray retval; + NDArray result (dv); + + bool err; + + for (int i = 0; i < nel; i++) + { + result (i) = gammainc (x, a(i), err); + + if (err) + goto done; + } + + retval = result; + + done: + + return retval; +} + +NDArray +gammainc (const NDArray& x, double a) +{ + dim_vector dv = x.dims (); + int nel = dv.numel (); + + NDArray retval; + NDArray result (dv); + + bool err; + + for (int i = 0; i < nel; i++) + { + result (i) = gammainc (x(i), a, err); + + if (err) + goto done; + } + + retval = result; + + done: + + return retval; +} + +NDArray +gammainc (const NDArray& x, const NDArray& a) +{ + dim_vector dv = x.dims (); + int nel = dv.numel (); + + NDArray retval; + NDArray result; + + if (dv == a.dims ()) + { + result.resize (dv); + + bool err; + + for (int i = 0; i < nel; i++) + { + result (i) = gammainc (x(i), a(i), err); + + if (err) + goto done; + } + + retval = result; + } + else + { + std::string x_str = dv.str (); + std::string a_str = a.dims ().str (); + + (*current_liboctave_error_handler) + ("gammainc: nonconformant arguments (arg 1 is %s, arg 2 is %s)", + x_str.c_str (), a_str. c_str ()); + } + + done: + + return retval; +} + /* ;;; Local Variables: *** ;;; mode: C++ ***
--- a/liboctave/lo-specfun.h +++ b/liboctave/lo-specfun.h @@ -24,10 +24,13 @@ #define octave_liboctave_specfun_h 1 #include "oct-cmplx.h" +#include "ArrayN.h" template <class T> class Array2; class Matrix; class ComplexMatrix; +class NDArray; +class ComplexNDArray; class RowVector; class ComplexColumnVector; class Range; @@ -145,6 +148,78 @@ besselh2 (const Matrix& alpha, const ComplexMatrix& x, bool scaled, Array2<int>& ierr); +extern ComplexNDArray +besselj (double alpha, const ComplexNDArray& x, bool scaled, + ArrayN<int>& ierr); + +extern ComplexNDArray +bessely (double alpha, const ComplexNDArray& x, bool scaled, + ArrayN<int>& ierr); + +extern ComplexNDArray +besseli (double alpha, const ComplexNDArray& x, bool scaled, + ArrayN<int>& ierr); + +extern ComplexNDArray +besselk (double alpha, const ComplexNDArray& x, bool scaled, + ArrayN<int>& ierr); + +extern ComplexNDArray +besselh1 (double alpha, const ComplexNDArray& x, bool scaled, + ArrayN<int>& ierr); + +extern ComplexNDArray +besselh2 (double alpha, const ComplexNDArray& x, bool scaled, + ArrayN<int>& ierr); + +extern ComplexNDArray +besselj (const NDArray& alpha, const Complex& x, bool scaled, + ArrayN<int>& ierr); + +extern ComplexNDArray +bessely (const NDArray& alpha, const Complex& x, bool scaled, + ArrayN<int>& ierr); + +extern ComplexNDArray +besseli (const NDArray& alpha, const Complex& x, bool scaled, + ArrayN<int>& ierr); + +extern ComplexNDArray +besselk (const NDArray& alpha, const Complex& x, bool scaled, + ArrayN<int>& ierr); + +extern ComplexNDArray +besselh1 (const NDArray& alpha, const Complex& x, bool scaled, + ArrayN<int>& ierr); + +extern ComplexNDArray +besselh2 (const NDArray& alpha, const Complex& x, bool scaled, + ArrayN<int>& ierr); + +extern ComplexNDArray +besselj (const NDArray& alpha, const ComplexNDArray& x, bool scaled, + ArrayN<int>& ierr); + +extern ComplexNDArray +bessely (const NDArray& alpha, const ComplexNDArray& x, bool scaled, + ArrayN<int>& ierr); + +extern ComplexNDArray +besseli (const NDArray& alpha, const ComplexNDArray& x, bool scaled, + ArrayN<int>& ierr); + +extern ComplexNDArray +besselk (const NDArray& alpha, const ComplexNDArray& x, bool scaled, + ArrayN<int>& ierr); + +extern ComplexNDArray +besselh1 (const NDArray& alpha, const ComplexNDArray& x, bool scaled, + ArrayN<int>& ierr); + +extern ComplexNDArray +besselh2 (const NDArray& alpha, const ComplexNDArray& x, bool scaled, + ArrayN<int>& ierr); + extern ComplexMatrix besselj (const RowVector& alpha, const ComplexColumnVector& x, bool scaled, Array2<int>& ierr); @@ -178,21 +253,40 @@ extern ComplexMatrix biry (const ComplexMatrix& z, bool deriv, bool scaled, Array2<int>& ierr); +extern ComplexNDArray +airy (const ComplexNDArray& z, bool deriv, bool scaled, ArrayN<int>& ierr); + +extern ComplexNDArray +biry (const ComplexNDArray& z, bool deriv, bool scaled, ArrayN<int>& ierr); + extern double betainc (double x, double a, double b); extern Matrix betainc (double x, double a, const Matrix& b); extern Matrix betainc (double x, const Matrix& a, double b); extern Matrix betainc (double x, const Matrix& a, const Matrix& b); +extern NDArray betainc (double x, double a, const NDArray& b); +extern NDArray betainc (double x, const NDArray& a, double b); +extern NDArray betainc (double x, const NDArray& a, const NDArray& b); + extern Matrix betainc (const Matrix& x, double a, double b); extern Matrix betainc (const Matrix& x, double a, const Matrix& b); extern Matrix betainc (const Matrix& x, const Matrix& a, double b); extern Matrix betainc (const Matrix& x, const Matrix& a, const Matrix& b); +extern NDArray betainc (const NDArray& x, double a, double b); +extern NDArray betainc (const NDArray& x, double a, const NDArray& b); +extern NDArray betainc (const NDArray& x, const NDArray& a, double b); +extern NDArray betainc (const NDArray& x, const NDArray& a, const NDArray& b); + extern double gammainc (double x, double a, bool& err); extern Matrix gammainc (double x, const Matrix& a); extern Matrix gammainc (const Matrix& x, double a); extern Matrix gammainc (const Matrix& x, const Matrix& a); +extern NDArray gammainc (double x, const NDArray& a); +extern NDArray gammainc (const NDArray& x, double a); +extern NDArray gammainc (const NDArray& x, const NDArray& a); + inline double gammainc (double x, double a) { bool err;
--- a/scripts/ChangeLog +++ b/scripts/ChangeLog @@ -1,3 +1,25 @@ +2004-04-02 David Bateman <dbateman@free.fr> + + * statistics/base/std.m: Allow optional args for type and dim. + * statistics/base/center.m statistics/base/meansq.m + statistics/base/moment.m statistics/base/range.m: Update for NDArrays. + * signal/fftshift.m: Fix dimensioning error. + + * statistics/base/std.m: Use repmat not ones(nr,1)*mean to allow + NDArrays. + + * general/mod.m, general/mod.m: Allow NDArrays with one scalar arg. + + * signal/fftshift.m: Update for NDArrays, allow optional dim arg. + + * specfun/erfinv.m, general/repmat.m: Update for NDArrays. + + * control/base/bode.m, control/base/lqg.m, control/system/ss2sys.m, + control/system/cellidx.m, control/system/dmr2d.m control/system/ss.m, + control/system/sysprune.m: Doc update for usage of cell arrays. + + * control/system/sysidx.m: Use cellidx and not listidx. + 2004-03-12 John W. Eaton <jwe@bevo.che.wisc.edu> * plot/__pltopt1__.m: Always add title clause to plot command with
--- a/scripts/control/base/bode.m +++ b/scripts/control/base/bode.m @@ -67,7 +67,7 @@ ## ## @strong{Example} ## @example -## bode(sys,[],"y_3",list("u_1","u_4"); +## bode(sys,[],"y_3", @{"u_1","u_4"@}); ## @end example ## @end table ## @strong{Outputs}
--- a/scripts/control/base/lqg.m +++ b/scripts/control/base/lqg.m @@ -40,7 +40,7 @@ ## @itemx r ## state, control weighting respectively. Control ARE is ## @item in_idx -## names or indices of controlled inputs (see @code{sysidx}, @code{listidx}) +## names or indices of controlled inputs (see @code{sysidx}, @code{cellidx}) ## ## default: last dim(R) inputs are assumed to be controlled inputs, all ## others are assumed to be noise inputs.
--- a/scripts/control/system/cellidx.m +++ b/scripts/control/system/cellidx.m @@ -17,7 +17,7 @@ ## Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111 USA. ## -*- texinfo -*- -## @deftypefn {Function File} {[@var{idxvec}, @var{errmsg}] =} listidx (@var{listvar}, @var{strlist}) +## @deftypefn {Function File} {[@var{idxvec}, @var{errmsg}] =} cellidx (@var{listvar}, @var{strlist}) ## Return indices of string entries in @var{listvar} that match strings ## in @var{strlist}. ## @@ -29,7 +29,7 @@ ## ## If @var{strlist} contains a string not in @var{listvar}, then ## an error message is returned in @var{errmsg}. If only one output -## argument is requested, then @var{listidx} prints @var{errmsg} to the +## argument is requested, then @var{cellidx} prints @var{errmsg} to the ## screen and exits with an error. ## @end deftypefn
--- a/scripts/control/system/dmr2d.m +++ b/scripts/control/system/dmr2d.m @@ -29,7 +29,7 @@ ## @code{dmr2d} exits with an error if @var{sys} is not discrete ## @item idx ## indices or names of states with sampling time -## @code{sysgettsam(@var{sys})} (may be empty); see @code{listidx} +## @code{sysgettsam(@var{sys})} (may be empty); see @code{cellidx} ## @item sprefix ## list of string prefixes of states with sampling time ## @code{sysgettsam(@var{sys})} (may be empty)
--- a/scripts/control/system/ss.m +++ b/scripts/control/system/ss.m @@ -46,18 +46,18 @@ ## see below for system partitioning ## ## @item stname -## list of strings of state signal names +## cell array of strings of state signal names ## ## default (@var{stname}=[] on input): @code{x_n} for continuous states, ## @code{xd_n} for discrete states ## ## @item inname -## list of strings of input signal names +## cell array of strings of input signal names ## ## default (@var{inname} = [] on input): @code{u_n} ## ## @item outname -## list of strings of input signal names +## cell array of strings of input signal names ## ## default (@var{outname} = [] on input): @code{y_n} ## @@ -140,7 +140,7 @@ ## octave:1> a = [1 2 3; 4 5 6; 7 8 10]; ## octave:2> b = [0 0 ; 0 1 ; 1 0]; ## octave:3> c = eye(3); -## octave:4> sys = ss(a,b,c,[],0,3,0,list("volts","amps","joules")); +## octave:4> sys = ss(a,b,c,[],0,3,0,@{"volts","amps","joules"@}); ## octave:5> sysout(sys); ## Input(s) ## 1: u_1
--- a/scripts/control/system/ss2sys.m +++ b/scripts/control/system/ss2sys.m @@ -46,18 +46,18 @@ ## see below for system partitioning ## ## @item stname -## list of strings of state signal names +## cell array of strings of state signal names ## ## default (@var{stname}=[] on input): @code{x_n} for continuous states, ## @code{xd_n} for discrete states ## ## @item inname -## list of strings of input signal names +## cell array of strings of input signal names ## ## default (@var{inname} = [] on input): @code{u_n} ## ## @item outname -## list of strings of input signal names +## cell array of strings of input signal names ## ## default (@var{outname} = [] on input): @code{y_n} ## @@ -140,7 +140,7 @@ ## octave:1> a = [1 2 3; 4 5 6; 7 8 10]; ## octave:2> b = [0 0 ; 0 1 ; 1 0]; ## octave:3> c = eye(3); -## octave:4> sys = ss(a,b,c,[],0,3,0,list("volts","amps","joules")); +## octave:4> sys = ss(a,b,c,[],0,3,0,@{"volts","amps","joules"@}); ## octave:5> sysout(sys); ## Input(s) ## 1: u_1
--- a/scripts/control/system/sysidx.m +++ b/scripts/control/system/sysidx.m @@ -38,13 +38,13 @@ endif ## extract correct set of signal names values - [idxvec, msg] = listidx (list ("in", "out", "st", "yd"), sigtype); + [idxvec, msg] = cellidx ({"in", "out", "st", "yd"}, sigtype); if (msg) error ("invalid sigtype = %s", sigtype); endif syssiglist = sysgetsignals (sys, sigtype); - [idxvec, msg] = listidx (syssiglist, signamelist); + [idxvec, msg] = cellidx (syssiglist, signamelist); if (length (msg)) error ("sysidx (sigtype = %s): %s", sigtype, strrep (msg, "strlist", "signamelist"));
--- a/scripts/control/system/sysprune.m +++ b/scripts/control/system/sysprune.m @@ -33,7 +33,7 @@ ## ## @example ## retsys = sysprune(Asys,[1:3,4],"u_1"); -## retsys = sysprune(Asys,list("tx","ty","tz"), 4); +## retsys = sysprune(Asys,@{"tx","ty","tz"@}, 4); ## @end example ## ## @end table
--- a/scripts/general/mod.m +++ b/scripts/general/mod.m @@ -44,7 +44,8 @@ usage ("r = mod (x, y)"); endif - if (any (size (x) != size (y)) && ! (isscalar (x) || isscalar (y))) + if (((ndims (x) != ndims (y)) || any (size (x) != size (y))) && + ! (isscalar (x) || isscalar (y))) error ("mod: argument sizes must agree"); endif
--- a/scripts/general/rem.m +++ b/scripts/general/rem.m @@ -39,7 +39,8 @@ usage ("rem (x, y)"); endif - if (any (size (x) != size (y)) && ! (isscalar (x) || isscalar (y))) + if (((ndims (x) != ndims (y)) || any (size (x) != size (y))) && + ! (isscalar (x) || isscalar (y))) error ("rem: argument sizes must agree"); endif
--- a/scripts/general/repmat.m +++ b/scripts/general/repmat.m @@ -34,21 +34,47 @@ usage ("repmat (a, m, n)"); endif - if (nargin == 2) + if (nargin == 3) + if (!isscalar (m) && !isscalar (n)) + error ("repmat: with 3 arguments m and n must be scalar"); + endif + idx = [m, n]; + else if (isscalar (m)) + idx = [m, m]; n = m; - elseif (isvector (m) && length (m) == 2) - n = m(2); - m = m(1); + elseif (isvector (m) && length (m) > 1) + # Ensure that we have a row vector + idx = m(:).'; else - error ("repmat: only builds 2D matrices") + error ("repmat: invalid dimensional argument"); endif endif - if (isstr (a)) - x = setstr (kron (ones (m, n), toascii (a))); + if (numel (a) == 1) + if (isstr (a)) + x = setstr (toascii (a) * ones (idx)); + else + x = a * ones(idx); + endif + elseif (ndims (a) == 2 && length (idx) < 3) + if (isstr (a)) + x = setstr (kron (ones (idx), toascii (a))); + else + x = kron (ones (idx), a); + endif else - x = kron (ones (m, n), a); + aidx = size(a); + if (length(aidx) > length(idx)) + idx = [idx, ones(1,length(aidx)-length(idx))]; + elseif (length(aidx) < length(idx)) + aidx = [aidx, ones(1,length(idx)-length(aidx))]; + endif + cidx = cell (); + for i=1:length(aidx) + cidx{i} = kron (ones (1, idx(i)), 1:aidx(i)); + endfor + x = a (cidx{:}); endif endfunction
--- a/scripts/signal/fftshift.m +++ b/scripts/signal/fftshift.m @@ -19,6 +19,7 @@ ## -*- texinfo -*- ## @deftypefn {Function File} {} fftshift (@var{v}) +## @deftypefnx {Function File} {} fftshift (@var{v}, @var{dim}) ## Perform a shift of the vector @var{v}, for use with the @code{fft} ## and @code{ifft} functions, in order the move the frequency 0 to the ## center of the vector or matrix. @@ -31,32 +32,55 @@ ## f = ((1:N) - ceil(N/2)) / N / Dt ## @end example ## -## If @var{v} is a matrix, the same holds for rows and columns. +## If @var{v} is a matrix, the same holds for rows and columns. If +## @var{v} is an array, then the same holds along each dimension. +## +## The optional @var{dim} argument can be used to limit the dimension +## along which the permutation occurs. ## @end deftypefn ## Author: Vincent Cautaerts <vincent@comf5.comm.eng.osaka-u.ac.jp> ## Created: July 1997 ## Adapted-By: jwe -function retval = fftshift (V) +function retval = fftshift (V, dim) retval = 0; - if (nargin != 1) - usage ("usage: fftshift (X)"); + if (nargin != 1 && nargin != 2) + usage ("usage: fftshift (X, dim)"); endif - if (isvector (V)) - x = length (V); - xx = ceil (x/2); - retval = V([xx+1:x, 1:xx]); - elseif (ismatrix (V)) - [x, y] = size (V); - xx = ceil (x/2); - yy = ceil (y/2); - retval = V([xx+1:x, 1:xx], [yy+1:y, 1:yy]); + if (nargin == 2) + if (!isscalar (dim)) + error ("fftshift: dimension must be an integer scalar"); + endif + nd = ndims (V); + sz = size (V); + sz2 = ceil (sz(dim) / 2); + idx = cell (); + for i=1:nd + idx {i} = 1:sz(i); + endfor + idx {dim} = [sz2+1:sz(dim), 1:sz2]; + retval = V (idx{:}); else - error ("fftshift: expecting vector or matrix argument"); + if (isvector (V)) + x = length (V); + xx = ceil (x/2); + retval = V([xx+1:x, 1:xx]); + elseif (ismatrix (V)) + nd = ndims (V); + sz = size (V); + sz2 = ceil (sz ./ 2); + idx = cell (); + for i=1:nd + idx{i} = [sz2(i)+1:sz(i), 1:sz2(i)]; + endfor + retval = V (idx{:}); + else + error ("fftshift: expecting vector or matrix argument"); + endif endif endfunction
--- a/scripts/specfun/erfinv.m +++ b/scripts/specfun/erfinv.m @@ -38,9 +38,11 @@ iterations = 0; - [m, n] = size (x); - x = reshape (x, m * n, 1); - y = zeros (m * n, 1); + sz = size (x); + nel = numel (x); + + x = reshape (x, nel, 1); + y = zeros (nel, 1); i = find ((x < -1) | (x > 1) | isnan(x)); if any (i) @@ -69,6 +71,6 @@ y(i) = z_new; endif - y = reshape (y, m, n); + y = reshape (y, sz); endfunction
--- a/scripts/statistics/base/center.m +++ b/scripts/statistics/base/center.m @@ -19,27 +19,40 @@ ## -*- texinfo -*- ## @deftypefn {Function File} {} center (@var{x}) +## @deftypefnx {Function File} {} center (@var{x}, @var{dim}) ## If @var{x} is a vector, subtract its mean. ## If @var{x} is a matrix, do the above for each column. +## If the optional argument @var{dim} is given, perform the above +## operation along this dimension ## @end deftypefn ## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> ## Description: Center by subtracting means -function retval = center (x) +function retval = center (x, varargin) - if (nargin != 1) + if (nargin != 1 && nargin != 2) usage ("center (x)"); endif if (isvector (x)) - retval = x - mean (x); + retval = x - mean (x, varargin{:}); elseif (ismatrix (x)) - retval = x - ones (rows (x), 1) * mean (x); + if nargin < 2 + dim = min (find (size (x) > 1)); + if isempty (dim), + dim=1; + endif; + else + dim = varargin {1}; + endif + sz = ones (1, ndims (x)); + sz (dim) = size (x, dim); + retval = x - repmat (mean (x, dim), sz); elseif (isempty (x)) retval = x; else error ("center: x must be a vector or a matrix"); endif -endfunction \ No newline at end of file +endfunction
--- a/scripts/statistics/base/meansq.m +++ b/scripts/statistics/base/meansq.m @@ -19,20 +19,22 @@ ## -*- texinfo -*- ## @deftypefn {Function File} {} meansq (@var{x}) +## @deftypefnx {Function File} {} meansq (@var{x}, @var{dim}) ## For vector arguments, return the mean square of the values. ## For matrix arguments, return a row vector contaning the mean square -## of each column. +## of each column. With the optional @var{dim} argument, returns the +## mean squared of the values along this dimension ## @end deftypefn ## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> ## Description: Compute mean square -function y = meansq (x) +function y = meansq (x, varargin) - if (nargin != 1) + if (nargin != 1 && nargin != 2) usage ("meansq (x)"); endif - y = mean (x.^2); + y = mean (x.^2, varargin{:}); endfunction
--- a/scripts/statistics/base/moment.m +++ b/scripts/statistics/base/moment.m @@ -18,7 +18,7 @@ ## 02111-1307, USA. ## -*- texinfo -*- -## @deftypefn {Function File} {} moment (@var{x}, @var{p}, @var{opt}) +## @deftypefn {Function File} {} moment (@var{x}, @var{p}, @var{opt},@var{dim}) ## If @var{x} is a vector, compute the @var{p}-th moment of @var{x}. ## ## If @var{x} is a matrix, return the row vector containing the @@ -34,6 +34,9 @@ ## ## @noindent ## computes the third central absolute moment of @var{x}. +## +## If the optional argument @var{dim} is supplied, work along dimension +## @var{dim}. ## @end deftypefn ## Can easily be made to work for continuous distributions (using quad) @@ -42,35 +45,73 @@ ## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> ## Description: Compute moments -function m = moment (x, p, opt) +function m = moment (x, p, opt1, opt2) - if ((nargin < 2) || (nargin > 3)) - usage ("moment (x, p, type"); + if ((nargin < 2) || (nargin > 4)) + usage ("moment (x, p, type, dim)"); endif - [nr, nc] = size (x); - if (nr == 0 || nc == 0) - error ("moment: x must not be empty"); - elseif (nr == 1) - x = reshape (x, nc, 1); - nr = nc; + need_dim = 0; + + if (nargin == 1) + opt = ""; + need_dim = 1; + elseif (nargin == 2) + if (isstr (opt1)) + opt = opt1; + need_dim = 1; + else + dim = opt1; + opt = ""; + endif + elseif (nargin == 3) + if (isstr (opt1)) + opt = opt1; + dim = opt2; + elseif (isstr (opt2)) + opt = opt2; + dim = opt1; + else + usage ("moment: expecting opt to be a string"); + endif + else + usage ("moment (x, p, dim, opt) or moment (x, p, dim, opt)"); endif - if (nargin == 3) - tmp = warn_str_to_num; - unwind_protect - warn_str_to_num = 0; - if any (opt == "c") - x = x - ones (nr, 1) * sum (x) / nr; - endif - if any (opt == "a") - x = abs (x); - endif - unwind_protect_cleanup - warn_str_to_num = tmp; - end_unwind_protect + sz = size(x); + n = sz (dim); + + if (need_dim) + t = find (size (x) != 1); + if (isempty (t)) + dim = 1; + else + dim = t(1); + endif + endif + + sz = size (x); + n = sz (dim); + + if (numels (x) < 1) + error ("moment: x must not be empty"); endif - m = sum(x .^ p) / nr; + tmp = warn_str_to_num; + unwind_protect + warn_str_to_num = 0; + if any (opt == "c") + rng = ones(1, length (sz)); + rng (dim) = sz (dim); + x = x - repmat (sum (x, dim), rng) / n; + endif + if any (opt == "a") + x = abs (x); + endif + unwind_protect_cleanup + warn_str_to_num = tmp; + end_unwind_protect + + m = sum(x .^ p, dim) / n; endfunction
--- a/scripts/statistics/base/std.m +++ b/scripts/statistics/base/std.m @@ -19,6 +19,8 @@ ## -*- texinfo -*- ## @deftypefn {Function File} {} std (@var{x}) +## @deftypefnx {Function File} {} std (@var{x}, @var{opt}) +## @deftypefnx {Function File} {} std (@var{x}, @var{opt}, @var{dim}) ## If @var{x} is a vector, compute the standard deviation of the elements ## of @var{x}. ## @iftex @@ -38,26 +40,50 @@ ## @end ifinfo ## If @var{x} is a matrix, compute the standard deviation for ## each column and return them in a row vector. +## +## The argument @var{opt} determines the type of normalization to use. Valid values +## are +## +## @table @asis +## @item 0: +## normalizes with N-1, provides the square root of best unbiased estimator of +## the variance [default] +## @item 1: +## normalizes with N, this provides the square root of the second moment around +## the mean +## @end table +## +## The third argument @var{dim} determines the dimension along which the standard +## deviation is calculated. ## @end deftypefn ## @seealso{mean and median} ## Author: jwe -function retval = std (a) +function retval = std (a, opt, dim) - if (nargin != 1) - usage ("std (a)"); + if (nargin < 1 || nargin > 3) + usage ("std (a, opt, dim)"); + endif + if nargin < 3 + dim = min(find(size(a)>1)); + if isempty(dim), dim=1; endif; + endif + if ((nargin < 2) || isempty(opt)) + opt = 0; endif - nr = rows (a); - nc = columns (a); - if (nc == 1 && nr == 1) - retval = 0; - elseif (nc == 1 || nr == 1) - n = length (a); - retval = sqrt (sumsq (a - mean (a)) / (n - 1)); - elseif (nr > 1 && nc > 0) - retval = sqrt (sumsq (a - ones (nr, 1) * mean (a)) / (nr - 1)); + sz = size(a); + if (sz (dim) == 1) + retval = zeros(sz); + elseif (numel (a) > 0) + rng = ones (1, length (sz)); + rng (dim) = sz (dim); + if (opt == 0) + retval = sqrt (sumsq (a - repmat(mean (a, dim), rng), dim) / (sz(dim) - 1)); + else + retval = sqrt (sumsq (a - repmat(mean (a, dim), rng), dim) / sz(dim)); + endif else error ("std: invalid matrix argument"); endif
--- a/src/ChangeLog +++ b/src/ChangeLog @@ -1,3 +1,26 @@ +2004-04-02 Quentin Spencer <qspencer@ieee.org> + + * parse.y: Use persistent instead of static in warnings messages. + +2004-04-02 John W. Eaton <jwe@bevo.che.wisc.edu> + + * pt-decl.cc (tree_static_command::do_init): Initialize to empty + matrix by default. + +2004-04-02 David Bateman <dbateman@free.fr> + + * ov-re-mat.cc (octave_matrix::convert_to_str_internal): + Return charNDArray. + * ov-bool-mat.cc (octave_bool_matrix::convert_to_str_internal): + Call array_value. + + * DLD-FUNCTIONS/besselj.cc, DLD-FUNCTIONS/betainc.cc, + DLD-FUNCTIONS/gammainc.cc, DLD-FUNCTIONS/minmax.cc, + DLD-FUNCTIONS/filter.cc: + Convert for NDArray, better Matlab compatibility. + + * load-save.cc (Fload): Better handling of non existent files. + 2004-03-19 John W. Eaton <jwe@bevo.che.wisc.edu> * ov-list.cc (octave_list::subsref): Correctly create return value.
--- a/src/DLD-FUNCTIONS/besselj.cc +++ b/src/DLD-FUNCTIONS/besselj.cc @@ -97,6 +97,24 @@ return retval; } +static inline NDArray +int_arrayN_to_array (const ArrayN<int>& a) +{ + dim_vector dv = a.dims (); + int nel = dv.numel (); + + NDArray retval (dv); + + for (int i = 0; i < nel; i++) + { + OCTAVE_QUIT; + + retval(i) = (double) (a(i)); + } + + return retval; +} + static void gripe_bessel_arg (const char *fn, const char *arg) { @@ -145,17 +163,17 @@ } else { - ComplexMatrix x = x_arg.complex_matrix_value (); + ComplexNDArray x = x_arg.complex_array_value (); if (! error_state) { - Array2<int> ierr; + ArrayN<int> ierr; octave_value result; DO_BESSEL (type, alpha, x, scaled, ierr, result); if (nargout > 1) - retval(1) = int_array2_to_matrix (ierr); + retval(1) = int_arrayN_to_array (ierr); retval(0) = result; } @@ -168,21 +186,28 @@ } else { - Matrix alpha = args(0).matrix_value (); + dim_vector dv0 = args(0).dims (); + dim_vector dv1 = args(1).dims (); + + bool args0_is_row_vector = (dv0 (1) == dv0.numel ()); + bool args1_is_col_vector = (dv1 (0) == dv1.numel ()); - if (! error_state) + if (args0_is_row_vector && args1_is_col_vector) { - if (x_arg.is_scalar_type ()) + RowVector ralpha = args(0).row_vector_value (); + + if (! error_state) { - Complex x = x_arg.complex_value (); + ComplexColumnVector cx = + x_arg.complex_column_vector_value (); if (! error_state) { Array2<int> ierr; octave_value result; - DO_BESSEL (type, alpha, x, scaled, ierr, result); - + DO_BESSEL (type, ralpha, cx, scaled, ierr, result); + if (nargout > 1) retval(1) = int_array2_to_matrix (ierr); @@ -192,45 +217,56 @@ gripe_bessel_arg (fn, "second"); } else - { - ComplexMatrix x = x_arg.complex_matrix_value (); - - if (! error_state) - { - if (alpha.rows () == 1 && x.columns () == 1) - { - RowVector ralpha = alpha.row (0); - ComplexColumnVector cx = x.column (0); + gripe_bessel_arg (fn, "first"); + } + else + { + NDArray alpha = args(0).array_value (); - Array2<int> ierr; - octave_value result; - - DO_BESSEL (type, ralpha, cx, scaled, ierr, result); + if (! error_state) + { + if (x_arg.is_scalar_type ()) + { + Complex x = x_arg.complex_value (); - if (nargout > 1) - retval(1) = int_array2_to_matrix (ierr); - - retval(0) = result; - } - else + if (! error_state) { - Array2<int> ierr; + ArrayN<int> ierr; octave_value result; DO_BESSEL (type, alpha, x, scaled, ierr, result); if (nargout > 1) - retval(1) = int_array2_to_matrix (ierr); + retval(1) = int_arrayN_to_array (ierr); retval(0) = result; } + else + gripe_bessel_arg (fn, "second"); } else - gripe_bessel_arg (fn, "second"); + { + ComplexNDArray x = x_arg.complex_array_value (); + + if (! error_state) + { + ArrayN<int> ierr; + octave_value result; + + DO_BESSEL (type, alpha, x, scaled, ierr, result); + + if (nargout > 1) + retval(1) = int_arrayN_to_array (ierr); + + retval(0) = result; + } + else + gripe_bessel_arg (fn, "second"); + } } + else + gripe_bessel_arg (fn, "first"); } - else - gripe_bessel_arg (fn, "first"); } } else @@ -422,7 +458,7 @@ int kind = 0; - ComplexMatrix z; + ComplexNDArray z; if (nargin > 1) { @@ -441,11 +477,11 @@ if (! error_state) { - z = args(nargin == 1 ? 0 : 1).complex_matrix_value (); + z = args(nargin == 1 ? 0 : 1).complex_array_value (); if (! error_state) { - Array2<int> ierr; + ArrayN<int> ierr; octave_value result; if (kind > 1) @@ -454,7 +490,7 @@ result = airy (z, kind == 1, scale, ierr); if (nargout > 1) - retval(1) = int_array2_to_matrix (ierr); + retval(1) = int_arrayN_to_array (ierr); retval(0) = result; }
--- a/src/DLD-FUNCTIONS/betainc.cc +++ b/src/DLD-FUNCTIONS/betainc.cc @@ -88,7 +88,7 @@ } else { - Matrix b = b_arg.matrix_value (); + NDArray b = b_arg.array_value (); if (! error_state) retval = betainc (x, a, b); @@ -97,7 +97,7 @@ } else { - Matrix a = a_arg.matrix_value (); + NDArray a = a_arg.array_value (); if (! error_state) { @@ -110,7 +110,7 @@ } else { - Matrix b = b_arg.matrix_value (); + NDArray b = b_arg.array_value (); if (! error_state) retval = betainc (x, a, b); @@ -120,7 +120,7 @@ } else { - Matrix x = x_arg.matrix_value (); + NDArray x = x_arg.array_value (); if (a_arg.is_scalar_type ()) { @@ -137,7 +137,7 @@ } else { - Matrix b = b_arg.matrix_value (); + NDArray b = b_arg.array_value (); if (! error_state) retval = betainc (x, a, b); @@ -146,7 +146,7 @@ } else { - Matrix a = a_arg.matrix_value (); + NDArray a = a_arg.array_value (); if (! error_state) { @@ -159,7 +159,7 @@ } else { - Matrix b = b_arg.matrix_value (); + NDArray b = b_arg.array_value (); if (! error_state) retval = betainc (x, a, b);
--- a/src/DLD-FUNCTIONS/filter.cc +++ b/src/DLD-FUNCTIONS/filter.cc @@ -39,34 +39,28 @@ #include "oct-obj.h" #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL) -extern MArray<double> -filter (MArray<double>&, MArray<double>&, MArray<double>&); +extern MArrayN<double> +filter (MArray<double>&, MArray<double>&, MArrayN<double>&, int dim); -extern MArray<Complex> -filter (MArray<Complex>&, MArray<Complex>&, MArray<Complex>&); +extern MArrayN<Complex> +filter (MArray<Complex>&, MArray<Complex>&, MArrayN<Complex>&, int dim); #endif template <class T> -MArray<T> -filter (MArray<T>& b, MArray<T>& a, MArray<T>& x, MArray<T>& si) +MArrayN<T> +filter (MArray<T>& b, MArray<T>& a, MArrayN<T>& x, MArrayN<T>& si, + int dim = 0) { - MArray<T> y; + MArrayN<T> y; int a_len = a.length (); int b_len = b.length (); - int x_len = x.length (); - - int si_len = si.length (); int ab_len = a_len > b_len ? a_len : b_len; b.resize (ab_len, 0.0); - - if (si.length () != ab_len - 1) - { - error ("filter: si must be a vector of length max (length (a), length (b)) - 1"); - return y; - } + if (a_len > 1) + a.resize (ab_len, 0.0); T norm = a (0); @@ -76,95 +70,183 @@ return y; } - y.resize (x_len, 0.0); + dim_vector x_dims = x.dims (); + if ((dim < 0) || (dim > x_dims.length ())) + { + error ("filter: filtering over invalid dimension"); + return y; + } + + int x_len = x_dims (dim); + + dim_vector si_dims = si.dims (); + int si_len = si_dims (0); + + if (si_len != ab_len - 1) + { + error ("filter: first dimension of si must be of length max (length (a), length (b)) - 1"); + return y; + } + + if (si_dims.length () != x_dims.length ()) + { + error ("filter: dimensionality of si and x must agree"); + return y; + } + + int si_dim = 0; + for (int i = 0; i < x_dims.length (); i++) + { + if (i == dim) + continue; + + if (si_dims (++si_dim) != x_dims (i)) + { + error ("filter: dimensionality of si and x must agree"); + return y; + } + } if (norm != 1.0) - b = b / norm; + { + a = a / norm; + b = b / norm; + } - if (a_len > 1) - { - a.resize (ab_len, 0.0); + if ((a_len <= 1) && (si_len <= 1)) + return b(0) * x; + + y.resize (x_dims, 0.0); + + int x_stride = 1; + for (int i = 0; i < dim; i++) + x_stride *= x_dims(i); - if (norm != 1.0) - a = a / norm; + int x_num = x_dims.numel () / x_len; + for (int num = 0; num < x_num; num++) + { + int x_offset; + if (x_stride == 1) + x_offset = num * x_len; + else + { + int x_offset2 = 0; + x_offset = num; + while (x_offset >= x_stride) + { + x_offset -= x_stride; + x_offset2++; + } + x_offset += x_offset2 * x_stride * x_len; + } + int si_offset = num * si_len; - for (int i = 0; i < x_len; i++) + if (a_len > 1) { - y (i) = si (0) + b (0) * x (i); - - if (si_len > 1) + for (int i = 0; i < x_len; i++) { - for (int j = 0; j < si_len - 1; j++) + int idx = i * x_stride + x_offset; + y (idx) = si (si_offset) + b (0) * x (idx); + + if (si_len > 1) { - OCTAVE_QUIT; + for (int j = 0; j < si_len - 1; j++) + { + OCTAVE_QUIT; - si (j) = si (j+1) - a (j+1) * y (i) - + b (j+1) * x (i); + si (j + si_offset) = si (j + 1 + si_offset) - + a (j+1) * y (idx) + b (j+1) * x (idx); + } + + si (si_len - 1 + si_offset) = b (si_len) * x (idx) + - a (si_len) * y (idx); } - - si (si_len-1) = b (si_len) * x (i) - - a (si_len) * y (i); + else + si (si_offset) = b (si_len) * x (idx) + - a (si_len) * y (idx); } - else - si (0) = b (si_len) * x (i) - - a (si_len) * y (i); + } + else if (si_len > 0) + { + for (int i = 0; i < x_len; i++) + { + int idx = i * x_stride + x_offset; + y (idx) = si (si_offset) + b (0) * x (idx); + + if (si_len > 1) + { + for (int j = 0; j < si_len - 1; j++) + { + OCTAVE_QUIT; + + si (j + si_offset) = si (j + 1 + si_offset) + + b (j+1) * x (idx); + } + + si (si_len - 1 + si_offset) = b (si_len) * x (idx); + } + else + si (si_offset) = b (1) * x (idx); + } } } - else if (si_len > 0) - { - for (int i = 0; i < x_len; i++) - { - y (i) = si (0) + b (0) * x (i); - - if (si_len > 1) - { - for (int j = 0; j < si_len - 1; j++) - { - OCTAVE_QUIT; - - si (j) = si (j+1) + b (j+1) * x (i); - } - - si (si_len-1) = b (si_len) * x (i); - } - else - si (0) = b (1) * x (i); - } - } - else - y = b (0) * x; - + return y; } #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL) -extern MArray<double> -filter (MArray<double>&, MArray<double>&, MArray<double>&, - MArray<double>&); +extern MArrayN<double> +filter (MArray<double>&, MArray<double>&, MArrayN<double>&, + MArrayN<double>&, int dim); -extern MArray<Complex> -filter (MArray<Complex>&, MArray<Complex>&, MArray<Complex>&, - MArray<Complex>&); +extern MArrayN<Complex> +filter (MArray<Complex>&, MArray<Complex>&, MArrayN<Complex>&, + MArrayN<Complex>&, int dim); #endif template <class T> -MArray<T> -filter (MArray<T>& b, MArray<T>& a, MArray<T>& x) +MArrayN<T> +filter (MArray<T>& b, MArray<T>& a, MArrayN<T>& x, int dim = -1) { + dim_vector x_dims = x.dims (); + + if (dim < 0) + { + // Find first non-singleton dimension + while ((dim < x_dims.length()) && (x_dims (dim) <= 1)) + dim++; + + // All dimensions singleton, pick first dimension + if (dim == x_dims.length ()) + dim = 0; + } + else + if (dim < 0 || dim > x_dims.length ()) + { + error ("filter: filtering over invalid dimension"); + return MArrayN<T> (); + } + int a_len = a.length (); int b_len = b.length (); int si_len = (a_len > b_len ? a_len : b_len) - 1; + dim_vector si_dims = x.dims (); + for (int i = dim; i > 0; i--) + si_dims (i) = si_dims (i-1); + si_dims (0) = si_len; + + MArrayN<T> si (si_dims, T (0.0)); - MArray<T> si (si_len, T (0.0)); - - return filter (b, a, x, si); + return filter (b, a, x, si, dim); } DEFUN_DLD (filter, args, nargout, "-*- texinfo -*-\n\ @deftypefn {Loadable Function} {y =} filter (@var{b}, @var{a}, @var{x})\n\ @deftypefnx {Loadable Function} {[@var{y}, @var{sf}] =} filter (@var{b}, @var{a}, @var{x}, @var{si})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{sf}] =} filter (@var{b}, @var{a}, @var{x}, [], @var{dim})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{sf}] =} filter (@var{b}, @var{a}, @var{x}, @var{si}, @var{dim})\n\ Return the solution to the following linear, time-invariant difference\n\ equation:\n\ @iftex\n\ @@ -194,7 +276,8 @@ $a \\in \\Re^{N-1}$, $b \\in \\Re^{M-1}$, and $x \\in \\Re^P$.\n\ @end tex\n\ @end iftex\n\ -An equivalent form of this equation is:\n\ +over the first non-singleton dimension of @var{x} or over @var{dim} if\n\ +supplied. An equivalent form of this equation is:\n\ @iftex\n\ @tex\n\ $$\n\ @@ -259,59 +342,91 @@ int nargin = args.length (); - if (nargin < 3 || nargin > 4) + if (nargin < 3 || nargin > 5) { print_usage ("filter"); return retval; } - const char *errmsg = "filter: arguments must be vectors"; + const char *errmsg = "filter: arguments a and b must be vectors"; + + int dim; + dim_vector x_dims = args(2).dims (); - bool x_is_row_vector = (args(2).rows () == 1); - - bool si_is_row_vector = (nargin == 4 && args(3).rows () == 1); + if (nargin == 5) + { + dim = args(4).nint_value() - 1; + if (dim < 0 || dim >= x_dims.length ()) + { + error ("filter: filtering over invalid dimension"); + return retval; + } + } + else + { + // Find first non-singleton dimension + dim = 0; + while ((dim < x_dims.length()) && (x_dims (dim) <= 1)) + dim++; + + // All dimensions singleton, pick first dimension + if (dim == x_dims.length ()) + dim = 0; + } if (args(0).is_complex_type () || args(1).is_complex_type () || args(2).is_complex_type () - || (nargin == 4 && args(3).is_complex_type ())) + || (nargin >= 4 && args(3).is_complex_type ())) { ComplexColumnVector b (args(0).complex_vector_value ()); ComplexColumnVector a (args(1).complex_vector_value ()); - ComplexColumnVector x (args(2).complex_vector_value ()); + + ComplexNDArray x (args(2).complex_array_value ()); if (! error_state) { - ComplexColumnVector si; + ComplexNDArray si; - if (nargin == 3) + if (nargin == 3 || args(3).is_empty ()) { int a_len = a.length (); int b_len = b.length (); int si_len = (a_len > b_len ? a_len : b_len) - 1; - si.resize (si_len, 0.0); + dim_vector si_dims = x.dims (); + for (int i = dim; i > 0; i--) + si_dims (i) = si_dims (i-1); + si_dims (0) = si_len; + + si.resize (si_dims, 0.0); } else - si = ComplexColumnVector (args(3).complex_vector_value ()); + { + dim_vector si_dims = args (3).dims (); + bool si_is_vector = true; + for (int i=0; i < si_dims.length (); i++) + if ((si_dims (i) != 1) && (si_dims (i) < si_dims.numel ())) + { + si_is_vector = false; + break; + } + + if (si_is_vector) + si = ComplexNDArray (args(3).complex_vector_value ()); + else + si = args(3).complex_array_value (); + } if (! error_state) { - ComplexColumnVector y (filter (b, a, x, si)); + ComplexNDArray y (filter (b, a, x, si, dim)); if (nargout == 2) - { - if (si_is_row_vector) - retval(1) = si.transpose (); - else - retval(1) = si; - } + retval(1) = si; - if (x_is_row_vector) - retval(0) = y.transpose (); - else - retval(0) = y; + retval(0) = y; } else error (errmsg); @@ -323,40 +438,52 @@ { ColumnVector b (args(0).vector_value ()); ColumnVector a (args(1).vector_value ()); - ColumnVector x (args(2).vector_value ()); + + NDArray x (args(2).array_value ()); if (! error_state) { - ColumnVector si; + NDArray si; - if (nargin == 3) + if (nargin == 3 || args(3).is_empty ()) { int a_len = a.length (); int b_len = b.length (); int si_len = (a_len > b_len ? a_len : b_len) - 1; - si.resize (si_len, 0.0); + dim_vector si_dims = x.dims (); + for (int i = dim; i > 0; i--) + si_dims (i) = si_dims (i-1); + si_dims (0) = si_len; + + si.resize (si_dims, 0.0); } else - si = ColumnVector (args(3).vector_value ()); + { + dim_vector si_dims = args (3).dims (); + bool si_is_vector = true; + for (int i=0; i < si_dims.length (); i++) + if ((si_dims (i) != 1) && (si_dims (i) < si_dims.numel ())) + { + si_is_vector = false; + break; + } + + if (si_is_vector) + si = NDArray (args(3).vector_value ()); + else + si = args(3).array_value (); + } if (! error_state) { - ColumnVector y (filter (b, a, x, si)); + NDArray y (filter (b, a, x, si, dim)); if (nargout == 2) - { - if (si_is_row_vector) - retval(1) = si.transpose (); - else - retval(1) = si; - } + retval(1) = si; - if (x_is_row_vector) - retval(0) = y.transpose (); - else - retval(0) = y; + retval(0) = y; } else error (errmsg); @@ -368,19 +495,19 @@ return retval; } -template MArray<double> -filter (MArray<double>&, MArray<double>&, MArray<double>&, - MArray<double>&); +template MArrayN<double> +filter (MArray<double>&, MArray<double>&, MArrayN<double>&, + MArrayN<double>&, int dim); -template MArray<double> -filter (MArray<double>&, MArray<double>&, MArray<double>&); +template MArrayN<double> +filter (MArray<double>&, MArray<double>&, MArrayN<double>&, int dim); -template MArray<Complex> -filter (MArray<Complex>&, MArray<Complex>&, MArray<Complex>&, - MArray<Complex>&); +template MArrayN<Complex> +filter (MArray<Complex>&, MArray<Complex>&, MArrayN<Complex>&, + MArrayN<Complex>&, int dim); -template MArray<Complex> -filter (MArray<Complex>&, MArray<Complex>&, MArray <Complex>&); +template MArrayN<Complex> +filter (MArray<Complex>&, MArray<Complex>&, MArrayN<Complex>&, int dim); /* ;;; Local Variables: ***
--- a/src/DLD-FUNCTIONS/gammainc.cc +++ b/src/DLD-FUNCTIONS/gammainc.cc @@ -86,7 +86,7 @@ } else { - Matrix a = a_arg.matrix_value (); + NDArray a = a_arg.array_value (); if (! error_state) retval = gammainc (x, a); @@ -95,7 +95,7 @@ } else { - Matrix x = x_arg.matrix_value (); + NDArray x = x_arg.array_value (); if (! error_state) { @@ -108,7 +108,7 @@ } else { - Matrix a = a_arg.matrix_value (); + NDArray a = a_arg.array_value (); if (! error_state) retval = gammainc (x, a);
--- a/src/DLD-FUNCTIONS/minmax.cc +++ b/src/DLD-FUNCTIONS/minmax.cc @@ -28,8 +28,8 @@ #include "lo-ieee.h" #include "lo-mappers.h" -#include "dMatrix.h" -#include "CMatrix.h" +#include "dNDArray.h" +#include "CNDArray.h" #include "quit.h" #include "defun-dld.h" @@ -37,13 +37,15 @@ #include "gripes.h" #include "oct-obj.h" +#include "ov-cx-mat.h" + #define MINMAX_BODY(FCN) \ \ octave_value_list retval; \ \ int nargin = args.length (); \ \ - if (nargin < 1 || nargin > 2 || nargout > 2) \ + if (nargin < 1 || nargin > 3 || nargout > 2) \ { \ print_usage (#FCN); \ return retval; \ @@ -51,9 +53,13 @@ \ octave_value arg1; \ octave_value arg2; \ + octave_value arg3; \ \ switch (nargin) \ { \ + case 3: \ + arg3 = args(2); \ + \ case 2: \ arg2 = args(1); \ \ @@ -66,97 +72,111 @@ break; \ } \ \ - if (nargin == 1 && (nargout == 1 || nargout == 0)) \ + int dim; \ + dim_vector dv = ((const octave_complex_matrix&) arg1) .dims (); \ + if (error_state) \ + { \ + gripe_wrong_type_arg (#FCN, arg1); \ + return retval; \ + } \ + \ + if (nargin == 3) \ + { \ + dim = arg3.nint_value () - 1; \ + if (dim < 0 || dim >= dv.length ()) \ + { \ + error ("%s: invalid dimension", #FCN); \ + return retval; \ + } \ + } \ + else \ + { \ + dim = 0; \ + while ((dim < dv.length ()) && (dv (dim) <= 1)) \ + dim++; \ + if (dim == dv.length ()) \ + dim = 0; \ + } \ + \ + bool single_arg = (nargin == 1) || arg2.is_empty(); \ + \ + if (single_arg) \ + { \ + dv(dim) = 1; \ + int n_dims = dv.length (); \ + for (int i = n_dims; i > 1; i--) \ + { \ + if (dv(i-1) == 1) \ + n_dims--; \ + else \ + break; \ + } \ + dv.resize (n_dims); \ + } \ + \ + if (single_arg && (nargout == 1 || nargout == 0)) \ { \ if (arg1.is_real_type ()) \ { \ - Matrix m = arg1.matrix_value (); \ + NDArray m = arg1.array_value (); \ \ if (! error_state) \ { \ - if (m.rows () == 1) \ - retval(0) = m.row_ ## FCN (); \ - else \ - { \ - if (m.rows () == 0 || m.columns () == 0) \ - retval(0) = Matrix (); \ - else \ - retval(0) = m.column_ ## FCN (); \ - } \ + NDArray n = m. FCN (dim); \ + n.resize (dv); \ + retval(0) = n; \ } \ } \ else if (arg1.is_complex_type ()) \ { \ - ComplexMatrix m = arg1.complex_matrix_value (); \ + ComplexNDArray m = arg1.complex_array_value (); \ \ if (! error_state) \ { \ - if (m.rows () == 1) \ - retval(0) = m.row_ ## FCN (); \ - else \ - { \ - if (m.rows () == 0 || m.columns () == 0) \ - retval(0) = Matrix (); \ - else \ - retval(0) = m.column_ ## FCN (); \ - } \ + ComplexNDArray n = m. FCN (dim); \ + n.resize (dv); \ + retval(0) = n; \ } \ } \ else \ gripe_wrong_type_arg (#FCN, arg1); \ } \ - else if (nargin == 1 && nargout == 2) \ + else if (single_arg && nargout == 2) \ { \ - Array<int> index; \ + ArrayN<int> index; \ \ if (arg1.is_real_type ()) \ { \ - Matrix m = arg1.matrix_value (); \ + NDArray m = arg1.array_value (); \ \ if (! error_state) \ { \ - retval.resize (2); \ - \ - if (m.rows () == 1) \ - retval(0) = m.row_ ## FCN (index); \ - else \ - { \ - if (m.rows () == 0 || m.columns () == 0) \ - retval(0) = Matrix (); \ - else \ - retval(0) = m.column_ ## FCN (index); \ - } \ + NDArray n = m. FCN (index, dim); \ + n.resize (dv); \ + retval(0) = n; \ } \ } \ else if (arg1.is_complex_type ()) \ { \ - ComplexMatrix m = arg1.complex_matrix_value (); \ + ComplexNDArray m = arg1.complex_array_value (); \ \ if (! error_state) \ { \ - retval.resize (2); \ - \ - if (m.rows () == 1) \ - retval(0) = m.row_ ## FCN (index); \ - else \ - { \ - if (m.rows () == 0 || m.columns () == 0) \ - retval(0) = Matrix (); \ - else \ - retval(0) = m.column_ ## FCN (index); \ - } \ + ComplexNDArray n = m. FCN (index, dim); \ + n.resize (dv); \ + retval(0) = n; \ } \ } \ else \ gripe_wrong_type_arg (#FCN, arg1); \ \ - int len = index.length (); \ + int len = index.numel (); \ \ if (len > 0) \ { \ double nan_val = lo_ieee_nan_value (); \ \ - RowVector idx (len); \ + NDArray idx (index.dims ()); \ \ for (int i = 0; i < len; i++) \ { \ @@ -169,9 +189,9 @@ retval(1) = idx; \ } \ else \ - retval(1) = Matrix (); \ + retval(1) = NDArray (); \ } \ - else if (nargin == 2) \ + else \ { \ int arg1_is_scalar = arg1.is_scalar_type (); \ int arg2_is_scalar = arg2.is_scalar_type (); \ @@ -184,10 +204,10 @@ if (arg1_is_complex || arg2_is_complex) \ { \ Complex c1 = arg1.complex_value (); \ - ComplexMatrix m2 = arg2.complex_matrix_value (); \ + ComplexNDArray m2 = arg2.complex_array_value (); \ if (! error_state) \ { \ - ComplexMatrix result = FCN (c1, m2); \ + ComplexNDArray result = FCN (c1, m2); \ if (! error_state) \ retval(0) = result; \ } \ @@ -195,11 +215,11 @@ else \ { \ double d1 = arg1.double_value (); \ - Matrix m2 = arg2.matrix_value (); \ + NDArray m2 = arg2.array_value (); \ \ if (! error_state) \ { \ - Matrix result = FCN (d1, m2); \ + NDArray result = FCN (d1, m2); \ if (! error_state) \ retval(0) = result; \ } \ @@ -209,24 +229,24 @@ { \ if (arg1_is_complex || arg2_is_complex) \ { \ - ComplexMatrix m1 = arg1.complex_matrix_value (); \ + ComplexNDArray m1 = arg1.complex_array_value (); \ \ if (! error_state) \ { \ Complex c2 = arg2.complex_value (); \ - ComplexMatrix result = FCN (m1, c2); \ + ComplexNDArray result = FCN (m1, c2); \ if (! error_state) \ retval(0) = result; \ } \ } \ else \ { \ - Matrix m1 = arg1.matrix_value (); \ + NDArray m1 = arg1.array_value (); \ \ if (! error_state) \ { \ double d2 = arg2.double_value (); \ - Matrix result = FCN (m1, d2); \ + NDArray result = FCN (m1, d2); \ if (! error_state) \ retval(0) = result; \ } \ @@ -236,15 +256,15 @@ { \ if (arg1_is_complex || arg2_is_complex) \ { \ - ComplexMatrix m1 = arg1.complex_matrix_value (); \ + ComplexNDArray m1 = arg1.complex_array_value (); \ \ if (! error_state) \ { \ - ComplexMatrix m2 = arg2.complex_matrix_value (); \ + ComplexNDArray m2 = arg2.complex_array_value (); \ \ if (! error_state) \ { \ - ComplexMatrix result = FCN (m1, m2); \ + ComplexNDArray result = FCN (m1, m2); \ if (! error_state) \ retval(0) = result; \ } \ @@ -252,15 +272,15 @@ } \ else \ { \ - Matrix m1 = arg1.matrix_value (); \ + NDArray m1 = arg1.array_value (); \ \ if (! error_state) \ { \ - Matrix m2 = arg2.matrix_value (); \ + NDArray m2 = arg2.array_value (); \ \ if (! error_state) \ { \ - Matrix result = FCN (m1, m2); \ + NDArray result = FCN (m1, m2); \ if (! error_state) \ retval(0) = result; \ } \ @@ -268,21 +288,18 @@ } \ } \ } \ - else \ - panic_impossible (); \ \ return retval DEFUN_DLD (min, args, nargout, "-*- texinfo -*-\n\ -@deftypefn {Mapping Function} {} min (@var{x}, @var{y})\n\ +@deftypefn {Mapping Function} {} min (@var{x}, @var{y}, @var{dim})\n\ @deftypefnx {Mapping Function} {[@var{w}, @var{iw}] =} min (@var{x})\n\ @cindex Utility Functions\n\ For a vector argument, return the minimum value. For a matrix\n\ argument, return the minimum value from each column, as a row\n\ -vector.\n\ -For two matrices (or a matrix and scalar),\n\ -return the pair-wise minimum.\n\ +vector, or over the dimension @var{dim} if defined. For two matrices\n\ +(or a matrix and scalar), return the pair-wise minimum.\n\ Thus,\n\ \n\ @example\n\ @@ -323,14 +340,13 @@ DEFUN_DLD (max, args, nargout, "-*- texinfo -*-\n\ -@deftypefn {Mapping Function} {} max (@var{x}, @var{y})\n\ +@deftypefn {Mapping Function} {} max (@var{x}, @var{y}, @var{dim})\n\ @deftypefnx {Mapping Function} {[@var{w}, @var{iw}] =} max (@var{x})\n\ @cindex Utility Functions\n\ For a vector argument, return the maximum value. For a matrix\n\ argument, return the maximum value from each column, as a row\n\ -vector.\n\ -For two matrices (or a matrix and scalar),\n\ -return the pair-wise maximum.\n\ +vector, or over the dimension @var{dim} if defined. For two matrices\n\ +(or a matrix and scalar), return the pair-wise maximum.\n\ Thus,\n\ \n\ @example\n\
--- a/src/load-save.cc +++ b/src/load-save.cc @@ -306,6 +306,13 @@ { load_save_format retval = LS_UNKNOWN; + // If the file doesn't exist do nothing + std::ifstream file_exist (fname.c_str ()); + if (file_exist) + file_exist.close (); + else + return LS_UNKNOWN; + #ifdef HAVE_HDF5 // check this before we open the file if (H5Fis_hdf5 (fname.c_str ()) > 0) @@ -722,19 +729,26 @@ { i++; - hdf5_ifstream hdf5_file (fname.c_str ()); - - if (hdf5_file.file_id >= 0) + // If the file doesn't exist do nothing + std::ifstream file (fname.c_str (), std::ios::in); + if (file) { - retval = do_load (hdf5_file, orig_fname, force, format, - flt_fmt, list_only, swap, verbose, - argv, i, argc, nargout); + file.close (); + + hdf5_ifstream hdf5_file (fname.c_str ()); - hdf5_file.close (); + if (hdf5_file.file_id >= 0) + { + retval = do_load (hdf5_file, orig_fname, force, format, + flt_fmt, list_only, swap, verbose, + argv, i, argc, nargout); + + hdf5_file.close (); + } + else + error ("load: couldn't open input file `%s'", + orig_fname.c_str ()); } - else - error ("load: couldn't open input file `%s'", - orig_fname.c_str ()); } else #endif /* HAVE_HDF5 */
--- a/src/ov-bool-mat.cc +++ b/src/ov-bool-mat.cc @@ -142,7 +142,7 @@ octave_value octave_bool_matrix::convert_to_str_internal (bool pad, bool force) const { - octave_value tmp = octave_value (matrix_value ()); + octave_value tmp = octave_value (array_value ()); return tmp.convert_to_str (pad, force); }
--- a/src/ov-re-mat.cc +++ b/src/ov-re-mat.cc @@ -186,66 +186,54 @@ octave_matrix::convert_to_str_internal (bool, bool) const { octave_value retval; + dim_vector dv = dims (); + int nel = dv.numel (); - int nr = matrix.rows (); - int nc = matrix.columns (); - - if (nr == 0 && nc == 0) + if (nel == 0) { char s = '\0'; retval = octave_value (&s); } else { - if (nr == 0 || nc == 0) - { - char s = '\0'; - retval = octave_value (&s); - } - else - { - charMatrix chm (nr, nc); + charNDArray chm (dv); - bool warned = false; + bool warned = false; + + for (int i = 0; i < nel; i++) + { + OCTAVE_QUIT; + + double d = matrix (i); - for (int j = 0; j < nc; j++) + if (xisnan (d)) { - for (int i = 0; i < nr; i++) + ::error ("invalid conversion from NaN to character"); + return retval; + } + else + { + int ival = NINT (d); + + if (ival < 0 || ival > UCHAR_MAX) { - OCTAVE_QUIT; + // XXX FIXME XXX -- is there something + // better we could do? - double d = matrix (i, j); + ival = 0; - if (xisnan (d)) + if (! warned) { - ::error ("invalid conversion from NaN to character"); - return retval; - } - else - { - int ival = NINT (d); - - if (ival < 0 || ival > UCHAR_MAX) - { - // XXX FIXME XXX -- is there something - // better we could do? - - ival = 0; - - if (! warned) - { - ::warning ("range error for conversion to character value"); - warned = true; - } - } - - chm (i, j) = static_cast<char> (ival); + ::warning ("range error for conversion to character value"); + warned = true; } } - } - retval = octave_value (chm, 1); + chm (i) = static_cast<char> (ival); + } } + + retval = octave_value (chm, 1); } return retval;
--- a/src/parse.y +++ b/src/parse.y @@ -2810,10 +2810,10 @@ else { if (reading_script_file) - warning ("ignoring static declaration near line %d of file `%s'", + warning ("ignoring persistent declaration near line %d of file `%s'", l, curr_fcn_file_full_name.c_str ()); else - warning ("ignoring static declaration near line %d", l); + warning ("ignoring persistent declaration near line %d", l); } break;
--- a/src/pt-decl.cc +++ b/src/pt-decl.cc @@ -147,13 +147,18 @@ { id->mark_as_static (); - tree_expression *expr = elt.expression (); + octave_lvalue ult = id->lvalue (); - if (expr) + if (ult.is_defined ()) { - octave_value init_val = expr->rvalue (); + tree_expression *expr = elt.expression (); + + octave_value init_val; - octave_lvalue ult = id->lvalue (); + if (expr) + init_val = expr->rvalue (); + else + init_val = Matrix (); ult.assign (octave_value::op_asn_eq, init_val); }
--- a/test/octave.test/arith/max-4.m +++ b/test/octave.test/arith/max-4.m @@ -1,1 +1,1 @@ -max (1, 2, 3) +max (1, 2, 3, 4)