Mercurial > octave-nkf
changeset 14970:0f156affb036
maint: periodic merge of default to jit
| author | Max Brister <max@2bass.com> |
|---|---|
| date | Mon, 25 Jun 2012 14:40:21 -0500 |
| parents | bbeef7b8ea2e (current diff) 3c5e6971064c (diff) |
| children | b23a98ca0e43 |
| files | |
| diffstat | 11 files changed, 445 insertions(+), 236 deletions(-) [+] |
line wrap: on
line diff
--- a/examples/embedded.cc Mon Jun 25 14:21:45 2012 -0500 +++ b/examples/embedded.cc Mon Jun 25 14:40:21 2012 -0500 @@ -13,30 +13,21 @@ octave_main (2, argv.c_str_vec(), 1); octave_idx_type n = 2; - Matrix a_matrix = Matrix (1, 2); + octave_value_list in; - std::cout << "GCD of ["; - for (octave_idx_type i = 0; i < n; i++) - { - a_matrix (i) = 5 * (i + 1); - if (i != 0) - std::cout << ", " << 5 * (i + 2); - else - std::cout << 5 * (i + 2); - } - std::cout << "] is "; - - octave_value_list in = octave_value (a_matrix); + for (octave_idx_type i = 0; i < n; i++) + in(i) = octave_value (5 * (i + 1)); + octave_value_list out = feval ("gcd", in, 1); + if (!error_state && out.length () > 0) - { - a_matrix = out(0).matrix_value (); - if (a_matrix.numel () == 1) - std::cout << a_matrix(0) << "\n"; - else - std::cout << "invalid\n"; - } + std::cout << "GCD of [" + << in(0).int_value () + << ", " + << in(1).int_value () + << "] is " << out(0).int_value () + << std::endl; else std::cout << "invalid\n";
--- a/liboctave/Sparse.cc Mon Jun 25 14:21:45 2012 -0500 +++ b/liboctave/Sparse.cc Mon Jun 25 14:40:21 2012 -0500 @@ -1240,15 +1240,31 @@ gripe_del_index_out_of_range (false, idx_j.extent (nc), nc); else if (idx_j.is_cont_range (nc, lb, ub)) { - const Sparse<T> tmp = *this; - octave_idx_type lbi = tmp.cidx(lb), ubi = tmp.cidx(ub), new_nz = nz - (ubi - lbi); - *this = Sparse<T> (nr, nc - (ub - lb), new_nz); - copy_or_memcpy (lbi, tmp.data (), data ()); - copy_or_memcpy (lbi, tmp.ridx (), ridx ()); - copy_or_memcpy (nz - ubi, tmp.data () + ubi, xdata () + lbi); - copy_or_memcpy (nz - ubi, tmp.ridx () + ubi, xridx () + lbi); - copy_or_memcpy (lb, tmp.cidx () + 1, cidx () + 1); - mx_inline_sub (nc - ub, xcidx () + 1, tmp.cidx () + ub + 1, ubi - lbi); + if (lb == 0 && ub == nc) + { + // Delete all rows and columns. + *this = Sparse<T> (nr, 0); + } + else if (nz == 0) + { + // No elements to preserve; adjust dimensions. + *this = Sparse<T> (nr, nc - (ub - lb)); + } + else + { + const Sparse<T> tmp = *this; + octave_idx_type lbi = tmp.cidx(lb), ubi = tmp.cidx(ub), + new_nz = nz - (ubi - lbi); + + *this = Sparse<T> (nr, nc - (ub - lb), new_nz); + copy_or_memcpy (lbi, tmp.data (), data ()); + copy_or_memcpy (lbi, tmp.ridx (), ridx ()); + copy_or_memcpy (nz - ubi, tmp.data () + ubi, xdata () + lbi); + copy_or_memcpy (nz - ubi, tmp.ridx () + ubi, xridx () + lbi); + copy_or_memcpy (lb, tmp.cidx () + 1, cidx () + 1); + mx_inline_sub (nc - ub, xcidx () + lb + 1, + tmp.cidx () + ub + 1, ubi - lbi); + } } else *this = index (idx_i, idx_j.complement (nc)); @@ -1261,24 +1277,40 @@ gripe_del_index_out_of_range (false, idx_i.extent (nr), nr); else if (idx_i.is_cont_range (nr, lb, ub)) { - // This is more memory-efficient than the approach below. - const Sparse<T> tmpl = index (idx_vector (0, lb), idx_j); - const Sparse<T> tmpu = index (idx_vector (ub, nr), idx_j); - *this = Sparse<T> (nr - (ub - lb), nc, tmpl.nnz () + tmpu.nnz ()); - for (octave_idx_type j = 0, k = 0; j < nc; j++) + if (lb == 0 && ub == nr) + { + // Delete all rows and columns. + *this = Sparse<T> (0, nc); + } + else if (nz == 0) { - for (octave_idx_type i = tmpl.cidx(j); i < tmpl.cidx(j+1); i++) + // No elements to preserve; adjust dimensions. + *this = Sparse<T> (nr - (ub - lb), nc); + } + else + { + // This is more memory-efficient than the approach below. + const Sparse<T> tmpl = index (idx_vector (0, lb), idx_j); + const Sparse<T> tmpu = index (idx_vector (ub, nr), idx_j); + *this = Sparse<T> (nr - (ub - lb), nc, + tmpl.nnz () + tmpu.nnz ()); + for (octave_idx_type j = 0, k = 0; j < nc; j++) { - xdata(k) = tmpl.data(i); - xridx(k++) = tmpl.ridx(i); + for (octave_idx_type i = tmpl.cidx(j); i < tmpl.cidx(j+1); + i++) + { + xdata(k) = tmpl.data(i); + xridx(k++) = tmpl.ridx(i); + } + for (octave_idx_type i = tmpu.cidx(j); i < tmpu.cidx(j+1); + i++) + { + xdata(k) = tmpu.data(i); + xridx(k++) = tmpu.ridx(i) + lb; + } + + xcidx(j+1) = k; } - for (octave_idx_type i = tmpu.cidx(j); i < tmpu.cidx(j+1); i++) - { - xdata(k) = tmpu.data(i); - xridx(k++) = tmpu.ridx(i) + lb; - } - - xcidx(j+1) = k; } } else @@ -1583,6 +1615,10 @@ { // It's actually vector indexing. The 1D index is specialized for that. retval = index (idx_i); + + // If nr == 1 then the vector indexing will return a column vector!! + if (nr == 1) + retval.transpose(); } else if (idx_i.is_scalar ()) { @@ -2736,6 +2772,18 @@ %!assert (speye (3,1)(3:-1:1), sparse ([0; 0; 1])) +## Test removing columns (bug #36656) + +%!test +%! s = sparse (magic (5)); +%! s(:,2:4) = []; +%! assert (s, sparse (magic (5)(:, [1,5]))); + +%!test +%! s = sparse([], [], [], 1, 1); +%! s(1,:) = []; +%! assert (s, sparse ([], [], [], 0, 1)); + */ template <class T>
--- a/liboctave/Sparse.h Mon Jun 25 14:21:45 2012 -0500 +++ b/liboctave/Sparse.h Mon Jun 25 14:40:21 2012 -0500 @@ -71,33 +71,34 @@ octave_idx_type ncols; octave_refcount<int> count; - SparseRep (void) : d (0), r (0), c (new octave_idx_type [1]), nzmx (0), nrows (0), - ncols (0), count (1) { c[0] = 0; } + SparseRep (void) + : d (0), r (0), c (new octave_idx_type [1]), nzmx (0), nrows (0), + ncols (0), count (1) + { + c[0] = 0; + } - SparseRep (octave_idx_type n) : d (0), r (0), c (new octave_idx_type [n+1]), nzmx (0), nrows (n), + SparseRep (octave_idx_type n) + : d (0), r (0), c (new octave_idx_type [n+1]), nzmx (0), nrows (n), ncols (n), count (1) { for (octave_idx_type i = 0; i < n + 1; i++) c[i] = 0; } - SparseRep (octave_idx_type nr, octave_idx_type nc) : d (0), r (0), c (new octave_idx_type [nc+1]), nzmx (0), - nrows (nr), ncols (nc), count (1) - { - for (octave_idx_type i = 0; i < nc + 1; i++) - c[i] = 0; - } - - SparseRep (octave_idx_type nr, octave_idx_type nc, octave_idx_type nz) : d (new T [nz]), - r (new octave_idx_type [nz]), c (new octave_idx_type [nc+1]), nzmx (nz), nrows (nr), + SparseRep (octave_idx_type nr, octave_idx_type nc, octave_idx_type nz = 0) + : d (new T [nz]), r (new octave_idx_type [nz]), + c (new octave_idx_type [nc+1]), nzmx (nz), nrows (nr), ncols (nc), count (1) { - for (octave_idx_type i = 0; i < nc + 1; i++) + c[nc] = nz; + for (octave_idx_type i = 0; i < nc; i++) c[i] = 0; } SparseRep (const SparseRep& a) - : d (new T [a.nzmx]), r (new octave_idx_type [a.nzmx]), c (new octave_idx_type [a.ncols + 1]), + : d (new T [a.nzmx]), r (new octave_idx_type [a.nzmx]), + c (new octave_idx_type [a.ncols + 1]), nzmx (a.nzmx), nrows (a.nrows), ncols (a.ncols), count (1) { octave_idx_type nz = a.nnz (); @@ -144,17 +145,17 @@ //-------------------------------------------------------------------- void make_unique (void) - { - if (rep->count > 1) - { - SparseRep *r = new SparseRep (*rep); + { + if (rep->count > 1) + { + SparseRep *r = new SparseRep (*rep); - if (--rep->count == 0) - delete rep; + if (--rep->count == 0) + delete rep; - rep = r; - } - } + rep = r; + } + } public: @@ -168,18 +169,18 @@ private: typename Sparse<T>::SparseRep *nil_rep (void) const - { - static typename Sparse<T>::SparseRep nr; - return &nr; - } + { + static typename Sparse<T>::SparseRep nr; + return &nr; + } public: Sparse (void) : rep (nil_rep ()), dimensions (dim_vector(0,0)) - { - rep->count++; - } + { + rep->count++; + } explicit Sparse (octave_idx_type n) : rep (new typename Sparse<T>::SparseRep (n)), @@ -193,7 +194,7 @@ Sparse (const dim_vector& dv, octave_idx_type nz) : rep (new typename Sparse<T>::SparseRep (dv(0), dv(1), nz)), - dimensions (dv) { } + dimensions (dv) { } Sparse (octave_idx_type nr, octave_idx_type nc, octave_idx_type nz) : rep (new typename Sparse<T>::SparseRep (nr, nc, nz)), @@ -207,20 +208,20 @@ template <class U> Sparse (const Sparse<U>& a) : rep (new typename Sparse<T>::SparseRep (a.rep->nrows, a.rep->ncols, a.rep->nzmx)), - dimensions (a.dimensions) - { - octave_idx_type nz = a.nnz (); - std::copy (a.rep->d, a.rep->d + nz, rep->d); - copy_or_memcpy (nz, a.rep->r, rep->r); - copy_or_memcpy (rep->ncols + 1, a.rep->c, rep->c); - } + dimensions (a.dimensions) + { + octave_idx_type nz = a.nnz (); + std::copy (a.rep->d, a.rep->d + nz, rep->d); + copy_or_memcpy (nz, a.rep->r, rep->r); + copy_or_memcpy (rep->ncols + 1, a.rep->c, rep->c); + } // No type conversion case. Sparse (const Sparse<T>& a) : rep (a.rep), dimensions (a.dimensions) - { - rep->count++; - } + { + rep->count++; + } public: @@ -249,9 +250,9 @@ // Querying the number of elements (incl. zeros) may overflow the index type, // so don't do it unless you really need it. octave_idx_type numel (void) const - { - return dimensions.safe_numel (); - } + { + return dimensions.safe_numel (); + } octave_idx_type nelem (void) const { return capacity (); } octave_idx_type length (void) const { return numel (); } @@ -265,18 +266,19 @@ octave_idx_type get_row_index (octave_idx_type k) { return ridx (k); } octave_idx_type get_col_index (octave_idx_type k) - { - octave_idx_type ret = 0; - while (cidx(ret+1) < k) - ret++; - return ret; - } + { + octave_idx_type ret = 0; + while (cidx(ret+1) < k) + ret++; + return ret; + } size_t byte_size (void) const - { - return (static_cast<size_t>(cols () + 1) * sizeof (octave_idx_type) - + static_cast<size_t> (capacity ()) * (sizeof (T) + sizeof (octave_idx_type))); - } + { + return (static_cast<size_t>(cols () + 1) * sizeof (octave_idx_type) + + static_cast<size_t> (capacity ()) + * (sizeof (T) + sizeof (octave_idx_type))); + } dim_vector dims (void) const { return dimensions; } @@ -296,145 +298,189 @@ // No checking, even for multiple references, ever. T& xelem (octave_idx_type n) - { - octave_idx_type i = n % rows (), j = n / rows(); - return xelem (i, j); - } + { + octave_idx_type i = n % rows (), j = n / rows(); + return xelem (i, j); + } T xelem (octave_idx_type n) const - { - octave_idx_type i = n % rows (), j = n / rows(); - return xelem (i, j); - } + { + octave_idx_type i = n % rows (), j = n / rows(); + return xelem (i, j); + } T& xelem (octave_idx_type i, octave_idx_type j) { return rep->elem (i, j); } - T xelem (octave_idx_type i, octave_idx_type j) const { return rep->celem (i, j); } + T xelem (octave_idx_type i, octave_idx_type j) const + { + return rep->celem (i, j); + } T& xelem (const Array<octave_idx_type>& ra_idx) - { return xelem (compute_index (ra_idx)); } + { return xelem (compute_index (ra_idx)); } T xelem (const Array<octave_idx_type>& ra_idx) const - { return xelem (compute_index (ra_idx)); } + { return xelem (compute_index (ra_idx)); } // FIXME -- would be nice to fix this so that we don't // unnecessarily force a copy, but that is not so easy, and I see no // clean way to do it. T& checkelem (octave_idx_type n) - { - if (n < 0 || n >= numel ()) - return range_error ("T& Sparse<T>::checkelem", n); - else - { - make_unique (); - return xelem (n); - } - } + { + if (n < 0 || n >= numel ()) + return range_error ("T& Sparse<T>::checkelem", n); + else + { + make_unique (); + return xelem (n); + } + } T& checkelem (octave_idx_type i, octave_idx_type j) - { - if (i < 0 || j < 0 || i >= dim1 () || j >= dim2 ()) - return range_error ("T& Sparse<T>::checkelem", i, j); - else - { - make_unique (); - return xelem (i, j); - } - } + { + if (i < 0 || j < 0 || i >= dim1 () || j >= dim2 ()) + return range_error ("T& Sparse<T>::checkelem", i, j); + else + { + make_unique (); + return xelem (i, j); + } + } T& checkelem (const Array<octave_idx_type>& ra_idx) - { - octave_idx_type i = compute_index (ra_idx); + { + octave_idx_type i = compute_index (ra_idx); - if (i < 0) - return range_error ("T& Sparse<T>::checkelem", ra_idx); - else - return elem (i); - } + if (i < 0) + return range_error ("T& Sparse<T>::checkelem", ra_idx); + else + return elem (i); + } T& elem (octave_idx_type n) - { - make_unique (); - return xelem (n); - } + { + make_unique (); + return xelem (n); + } T& elem (octave_idx_type i, octave_idx_type j) - { - make_unique (); - return xelem (i, j); - } + { + make_unique (); + return xelem (i, j); + } T& elem (const Array<octave_idx_type>& ra_idx) - { return Sparse<T>::elem (compute_index (ra_idx)); } + { return Sparse<T>::elem (compute_index (ra_idx)); } #if defined (BOUNDS_CHECKING) - T& operator () (octave_idx_type n) { return checkelem (n); } - T& operator () (octave_idx_type i, octave_idx_type j) { return checkelem (i, j); } - T& operator () (const Array<octave_idx_type>& ra_idx) { return checkelem (ra_idx); } + T& operator () (octave_idx_type n) + { + return checkelem (n); + } + + T& operator () (octave_idx_type i, octave_idx_type j) + { + return checkelem (i, j); + } + + T& operator () (const Array<octave_idx_type>& ra_idx) + { + return checkelem (ra_idx); + } + #else - T& operator () (octave_idx_type n) { return elem (n); } - T& operator () (octave_idx_type i, octave_idx_type j) { return elem (i, j); } - T& operator () (const Array<octave_idx_type>& ra_idx) { return elem (ra_idx); } + T& operator () (octave_idx_type n) + { + return elem (n); + } + + T& operator () (octave_idx_type i, octave_idx_type j) + { + return elem (i, j); + } + + T& operator () (const Array<octave_idx_type>& ra_idx) + { + return elem (ra_idx); + } + #endif T checkelem (octave_idx_type n) const - { - if (n < 0 || n >= numel ()) - return range_error ("T Sparse<T>::checkelem", n); - else - return xelem (n); - } + { + if (n < 0 || n >= numel ()) + return range_error ("T Sparse<T>::checkelem", n); + else + return xelem (n); + } T checkelem (octave_idx_type i, octave_idx_type j) const - { - if (i < 0 || j < 0 || i >= dim1 () || j >= dim2 ()) - return range_error ("T Sparse<T>::checkelem", i, j); - else - return xelem (i, j); - } + { + if (i < 0 || j < 0 || i >= dim1 () || j >= dim2 ()) + return range_error ("T Sparse<T>::checkelem", i, j); + else + return xelem (i, j); + } T checkelem (const Array<octave_idx_type>& ra_idx) const - { - octave_idx_type i = compute_index (ra_idx); + { + octave_idx_type i = compute_index (ra_idx); - if (i < 0) - return range_error ("T Sparse<T>::checkelem", ra_idx); - else - return Sparse<T>::elem (i); - } + if (i < 0) + return range_error ("T Sparse<T>::checkelem", ra_idx); + else + return Sparse<T>::elem (i); + } T elem (octave_idx_type n) const { return xelem (n); } T elem (octave_idx_type i, octave_idx_type j) const { return xelem (i, j); } T elem (const Array<octave_idx_type>& ra_idx) const - { return Sparse<T>::elem (compute_index (ra_idx)); } + { return Sparse<T>::elem (compute_index (ra_idx)); } #if defined (BOUNDS_CHECKING) T operator () (octave_idx_type n) const { return checkelem (n); } - T operator () (octave_idx_type i, octave_idx_type j) const { return checkelem (i, j); } - T operator () (const Array<octave_idx_type>& ra_idx) const { return checkelem (ra_idx); } + T operator () (octave_idx_type i, octave_idx_type j) const + { + return checkelem (i, j); + } + + T operator () (const Array<octave_idx_type>& ra_idx) const + { + return checkelem (ra_idx); + } + #else T operator () (octave_idx_type n) const { return elem (n); } - T operator () (octave_idx_type i, octave_idx_type j) const { return elem (i, j); } - T operator () (const Array<octave_idx_type>& ra_idx) const { return elem (ra_idx); } + T operator () (octave_idx_type i, octave_idx_type j) const + { + return elem (i, j); + } + + T operator () (const Array<octave_idx_type>& ra_idx) const + { + return elem (ra_idx); + } #endif Sparse<T> maybe_compress (bool remove_zeros = false) - { - if (remove_zeros) - make_unique (); // Needs to unshare because elements are removed. + { + if (remove_zeros) + make_unique (); // Needs to unshare because elements are removed. - rep->maybe_compress (remove_zeros); - return (*this); - } + rep->maybe_compress (remove_zeros); + return (*this); + } Sparse<T> reshape (const dim_vector& new_dims) const; Sparse<T> permute (const Array<octave_idx_type>& vec, bool inv = false) const; Sparse<T> ipermute (const Array<octave_idx_type>& vec) const - { return permute (vec, true); } + { + return permute (vec, true); + } void resize1 (octave_idx_type n); @@ -443,11 +489,11 @@ void resize (const dim_vector& dv); void change_capacity (octave_idx_type nz) - { - if (nz < nnz ()) - make_unique (); // Unshare now because elements will be truncated. - rep->change_length (nz); - } + { + if (nz < nnz ()) + make_unique (); // Unshare now because elements will be truncated. + rep->change_length (nz); + } Sparse<T>& insert (const Sparse<T>& a, octave_idx_type r, octave_idx_type c); Sparse<T>& insert (const Sparse<T>& a, const Array<octave_idx_type>& idx); @@ -468,7 +514,11 @@ T* data (void) const { return rep->d; } octave_idx_type* ridx (void) { make_unique (); return rep->r; } - octave_idx_type& ridx (octave_idx_type i) { make_unique (); return rep->ridx (i); } + octave_idx_type& ridx (octave_idx_type i) + { + make_unique (); return rep->ridx (i); + } + octave_idx_type* xridx (void) { return rep->r; } octave_idx_type& xridx (octave_idx_type i) { return rep->ridx (i); } @@ -477,7 +527,11 @@ octave_idx_type* ridx (void) const { return rep->r; } octave_idx_type* cidx (void) { make_unique (); return rep->c; } - octave_idx_type& cidx (octave_idx_type i) { make_unique (); return rep->cidx (i); } + octave_idx_type& cidx (octave_idx_type i) + { + make_unique (); return rep->cidx (i); + } + octave_idx_type* xcidx (void) { return rep->c; } octave_idx_type& xcidx (octave_idx_type i) { return rep->cidx (i); } @@ -495,7 +549,8 @@ Sparse<T> index (const idx_vector& i, bool resize_ok = false) const; - Sparse<T> index (const idx_vector& i, const idx_vector& j, bool resize_ok = false) const; + Sparse<T> index (const idx_vector& i, const idx_vector& j, + bool resize_ok = false) const; void assign (const idx_vector& i, const Sparse<T>& rhs); @@ -507,13 +562,19 @@ // You should not use them anywhere else. void *mex_get_data (void) const { return const_cast<T *> (data ()); } - octave_idx_type *mex_get_ir (void) const { return const_cast<octave_idx_type *> (ridx ()); } + octave_idx_type *mex_get_ir (void) const + { + return const_cast<octave_idx_type *> (ridx ()); + } - octave_idx_type *mex_get_jc (void) const { return const_cast<octave_idx_type *> (cidx ()); } + octave_idx_type *mex_get_jc (void) const + { + return const_cast<octave_idx_type *> (cidx ()); + } Sparse<T> sort (octave_idx_type dim = 0, sortmode mode = ASCENDING) const; Sparse<T> sort (Array<octave_idx_type> &sidx, octave_idx_type dim = 0, - sortmode mode = ASCENDING) const; + sortmode mode = ASCENDING) const; Sparse<T> diag (octave_idx_type k = 0) const;
--- a/liboctave/lo-specfun.cc Mon Jun 25 14:21:45 2012 -0500 +++ b/liboctave/lo-specfun.cc Mon Jun 25 14:40:21 2012 -0500 @@ -3152,7 +3152,7 @@ // This algorithm is due to P. J. Acklam. // See http://home.online.no/~pjacklam/notes/invnorm/ // The rational approximation has relative accuracy 1.15e-9 in the whole region. -// For doubles, it is refined by a single step of Higham's 3rd order method. +// For doubles, it is refined by a single step of Halley's 3rd order method. // For single precision, the accuracy is already OK, so we skip it to get // faster evaluation. @@ -3175,7 +3175,7 @@ { 7.784695709041462e-03, 3.224671290700398e-01, 2.445134137142996e+00, 3.754408661907416e+00 }; - static const double spi2 = 8.862269254527579e-01; // sqrt(pi)/2. + static const double spi2 = 8.862269254527579e-01; // sqrt(pi)/2. static const double pbreak = 0.95150; double ax = fabs (x), y; @@ -3204,7 +3204,7 @@ if (refine) { // One iteration of Halley's method gives full precision. - double u = (erf(y) - x) * spi2 * exp (y*y); + double u = (erf (y) - x) * spi2 * exp (y*y); y -= u / (1 + y*u); } @@ -3221,6 +3221,10 @@ return do_erfinv (x, false); } +// The algorthim for erfcinv is an adaptation of the erfinv algorithm above +// from P. J. Acklam. It has been modified to run over the different input +// domain of erfcinv. See the notes for erfinv for an explanation. + static double do_erfcinv (double x, bool refine) { // Coefficients of rational approximation. @@ -3241,12 +3245,12 @@ 2.445134137142996e+00, 3.754408661907416e+00 }; static const double spi2 = 8.862269254527579e-01; // sqrt(pi)/2. - static const double pi = 3.14159265358979323846; - static const double pbreak = 0.95150; + static const double pbreak_lo = 0.04850; // 1-pbreak + static const double pbreak_hi = 1.95150; // 1+pbreak double y; // Select case. - if (x <= 1+pbreak && x >= 1-pbreak) + if (x >= pbreak_lo && x <= pbreak_hi) { // Middle region. const double q = 0.5*(1-x), r = q*q; @@ -3254,15 +3258,15 @@ const double yd = ((((b[0]*r + b[1])*r + b[2])*r + b[3])*r + b[4])*r + 1.0; y = yn / yd; } - else if (x < 2.0 && x > 0.0) + else if (x > 0.0 && x < 2.0) { // Tail region. const double q = x < 1 ? sqrt (-2*log (0.5*x)) : sqrt (-2*log (0.5*(2-x))); const double yn = ((((c[0]*q + c[1])*q + c[2])*q + c[3])*q + c[4])*q + c[5]; const double yd = (((d[0]*q + d[1])*q + d[2])*q + d[3])*q + 1.0; y = yn / yd; - if (x < 1-pbreak) - y *= -1; + if (x < pbreak_lo) + y = -y; } else if (x == 0.0) return octave_Inf;
--- a/scripts/io/strread.m Mon Jun 25 14:21:45 2012 -0500 +++ b/scripts/io/strread.m Mon Jun 25 14:40:21 2012 -0500 @@ -154,7 +154,21 @@ ## ## @end table ## -## @seealso{textscan, textread, load, dlmread, fscanf} +## When the number of words in @var{str} doesn't match an exact multiple +## of the number of format conversion specifiers, strread's behavior +## depends on the last character of @var{str}: +## +## @table @asis +## @item last character = "\n" +## Data columns are padded with empty fields or Nan so that all columns +## have equal length +## +## @item last character is not "\n" +## Data columns are not padded; strread returns columns of unequal length +## +## @end table +## +# @seealso{textscan, textread, load, dlmread, fscanf} ## @end deftypefn function varargout = strread (str, format = "%f", varargin) @@ -445,7 +459,8 @@ ## Alternative below goes by simply parsing a first grab of words ## and counting words until the fmt_words array is exhausted: iwrd = 1; iwrdp = 0; iwrdl = length (words{iwrd}); - for ii = 1:numel (fmt_words) + ii = 1; + while ii <= numel (fmt_words) nxt_wrd = 0; @@ -507,12 +522,16 @@ if (nxt_wrd) ++iwrd; iwrdp = 0; - if (ii < numel (fmt_words)) + if (iwrd > numel (words)) + ## Apparently EOF; assume incomplete row already at L.1 of data + ii = numel (fmt_words); + elseif (ii < numel (fmt_words)) iwrdl = length (words{iwrd}); endif endif + ++ii; - endfor + endwhile ## Done words_period = max (iwrd - 1, 1); num_lines = ceil (num_words / words_period); @@ -896,6 +915,16 @@ %! assert (c, [0.94; 0.87], 0.01) %!test +%! [a, b] = strread (['Empty 1' char(10)], 'Empty%s %f'); +%! assert (a{1}, '1'); +%! assert (b, NaN); + +%!test +%! [a, b] = strread (['Empty' char(10)], 'Empty%f %f'); +%! assert (a, NaN); +%! assert (b, NaN); + +%!test %! # Bug #35999 %! [a, b, c] = strread ("", "%f"); %! assert (isempty (a));
--- a/scripts/io/textscan.m Mon Jun 25 14:21:45 2012 -0500 +++ b/scripts/io/textscan.m Mon Jun 25 14:40:21 2012 -0500 @@ -25,13 +25,13 @@ ## @deftypefnx {Function File} {[@var{C}, @var{position}] =} textscan (@var{fid}, @dots{}) ## Read data from a text file or string. ## -## The file associated with @var{fid} is read and parsed according to -## @var{format}. The function behaves like @code{strread} except it works by -## parsing a file instead of a string. See the documentation of -## @code{strread} for details. +## The string @var{str} or file associated with @var{fid} is read from and +## parsed according to @var{format}. The function behaves like @code{strread} +## except it can also read from file instead of a string. See the documentation +## of @code{strread} for details. ## -## In addition to the options supported by -## @code{strread}, this function supports a few more: +## In addition to the options supported by @code{strread}, this function +## supports a few more: ## ## @itemize ## @item "collectoutput": @@ -50,16 +50,19 @@ ## @item "returnonerror": ## If set to numerical 1 or true (default), return normally when read errors ## have been encountered. If set to 0 or false, return an error and no data. +## As the string or file is read by columns rather than by rows, and because +## textscan is fairly forgiving as regards read errors, setting this option +## may have little or no actual effect. ## @end itemize ## ## When reading from a character string, optional input argument @var{n} ## specifies the number of times @var{format} should be used (i.e., to limit ## the amount of data read). -## When reading fro file, @var{n} specifies the number of data lines to read; +## When reading from file, @var{n} specifies the number of data lines to read; ## in this sense it differs slightly from the format repeat count in strread. ## -## The output @var{C} is a cell array whose length is given by the number -## of format specifiers. +## The output @var{C} is a cell array whose second dimension is determined +## by the number of format specifiers. ## ## The second output, @var{position}, provides the position, in characters, ## from the beginning of the file. @@ -80,14 +83,18 @@ format = "%f"; endif + if (! ischar (format)) + error ("textscan: FORMAT must be a string"); + endif + + ## Determine the number of data fields & initialize output array + num_fields = numel (strfind (format, "%")) - numel (strfind (format, "%*")); + C = cell (1, num_fields); + if (! (isa (fid, "double") && fid > 0) && ! ischar (fid)) error ("textscan: first argument must be a file id or character string"); endif - if (! ischar (format)) - error ("textscan: FORMAT must be a string"); - endif - args = varargin; if (nargin > 2 && isnumeric (args{1})) nlines = args{1}; @@ -96,7 +103,6 @@ endif if (nlines < 1) printf ("textscan: N = 0, no data read\n"); - C = []; return endif @@ -174,7 +180,6 @@ ## Check for empty result if (isempty (str)) warning ("textscan: no data read"); - C = []; return; endif @@ -249,10 +254,18 @@ ## Determine the number of data fields num_fields = numel (strfind (format, "%")) - numel (strfind (format, "%*")); - ## Strip trailing EOL to avoid returning stray missing values (f. strread) - if (strcmp (str(end-length (eol_char) + 1 : end), eol_char)); - str(end-length (eol_char) + 1 : end) = ""; - endif + ## Strip trailing EOL to avoid returning stray missing values (f. strread). + ## However, in case of CollectOutput request, presence of EOL is required + eol_at_end = strcmp (str(end-length (eol_char) + 1 : end), eol_char); + if (collop) + if (! eol_at_end) + str(end+1 : end+length (eol_char)) = eol_char; + endif + elseif (eol_at_end) + str(end-length (eol_char) + 1 : end) = ""; + ## A corner case: str may now be empty.... + if (isempty (str)); return; endif + endif ## Call strread to make it do the real work C = cell (1, num_fields); @@ -316,14 +329,14 @@ %! assert (b(1,:)', c{1}, 1e-5); %! assert (b(2,:)', c{2}, 1e-5); -#%!test -#%! str = "13, 72, NA, str1, 25\r\n// Middle line\r\n36, na, 05, str3, 6"; -#%! a = textscan (str, "%d %n %f %s %n", "delimiter", ",","treatAsEmpty", {"NA", "na"},"commentStyle", "//"); -#%! assert (a{1}, int32([13; 36])); -#%! assert (a{2}, [72; NaN]); -#%! assert (a{3}, [NaN; 5]); -#%! assert (a{4}, {"str1"; "str3"}); -#%! assert (a{5}, [25; 6]); +%!test +%! str = "13, 72, NA, str1, 25\r\n// Middle line\r\n36, na, 05, str3, 6"; +%! a = textscan (str, "%d %n %f %s %n", "delimiter", ",","treatAsEmpty", {"NA", "na"},"commentStyle", "//"); +%! assert (a{1}, int32([13; 36])); +%! assert (a{2}, [72; NaN]); +%! assert (a{3}, [NaN; 5]); +%! assert (a{4}, {"str1"; "str3"}); +%! assert (a{5}, [25; 6]); %!test %! str = "Km:10 = hhhBjjj miles16hour\r\n"; @@ -362,6 +375,21 @@ %! assert (size(c{3}), [10, 2]); %! assert (size(c{2}), [10, 2]); +%!test +%% CollectOutput test with uneven column length files +%! b = [10:10:100]; +%! b = [b; 8*b/5; 8*b*1000/5]; +%! str = sprintf ("%g miles/hr = %g (%g) kilometers (meters)/hr\n", b); +%! str = [str "110 miles/hr"]; +%! fmt = "%f miles%s %s %f (%f) kilometers %*s"; +%! c = textscan (str, fmt, "collectoutput", 1); +%! assert (size(c{1}), [11, 1]); +%! assert (size(c{3}), [11, 2]); +%! assert (size(c{2}), [11, 2]); +%! assert (c{3}(end), NaN); +%! assert (c{2}{11, 1}, "/hr"); +%! assert (isempty (c{2}{11, 2}), true); + %% Test input validation %!error textscan () %!error textscan (single (4)) @@ -370,3 +398,7 @@ %!error <cannot provide position information> [C, pos] = textscan ("Hello World") %!error <character value required> textscan ("Hello World", '%s', 'EndOfLine', 3) +%! Test incomplete first data line +%! R = textscan (['Empty1' char(10)], 'Empty%d %f'); +%! assert (R{1}, int32(1)); +%! assert (isempty(R{2}), true);
--- a/scripts/optimization/lsqnonneg.m Mon Jun 25 14:21:45 2012 -0500 +++ b/scripts/optimization/lsqnonneg.m Mon Jun 25 14:40:21 2012 -0500 @@ -21,6 +21,7 @@ ## -*- texinfo -*- ## @deftypefn {Function File} {@var{x} =} lsqnonneg (@var{c}, @var{d}) ## @deftypefnx {Function File} {@var{x} =} lsqnonneg (@var{c}, @var{d}, @var{x0}) +## @deftypefnx {Function File} {@var{x} =} lsqnonneg (@var{c}, @var{d}, @var{x0}, @var{options}) ## @deftypefnx {Function File} {[@var{x}, @var{resnorm}] =} lsqnonneg (@dots{}) ## @deftypefnx {Function File} {[@var{x}, @var{resnorm}, @var{residual}] =} lsqnonneg (@dots{}) ## @deftypefnx {Function File} {[@var{x}, @var{resnorm}, @var{residual}, @var{exitflag}] =} lsqnonneg (@dots{}) @@ -29,6 +30,9 @@ ## Minimize @code{norm (@var{c}*@var{x} - d)} subject to ## @code{@var{x} >= 0}. @var{c} and @var{d} must be real. @var{x0} is an ## optional initial guess for @var{x}. +## Currently, @code{lsqnonneg} +## recognizes these options: @code{"MaxIter"}, @code{"TolX"}. +## For a description of these options, see @ref{doc-optimset,,optimset}. ## ## Outputs: ## @@ -147,7 +151,8 @@ ## compute the gradient. w = c'*(d - c*x); w(p) = []; - if (! any (w > 0)) + tolx = optimget (options, "TolX", 10*eps*norm (c, 1)*length (c)); + if (! any (w > tolx)) if (useqr) ## verify the solution achieved using qr updating. ## in the best case, this should only take a single step.
--- a/scripts/plot/plotyy.m Mon Jun 25 14:21:45 2012 -0500 +++ b/scripts/plot/plotyy.m Mon Jun 25 14:40:21 2012 -0500 @@ -91,6 +91,7 @@ ax = ax(1:2); elseif (length (ax) == 1) ax(2) = axes (); + set (ax(2), "nextplot", get (ax(1), "nextplot")) elseif (isempty (ax)) ax(1) = axes (); ax(2) = axes (); @@ -156,6 +157,7 @@ axes (ax(2)); else ax(2) = axes (); + set (ax(2), "nextplot", get (ax(1), "nextplot")) endif newplot (); @@ -172,7 +174,6 @@ set (ax(2), "yaxislocation", "right"); set (ax(2), "ycolor", getcolor (h2(1))); - if (strcmp (get(ax(1), "activepositionproperty"), "position")) set (ax(2), "position", get (ax(1), "position")); else @@ -208,6 +209,8 @@ addlistener (ax(2), "plotboxaspectratio", {@update_position, ax(1)}); addlistener (ax(1), "plotboxaspectratiomode", {@update_position, ax(2)}); addlistener (ax(2), "plotboxaspectratiomode", {@update_position, ax(1)}); + addlistener (ax(1), "nextplot", {@update_nextplot, ax(2)}); + addlistener (ax(2), "nextplot", {@update_nextplot, ax(1)}); ## Store the axes handles for the sister axes. if (ishandle (ax(1)) && ! isprop (ax(1), "__plotyy_axes__")) @@ -261,10 +264,24 @@ %! clf; %! x = linspace (-1, 1, 201); %! hax = plotyy (x, sin (pi*x), x, cos (pi*x)); -%! ylabel ('Blue on the Left'); +%! ylabel (hax(1), 'Blue on the Left'); %! ylabel (hax(2), 'Green on the Right'); %! xlabel ('xlabel'); +%!demo +%! clf +%! hold on +%! t = (0:0.1:9); +%! x = sin (t); +%! y = 5 * cos (t); +%! [hax, h1, h2] = plotyy (t, x, t, y); +%! [~, h3, h4] = plotyy (t+1, x, t+1, y); +%! set ([h3, h4], "linestyle", "--") +%! xlabel (hax(1), 'xlabel') +%! title (hax(2), 'title') +%! ylabel (hax(1), 'Left axis is Blue') +%! ylabel (hax(2), 'Right axis is Green') + function deleteplotyy (h, d, ax2, t2) if (ishandle (ax2) && strcmp (get (ax2, "type"), "axes") && (isempty (gcbf()) || strcmp (get (gcbf(), "beingdeleted"),"off")) @@ -274,6 +291,19 @@ endif endfunction +function update_nextplot (h, d, ax2) + persistent recursion = false; + prop = "nextplot"; + if (! recursion) + unwind_protect + recursion = true; + set (ax2, prop, get (h, prop)); + unwind_protect_cleanup + recursion = false; + end_unwind_protect + endif +endfunction + function update_position (h, d, ax2) persistent recursion = false;
--- a/scripts/plot/private/__scatter__.m Mon Jun 25 14:21:45 2012 -0500 +++ b/scripts/plot/private/__scatter__.m Mon Jun 25 14:40:21 2012 -0500 @@ -268,7 +268,8 @@ ## Does gnuplot only support triangles with different vertex colors ? ## TODO - Verify gnuplot can only support one color. If RGB triplets ## can be assigned to each vertex, then fix __go_draw_axe__.m - gnuplot_hack = numel (x) > 1 && strcmp (toolkit, "gnuplot"); + gnuplot_hack = (numel (x) > 1 && size(c, 2) == 3 + && strcmp (toolkit, "gnuplot")); if (ischar (c) || ! isflat || gnuplot_hack) if (filled) h = __go_patch__ (hg, "xdata", x, "ydata", y, "zdata", z,
--- a/scripts/plot/scatter.m Mon Jun 25 14:21:45 2012 -0500 +++ b/scripts/plot/scatter.m Mon Jun 25 14:40:21 2012 -0500 @@ -92,6 +92,14 @@ %! clf; %! x = randn (100, 1); %! y = randn (100, 1); +%! c = x .* y; +%! scatter (x, y, 20, c, "filled"); +%! title ('Scatter with colors'); + +%!demo +%! clf; +%! x = randn (100, 1); +%! y = randn (100, 1); %! scatter (x, y, [], sqrt (x.^2 + y.^2)); %! title ('Scatter plot with bubble color determined by distance from origin');
--- a/src/graphics.cc Mon Jun 25 14:21:45 2012 -0500 +++ b/src/graphics.cc Mon Jun 25 14:40:21 2012 -0500 @@ -4284,7 +4284,7 @@ xticklabelmode = "auto"; yticklabelmode = "auto"; zticklabelmode = "auto"; - color = "none"; + color = color_values ("white"); xcolor = color_values ("black"); ycolor = color_values ("black"); zcolor = color_values ("black");