bugGNU Octave - Bugs: bug #61863, interp2 gives wrong results for...

 
 

bug #61863: interp2 gives wrong results for complex matrices

Submitter:  Christof Kaufmann <christofkaufmann>
Submitted:  Thu 20 Jan 2022 01:55:17 AM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Fixed Assigned to:  None
Originator Name:  Open/Closed:  * Closed
Release:  * dev Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Mon 31 Jan 2022 11:29:01 AM UTC, comment #32: 

Thanks!

Christof Kaufmann <christofkaufmann>
Mon 31 Jan 2022 03:09:58 AM UTC, comment #31: 

Verified changeset.  Marking as Fixed and closing report.

Rik <rik5>
Group administrator
Sun 30 Jan 2022 12:00:04 PM UTC, comment #30: 

Pushed the patch for the C++ code here:
https://hg.savannah.gnu.org/hgweb/octave/rev/4b367bf5eb16

Markus Mützel <mmuetzel>
Group administrator
Thu 27 Jan 2022 08:27:03 PM UTC, comment #29: 

I rebased the proper patch onto default, since the m-file patches have been applied to stable and stable got merged into default.

This patch does the following.

  • Fix _pchip_deriv_ for complex inputs and adds some tests to it.
  • Remove the workaround in interp2.
  • Remove the workaround in pchip.


(file #52760)

Christof Kaufmann <christofkaufmann>
Thu 27 Jan 2022 07:26:13 PM UTC, comment #28: 

Expected failures and Best practices section added at teh end of https://wiki.octave.org/Tests#Writing_tests

The latter is a draft of current thoughts and could probably use a seasoned developer to add/update/trim/correct.

Nicholas Jankowski <nrjank>
Group Member
Thu 27 Jan 2022 05:12:58 PM UTC, comment #27: 

Pushed a patch describing the `*` with bug numbers in the manual here:
https://hg.savannah.gnu.org/hgweb/octave/rev/c3494f408895

I hope the language is ok.

If you'd like to document best practice for Octave's sources in the Wiki, please go ahead.

Markus Mützel <mmuetzel>
Group administrator
Thu 27 Jan 2022 04:33:45 PM UTC, comment #26: 

ignore comment #24 - a grep shows no signifciant usage of xtest in core octave, and significant usage of %!test <#> and <*#>

Nicholas Jankowski <nrjank>
Group Member
Thu 27 Jan 2022 04:29:10 PM UTC, comment #25: 

at the risk of digressing further in this bug thread -

@markus - since contribution guidelines were moved out of the manual and now only exist on the wiki, that's where I had typcially referred to for BIST guidelines.  The guidance seems to be a bit bifurcated, and I notice neither mentions anything about <*bugid>. In the doc page the xtest does mention 'xtest <number>' is interpreted as a bug id. neither xtest or test <id> appear on the wiki. I agree the docs should capture all syntax related issues. whether suggested usage lives in either or both place is I guess up for some discussion.

Nicholas Jankowski <nrjank>
Group Member
Thu 27 Jan 2022 04:12:40 PM UTC, comment #24: 

Worth creating a separate bug report to replace all uses of xtest with test <bugid> to ensure better match of tests to bug reports?

Nicholas Jankowski <nrjank>
Group Member
Thu 27 Jan 2022 03:54:35 PM UTC, comment #23: 

The behavior of "test <bug-id>" and "xtest <bug-id>" is the same.  The only real need for xtest is to mark a test as a "known failure" when there is no associated bug report.  But we should try to avoid that and instead create bug reports for all known failures and just use "test <bug-id>" instead of "xtest".  That way, we have some documentation about why a test is known to fail.

For Octave, there is a script that checks the status of bug reports and adds the "*" to the tag.  The Makefile target "update-bug-status" will execute that script for all source files that contain tests that are tagged with bug numbers.

Also, please try to use the "Matlab Compatibility" item group only when Octave behavior is arguably reasonable but different from Matlab.  Otherwise, use an item group that most closely describes the actual problem.  For example, if the problem is that the feature is missing in Octave, then use "feature request" or if the bug is about incorrect behavior, use "incorrect result".

John W. Eaton <jwe>
Group administrator
Thu 27 Jan 2022 02:33:06 PM UTC, comment #22: 

I updated the last patch accordingly.

(file #52750)

Christof Kaufmann <christofkaufmann>
Thu 27 Jan 2022 08:49:20 AM UTC, comment #21: 

Imho, the documentation for the `test` syntax should better be done in the manual, not in the Wiki. Its syntax might change (and has changed) between versions of Octave and that might get messy to describe in the Wiki.

See:
https://octave.org/doc/v6.4.0/Test-Functions.html#Test-Functions

Afaict, it's currently not mentioned how the asterisk `*` affects the interpretation of the test result.

Markus Mützel <mmuetzel>
Group administrator
Thu 27 Jan 2022 08:04:10 AM UTC, comment #20: 

@Nik: What should be the difference between `xtest` with bug number and `test` with bug number (both without `*`)?

Markus Mützel <mmuetzel>
Group administrator
Thu 27 Jan 2022 03:00:53 AM UTC, comment #19: 

I think there was a discussion pointing out ho that worked, but could someone add the details to
https://wiki.octave.org/Tests

I assume xtests should have a bug # too?  Then a fixed bug should just change from xtest to test, adding  * to the bug number?

Nicholas Jankowski <nrjank>
Group Member
Thu 27 Jan 2022 12:08:50 AM UTC, comment #18: 

The '*' character indicates that the bug has been fixed.  If your patch fixes the problem, and then adds a BIST test to verify the fix, then it is appropriate to write


%!assert <*61863>



Rik <rik5>
Group administrator
Wed 26 Jan 2022 11:47:34 PM UTC, comment #17: 

OK, I added the bug number in the same style Markus added it to my patch from comment #2, i. e.

%!assert <61903> ...

However, in other tests I often found it prepended with a '*', like

%!assert <*61903> ...

What is this for?

(file #52745)

Christof Kaufmann <christofkaufmann>
Wed 26 Jan 2022 03:56:58 PM UTC, comment #16: 

@Christof: Yes, adding the bug numbers would be great!

Rik <rik5>
Group administrator
Wed 26 Jan 2022 01:57:53 PM UTC, comment #15: 

I also fixed `interpn()` and `__lin_interpn__()` in related bug #61907 now.

I notices, I forgot the bug number in the `assert`s in `interp1()` in the last patch here in comment #12. Do you want me to add them and update the patch?

Christof Kaufmann <christofkaufmann>
Tue 25 Jan 2022 05:20:59 PM UTC, comment #14: 

@Markus: I agree.  There shouldn't be a problem, but let's put the m-file fix on stable and the C++ fix on the development branch.

Rik <rik5>
Group administrator
Tue 25 Jan 2022 04:55:55 PM UTC, comment #13: 

I agree that fixing `__pchip_deriv__` is the "real" fix for this issue. (See also the FIXME note I added in the patch for `interp1`.)
However, I'm worried (maybe without reason) that this change might be too big for the stable branch (not too long before a major release).
Scrolling through the patch, I don't see any major issues. But one never knows...
Maybe, we could take the "saver" change in the .m file to fix the bug on the stable branch and then apply the more general fix on default.

@Rik: What's your opinion on this?

Markus Mützel <mmuetzel>
Group administrator
Tue 25 Jan 2022 04:35:44 PM UTC, comment #12: 

I set up the tools and prepared the attached patch to fix `interp1()` and `interp2()` for complex inputs, since both functions are affected by modifying `__pchip_deriv__()`. I will fix `interpn()` in bug #61907.

Regarding the handling of the complex matrices in `__pchip_deriv__()`, I would like to hear your opinion. To avoid copies, I used `reinterpret_cast<T*>(mat.data ())` for the real part and `reinterpret_cast<T*>(mat.data ()) + 1` for the imaginary part. This seems a bit hacky, but the standard covers it:
https://en.cppreference.com/w/cpp/numeric/complex#Array-oriented_access

The alternative way would be to split up `ymat` into `ymat_real` and `ymat_imag` and similarly for `dmat`. The splitting required a copy. So, if you prefer this, let me know.

(file #52729)

Christof Kaufmann <christofkaufmann>
Tue 25 Jan 2022 09:01:32 AM UTC, comment #11: 

Regarding the performance of `complex()`: In my patch it was `NA + 1i*NA`, but I think it doesn't matter, since it is just a single scalar.

Yes, I agree, that it is reasonable to fix it in the C++ code. I just have not set up the tools for compilation and have not looked at the C++ API. There is even one more file, that has to be considered. Here is an overview:

  • `interp1()` uses `pchip()` and pchip()` uses `__pchip_deriv__()`
  • `interp2()` uses `__pchip_deriv__()` directly
  • `interpn()` might use `__pchip_deriv__()` in future, but "pchip" is not yet implemented.
  • `interpn()` uses `__lin_interpn__()`, which also ignores the imaginary part.

So libinterp/corefcn/__lin_interpn__.cc and libinterp/corefcn/__pchip_deriv__.cc would have to be fixed.

Christof Kaufmann <christofkaufmann>
Mon 24 Jan 2022 10:03:55 PM UTC, comment #10: 

Since all of the various interpX functions call _pchip_deriv_(), wouldn't it be better to change that C++ routine to handle complex numbers rather than solving it three times in interp1.m, interp2.m, interpn.m?

Rik <rik5>
Group administrator
Mon 24 Jan 2022 10:02:21 PM UTC, comment #9: 

On a related note, see bug #61912.  The performance of complex() is slow for some reason.  For performance, it would be better to write


extrap = NA + NA*i;


rather than


extrap = complex (NA, NA);


Rik <rik5>
Group administrator
Mon 24 Jan 2022 11:08:30 AM UTC, comment #8: 

This bug report is about interp2, so I opened new bug reports regarding complex matrices for

These contain patches already.

Christof Kaufmann <christofkaufmann>
Sun 23 Jan 2022 03:34:55 PM UTC, comment #7: 

@Christof: In general, we try to use one report for one issue. If there are follow up changes needed to fix this issue you can post them here. (Or open a new report. Whatever you prefer.) If it is fixing a different issue, please open a new report.

Markus Mützel <mmuetzel>
Group administrator
Sun 23 Jan 2022 03:32:59 PM UTC, comment #6: 

Re comment #5: One was fixing a bug, the other was just style.

Markus Mützel <mmuetzel>
Group administrator
Sun 23 Jan 2022 12:55:57 PM UTC, comment #5: 

Why was one fix applied to stable and another to default?

Anonymous
Sun 23 Jan 2022 12:22:04 PM UTC, comment #4: 

Thanks! BTW: The style changes were necessary, when I had the `i` instead of `1i`.

When I prepare similar patches for interp functions (like `interpn`), should I open a new bug report or add it here?

Christof Kaufmann <christofkaufmann>
Sun 23 Jan 2022 11:31:56 AM UTC, comment #3: 

Thanks again.

I split the patch into two parts. One with the bug fixes and another one with the style change.
Pushed the part with the bug fixes to the stable branch here:
https://hg.savannah.gnu.org/hgweb/octave/rev/27566803b3b1

And the part with the style change here:
https://hg.savannah.gnu.org/hgweb/octave/rev/9b8c74d85fa7

Looking at this a little more closely, I'm wondering if it would make sense to move the respective change into `__pchip_deriv__.cc`. (Added a FIXME note for that.)

Marking as ready for test.

Markus Mützel <mmuetzel>
Group administrator
Sat 22 Jan 2022 10:11:24 AM UTC, comment #2: 

OK, replaced all `i`s by `1i`s. Modified patch attached.

Regarding the contributors list, if you like, you can add me. Yes, Christof Kaufmann is my real name.

(file #52706)

Christof Kaufmann <christofkaufmann>
Fri 21 Jan 2022 03:54:12 PM UTC, comment #1: 

Looks good to me. Just some minor remark:
It's generally a good idea to write `1i` instead of `i` when you mean the imaginary unit `sqrt(-1)`. That way, you eliminate the risk if lateron a (loop) variable `i` (which imho should also be avoided for the same reason) is used somewhere else in that function.

Would you like to be added to the list of contributors that appears close to the front of the manual? Is Christof Kaufmann your real name?


Markus Mützel <mmuetzel>
Group administrator
Thu 20 Jan 2022 01:55:17 AM UTC, original submission:  

Interpolating complex matrices gives a default extrapolation value of NA + 0i instead of NA + NAi. Example:

> interp2 (i*eye(2), 0:1, 1)
ans =

    NA +   0i     0 +   1i


Additionally "pchip" gives wrong results for the imaginary part:

> interp2(i*eye(2), 1.25, 1, "pchip")
ans =       0 + 0.8438i

while it should give:

interp2(eye(2), 1.25, 1, "pchip") * i
ans =       0 + 0.7500i


The appended patch fixes interp2 for complex matrices.

Christof Kaufmann <christofkaufmann>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #52760:  interp1_interp2_complex_fix4.patch added by christofkaufmann (12KiB - text/x-patch - Fix interp1 and interp2 on current default branch - v4)
file #52750:  interp1_interp2_complex_fix3.patch added by christofkaufmann (19KiB - text/x-patch - Fix interp1 and interp2 on current default branch - v2 (added asterisk to bug numbers))
file #52745:  interp1_interp2_complex_fix2.patch added by christofkaufmann (19KiB - text/x-patch - Fix interp1 and interp2 on current default branch - v2 (added bug number))
file #52729:  interp1_interp2_complex_fix.patch added by christofkaufmann (18KiB - text/x-patch - Fix interp1 and interp2 on current default branch - v1)
file #52706:  interp2_fix_complex2.diff added by christofkaufmann (2KiB - text/x-patch - Fix interp2 for complex matrices - v2 (Use 1i instead of i))
file #52692:  interp2_fix_complex1.diff added by christofkaufmann (2KiB - text/x-patch - Fix interp2 for complex matrices - v1)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by jwe (Posted a comment)
  • -email is unavailable- added by nrjank (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by mmuetzel (Posted a comment)
  • -email is unavailable- added by christofkaufmann (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 11 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2022-01-31 rik5 StatusReady For Test Fixed
        Open/ClosedOpen Closed
    2022-01-27 christofkaufmann Attached File- Added interp1_interp2_complex_fix4.patch, #52760
    2022-01-27 jwe Item GroupMatlab Compatibility Incorrect Result
    2022-01-27 christofkaufmann Attached File- Added interp1_interp2_complex_fix3.patch, #52750
    2022-01-26 christofkaufmann Attached File- Added interp1_interp2_complex_fix2.patch, #52745
    2022-01-25 christofkaufmann Attached File- Added interp1_interp2_complex_fix.patch, #52729
    2022-01-23 mmuetzel StatusPatch Submitted Ready For Test
    2022-01-22 christofkaufmann Attached File- Added interp2_fix_complex2.diff, #52706
    2022-01-21 mmuetzel StatusNone Patch Submitted
    2022-01-20 christofkaufmann Attached File- Added interp2_fix_complex1.diff, #52692

    Back to the top

    Powered by Savane 3.13-02a9.
    Corresponding source code