Pyminuit Richardson lucy implementation

[GoogleCodeExporter]

I am trying to use Pyminuit to iterate the richardson lucy(RL) algorithm as a 
possible method of image deconvolution. However I am getting a few errors 
through migrad(). I have tried two RL lambda functions, one using a summation 
method and the other a convolution method. Both can be found here, 
http://en.wikipedia.org/wiki/Richardson%E2%80%93Lucy_deconvolution. Here is my 
code snippet. 

    def Richardsonlucy_deconvolve_min(image_data='/Users/krishnamooroogen/Documents/physics/MSSL_projectfiles/VSOdata/AIA/2011_06_07/aia_lev1_211a_2011_06_07t06_27_36_63z_image_lev1_fits.fits' ,psf='/Users/krishnamooroogen/Documents/physics/MSSL_projectfiles/psf/psf_211_AIA.fits'):
    #opens image data and psf
    #creates inverse psf
    #sets lambda function (ricahrdson lucy)
    #sets start value
    #iterates gradient

        image_data = pyfits.open(image_data)[0].data
        psf = pyfits.open(psf)[0].data
        psf_hat = psf[::-1]

        #RL = lambda latent: latent*(convolute((image_data/(convolute(latent,psf))),psf_hat))
        RL = lambda latent: latent*(((image_data/(latent*psf).sum())*psf).sum())

        m = minuit.Minuit(RL,latent = numpy.array([numpy.array([1.0]*len(image_data)),numpy.array([1.0]*len(image_data))]*int(len(image_data)/2)))
        m.tol = 1000000000
        m.printMode = 1
        m.migrad()
        m.hesse()
        latent = m.values['latent']

The convolute function is here, 

    def convolute(input,input2):

    fourier_space_convolution = (numpy.fft.rfft2(input,input.shape))*(numpy.fft.rfft2(input2,input2.shape))
    convert_realspace_convolution = numpy.fft.fftshift(numpy.fft.irfft2(fourier_space_convolution,fourier_space_convolution.shape))
    return convert_realspace_convolution

When using the summation method of RL I get the following error, 

    --> 359     Richardsonlucy_deconvolve_min()

    /Users/krishnamooroogen/Documents/Modules/MSSL_solar/solar_usr.pyc in Richardsonlucy_deconvolve_min(image_data, psf)
    276         m.tol = 1000000000
    277         m.printMode = 1
    --> 278         m.migrad()
    279         m.hesse()
    280         latent = m.values['latent']

    /Users/krishnamooroogen/Documents/Modules/MSSL_solar/solar_usr.pyc in <lambda>(latent)
    271 
    272         #RL = lambda latent: latent*  (convolute((image_data/(convolute(latent,psf))),psf_hat))
    --> 273         RL = lambda latent: latent*(((image_data/(latent*psf).sum())*psf).sum())
    274 
    275         m = minuit.Minuit(RL,latent = numpy.array([numpy.array([1.0]*len(image_data)),numpy.array([1.0]*len(image_data))]*int(len(image_data)/2)))

    TypeError: only length-1 arrays can be converted to Python scalars 

The error when using the convolute method, 

    AttributeError: 'float' object has no attribute 'shape' 


However outside of minuit, both the convolute function and the RL summation 
method for the latent image estimation work fine. Can someone shed some light, 
is this a Minuit specific problem? Can it be resolved? 

Original issue reported on code.google.com by kpmooroo...@gmail.com on 21 Nov 2013 at 1:00

[GoogleCodeExporter]

(It has been a while since I've worked on this.)  Is the output of the 
objective function non-scalar, such as [1.0, 2.0, 3.0]?  Minuit searches for a 
minimum in a f: *float -> float function (potentially many input floats, one 
output float).  If you're outputting a vector, then you probably want to output 
the norm of that vector.

If that's not the problem, perhaps the scalar you're outputting is a Numpy 
number.  Scalars constructed by Numpy are often `<type 'numpy.float64'>` rather 
than `<type 'float'>` (try `type(myNumber)`.).  PyMinuit uses 
`PyFloat_AsDouble(pyObject)` to read in numbers and `PyFloat_FromDouble(x)` to 
write them out--- these functions are not explicitly Numpy-aware and might not 
work.  The `TypeError` (which complains about not being able to produce a 
Python scalar from Numpy) and the `AttributeError` (which complains about not 
being able to treat a Python scalar as a Numpy object) sound suspiciously like 
Numpy incompatibilities.  If so, then the solution is simple: cast the objects 
with `float(npObject)` and `numpy.array(pyObject)` or 
`numpy.float64(pyObject)`, depending on whether you want a one-element array or 
a Numpy number.

(Note that when PyMinuit was written, Numpy didn't exist or it wasn't mature--- 
I remember that NumArray and Numeric were still being used.  PyMinuit doesn't 
do anything explicitly with large arrays (though you could feed it a large 
array, one element at a time), so there's no reason for it to be Numpy-aware.  
It works with Python floats.)

If this doesn't work out for you, consider the SciPy optimization functions 
(http://docs.scipy.org/doc/scipy/reference/optimize.html), which also have 
general non-local minimizers.  The SciPy optimizers look to me like basic 
implementations of well-known optimization algorithms (with names and 
references to the literature).  Minuit is a kitchen sink of algorithms that 
switches modes when it senses that the search might be diverging, and it has a 
long history of special cases built in to handle problems that people have 
encountered in the past.  That can be good or bad, depending on your 
application (good because it can recover from common problems, bad because its 
internal workings are less transparent).

Original comment by jpivar...@gmail.com on 21 Nov 2013 at 6:25

[GoogleCodeExporter]

I think you're right, it makes sense for it to be a numpy compatibility 
problem. The data sent in are of type numpy.array but the data in the array are 
floats. Can you clarify what you mean by `float(npObject)` and 
`numpy.array(pyObject)` or `numpy.float64(pyObject). Thanks a lot for your 
help. 

Original comment by kpmooroo...@gmail.com on 22 Nov 2013 at 12:22

[GoogleCodeExporter]

I mean convert the Numpy scalar to a Python float using the float() function, 
and convert the other way with numpy.array or numpy.float64.

Original comment by jpivar...@gmail.com on 22 Nov 2013 at 1:30

[GoogleCodeExporter]

Thanks, I managed to get the second implementation to work(the one without the 
convolution) with your suggestions. Though this led to the minuit convarience 
matrix error which i believe is a start value issue? 

Original comment by kpmooroo...@gmail.com on 22 Nov 2013 at 7:16