Issue with 2D arrays and RM-CLEAN

[AlecThomson]

For 2D arrays (e.g. HEALPIX), I think there is a bug currently present in RM-CLEAN. Specifically, in util_RM.py -> do_rmclean_hogbom, lines 504 and 505. Currently, they read:

elif nDims==2:
        dirtyFDF = np.reshape(list(dirtyFDF.shape[:2])+[1])
        RMSFArr = np.reshape(list(RMSFArr.shape[:2])+[1])

This is syntactically incorrect for np.reshape. I think they should be:

elif nDims==2:
        dirtyFDF = np.reshape(dirtyFDF, list(dirtyFDF.shape[:2])+[1])
        RMSFArr = np.reshape(RMSFArr, list(RMSFArr.shape[:2])+[1])

Following this, in the main loop (i.e. lines 546+) there is an issue with the indexing of 2D arrays. I found that for a 2D array of shape (N_phi_Chan, N_pix), as produced by RM-Synthesis, the indexing variables xi and yi were not in the correct order for lines 586-588, specifically:

# Restore the CC * a Gaussian to the cleaned FDF
  cleanFDF[:, yi, xi] += \
  gauss1D(CC, phiPeak, fwhmRMSFArr[yi, xi])(phiArr_radm2)

That is, for fwhmRMSFArr the order was flipped. I was able to temporarily sidestep this issue by replacing this line with:

# Restore the CC * a Gaussian to the cleaned FDF
  cleanFDF[:, yi, xi] += \
  gauss1D(CC, phiPeak, fwhmRMSFArr.T[yi, xi])(phiArr_radm2)

But, this seems like it might be missing the underlying issue.

[Cameron-Van-Eck]

Hi Alec: Thanks a lot for testing this out; I don’t think anyone’s run HEALPIX files through RM-tools yet. Also thanks a million for narrowing down the specific problems in the code. I don’t have any HEALPIX data to test the bug and fix on. Could you send me some? That’ll make it a lot easier to confirm that things are fixed. If the data’s too large to send easily, just a few channels would be enough. Thanks, Cameron

…

On Mar 18, 2019, at 11:03 PM, Alec Thomson ***@***.***> wrote: For 2D arrays (e.g. HEALPIX), I think there is a bug currently present in RM-CLEAN. Specifically, in util_RM.py -> do_rmclean_hogbom, lines 504 and 505. Currently, they read: elif nDims==2: dirtyFDF = np.reshape(list(dirtyFDF.shape[:2])+[1]) RMSFArr = np.reshape(list(RMSFArr.shape[:2])+[1]) This is syntactically incorrect for np.reshape. I think they should be: elif nDims==2: dirtyFDF = np.reshape(dirtyFDF, list(dirtyFDF.shape[:2])+[1]) RMSFArr = np.reshape(RMSFArr, list(RMSFArr.shape[:2])+[1]) Following this, in the main loop (i.e. lines 546+) there is an issue with the indexing of 2D arrays. I found that for a 2D array of shape (N_phi_Chan, N_pix), as produced by RM-Synthesis, the indexing variables xi and yi were not in the correct order for lines 586-588, specifically: # Restore the CC * a Gaussian to the cleaned FDF cleanFDF[:, yi, xi] += \ gauss1D(CC, phiPeak, fwhmRMSFArr[yi, xi])(phiArr_radm2) That is, for fwhmRMSFArr the order was flipped. I was able to temporarily sidestep this issue by replacing this line with: # Restore the CC * a Gaussian to the cleaned FDF cleanFDF[:, yi, xi] += \ gauss1D(CC, phiPeak, fwhmRMSFArr.T[yi, xi])(phiArr_radm2) But, this seems like it might be missing the underlying issue. — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#4>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AIipy3cWw7TnEUwj_gb5KKPqovFX5ifGks5vYFONgaJpZM4b7Lls>.

[AlecThomson]

Hi Cameron,

Here's a script for making some dummy HEALpix data. This script saves out both cubes (2-D array) and maps (per frequency). The latter is the current HEALpix standard, but more cumbersome for broadband data sets.

I encountered the above issue when running the 2-D array formatted data. This data can be generated from either reading in a FITS cube or stacking the arrays from a bunch of maps.

There are downsides to each approach here. The FITS cube is not a HEALpix standard, but stores the data more efficiently. As HEALpix currently doesn't natively support cubes, saving individual maps preserves the standard, but is a pain for broadband data.

Cheers,
Alec

[Cameron-Van-Eck]

Hi Alec:

So I found that the problem with the reshape commands was already fixed in the CIRADA branch (based on the line numbers, it looks like you've been using the master branch?).

I'm not seeing any problems with indices in the RMSynth or RMClean codes when I run them. The outputs look as I expect them to.

I've spent the last few days redoing a lot of the input/output code to have more flexibility with the number of dimensions and also determining which axis is frequency (previously it was assuming the last axis was frequency, which wasn't the case for data produced by the HEALPIX dummy script?).

Can you try the most recent version (of the CIRADA branch) and see if you still get the same problems?

Thanks,
Cameron

[AlecThomson]

Hi Cameron,

Whoops, yes sorry the numbers were from the master branch, but I think the bug is still present in CIRADA on lines 509/510.

I just tested on the dummy cubes I produces and RMsyth worked great. I'll note that the frequency first part was purely my choice for speed in the other processing I was doing and not a standard at all. For example GMIMS has frequency last. But, other codes such as cuFFs requires frequency first for speed, so I think the flexibility for both is great.

I now get the following error when attempting to run RMclean on the output of synth:

Traceback (most recent call last):
  File "/home/athomson/bin/RM-tools/RMtools_3D/do_RMclean_3D.py", line 214, in <module>
    main()
  File "/home/athomson/bin/RM-tools/RMtools_3D/do_RMclean_3D.py", line 99, in main
    write_separate_FDF=args.write_separate_FDF)
  File "/home/athomson/bin/RM-tools/RMtools_3D/do_RMclean_3D.py", line 140, in run_rmclean
    doPlots          = False)
  File "/home/athomson/bin/RM-tools/RMtools_3D/RMutils/util_RM.py", line 472, in do_rmclean_hogbom
    nPhi = phiArr_radm2.shape[0]
AttributeError: 'list' object has no attribute 'shape'

[Cameron-Van-Eck]

Hi Alec:

Aha, I forgot to update RMclean_3D when I updated RMSynth_3D. I've done that today. I found the exactly problems you reported, and they should now be fixed. Rather than using the transpose on fwhmRMSFArr, I redefined it to have the expected shape (for 2D input data, it appends an additional (degenerate) dimension; but it was putting it on the wrong axis and thus breaking things).

So it should be working now? I think it's working for the (FREQ x HEALPIX) data, and I don't think I've broken anything for (POSITION x POSITION x FREQ) input cubes. I'm not entirely sure if my rapidly spliced-together new read-in code will work as intended for other possible configurations, but that's a problem for another day.

Cheers,
Cameron