SLIT V1.0

Examples/Test_SLIT.py

@@ -97,7 +97,6 @@ def SIE(x0,y0,n1,n2,b,beta,q,xc,theta):
                                            -sourcesl[np.where(real_source!=0)])**2
                                            /real_source[np.where(real_source!=0)]**2)/(np.size(
                                            np.where(real_source!=0))/2.)
-
 image_chi2 = np.std(Image-Imsl)**2/sigma**2
 print('Residuals in source space', source_error)
 print('Residuals in image space',image_chi2)

Examples/Test_SLIT_MCA.py

@@ -63,15 +63,12 @@ def SIE(x0,y0,n1,n2,b,beta,q,xc,theta):
 
 #Generation of lensed source
 I2 = slit.Lens.source_to_image(newsource, nt1 ,nt2 , Fkappa)
-
+HI2 = scp.fftconvolve(I2, PSF.astype(float), mode = 'same')
 
 #Noise levels
 SNR = 500
 sigma = np.sqrt(np.sum(I2**2)/SNR/(nt1*nt2*size**2))
 
-#Convolution by the PSF and generation of the final image
-I2 = scp.fftconvolve(I2, PSF, mode = 'same')
-
 #Convolution of the observed image
 simu = scp.fftconvolve(gal0.astype(float)+I2, PSF.astype(float), mode = 'same')
 #Sotring the convolved lens light profile:
@@ -94,7 +91,7 @@ def SIE(x0,y0,n1,n2,b,beta,q,xc,theta):
 start = time.clock()
 
 #Running SLIT_MCA
-sourcesl, Imsl, G = slit.SLIT_MCA(Image, Fkappa, kmax, niter,riter, size,PSF, PSFconj, levels = levels)
+S, FS, G = slit.SLIT_MCA(Image, Fkappa, kmax, niter,riter, size,PSF, PSFconj, levels = levels)
 
 #Stop clock
 elapsed = (time.clock()-start)
@@ -104,11 +101,11 @@ def SIE(x0,y0,n1,n2,b,beta,q,xc,theta):
 real_source = newsource
 
 source_error = np.sum(np.abs(real_source[np.where(real_source!=0)]
-                                           -sourcesl[np.where(real_source!=0)])**2
+                                           -S[np.where(real_source!=0)])**2
                                            /real_source[np.where(real_source!=0)]**2)/(np.size(
                                            np.where(real_source!=0))/2.)
 
-image_chi2 = np.std(Image-Imsl-G)**2/sigma**2
+image_chi2 = np.std(Image-FS-G)**2/sigma**2
 print('Residuals in source space', source_error)
 print('Residuals in image space',image_chi2)
 
@@ -117,7 +114,7 @@ def SIE(x0,y0,n1,n2,b,beta,q,xc,theta):
 
     plt.figure(0)
     plt.title('Source from SLIT')
-    plt.imshow((sourcesl), vmin = np.min(real_source), vmax = np.max(real_source), cmap = cm.gist_stern, interpolation = 'nearest')
+    plt.imshow((S), vmin = np.min(real_source), vmax = np.max(real_source), cmap = cm.gist_stern, interpolation = 'nearest')
     plt.axis('off')
     plt.colorbar()
 
@@ -129,26 +126,26 @@ def SIE(x0,y0,n1,n2,b,beta,q,xc,theta):
 
     plt.figure(2)
     plt.title('relative difference')
-    diff = (real_source-sourcesl)
+    diff = (real_source-S)
     plt.imshow((np.abs(diff)), cmap = cm.gist_stern, interpolation = 'nearest')
     plt.axis('off')
     plt.colorbar()
     ####Lensed source
     plt.figure(3)
     plt.title('Original lensed galaxy')
-    plt.imshow(I2, cmap = cm.gist_stern, interpolation = 'nearest')
+    plt.imshow(HI2, cmap = cm.gist_stern, interpolation = 'nearest')
     plt.axis('off')
     plt.colorbar()
 
     plt.figure(4)
     plt.title('reconstructed lensed source')
-    plt.imshow((Imsl), vmin = np.min(I2), vmax = np.max(I2), cmap = cm.gist_stern, interpolation = 'nearest')
+    plt.imshow((FS), vmin = np.min(I2), vmax = np.max(I2), cmap = cm.gist_stern, interpolation = 'nearest')
     plt.axis('off')
     plt.colorbar()
 
     plt.figure(5)
     plt.title('error on the source in image plane')
-    plt.imshow((I2-Imsl), cmap = cm.gist_stern, interpolation = 'nearest')
+    plt.imshow((HI2-FS), cmap = cm.gist_stern, interpolation = 'nearest')
     plt.axis('off')
     plt.colorbar()
     ###Galaxy
@@ -178,13 +175,13 @@ def SIE(x0,y0,n1,n2,b,beta,q,xc,theta):
 
     plt.figure(9)
     plt.title('Reconstructed image')
-    plt.imshow(Imsl+G, cmap = cm.gist_stern, interpolation = 'nearest')
+    plt.imshow(FS+G, cmap = cm.gist_stern, interpolation = 'nearest')
     plt.axis('off')
     plt.colorbar()
 
     plt.figure(10)
     plt.title('difference with reconstructed image')
-    plt.imshow(Image-Imsl-G,cmap = cm.gist_stern, interpolation = 'nearest', vmin = -5*sigma, vmax = 5*sigma)#slit.fft_convolve(Im,PSF)
+    plt.imshow(Image-FS-G,cmap = cm.gist_stern, interpolation = 'nearest', vmin = -5*sigma, vmax = 5*sigma)#slit.fft_convolve(Im,PSF)
     plt.axis('off')
     plt.colorbar()
 

Examples/Test_SLIT_forUsers.py

@@ -0,0 +1,102 @@
+import pyfits as pf
+import matplotlib.pyplot as plt 
+import numpy as np
+import matplotlib.cm as cm
+from scipy import signal as scp
+import SLIT 
+import time
+from scipy import signal as scp
+import warnings
+warnings.simplefilter("ignore")
+
+#Example of a run of the SLIT algorithm on simulated images. 
+#Here the first part of the file shows how simulations are generated.
+#For users intereseted only in seeing the code run, have a look at the running SLIT section.
+#The command line that builds the Fkappa operator is also of outmost importance.
+
+Image = '''Input 2D image of the lens to invert '''
+nt1,nt2  = np.shape(Image)
+###############################Mass profile###############################
+x0,y0 = '''Input the center of mass of the lens with regard to coodinates in Image '''
+kappa = '''Input dimensionless pixelated mass density profile here '''
+size = '''Input the desired size of the output with regard to Image. Chosing 1 will result in a source with the same number of pixels as in Image. '''
+#Mapping between lens and source IMPORTANT
+Fkappa = SLIT.Lens.F(kappa, nt1,nt2, size,x0,y0)
+
+PSF = '''Input the PSF for Image here'''
+PSFconj = '''Input the conjugate of the PSF here given by
+                    np.real(np.fft.ifft2(np.conjugate(np.fft.fft2(PSF0[:-1,:-1])))), but be carefull that the result is still centered '''
+
+################################Running SLIT############################
+#Parameters
+kmax = 5
+niter =50
+
+#Start clock
+start = time.clock()
+
+#Running SLIT
+S, FS = SLIT.SLIT(Image, Fkappa, kmax, niter, size, PSF, PSFconj)
+
+#Stop clock
+elapsed = (time.clock()-start)
+print('execution time:', elapsed, 'seconds')
+
+#Reconstruction goodness
+real_source = newsource
+source_error = np.sum(np.abs(real_source[np.where(real_source!=0)]
+                                           -S[np.where(real_source!=0)])**2
+                                           /real_source[np.where(real_source!=0)]**2)/(np.size(
+                                           np.where(real_source!=0))/2.)
+image_chi2 = np.std(Image-FS)**2/sigma**2
+print('Residuals in source space', source_error)
+print('Residuals in image space',image_chi2)
+
+#Display of results
+for i in [1]:
+    plt.figure(2)
+ #   plt.suptitle('FISTA: error per pixel on the source: '+str(source_error)+' image chi2:'+str(image_chi2))
+ #   plt.subplot(2,3,1)
+    plt.title('Source from SLIT')
+    plt.imshow((S), vmin = np.min(real_source), vmax = np.max(real_source),cmap = cm.gist_stern, interpolation = 'nearest')
+    plt.axis('off')
+    plt.colorbar()
+#    plt.subplot(2,3,2)
+    plt.figure(3)
+    plt.title('Original source')
+    plt.imshow(real_source, cmap = cm.gist_stern, interpolation = 'nearest')
+    plt.axis('off')
+    plt.colorbar()
+ #   plt.subplot(2,3,3)
+    plt.figure(4)
+    plt.title('Lensed source')
+    plt.imshow(Image, cmap = cm.gist_stern, interpolation = 'nearest')
+    plt.axis('off')
+    plt.colorbar()
+    plt.figure(41)
+    plt.title('Reconstructed lensed source')
+    plt.imshow(FS, vmin = np.min(Image), vmax = np.max(Image), cmap = cm.gist_stern, interpolation = 'nearest')
+    plt.axis('off')
+    plt.colorbar()
+ #   plt.subplot(2,3,4)
+    plt.figure(5)
+    plt.title('relative difference')
+    diff = (real_source-S)/real_source
+    diff[np.where(real_source==0)] = 0
+    diff[np.where(diff>1)]= np.log(0.)
+    plt.imshow(np.abs(diff), vmax = 1., vmin = 0., cmap = cm.gist_stern, interpolation = 'nearest' )
+    plt.axis('off')
+    plt.colorbar()
+ #   plt.subplot(2,3,5)
+    plt.figure(6)
+    plt.title('difference with reconstructed lensed image')
+    plt.imshow(Image-FS, vmin = -5*sigma, vmax = 5*sigma, cmap = cm.gist_stern, interpolation = 'nearest')
+    plt.axis('off')
+    plt.colorbar()
+ #   plt.subplot(2,3,6)
+    plt.figure(7)
+    plt.title('difference with true source')
+    plt.imshow((np.abs(real_source-S)), cmap = cm.gist_stern, interpolation = 'nearest')
+    plt.axis('off')
+    plt.colorbar()
+plt.show()

README.md

@@ -1,2 +1,47 @@
-# SLIT
-SLIT is a python code that provides a method for invertion of lensed images in the frame of strong gravitational lensing
+When using this algorithm for publications, please acknowledge the work presented in:
+http://arxiv.org/abs/17soon
+
+========================================================
+INSTALLATION: 
+- Download and unzip the SLIT.zip file
+- In command line, open the SLIT folder with "cd SLIT"
+- Install the python package:
+"python setup.py install"
+
+========================================================
+TESTS:
+In SLIT, open the "Examples" folder and execute either of the examples:
+"python Test_SLIT.py" to test the reconstruction of a deblended lensed source.
+"python Test_SLIT_HR.py" to test the reconstruction of a deblended lensed source at higher resolution.
+"python Test_SLIT_MCA.py" to test the joint separation of lens and source light profiles and reconstruction of the source.
+"python Test_sparsity.py" to reproduce the non-linear approximation error plot from the paper to show the sparsity of lensed and unlensed galaxies.
+
+Inputs provided in the "Files" folder:
+"source.fits": the source used to generate simulations in Test_SLIT.py and Test_SLIT_MCA.py.
+"Source_HR.fits": the source used to generate simulations in Test_SLIT_HR.py.
+"Galaxy.fits": the lens galaxy light profile used in Test_SLIT_MCA.py.
+"kappa.fits": the lens mass density profile used in all simulations.
+"PSF.fits": the TinyTim PSF used in all simulations.
+"Noise_levels_SLIT_HR.fits": noise levels used in "Test_SLIT_HR" saved as a fits cube. They are usually computed in the algorithm, but we provide them as inputs to save time.
+"Noise_levels_SLIT_MCA.fits": noise levels used in "Test_SLIT_MCA" saved as a fits cube. They are usually computed in the algorithm, but we provide them as inputs to save time.
+"Noise_levels_SLIT.fits": noise levels used in "Test_SLIT_HR" saved as a fits cube. They are usually computed in the algorithm, but we provide them as inputs to save time.
+
+Results are displayed in pyplot windows.
+
+Users may change the input parameters via the same files and are encouraged to use them for further applications.
+
+========================================================
+HOW TO USE:
+We provide a Test_for_user.py file, that users may fill in with their own inputs and where the commands are a bit more detailed than in the other Test files.
+
+SLIT needs very little amount of inputs. It requires the input image along with lens mass profile and a PSF. The user then has to chose a maximum number of iterations after which the algorithm will stop if not converged and a image size ratio to the input image to set the resolution of the reconstructed source.
+
+Please do not hesitate emailing me for support via github, or at: [email protected].
+
+
+
+
+Contact GitHub API Training Shop Blog About 
+
+© 2017 GitHub, Inc. Terms Privacy Security Status Help 
+

SLIT/SLIT.py

@@ -67,6 +67,7 @@ def SLIT(img, Fkappa, kmax, niter, size, PSF,  PSFconj, levels = [0], mask = [0]
 
     #Noise simulations to estimate noise levels in source plane
     if np.sum(levels)==0:
+        print('Calculating noise levels')
         levels = simulate_noise(nt1,nt2, size, Fkappa, lensed, PSFconj)
         #Saves levels
         hdus = pf.PrimaryHDU(levels)
@@ -203,6 +204,7 @@ def SLIT_MCA(img, Fkappa, kmax, niter, riter, size,PSF, PSFconj, levels = [0], m
     levelg = level(nt1,nt1)
     #Noise simulations to estimate noise levels in source plane
     if np.sum(levels)==0:
+        print('Calculating noise levels')
         levels = simulate_noise(nt1,nt2, size, Fkappa, lensed, PSFconj)
         #Saves levels
         hdus = pf.PrimaryHDU(levels)
@@ -236,8 +238,6 @@ def star_inv(x):
 
 
     #Initialisations
-    RS = Lens.image_to_source(img, size, Fkappa, lensed = lensed)
-    RG = np.copy(img)
     if np.sum(Ginit)==0:
         G = np.random.randn(nt1,nt2)*sigma0
     else:
@@ -315,8 +315,7 @@ def star_inv(x):
     G = mw.iuwt(alphag)
     S[np.where(S<0)] = 0
     FS = Lens.source_to_image(S, nt1, nt2,Fkappa)
-    if np.sum(PSF)!=0:
-        FS = scp.fftconvolve(FS,PSF,mode = 'same')
+    FS = scp.fftconvolve(FS,PSF,mode = 'same')
 
     return S, FS,G
 

build/lib/SLIT/SLIT.py

@@ -67,6 +67,7 @@ def SLIT(img, Fkappa, kmax, niter, size, PSF,  PSFconj, levels = [0], mask = [0]
 
     #Noise simulations to estimate noise levels in source plane
     if np.sum(levels)==0:
+        print('Calculating noise levels')
         levels = simulate_noise(nt1,nt2, size, Fkappa, lensed, PSFconj)
         #Saves levels
         hdus = pf.PrimaryHDU(levels)
@@ -203,6 +204,7 @@ def SLIT_MCA(img, Fkappa, kmax, niter, riter, size,PSF, PSFconj, levels = [0], m
     levelg = level(nt1,nt1)
     #Noise simulations to estimate noise levels in source plane
     if np.sum(levels)==0:
+        print('Calculating noise levels')
         levels = simulate_noise(nt1,nt2, size, Fkappa, lensed, PSFconj)
         #Saves levels
         hdus = pf.PrimaryHDU(levels)
@@ -236,8 +238,6 @@ def star_inv(x):
 
 
     #Initialisations
-    RS = Lens.image_to_source(img, size, Fkappa, lensed = lensed)
-    RG = np.copy(img)
     if np.sum(Ginit)==0:
         G = np.random.randn(nt1,nt2)*sigma0
     else:
@@ -315,8 +315,7 @@ def star_inv(x):
     G = mw.iuwt(alphag)
     S[np.where(S<0)] = 0
     FS = Lens.source_to_image(S, nt1, nt2,Fkappa)
-    if np.sum(PSF)!=0:
-        FS = scp.fftconvolve(FS,PSF,mode = 'same')
+    FS = scp.fftconvolve(FS,PSF,mode = 'same')
 
     return S, FS,G