bugGNU Octave - Bugs: bug #64733, IFFT returns wrong data for...

 
 

bug #64733: IFFT returns wrong data for dNDArray input

Submitter:  Hg200 <hg200>
Submitted:  Sun 01 Oct 2023 02:04:27 PM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Fixed Assigned to:  None
Originator Name:  hg200 Open/Closed:  * Closed
Release:  * 8.3.0 Operating System:  * Any
Fixed Release:  8.4.0 Planned Release:  8.4.0
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Wed 04 Oct 2023 05:23:32 PM UTC, comment #13: 

Thanks for testing.

Closing as fixed.

Markus Mützel <mmuetzel>
Group administrator
Tue 03 Oct 2023 06:55:18 PM UTC, comment #12: 

Thank you. I have tested the patch and i haven't found any problems. From my side the issue can be closed.

Hg200 <hg200>
Tue 03 Oct 2023 12:46:45 PM UTC, comment #11: 

Thanks very much.

I pushed the change to the stable branch:
https://hg.savannah.gnu.org/hgweb/octave/rev/903c9100178b

Marking as ready for test.

Markus Mützel <mmuetzel>
Group administrator
Tue 03 Oct 2023 08:56:36 AM UTC, comment #10: 

The extra space for alignment (+32) is also divided by 2. Attaching a new patch that should avoid this (untested). What do you think?

(file #55194)

Hg200 <hg200>
Mon 02 Oct 2023 08:15:04 PM UTC, comment #9: 

I haven't checked if i am right, but i feel

OCTAVE_SCOPED_BUFFER (Complex, itmp,(nn * howmany + 32) / ((! in_place) + 1))

+

ifft(fft(rand (3, 5, 2), [], 2), [], 2)

then nn=5 and hownmany=3, so 3*5/2 is rounding towards zero, so buffer reserved for 14 doubles (after pointercast), but size must be 15. Need to verify this.

Hg200 <hg200>
Mon 02 Oct 2023 06:38:59 PM UTC, comment #8: 

Oops. Couple of copy-paste errors. Please, ignore the patch from comment #7.

Updated version here.


(file #55189)

Markus Mützel <mmuetzel>
Group administrator
Mon 02 Oct 2023 06:34:13 PM UTC, comment #7: 

Thank you for pointing that out.

What about the attached patch for that case?

(file #55188)

Markus Mützel <mmuetzel>
Group administrator
Mon 02 Oct 2023 06:02:52 PM UTC, comment #6: 

Okay i have just finished testing the first patch "bug64733-fft-planner-inplace.patch". I haven't checked the second patch yet. I will do this afterwards, but here is an intermediate result.

Seems everything to be correct. There is only one question for the case "R2C + MEASURE + Inplace Transform". The FFTW documentation is not very specific for R2C but it states that inplace transform is possible for R2C. You run the planner with the input pointing to "itmp" and the output pointing to "otmp". However, since this is an Inplace Transform case, I suspect "in" and "out" should point to the same address. In this case, it probably needs to be the address for the complex output buffer array "otmp", since the data will inflate (get larger due to the conversion from real to complex) during the transform. I had to modify some Octave code to test this case, but it seems to work. The code is executed correctly and it gives the correct result. Any comments on this?

OCTAVE_LOCAL_BUFFER (Complex, otmp, nn * howmany + 32);
otmp = reinterpret_cast<Complex *>
       (((reinterpret_cast<std::ptrdiff_t> (otmp) + 15) & ~ 0xF) +
       ((reinterpret_cast<std::ptrdiff_t> (out)) & 0xF));
*cur_plan_p
     = fftw_plan_many_dft_r2c (rank, tmp, howmany, reinterpret_cast<double *> (otmp),
                               nullptr, stride, dist,
                               reinterpret_cast<fftw_complex *> (otmp),
                               nullptr, stride, dist, plan_flags);

Maybe it is also helpful to rename "plan_destroys_in" to "plan_destroys_arrays" or s.th. similar.

Hg200 <hg200>
Mon 02 Oct 2023 05:50:41 PM UTC, comment #5: 

Attached is a patch with the changes that reduce the code duplication in that area. I'm not sure if it can go on the stable branch. It adds new macros in a public header. But it doesn't remove or modify existing ones. So, it could still be ok.

It also fixes the issue from comment #0 for me.
Actually, the patch from comment #4 contains superfluous parts because for *_r2c operations, the output buffer can never be the input buffer (because it wouldn't be large enough).

Still not much more testing than running the test suite with it.


(file #55187)

Markus Mützel <mmuetzel>
Group administrator
Mon 02 Oct 2023 01:26:27 PM UTC, comment #4: 

I'm not sure what the intensions of the original author were. I can only guess.
It is true that for the "end usage" in Octave, there is probably always a allocation for the result that is separate from the input. I'm not sure what is happening for inplace operations like `x = ifft(x);`.
Technically, the functions in liboctave could also be used by code outside the Octave interpreter. In that case, in-place operation might also be of advantage. (I don't know how far our support for usage outside Octave goes though.)

I agree that changing how we use this function in liboctinterp could also work around this issue. The issue in liboctave would still remain though (if it is an issue).

I was considering this change for the stable branch (because it should have low risk of breaking things and wouldn't change the API).
I was thinking about introducing new preprocessor macros similar to OCTAVE_LOCAL_BUFFER that would allow extending the scope of a buffer. If that works, it could allow reducing the current code duplication in this area. That change would need to go to the default branch though.

I haven't fully explored how these functions are used though. So, please let me know if I'm missing something.

Markus Mützel <mmuetzel>
Group administrator
Sun 01 Oct 2023 10:36:57 PM UTC, comment #3: 

Thank you very much for your fast response. I have read some FFTW documentation. I think your proposal goes in the direction of what the oct-fftw.cc author originally intended with the "plan_destroys_in" variable: If the FFTW planner is in "ESTIMATE", then it creates a default plan. "in" and "out" arrays remains unchanged. In this case the "else" branch should be executed. Where as in the "MEASURE" state, the planner performs test transformations to improve its plan. Then the "in" and "out" arrays are destroyed. In this case, the "if" branch should be executed. But only the handling of the input array is implemented. I also agree that the code becomes convoluted now.

Comment: If the Octave command "y=ifft(x)" is executed, then finally two different memory spaces must exist. One for the input array "x" and a separate one for the ouptut array "y". So what is the benefit of an in place transformation, if in the end a copy action of the result to "y" has to be performed anyway? You can do this copy action directly with an out of place transformation. Are there any other disadvantages if we would change everything to out of place transformations?

Side note: Actually "plan_destroys_in" exactly tells us what the code is doing. Translation: "This var indicates when the planner destroys its input array". Nothing is said about output array ;-)

Hg200 <hg200>
Sun 01 Oct 2023 07:35:57 PM UTC, comment #2: 

Please, ignore my comment #1. I didn't understand the issue correctly.

Attached is a patch that might avoid this issue by creating a buffer for the output if necessary.
I don't really like how much code duplication there is. That might be avoidable if we had a variant of OCTAVE_LOCAL_BUFFER for which it would be possible to select a wider scope. But that's a different issue entirely.

With those changes, your examples are returning the correct results on the first invocation for me.
The test suite still passes here. But I didn't do any extensive testing yet.


(file #55182)

Markus Mützel <mmuetzel>
Group administrator
Sun 01 Oct 2023 03:18:20 PM UTC, comment #1: 

Thank you for the report and the analysis.

This might be a very old bug going back more than 15 years...

Would initializing `plan_destroys_in` with `in == out` (instead of setting it to `true`) do the right thing? I haven't tested that though.

Markus Mützel <mmuetzel>
Group administrator
Sun 01 Oct 2023 02:04:27 PM UTC, original submission:  

Start a new Octave session. I have used octave-8.3.0. Then switch the planner to "measure" and perform an inverse transformation with a real-valued input (dNDArray). The return value of the ifft command is zero if the planner performs a measurement in the current matrix dimension. if "ifft" is called again, then the correct result appears. if "ifft" is called again with a new NDArray size, then zero appears again, and so on.

octave:1> fftw ('planner','measure');
octave:2> ifft ([2, 4, 6, 8])
ans =

   0   0   0   0

octave:3> ifft ([2, 4, 6, 8])
ans =

   5 + 0i  -1 - 1i  -1 + 0i  -1 + 1i

octave:4> ifft ([2, 4, 6, 8; 10, 12, 14, 16])
ans =

   0   0   0   0
   0   0   0   0

octave:5> ifft ([2, 4, 6, 8; 10, 12, 14, 16])
ans =

    6    8   10   12
   -4   -4   -4   -4

My interpretation: The NDArray datatype uses an inplace transform where as the CNDarray data type (which is not affected) uses an out-of-place. At the other side, the planner input points to fresh reallocated memory ("itmp"), but its output "out" points to the input data, that is to be transformed (because of the in place transform). When the planner performs a measurement, it overwrites the input data and thus the zero sequence is obtained. This happens here:

https://hg.savannah.gnu.org/hgweb/octave/file/43d010974a89/liboctave/numeric/oct-fftw.cc#l256

One way to solve this would be to change from "in place" to "out of place" in dNDArray.cc (and in the associated short data type). That solves the problem for me. But my question is if anyone has a more elegant solution?

Hg200 <hg200>

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #55189:  bug64733-fft-planner-inplace-v3a.patch added by mmuetzel (11KiB - application/octet-stream)
file #55188:  bug64733-fft-planner-inplace-v3.patch added by mmuetzel (11KiB - application/octet-stream)
file #55187:  bug64733-fft-planner-inplace-v2.patch added by mmuetzel (11KiB - application/octet-stream)
file #55182:  bug64733-fft-planner-inplace.patch added by mmuetzel (7KiB - application/octet-stream)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by mmuetzel (Posted a comment)
  • -email is unavailable- added by hg200 (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 12 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2023-10-04 mmuetzel StatusReady For Test Fixed
        Open/ClosedOpen Closed
    2023-10-03 mmuetzel StatusPatch Submitted Ready For Test
        Fixed ReleaseNone 8.4.0
        Planned ReleaseNone 8.4.0
    2023-10-03 hg200 Attached File- Added bug64733-fft-planner-inplace-v3b.patch, #55194
    2023-10-02 mmuetzel Attached File- Added bug64733-fft-planner-inplace-v3a.patch, #55189
    2023-10-02 mmuetzel Attached File- Added bug64733-fft-planner-inplace-v3.patch, #55188
    2023-10-02 mmuetzel Attached File- Added bug64733-fft-planner-inplace-v2.patch, #55187
    2023-10-01 mmuetzel Attached File- Added bug64733-fft-planner-inplace.patch, #55182
        StatusNone Patch Submitted
        Operating SystemGNU/Linux Any

    Back to the top

    Powered by Savane 3.13-4b48.
    Corresponding source code