SKA-ScienceDataProcessor / algorithm-reference-library

Reference code for simplified aperture synthesis imaging
Apache License 2.0
26 stars 21 forks source link

MemoryError in predict_skycomponent_visibility #55

Closed rstofi closed 5 years ago

rstofi commented 5 years ago

I believe there is an indexing error in predict_skycomponent_visibility() because I get the following error:

File "/home/krozgonyi/algorithm-reference-library/processing_components/imaging/base.py", line 221, in predict_skycomponent_visibility
    vis.data['vis'][...] += comp[:,:] * phasor[:, numpy.newaxis]
MemoryError

The arrays vis.data['vis'][...] and comp[:,:] both have a shape of (N,N_pol) However, the array phasor[:,numpy.newaxis] have shape of (N,1,1) Thus, the multiplication will yield an array with shape of (N,N,N_pol)

The resulting ndarray is too big even for small vis objects, so the MemoryError occurs before a Value Error of shape inconsistency raise.

Getting rid of the newaxis solves the problem:

vis.data['vis'][...] += comp[:,:] * phasor[:]

timcornwell commented 5 years ago

I cannot reproduce this which is probably an oversight in my test code. There is a branch issue_55 with a slight change to predict_skycomponent_visibility and a separate test case trying to isolate this problem. Can you have a look at the test program processing_components/test_predict_skycomponent_visibility.py and compare it to your code?

Thanks, Tim

rstofi commented 5 years ago

I create visibilities from an MS, which is linearly polarized, and it only have one channel. Might be how a single channel handled... I can e-mail you an example code and the MS if you want to have a look!

timcornwell commented 5 years ago

Sure - send me the MS and your code.

rstofi commented 5 years ago

Actually, I cut my code to the minimum and put the MS here. So everyone can test it!

The code:

#The base code produce the MemoryError
#=== Imports ===
import os,numpy
from wrappers.arlexecute.visibility.base import create_visibility_from_ms
from wrappers.serial.imaging.base import predict_skycomponent_visibility
from wrappers.serial.skycomponent.operations import create_skycomponent

#=== MAIN ===
MSname = os.getcwd() + '/memory_error_example.ms'
vis = create_visibility_from_ms(MSname)[0]
skycomponent_list = [create_skycomponent(flux=numpy.array([numpy.ones(4)]),frequency=numpy.array([vis.frequency[0]]), 
    direction=vis.phasecentre, polarisation_frame=vis.polarisation_frame)]

vis = predict_skycomponent_visibility(vis,skycomponent_list)

I got the MemoryError with this setup.

timcornwell commented 5 years ago

I found the problem: the MSreader was creating an array of SkyCoords rather than a SkyCoord as the phase centre. This propagated to the calculation of the lmm coordinates.

The fix is in the branch issue_55. Once the build completes, I will commit to the master (a couple of hours).

timcornwell commented 5 years ago

Now in master