diff --git a/Examples/Test_sparsity.py b/Examples/Test_sparsity.py index 7284ff6..77d0514 100644 --- a/Examples/Test_sparsity.py +++ b/Examples/Test_sparsity.py @@ -27,17 +27,20 @@ Fkappa = Lens.F(kappa, nt1,nt2, size,nt1/2.,nt2/2.) lensed = slit.lens_one(Fkappa, nt1,nt2, size) +#Levels for normalisation +lev = slit.level(nt1,nt1) + #Starlet transforms of the lens and source in their respective planes -wG = mw.wave_transform(G, lvl = 6) -wS = mw.wave_transform(S, lvl = 6) +wG = mw.wave_transform(G, lvl = 6)/lev +wS = mw.wave_transform(S, lvl = 6)/lev #Lensed source FS = Lens.source_to_image(S, nt1, nt2,Fkappa) #Unlensed lens FG = Lens.image_to_source(G, size, Fkappa, lensed=lensed) #Starlet transform of the unlensed lens -wFG = mw.wave_transform(FG, 6) +wFG = mw.wave_transform(FG, 6)/lev #Starlet transform of the lensed -wFS = mw.wave_transform(FS, 6) +wFS = mw.wave_transform(FS, 6)/lev def mk_sort(X): Y = np.sort(np.resize(np.abs(X), X.size)) @@ -72,27 +75,35 @@ def error_rec_from(X, p, wave = 0): error_wG[i] = error_rec_from(wG, i, wave = 1) error_wFG[i] = error_rec_from(wFG, i, wave = 1) - +print('NLA on the source at 10%: ',error_wS[100]/np.max(error_wS)) +print('NLA on the lens at 10%: ', error_wG[100]/np.max(error_wG)) +print('NLA on the lensed source at 10%: ', error_wFS[100]/np.max(error_wFS)) +print('NLA on the delensed lens at 10%: ', error_wFG[100]/np.max(error_wFG)) #Display -plt.plot(np.linspace(0,100, 1000), error_S/np.max(error_S), '--r', label = 'Source in direct space', linewidth = 3) -plt.plot(np.linspace(0,100, 1000), error_G/np.max(error_G), '--b', label = 'Lens in direct space', linewidth = 3) +plt.figure(1) plt.plot(np.linspace(0,100, 1000), error_wS/np.max(error_wS), 'r', label = 'Source in starlet space', linewidth = 3) +plt.plot(np.linspace(0,100, 1000), error_wFG/np.max(error_wFG), 'c', label = 'Lens in source plane in starlet space', linewidth = 3) + +plt.xlabel('percentage of coefficients used in reconstruction', fontsize=25) +plt.ylabel('Error on reconstruction', fontsize=25) +plt.title('Non-linear approximation error in source plane', fontsize=25) +plt.legend(fontsize = 25) +a = plt.axes([0.4, 0.2, 0.45, 0.4]) +plt.semilogy(np.linspace(0,100, 1000), (error_wFG/np.max(error_wFG)), 'c', linewidth = 3) +plt.semilogy(np.linspace(0,100, 1000), error_wS/np.max(error_wS), 'r', linewidth = 3) +plt.xlim(20,100) + +plt.figure(2) plt.plot(np.linspace(0,100, 1000), error_wG/np.max(error_wG), 'b', label = 'Galaxy in starlet space', linewidth = 3) plt.plot(np.linspace(0,100, 1000), error_wFS/np.max(error_wFS), 'm', label = 'Lensed source in starlet space', linewidth = 3) -plt.plot(np.linspace(0,100, 1000), error_wFG/np.max(error_wFG), 'c', label = 'Unlensed lens in starlet space', linewidth = 3) + plt.xlabel('percentage of coefficients used in reconstruction', fontsize=25) plt.ylabel('Error on reconstruction', fontsize=25) -plt.title('Non-linear approximation error', fontsize=25) +plt.title('Non-linear approximation error in lens plane', fontsize=25) plt.legend(fontsize = 25) - a = plt.axes([0.4, 0.2, 0.45, 0.4]) - -plt.semilogy(np.linspace(0,100, 1000), (error_wS/np.max(error_wS)), 'r', linewidth = 3) -plt.semilogy(np.linspace(0,100, 1000), (error_wG/np.max(error_wG)), 'b', linewidth = 3) plt.semilogy(np.linspace(0,100, 1000), (error_wFS/np.max(error_wFS)), 'm', linewidth = 3) -plt.semilogy(np.linspace(0,100, 1000), (error_wFG/np.max(error_wFG)), 'c', linewidth = 3) -plt.semilogy(np.linspace(0,100, 1000), error_S/np.max(error_S), '--r', linewidth = 3) -plt.semilogy(np.linspace(0,100, 1000), error_G/np.max(error_G), '--b', linewidth = 3) +plt.semilogy(np.linspace(0,100, 1000), error_wG/np.max(error_wG), 'b', linewidth = 3) plt.xlim(20,100) plt.show() diff --git a/SLIT/Lens.py b/SLIT/Lens.py index 74bcfae..07640f5 100644 --- a/SLIT/Lens.py +++ b/SLIT/Lens.py @@ -96,6 +96,7 @@ def F(kappa, nt1,nt2, size, x0, y0, extra=100,local = False, alpha_file = 'none' F = [] + F2 = [] rec = [] for i in range(np.size(xb)): #Deflection of photons emitted in xb[i],yb[i] @@ -105,23 +106,31 @@ def F(kappa, nt1,nt2, size, x0, y0, extra=100,local = False, alpha_file = 'none' xprox = np.int_(np.abs((xa2d-theta_x)*2)) yprox = np.int_(np.abs((ya2d-theta_y)*2)) - if np.min(xprox+yprox) <3: - loc = np.where((xprox+yprox)==np.min(xprox+yprox))# - else: - loc = [] + loc = np.where((xprox+yprox)==0)# + if (np.size(loc)==0): F.append([0]) else: F.append(np.int_(np.array(loc))) - return F + + if np.min(xprox+yprox) <3: + loc2 = np.where((xprox+yprox)==np.min(xprox+yprox))# + else: + loc2 = [] + if (np.size(loc2)==0): + + F2.append([0]) + else: + F2.append(np.int_(np.array(loc2))) + return F,F2 def source_to_image(Source, nt1,nt2, theta, ones = 1): # Source: Image of the source in the source plane # n1,n2: size in pixels of the postage stamp in image plane # F: the coordinates of the lens mapping - F = (theta) + F = (theta[0]) nb1,nb2 = np.shape(Source) if ones == 1: @@ -152,7 +161,7 @@ def image_to_source(Image, size,beta,lensed = 0): # nsize1,nsize2: size of the postagestamp in source plane # F: lens mapping matrix - F = (beta) + F = (beta[0]) nt1,nt2 = np.shape(Image) nb1 = nt1*size nb2 = nt2*size @@ -171,10 +180,27 @@ def image_to_source(Image, size,beta,lensed = 0): return Source +def image_to_source_bound(Image, size,beta,lensed = 0): + # Image: postagestamp of the observed image + # nsize1,nsize2: size of the postagestamp in source plane + # F: lens mapping matrix - - - - - + F = (beta[1]) + nt1,nt2 = np.shape(Image) + nb1 = nt1*size + nb2 = nt2*size + Source = np.zeros((nb1,nb2)) + xb,yb = np.where(Source == 0) + N = np.size(xb) + for k in range(N): + pos = F[k] + if np.size(np.shape(pos)) > 1: + if np.sum(lensed) !=0: + Source[xb[k],yb[k]] += np.sum(Image[np.array(pos[0][:]), + np.array(pos[1][:])])/np.max([1,np.size(pos[0][:])]) + else: + Source[xb[k],yb[k]] += np.sum(Image[np.array(pos[0][:]), + np.array(pos[1][:])]) + + return Source diff --git a/SLIT/SLIT.py b/SLIT/SLIT.py index 74c7b55..04afeea 100644 --- a/SLIT/SLIT.py +++ b/SLIT/SLIT.py @@ -44,7 +44,7 @@ def SLIT(img, Fkappa, kmax, niter, size, PSF, PSFconj, levels = [0], mask = [0] ##EXAMPLE: ## S,FS = SLIT(img, Fkappa, 5, 100, 1, PSF, PSFconj) - + PSFconj = PSF.T nt1,nt2 = np.shape(img) #Size of the source @@ -181,6 +181,7 @@ def SLIT_MCA(img, Fkappa, kmax, niter, riter, size,PSF, PSFconj, levels = [0], m ##EXAMPLE: ## S,FS = SLIT(img, Fkappa, 5, 100, 1, PSF, PSFconj) + PSFconj = PSF.T #Shape of the image nt1,nt2 = np.shape(img) #Initialisation of the source @@ -198,8 +199,7 @@ def SLIT_MCA(img, Fkappa, kmax, niter, riter, size,PSF, PSFconj, levels = [0], m #Mapping of an all-at-one image lensed = lens_one(Fkappa, nt1,nt2, size) #Limits of the image plane in source plane - bound = lensed/lensed - bound[lensed == 0] = 0 + bound = mk_bound(Fkappa, nt1,nt2, size) #Noise levels in image plane in starlet space levelg = level(nt1,nt1) #Noise simulations to estimate noise levels in source plane @@ -420,6 +420,23 @@ def lens_one(Fkappa, nt1,nt2,size): lensed = Lens.image_to_source(dirac, size,Fkappa,lensed = [0]) return lensed +def mk_bound(Fkappa, nt1,nt2,size): + ##DESCRIPTION: + ## Function that returns the support of the lens image in source plane. + ## + ##INPUTS: + ## -Fkappa: the mapping between source and image planes + ## -nt1,nt2: the shape of the image. + ## -size: the factor that scales the shape of the source relative to the shape of the image + ## + ##OUTPUTS: + ## -lensed: the projection to source plane of an all at aone image. + dirac = np.ones((nt1,nt2)) + lensed = Lens.image_to_source_bound(dirac, size,Fkappa,lensed = [0]) + bound = lensed/lensed + bound[lensed==0]=0 + return bound + def MAD(x,n=3): ##DESCRIPTION: @@ -455,7 +472,7 @@ def ST(alpha, k, levels, sigma): ## -alpha: The thresholded coefficients. lvl, n1,n2 = np.shape(alpha) th = np.ones((lvl,n1,n2))*k - th[0,:,:] = th[0,:,:]+1 + th[0,:,:] = th[0,:,:]+3 th[-1,:,:] = 0 @@ -490,14 +507,8 @@ def simulate_noise(nt1,nt2, size, Fkappa, lensed, PSFconj): nb1,nb2 = nt1*size, nt2*size lvl = np.int(np.log2(nb1)) w_levels = np.zeros((lvl,nb1,nb2)) - - dirac = np.zeros((nb1,nb2)) - dirac[nb1/2.,nb2/2.] = 1. - wave_dirac = mw.wave_transform(dirac, lvl, newwave = 1) - lens = 0 lvl = np.int_(np.log2(nb1)) - nn = 2**(lvl+2) storage = np.zeros((lvl, nb1,nb2,n)) for i in range(np.int(n)): noise = np.random.randn(nt1,nt2) diff --git a/build/lib/SLIT/Lens.py b/build/lib/SLIT/Lens.py index 74bcfae..07640f5 100644 --- a/build/lib/SLIT/Lens.py +++ b/build/lib/SLIT/Lens.py @@ -96,6 +96,7 @@ def F(kappa, nt1,nt2, size, x0, y0, extra=100,local = False, alpha_file = 'none' F = [] + F2 = [] rec = [] for i in range(np.size(xb)): #Deflection of photons emitted in xb[i],yb[i] @@ -105,23 +106,31 @@ def F(kappa, nt1,nt2, size, x0, y0, extra=100,local = False, alpha_file = 'none' xprox = np.int_(np.abs((xa2d-theta_x)*2)) yprox = np.int_(np.abs((ya2d-theta_y)*2)) - if np.min(xprox+yprox) <3: - loc = np.where((xprox+yprox)==np.min(xprox+yprox))# - else: - loc = [] + loc = np.where((xprox+yprox)==0)# + if (np.size(loc)==0): F.append([0]) else: F.append(np.int_(np.array(loc))) - return F + + if np.min(xprox+yprox) <3: + loc2 = np.where((xprox+yprox)==np.min(xprox+yprox))# + else: + loc2 = [] + if (np.size(loc2)==0): + + F2.append([0]) + else: + F2.append(np.int_(np.array(loc2))) + return F,F2 def source_to_image(Source, nt1,nt2, theta, ones = 1): # Source: Image of the source in the source plane # n1,n2: size in pixels of the postage stamp in image plane # F: the coordinates of the lens mapping - F = (theta) + F = (theta[0]) nb1,nb2 = np.shape(Source) if ones == 1: @@ -152,7 +161,7 @@ def image_to_source(Image, size,beta,lensed = 0): # nsize1,nsize2: size of the postagestamp in source plane # F: lens mapping matrix - F = (beta) + F = (beta[0]) nt1,nt2 = np.shape(Image) nb1 = nt1*size nb2 = nt2*size @@ -171,10 +180,27 @@ def image_to_source(Image, size,beta,lensed = 0): return Source +def image_to_source_bound(Image, size,beta,lensed = 0): + # Image: postagestamp of the observed image + # nsize1,nsize2: size of the postagestamp in source plane + # F: lens mapping matrix - - - - - + F = (beta[1]) + nt1,nt2 = np.shape(Image) + nb1 = nt1*size + nb2 = nt2*size + Source = np.zeros((nb1,nb2)) + xb,yb = np.where(Source == 0) + N = np.size(xb) + for k in range(N): + pos = F[k] + if np.size(np.shape(pos)) > 1: + if np.sum(lensed) !=0: + Source[xb[k],yb[k]] += np.sum(Image[np.array(pos[0][:]), + np.array(pos[1][:])])/np.max([1,np.size(pos[0][:])]) + else: + Source[xb[k],yb[k]] += np.sum(Image[np.array(pos[0][:]), + np.array(pos[1][:])]) + + return Source diff --git a/build/lib/SLIT/SLIT.py b/build/lib/SLIT/SLIT.py index 74c7b55..2a0f931 100644 --- a/build/lib/SLIT/SLIT.py +++ b/build/lib/SLIT/SLIT.py @@ -44,7 +44,7 @@ def SLIT(img, Fkappa, kmax, niter, size, PSF, PSFconj, levels = [0], mask = [0] ##EXAMPLE: ## S,FS = SLIT(img, Fkappa, 5, 100, 1, PSF, PSFconj) - + PSFconj = PSF.T nt1,nt2 = np.shape(img) #Size of the source @@ -181,6 +181,7 @@ def SLIT_MCA(img, Fkappa, kmax, niter, riter, size,PSF, PSFconj, levels = [0], m ##EXAMPLE: ## S,FS = SLIT(img, Fkappa, 5, 100, 1, PSF, PSFconj) + PSFconj = PSF.T #Shape of the image nt1,nt2 = np.shape(img) #Initialisation of the source @@ -198,8 +199,7 @@ def SLIT_MCA(img, Fkappa, kmax, niter, riter, size,PSF, PSFconj, levels = [0], m #Mapping of an all-at-one image lensed = lens_one(Fkappa, nt1,nt2, size) #Limits of the image plane in source plane - bound = lensed/lensed - bound[lensed == 0] = 0 + bound = mk_bound(Fkappa, nt1,nt2, size) #Noise levels in image plane in starlet space levelg = level(nt1,nt1) #Noise simulations to estimate noise levels in source plane @@ -420,6 +420,23 @@ def lens_one(Fkappa, nt1,nt2,size): lensed = Lens.image_to_source(dirac, size,Fkappa,lensed = [0]) return lensed +def mk_bound(Fkappa, nt1,nt2,size): + ##DESCRIPTION: + ## Function that returns the support of the lens image in source plane. + ## + ##INPUTS: + ## -Fkappa: the mapping between source and image planes + ## -nt1,nt2: the shape of the image. + ## -size: the factor that scales the shape of the source relative to the shape of the image + ## + ##OUTPUTS: + ## -lensed: the projection to source plane of an all at aone image. + dirac = np.ones((nt1,nt2)) + lensed = Lens.image_to_source_bound(dirac, size,Fkappa,lensed = [0]) + bound = lensed/lensed + bound[lensed==0]=0 + return bound + def MAD(x,n=3): ##DESCRIPTION: @@ -490,14 +507,8 @@ def simulate_noise(nt1,nt2, size, Fkappa, lensed, PSFconj): nb1,nb2 = nt1*size, nt2*size lvl = np.int(np.log2(nb1)) w_levels = np.zeros((lvl,nb1,nb2)) - - dirac = np.zeros((nb1,nb2)) - dirac[nb1/2.,nb2/2.] = 1. - wave_dirac = mw.wave_transform(dirac, lvl, newwave = 1) - lens = 0 lvl = np.int_(np.log2(nb1)) - nn = 2**(lvl+2) storage = np.zeros((lvl, nb1,nb2,n)) for i in range(np.int(n)): noise = np.random.randn(nt1,nt2) diff --git a/dist/SLIT-0.1-py2.7.egg b/dist/SLIT-0.1-py2.7.egg index 1621ea2..a713e85 100644 Binary files a/dist/SLIT-0.1-py2.7.egg and b/dist/SLIT-0.1-py2.7.egg differ