Mercurial > octave
changeset 30678:79d6280fb00a stable
interp2: Fix issues with complex input (bug #61903).
* scripts/polynomial/pchip.m: Fix interpolation of imaginary part of complex
numbers.
* scripts/general/interp1.m: Fix value for extrapolation of complex values.
Add BISTs.
| author | Christof Kaufmann <christofkaufmann.dev@gmail.com> |
|---|---|
| date | Mon, 24 Jan 2022 12:52:29 +0100 |
| parents | c3494f408895 |
| children | f9a922cdb694 7e5e77ef09d7 |
| files | scripts/general/interp1.m scripts/polynomial/pchip.m |
| diffstat | 2 files changed, 56 insertions(+), 1 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/general/interp1.m Thu Jan 27 18:06:01 2022 +0100 +++ b/scripts/general/interp1.m Mon Jan 24 12:52:29 2022 +0100 @@ -123,7 +123,7 @@ endif method = "linear"; - extrap = NA; + extrap = []; xi = []; ispp = false; have_xi = false; @@ -167,6 +167,14 @@ endif endif + if (isempty (extrap)) + if (iscomplex (y)) + extrap = NA + 1i*NA; + else + extrap = NA; + endif + endif + ## reshape matrices for convenience x = x(:); nx = rows (x); @@ -503,6 +511,9 @@ %!assert (isempty (interp1 (xp,yp,[],style))) %!assert (interp1 (xp,[yp',yp'],xi(:),style),... %! [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)]) +%!assert <*61903> ... +%! (interp1 (xp, yp + 1i * yp.^2, xi, style),... +%! interp1 (xp,yp,xi,style) + 1i * interp1 (xp,yp.^2,xi,style)) %!assert (interp1 (xp,yp,xi,style),... %! interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps) %!assert (ppval (interp1 (xp,yp,style,"pp"),xi), @@ -520,6 +531,9 @@ %!assert (isempty (interp1 (xp,yp,[],style))) %!assert (interp1 (xp,[yp',yp'],xi(:),style),... %! [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)]) +%!assert <*61903> ... +%! (interp1 (xp, yp + 1i * yp.^2, xi, style),... +%! interp1 (xp,yp,xi,style) + 1i * interp1 (xp,yp.^2,xi,style)) %!assert (interp1 (xp,yp,xi,style),... %! interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps) %!assert (ppval (interp1 (xp,yp,style,"pp"),xi), @@ -540,6 +554,9 @@ %!assert (isempty (interp1 (xp,yp,[],style))) %!assert (interp1 (xp,[yp',yp'],xi(:),style),... %! [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)]) +%!assert <*61903> ... +%! (interp1 (xp, yp + 1i * yp.^2, xi, style),... +%! interp1 (xp,yp,xi,style) + 1i * interp1 (xp,yp.^2,xi,style)) ## This test is expected to fail, so commented out. ## "previous" and "next" options are not symmetric w.r.t to flipping xp,yp #%!assert (interp1 (xp,yp,xi,style),... @@ -559,6 +576,9 @@ %!assert (isempty (interp1 (xp,yp,[],style))) %!assert (interp1 (xp,[yp',yp'],xi(:),style),... %! [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)]) +%!assert <*61903> ... +%! (interp1 (xp, yp + 1i * yp.^2, xi, style),... +%! interp1 (xp,yp,xi,style) + 1i * interp1 (xp,yp.^2,xi,style)) #%!assert (interp1 (xp,yp,xi,style),... #%! interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps) %!assert (ppval (interp1 (xp,yp,style,"pp"),xi), @@ -579,6 +599,9 @@ %!assert (isempty (interp1 (xp,yp,[],style))) %!assert (interp1 (xp,[yp',yp'],xi(:),style),... %! [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)]) +%!assert <*61903> ... +%! (interp1 (xp, yp + 1i * yp.^2, xi, style),... +%! interp1 (xp,yp,xi,style) + 1i * interp1 (xp,yp.^2,xi,style)) #%!assert (interp1 (xp,yp,xi,style),... #%! interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps) %!assert (ppval (interp1 (xp,yp,style,"pp"),xi), @@ -596,6 +619,9 @@ %!assert (isempty (interp1 (xp,yp,[],style))) %!assert (interp1 (xp,[yp',yp'],xi(:),style),... %! [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)]) +%!assert <*61903> ... +%! (interp1 (xp, yp + 1i * yp.^2, xi, style),... +%! interp1 (xp,yp,xi,style) + 1i * interp1 (xp,yp.^2,xi,style)) #%!assert (interp1 (xp,yp,xi,style),... #%! interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps) %!assert (ppval (interp1 (xp,yp,style,"pp"),xi), @@ -616,6 +642,9 @@ %!assert (isempty (interp1 (xp,yp,[],style))) %!assert (interp1 (xp,[yp',yp'],xi(:),style),... %! [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)]) +%!assert <*61903> ... +%! (interp1 (xp, yp + 1i * yp.^2, xi, style),... +%! interp1 (xp,yp,xi,style) + 1i * interp1 (xp,yp.^2,xi,style)) %!assert (interp1 (xp,yp,xi,style),... %! interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps) %!assert (ppval (interp1 (xp,yp,style,"pp"),xi), @@ -633,6 +662,9 @@ %!assert (isempty (interp1 (xp,yp,[],style))) %!assert (interp1 (xp,[yp',yp'],xi(:),style),... %! [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)]) +%!assert <*61903> ... +%! (interp1 (xp, yp + 1i * yp.^2, xi, style),... +%! interp1 (xp,yp,xi,style) + 1i * interp1 (xp,yp.^2,xi,style)) %!assert (interp1 (xp,yp,xi,style),... %! interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps) %!assert (ppval (interp1 (xp,yp,style,"pp"),xi), @@ -655,6 +687,9 @@ %!assert (isempty (interp1 (xp,yp,[],style))) %!assert (interp1 (xp,[yp',yp'],xi(:),style),... %! [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)]) +%!assert <*61903> ... +%! (interp1 (xp, yp + 1i * yp.^2, xi, style),... +%! interp1 (xp,yp,xi,style) + 1i * interp1 (xp,yp.^2,xi,style)) %!assert (interp1 (xp,yp,xi,style),... %! interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps) %!assert (ppval (interp1 (xp,yp,style,"pp"),xi), @@ -672,6 +707,9 @@ %!assert (isempty (interp1 (xp,yp,[],style))) %!assert (interp1 (xp,[yp',yp'],xi(:),style),... %! [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)]) +%!assert <*61903> ... +%! (interp1 (xp, yp + 1i * yp.^2, xi, style),... +%! interp1 (xp,yp,xi,style) + 1i * interp1 (xp,yp.^2,xi,style)) %!assert (interp1 (xp,yp,xi,style),... %! interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps) %!assert (ppval (interp1 (xp,yp,style,"pp"),xi), @@ -692,6 +730,9 @@ %!assert (isempty (interp1 (xp,yp,[],style))) %!assert (interp1 (xp,[yp',yp'],xi(:),style),... %! [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)]) +%!assert <*61903> ... +%! (interp1 (xp, yp + 1i * yp.^2, xi, style),... +%! interp1 (xp,yp,xi,style) + 1i * interp1 (xp,yp.^2,xi,style)) %!assert (interp1 (xp,yp,xi,style),... %! interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps) %!assert (ppval (interp1 (xp,yp,style,"pp"),xi), @@ -709,6 +750,9 @@ %!assert (isempty (interp1 (xp,yp,[],style))) %!assert (interp1 (xp,[yp',yp'],xi(:),style),... %! [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)]) +%!assert <*61903> ... +%! (interp1 (xp, yp + 1i * yp.^2, xi, style),... +%! interp1 (xp,yp,xi,style) + 1i * interp1 (xp,yp.^2,xi,style)) %!assert (interp1 (xp,yp,xi,style),... %! interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps) %!assert (ppval (interp1 (xp,yp,style,"pp"),xi), @@ -729,6 +773,9 @@ %!assert (isempty (interp1 (xp,yp,[],style))) %!assert (interp1 (xp,[yp',yp'],xi(:),style),... %! [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)]) +%!assert <*61903> ... +%! (interp1 (xp, yp + 1i * yp.^2, xi, style),... +%! interp1 (xp,yp,xi,style) + 1i * interp1 (xp,yp.^2,xi,style)) %!assert (interp1 (xp,yp,xi,style),... %! interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps) %!assert (ppval (interp1 (xp,yp,style,"pp"),xi), @@ -746,6 +793,9 @@ %!assert (isempty (interp1 (xp,yp,[],style))) %!assert (interp1 (xp,[yp',yp'],xi(:),style),... %! [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)]) +%!assert <*61903> ... +%! (interp1 (xp, yp + 1i * yp.^2, xi, style),... +%! interp1 (xp,yp,xi,style) + 1i * interp1 (xp,yp.^2,xi,style)) %!assert (interp1 (xp,yp,xi,style),... %! interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps) %!assert (ppval (interp1 (xp,yp,style,"pp"),xi),
--- a/scripts/polynomial/pchip.m Thu Jan 27 18:06:01 2022 +0100 +++ b/scripts/polynomial/pchip.m Mon Jan 24 12:52:29 2022 +0100 @@ -111,6 +111,11 @@ ## Compute derivatives. d = __pchip_deriv__ (x, y, 2); + if (iscomplex (y)) + ## __pchip_deriv__ ignores imaginary part. Do it again for imag part. + ## FIXME: Adapt __pchip_deriv__ to correctly handle complex input. + d += 1i * __pchip_deriv__ (x, imag (y), 2); + endif d1 = d(:, 1:n-1); d2 = d(:, 2:n);