Mercurial > octave
changeset 31927:fef004fa8575
Improve Matlab compatibility of linspace for Inf endpoints.
Return an N-element vector if both START and END are the same.
When N is odd and [START, END] are [-Inf, Inf] return 0 for middle return value.
For all other elements besides START and END return NaN when one endpoint is
not finite.
* data.cc (Flinspace): Update and add new BIST tests for linspace.
* CRowVector.cc, dRowVector.cc, fCRowVector.cc, fRowVector.cc (linspace):
Add special case in input validation for START==END and return an N-element
vector of value END. Check whether variable delta is Inf value and change
delta to NaN if it is. For complex calculations, a separate isinf() check
is done for real and imaginary parts of delta. For middle point when N
is odd, return 0 if START == -END, NaN if delta was NaN, or otherwise
(START + END) / 2.
* logspace.m: Update BIST tests.
| author | Rik <rik@octave.org> |
|---|---|
| date | Fri, 24 Mar 2023 14:42:24 -0700 |
| parents | cd5be42a81ab |
| children | a8473e7db0f4 69c15474d829 |
| files | 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 | 6 files changed, 130 insertions(+), 20 deletions(-) [+] |
line wrap: on
line diff
--- a/libinterp/corefcn/data.cc Fri Mar 24 16:08:15 2023 -0400 +++ b/libinterp/corefcn/data.cc Fri Mar 24 14:42:24 2023 -0700 @@ -5745,27 +5745,41 @@ %!assert (numel (linspace (0, 1, 2-eps)), 1) %!assert (linspace (10, 20, 2.1), [10 20]) %!assert (linspace (10, 20, 2.9), [10 20]) -%!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, 0, Inf]) +## Octave prefers to return NaN which indicates failure of algorithm. +%!assert (linspace (-Inf, Inf, 4), [-Inf, NaN, NaN, Inf]) +%!assert (linspace (-Inf, 0, 3), [-Inf, NaN, 0]) +%!assert (linspace (-Inf, 0, 4), [-Inf, NaN, NaN, 0]) %!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 (-Inf - 1i, Inf + 1i, 3), [-Inf - 1i, 0 + 0i, Inf + 1i]) +%!assert (linspace (-Inf - 1i, Inf + 2i, 3), [-Inf - 1i, NaN + 0.5i, Inf + 2i]) +%!assert (linspace (-Inf - 3i, Inf + 0i, 4), +%! [-Inf - 3i, NaN - 2i, NaN - 1i, Inf + 0i]) +%!assert (linspace (complex (-1, -Inf), complex (1, Inf), 3), +%! [complex(-1, -Inf), 0 + 0i, complex(1, Inf)]) +%!assert (linspace (complex (-1, -Inf), complex (2, Inf), 3), +%! [complex(-1, -Inf), complex(0.5, NaN), complex(2, Inf)]) +%!assert (linspace (complex (-3, -Inf), complex (0, Inf), 4), +%! [complex(-3, -Inf), complex(-2, NaN), complex(-1, NaN), complex(0, 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. +## as Octave prefers to return NaN for some of these conditions to +## better reflect that the algorithm has failed. 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) +##%!assert (1 ./ linspace (-0, 0, 4), [-Inf, Inf, Inf, Inf]) + +## Test input validation +%!error <Invalid call> linspace () +%!error <Invalid call> linspace (1, 2, 3, 4) %!error <N must be a scalar> linspace (1, 2, [3, 4]) %!error <START, END must be scalars or vectors> linspace (ones (2,2), 2, 3) %!error <START, END must be scalars or vectors> linspace (2, ones (2,2), 3)
--- a/liboctave/array/CRowVector.cc Fri Mar 24 16:08:15 2023 -0400 +++ b/liboctave/array/CRowVector.cc Fri Mar 24 14:42:24 2023 -0700 @@ -436,6 +436,11 @@ retval.resize (1, x2); return retval; } + else if (x1 == x2) + { + retval.resize (n_in, 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. @@ -449,7 +454,17 @@ retval.xelem (n-1) = x2; // Construct linspace symmetrically from both ends. + bool isnan_delta = false; Complex delta = (x2 - x1) / (n - 1.0); + if (octave::math::isinf (delta)) + { + if (octave::math::isinf (delta.real ())) + delta.real (octave::numeric_limits<double>::NaN ()); + if (octave::math::isinf (delta.imag ())) + delta.imag (octave::numeric_limits<double>::NaN ()); + isnan_delta = true; + } + unsigned_octave_idx_type n2 = n/2; for (unsigned_octave_idx_type i = 1; i < n2; i++) { @@ -457,7 +472,22 @@ retval.xelem (n-1-i) = x2 - static_cast<double> (i)*delta; } if (n % 2 == 1) // Middle element if number of elements is odd. - retval.xelem (n2) = (x1 == -x2 ? 0 : (x1 + x2) / 2.0); + { + if (x1 == -x2) + retval.xelem (n2) = 0; + else + { + Complex c = (x1 + x2) / 2.0; + if (isnan_delta) + { + if (octave::math::isnan (delta.real ())) + c.real (octave::numeric_limits<double>::NaN ()); + if (octave::math::isnan (delta.imag ())) + c.imag (octave::numeric_limits<double>::NaN ()); + } + retval.xelem (n2) = c; + } + } return retval; }
--- a/liboctave/array/dRowVector.cc Fri Mar 24 16:08:15 2023 -0400 +++ b/liboctave/array/dRowVector.cc Fri Mar 24 14:42:24 2023 -0700 @@ -279,6 +279,11 @@ retval.resize (1, x2); return retval; } + else if (x1 == x2) + { + retval.resize (n_in, 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. @@ -292,7 +297,14 @@ retval.xelem (n-1) = x2; // Construct linspace symmetrically from both ends. + bool isnan_delta = false; double delta = (x2 - x1) / (n - 1); + if (octave::math::isinf (delta)) + { + delta = octave::numeric_limits<double>::NaN (); + isnan_delta = true; + } + unsigned_octave_idx_type n2 = n/2; for (unsigned_octave_idx_type i = 1; i < n2; i++) { @@ -300,7 +312,14 @@ retval.xelem (n-1-i) = x2 - i*delta; } if (n % 2 == 1) // Middle element if number of elements is odd. - retval.xelem (n2) = (x1 == -x2 ? 0 : (x1 + x2) / 2); + { + if (x1 == -x2) + retval.xelem (n2) = 0; + else if (isnan_delta) + retval.xelem (n2) = octave::numeric_limits<double>::NaN (); + else + retval.xelem (n2) = (x1 + x2) / 2; + } return retval; }
--- a/liboctave/array/fCRowVector.cc Fri Mar 24 16:08:15 2023 -0400 +++ b/liboctave/array/fCRowVector.cc Fri Mar 24 14:42:24 2023 -0700 @@ -438,6 +438,11 @@ retval.resize (1, x2); return retval; } + else if (x1 == x2) + { + retval.resize (n_in, 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. @@ -451,7 +456,17 @@ retval.xelem (n-1) = x2; // Construct linspace symmetrically from both ends. + bool isnan_delta = false; FloatComplex delta = (x2 - x1) / (n - 1.0f); + if (octave::math::isinf (delta)) + { + if (octave::math::isinf (delta.real ())) + delta.real (octave::numeric_limits<float>::NaN ()); + if (octave::math::isinf (delta.imag ())) + delta.imag (octave::numeric_limits<float>::NaN ()); + isnan_delta = true; + } + unsigned_octave_idx_type n2 = n/2; for (unsigned_octave_idx_type i = 1; i < n2; i++) { @@ -459,7 +474,22 @@ retval.xelem (n-1-i) = x2 - static_cast<float> (i)*delta; } if (n % 2 == 1) // Middle element if number of elements is odd. - retval.xelem (n2) = (x1 == -x2 ? 0 : (x1 + x2) / 2.0f); + { + if (x1 == -x2) + retval.xelem (n2) = 0; + else + { + FloatComplex c = (x1 + x2) / 2.0f; + if (isnan_delta) + { + if (octave::math::isnan (delta.real ())) + c.real (octave::numeric_limits<float>::NaN ()); + if (octave::math::isnan (delta.imag ())) + c.imag (octave::numeric_limits<float>::NaN ()); + } + retval.xelem (n2) = c; + } + } return retval; }
--- a/liboctave/array/fRowVector.cc Fri Mar 24 16:08:15 2023 -0400 +++ b/liboctave/array/fRowVector.cc Fri Mar 24 14:42:24 2023 -0700 @@ -279,6 +279,11 @@ retval.resize (1, x2); return retval; } + else if (x1 == x2) + { + retval.resize (n_in, 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. @@ -292,7 +297,14 @@ retval.xelem (n-1) = x2; // Construct linspace symmetrically from both ends. + bool isnan_delta = false; float delta = (x2 - x1) / (n - 1); + if (octave::math::isinf (delta)) + { + delta = octave::numeric_limits<float>::NaN (); + isnan_delta = true; + } + unsigned_octave_idx_type n2 = n/2; for (unsigned_octave_idx_type i = 1; i < n2; i++) { @@ -300,7 +312,14 @@ retval.xelem (n-1-i) = x2 - i*delta; } if (n % 2 == 1) // Middle element if number of elements is odd. - retval.xelem (n2) = (x1 == -x2 ? 0 : (x1 + x2) / 2); + { + if (x1 == -x2) + retval.xelem (n2) = 0; + else if (isnan_delta) + retval.xelem (n2) = octave::numeric_limits<float>::NaN (); + else + retval.xelem (n2) = (x1 + x2) / 2; + } return retval; }
--- a/scripts/general/logspace.m Fri Mar 24 16:08:15 2023 -0400 +++ b/scripts/general/logspace.m Fri Mar 24 14:42:24 2023 -0700 @@ -130,15 +130,13 @@ %!testif ; ismac () <55538> %! 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]) +## Octave prefers to return NaN which indicates failure of algorithm. +## Tests can be re-instated if full Matlab-compatibility is coded. +%!#assert (logspace (0, Inf, 3), [1, Inf, Inf]) +%!#assert (logspace (0, -Inf, 3), [1, 0, 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]) +%!assert (logspace (-Inf, 0, 3), [0, NaN, 1]) +%!assert (logspace (Inf, 0, 3), [Inf, NaN, 1]) ## Test input validation %!error <Invalid call> logspace ()