Mercurial > octave
changeset 27435:a8a5d2e8807f
Produce symmetric sequences for linspace (bug #56659).
* NEWS: Announce change.
* data.cc (Flinspace): Add new BIST tests.
* CRowVector.cc, dRowVector.cc, fCRowVector.cc, fRowVector.cc (linspace):
Change algorithm to construct sequence starting from both endpoints and
meeting in the middle to create symmetric sequences.
* logspace.m: Update BIST tests.
| author | Rik <rik@octave.org> |
|---|---|
| date | Fri, 20 Sep 2019 09:03:14 -0700 |
| parents | 57f5c5768eb3 |
| children | ece17410605e |
| files | NEWS libinterp/corefcn/data.cc liboctave/array/CRowVector.cc liboctave/array/dRowVector.cc liboctave/array/fCRowVector.cc liboctave/array/fRowVector.cc scripts/general/logspace.m |
| diffstat | 7 files changed, 160 insertions(+), 47 deletions(-) [+] |
line wrap: on
line diff
--- a/NEWS Sat Sep 21 17:17:02 2019 -0700 +++ b/NEWS Fri Sep 20 09:03:14 2019 -0700 @@ -13,6 +13,10 @@ accept a new sorting option `"stable"` which will return output values in the same order as the input, rather than in ascending order. +- The `linspace` function now produces symmetrical sequences when the + endpoints are symmetric. This is more intuitive and also compatible + with recent changes made in Matlab R2019b. + #### Graphics backend - Graphic primitives now accept a color property value of `"none"`
--- a/libinterp/corefcn/data.cc Sat Sep 21 17:17:02 2019 -0700 +++ b/libinterp/corefcn/data.cc Fri Sep 20 09:03:14 2019 -0700 @@ -5287,6 +5287,31 @@ %!assert (class (linspace (single (1), 2)), "single") %!assert (class (linspace (1, single (2))), "single") +## Test symmetry +%!test <*56659> +%! x = linspace (-1, 1, 10); +%! assert (all (x == -fliplr (x))); +%! x = linspace (-1, 1, 11); +%! assert (all (x == -fliplr (x))); + +%!test <*56659> +%! x = linspace (-1-1i, 1+1i, 10); +%! assert (all (x == -fliplr (x))); +%! x = linspace (-1-1i, 1+1i, 11); +%! assert (all (x == -fliplr (x))); + +%!test <*56659> +%! x = linspace (single (-1), 1, 10); +%! assert (all (x == -fliplr (x))); +%! x = linspace (single (-1), 1, 11); +%! assert (all (x == -fliplr (x))); + +%!test <*56659> +%! x = linspace (single (-1-1i), 1+1i, 10); +%! assert (all (x == -fliplr (x))); +%! x = linspace (single (-1-1i), 1+1i, 11); +%! assert (all (x == -fliplr (x))); + ## Test obscure Matlab compatibility options %!assert (linspace (0, 1, []), 1) %!assert (linspace (10, 20, 2), [10 20]) @@ -5300,14 +5325,21 @@ %!assert (1 ./ linspace (-0, 0, 4), [-Inf, Inf, Inf, Inf]) %!assert (linspace (Inf, Inf, 3), [Inf, Inf, Inf]) %!assert (linspace (-Inf, -Inf, 3), [-Inf, -Inf, -Inf]) -%!assert (linspace (-Inf, Inf, 3), [-Inf, NaN, Inf]) +%!assert (linspace (-Inf, Inf, 3), [-Inf, 0, Inf]) %!assert (linspace (Inf + 1i, Inf + 1i, 3), [Inf + 1i, Inf + 1i, Inf + 1i]) %!assert (linspace (-Inf + 1i, Inf + 1i, 3), [-Inf + 1i, NaN + 1i, Inf + 1i]) -%!assert (linspace (0, Inf, 3), [0, Inf, Inf]) -%!assert (linspace (0, -Inf, 3), [0, -Inf, -Inf]) -%!assert (linspace (-Inf, 0, 3), [-Inf, NaN, 0]) -%!assert (linspace (Inf, 0, 3), [Inf, NaN, 0]) -%!assert (linspace (Inf, -Inf, 3), [Inf, NaN, -Inf]) + +## FIXME: Octave is not fully Matlab-compatible for some combinations of +## Inf/-Inf endpoints. See bug #56933. This was dubbed "Won't Fix" +## so these tests have been removed from the test suite by commenting +## them out. If the behavior in the future is made compatible these +## tests can be re-instated. +##%!assert <56933> (linspace (-Inf, Inf, 4), [-Inf, -Inf, Inf, Inf]) +##%!assert <56933> (linspace (-Inf, Inf, 5), [-Inf, -Inf, 0, Inf, Inf]) +##%!assert <56933> (linspace (0, Inf, 4), [0, Inf, Inf, Inf]) +##%!assert <56933> (linspace (0, -Inf, 4), [0, -Inf, -Inf, -Inf]) +##%!assert <56933> (linspace (-Inf, 0, 4), [-Inf, NaN, NaN, 0]) +##%!assert <56933> (linspace (Inf, 0, 4), [Inf, NaN, NaN, 0]) %!error linspace () %!error linspace (1, 2, 3, 4)
--- a/liboctave/array/CRowVector.cc Sat Sep 21 17:17:02 2019 -0700 +++ b/liboctave/array/CRowVector.cc Fri Sep 20 09:03:14 2019 -0700 @@ -26,6 +26,7 @@ #include <istream> #include <ostream> +#include <type_traits> #include "Array-util.h" #include "lo-blas-proto.h" @@ -421,22 +422,39 @@ // other operations ComplexRowVector -linspace (const Complex& x1, const Complex& x2, octave_idx_type n) +linspace (const Complex& x1, const Complex& x2, octave_idx_type n_in) { NoAlias<ComplexRowVector> retval; - if (n < 1) + if (n_in < 1) return retval; - else - retval.clear (n); + else if (n_in == 1) + { + retval.resize (1, x2); + return retval; + } + // Use unsigned type (guaranteed n_in > 1 at this point) so that divisions + // by 2 can be replaced by compiler with shift right instructions. + typedef std::make_unsigned<octave_idx_type>::type unsigned_octave_idx_type; + + unsigned_octave_idx_type n = n_in; + + // Set endpoints, rather than calculate, for maximum accuracy. + retval.clear (n); retval(0) = x1; + retval(n-1) = x2; - Complex delta = (x1 == x2) ? 0 : (x2 - x1) / (n - 1.0); - for (octave_idx_type i = 1; i < n-1; i++) - retval(i) = x1 + static_cast<double> (i)*delta; - - retval(n-1) = x2; + // Construct linspace symmetrically from both ends. + Complex delta = (x2 - x1) / (n - 1.0); + unsigned_octave_idx_type n2 = n/2; + for (unsigned_octave_idx_type i = 1; i < n2; i++) + { + retval(i) = x1 + static_cast<double> (i)*delta; + retval(n-1-i) = x2 - static_cast<double> (i)*delta; + } + if (n % 2 == 1) // Middle element if number of elements is odd. + retval(n2) = (x1 == -x2 ? 0 : (x1 + x2) / 2.0); return retval; }
--- a/liboctave/array/dRowVector.cc Sat Sep 21 17:17:02 2019 -0700 +++ b/liboctave/array/dRowVector.cc Fri Sep 20 09:03:14 2019 -0700 @@ -26,6 +26,7 @@ #include <istream> #include <ostream> +#include <type_traits> #include "Array-util.h" #include "lo-blas-proto.h" @@ -264,22 +265,39 @@ // other operations RowVector -linspace (double x1, double x2, octave_idx_type n) +linspace (double x1, double x2, octave_idx_type n_in) { NoAlias<RowVector> retval; - if (n < 1) + if (n_in < 1) return retval; - else - retval.clear (n); + else if (n_in == 1) + { + retval.resize (1, x2); + return retval; + } + // Use unsigned type (guaranteed n_in > 1 at this point) so that divisions + // by 2 can be replaced by compiler with shift right instructions. + typedef std::make_unsigned<octave_idx_type>::type unsigned_octave_idx_type; + + unsigned_octave_idx_type n = n_in; + + // Set endpoints, rather than calculate, for maximum accuracy. + retval.clear (n); retval(0) = x1; + retval(n-1) = x2; - double delta = (x1 == x2) ? 0 : ((x2 - x1) / (n - 1)); - for (octave_idx_type i = 1; i < n-1; i++) - retval(i) = x1 + i*delta; - - retval(n-1) = x2; + // Construct linspace symmetrically from both ends. + double delta = (x2 - x1) / (n - 1); + unsigned_octave_idx_type n2 = n/2; + for (unsigned_octave_idx_type i = 1; i < n2; i++) + { + retval(i) = x1 + i*delta; + retval(n-1-i) = x2 - i*delta; + } + if (n % 2 == 1) // Middle element if number of elements is odd. + retval(n2) = (x1 == -x2 ? 0 : (x1 + x2) / 2); return retval; }
--- a/liboctave/array/fCRowVector.cc Sat Sep 21 17:17:02 2019 -0700 +++ b/liboctave/array/fCRowVector.cc Fri Sep 20 09:03:14 2019 -0700 @@ -26,6 +26,7 @@ #include <istream> #include <ostream> +#include <type_traits> #include "Array-util.h" #include "lo-blas-proto.h" @@ -423,22 +424,39 @@ // other operations FloatComplexRowVector -linspace (const FloatComplex& x1, const FloatComplex& x2, octave_idx_type n) +linspace (const FloatComplex& x1, const FloatComplex& x2, octave_idx_type n_in) { NoAlias<FloatComplexRowVector> retval; - if (n < 1) + if (n_in < 1) return retval; - else - retval.clear (n); + else if (n_in == 1) + { + retval.resize (1, x2); + return retval; + } + // Use unsigned type (guaranteed n_in > 1 at this point) so that divisions + // by 2 can be replaced by compiler with shift right instructions. + typedef std::make_unsigned<octave_idx_type>::type unsigned_octave_idx_type; + + unsigned_octave_idx_type n = n_in; + + // Set endpoints, rather than calculate, for maximum accuracy. + retval.clear (n); retval(0) = x1; + retval(n-1) = x2; - FloatComplex delta = (x1 == x2) ? 0 : (x2 - x1) / (n - 1.0f); - for (octave_idx_type i = 1; i < n-1; i++) - retval(i) = x1 + static_cast<float> (i)*delta; - - retval(n-1) = x2; + // Construct linspace symmetrically from both ends. + FloatComplex delta = (x2 - x1) / (n - 1.0f); + unsigned_octave_idx_type n2 = n/2; + for (unsigned_octave_idx_type i = 1; i < n2; i++) + { + retval(i) = x1 + static_cast<float> (i)*delta; + retval(n-1-i) = x2 - static_cast<float> (i)*delta; + } + if (n % 2 == 1) // Middle element if number of elements is odd. + retval(n2) = (x1 == -x2 ? 0 : (x1 + x2) / 2.0f); return retval; }
--- a/liboctave/array/fRowVector.cc Sat Sep 21 17:17:02 2019 -0700 +++ b/liboctave/array/fRowVector.cc Fri Sep 20 09:03:14 2019 -0700 @@ -26,6 +26,7 @@ #include <istream> #include <ostream> +#include <type_traits> #include "Array-util.h" #include "lo-blas-proto.h" @@ -264,22 +265,39 @@ // other operations FloatRowVector -linspace (float x1, float x2, octave_idx_type n) +linspace (float x1, float x2, octave_idx_type n_in) { NoAlias<FloatRowVector> retval; - if (n < 1) + if (n_in < 1) return retval; - else - retval.clear (n); + else if (n_in == 1) + { + retval.resize (1, x2); + return retval; + } + // Use unsigned type (guaranteed n_in > 1 at this point) so that divisions + // by 2 can be replaced by compiler with shift right instructions. + typedef std::make_unsigned<octave_idx_type>::type unsigned_octave_idx_type; + + unsigned_octave_idx_type n = n_in; + + // Set endpoints, rather than calculate, for maximum accuracy. + retval.clear (n); retval(0) = x1; + retval(n-1) = x2; - float delta = (x1 == x2) ? 0 : (x2 - x1) / (n - 1); - for (octave_idx_type i = 1; i < n-1; i++) - retval(i) = x1 + i*delta; - - retval(n-1) = x2; + // Construct linspace symmetrically from both ends. + float delta = (x2 - x1) / (n - 1); + unsigned_octave_idx_type n2 = n/2; + for (unsigned_octave_idx_type i = 1; i < n2; i++) + { + retval(i) = x1 + i*delta; + retval(n-1-i) = x2 - i*delta; + } + if (n % 2 == 1) // Middle element if number of elements is odd. + retval(n2) = (x1 == -x2 ? 0 : (x1 + x2) / 2); return retval; }
--- a/scripts/general/logspace.m Sat Sep 21 17:17:02 2019 -0700 +++ b/scripts/general/logspace.m Fri Sep 20 09:03:14 2019 -0700 @@ -97,16 +97,21 @@ %! assert (size (x2) == [1, 10] && abs (x2(1) - 10) < eps && abs (x2(10) - 100) < eps); %! assert (size (x3) == [1, 10] && abs (x3(1) - 10) < eps && abs (x3(10) - 0.01) < eps); %! assert (size (x4) == [1, 10] && abs (x4(1) - 10) < eps && abs (x4(10) - pi) < sqrt (eps)); + +## Edge cases %!assert (logspace (Inf, Inf, 3), [Inf, Inf, Inf]) -%!assert (logspace (-Inf, -Inf, 3), [0, 0, 0]) -%!assert (logspace (-Inf, Inf, 3), [0, NaN, Inf]) +%!assert (logspace (-Inf, Inf, 3), [0, 1, Inf]) %!assert (logspace (Inf + 1i, Inf + 1i, 3), repmat (complex (-Inf,Inf), [1, 3])) %!assert (logspace (-Inf + 1i, Inf + 1i, 3), [0, NaN + NaN * 1i, complex(-Inf, Inf)]) %!assert (logspace (0, Inf, 3), [1, Inf, Inf]) %!assert (logspace (0, -Inf, 3), [1, 0, 0]) -%!assert (logspace (-Inf, 0, 3), [0, NaN, 1]) -%!assert (logspace (Inf, 0, 3), [Inf, NaN, 1]) -%!assert (logspace (Inf, -Inf, 3), [Inf, NaN, 0]) +%!assert (logspace (Inf, -Inf, 3), [Inf, 1, 0]) + +## FIXME: These are bizarre corner cases for Matlab compatibility. See +## bug #56933. This is marked as "Won't Fix", but if linspace is updated at +## some point then these tests can be re-instated. +##%!assert (logspace (-Inf, 0, 3), [0, NaN, 1]) +##%!assert (logspace (Inf, 0, 3), [Inf, NaN, 1]) ## Test input validation %!error logspace ()