Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

dev-aym branch #1

Open
wants to merge 17 commits into
base: master
Choose a base branch
from
Open
Changes from 1 commit
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
Add support of PySAP starlet (1st gen) transform
  • Loading branch information
aymgal committed Mar 16, 2019
commit ba22bbe7d08c2d492c43112a3c4f7d508281e073
16 changes: 8 additions & 8 deletions SLIT/Lens.py
Original file line number Diff line number Diff line change
@@ -151,8 +151,8 @@ def F(kappa, nt1,nt2, size, extra=100, x_shear = 0, y_shear = 0, alpha_x_in = [-
na1,na2 = np.shape(alpha_x)
xa,ya = np.where(np.zeros((na1,na2)) == 0)

nb1=nt1*size
nb2=nt2*size
nb1 = int(nt1*size)
nb2 = int(nt2*size)
xb, yb = np.where(np.zeros((nb1,nb2))==0)

#Scaling of the source grid
@@ -224,10 +224,10 @@ def image_to_source(Image, size,beta,lensed = 0, square = 0):
# nsize1,nsize2: size of the postagestamp in source plane
# F: lens mapping matrix

F = (beta)
F = beta
nt1,nt2 = np.shape(Image)
nb1 = (nt1*size)
nb2 = (nt2*size)
nb1 = int(nt1*size)
nb2 = int(nt2*size)

Source = np.zeros((nb1,nb2))
xb,yb = np.where(Source == 0)
@@ -255,10 +255,10 @@ def image_to_source_bound(Image, size,beta,lensed = 0):
# nsize1,nsize2: size of the postagestamp in source plane
# F: lens mapping matrix

F = (beta)
F = beta
nt1,nt2 = np.shape(Image)
nb1 = nt1*size
nb2 = nt2*size
nb1 = int(nt1*size)
nb2 = int(nt2*size)
Source = np.zeros((nb1,nb2))
xb,yb = np.where(Source == 0)
N = np.size(xb)
33 changes: 18 additions & 15 deletions SLIT/Solve.py
Original file line number Diff line number Diff line change
@@ -51,7 +51,7 @@ def SLIT(Y, Fkappa, kmax, niter, size, PSF, PSFconj, S0 = [0], levels = [0], sch
n1,n2 = np.shape(Y)
PSFconj = PSF.T
#Size of the source
ns1,ns2 = n1*size, n2*size
ns1,ns2 = int(n1*size), int(n2*size)
#Number of starlet scales in source plane
if lvl ==0:
lvl = np.int(np.log2(ns2))
@@ -97,9 +97,10 @@ def PSF_apply(i):
def PSFT_apply(ii):
return scp.fftconvolve(ii,PSFconj,mode = 'same')
def transform(x):
return tools.wave_transform(x, lvl, newwave = 1)
coeffs, _ = tools.wave_transform(x, lvl, newwave=1)
return coeffs
def inverse(x):
return tools.iuwt(x)
return tools.iuwt(x, newwave=1)

#Forward operator
def F_op(X):
@@ -324,8 +325,8 @@ def SLIT_MCA(Y, Fkappa, kmax, niter, riter, size,PSF, PSFconj, lvlg = 0, lvls =
#Shape of the image
n1,n2 = np.shape(Y)
#Initialisation of the source
ns1= n1*size
ns2 = n2*size
ns1 = int(n1*size)
ns2 = int(n2*size)
PSFconj = PSF.T
#Number of starlet scales in source and image planes
if lvlg ==0:
@@ -375,9 +376,10 @@ def PSF_apply(i):
def PSFT_apply(ii):
return scp.fftconvolve(ii,PSFconj,mode = 'same')
def transform(x):
return tools.wave_transform(x, lvlg)
coeffs, _ = tools.wave_transform(x, lvlg, newwave=1)
return coeffs
def inverse(x):
return tools.iuwt(x)
return tools.iuwt(x, newwave=1)

def FWS_op(X):
return PSF_apply(F_apply(inverse(X)))
@@ -743,10 +745,11 @@ def PSFT_apply(ii):
return scp.fftconvolve(ii, PSFconj, mode='same')

def transform(x):
return tools.wave_transform(x, lvlg)
coeffs, _ = tools.wave_transform(x, lvlg, newwave=1)
return coeffs

def inverse(x):
return tools.iuwt(x)
return tools.iuwt(x, newwave=1)

def FWS_op(X):
return Down(PSF_apply(F_apply(inverse(X))))
@@ -1039,7 +1042,7 @@ def plot_cube(cube):


def level_source_HR(n1,n2, size, sigma,PSFT, Lens_op2, Up, lvl):
ns1,ns2 = n1*size, n2*size
ns1,ns2 = int(n1*size), int(n2*size)
ones = np.ones((n1,n2))
noise = Up(ones*sigma)*(size)**2
Hnoise = np.sqrt(scp.fftconvolve(noise**2, PSFT**2, mode = 'same'))#noise*np.sqrt(np.sum(PSFT**2))##
@@ -1050,7 +1053,7 @@ def level_source_HR(n1,n2, size, sigma,PSFT, Lens_op2, Up, lvl):
dirac = np.zeros((ns1,ns2))
dirac[int(ns1/2),int(ns2/2)] = 1
print(dirac.shape, ns1,ns2)
wave_dirac = tools.wave_transform(dirac, lvl)
wave_dirac, _ = tools.wave_transform(dirac, lvl)
levels = np.zeros(wave_dirac.shape)
for i in range(lvl):
if np.size(noise.shape) > 2:
@@ -1066,7 +1069,7 @@ def level_source_HR(n1,n2, size, sigma,PSFT, Lens_op2, Up, lvl):


def level_source(n1,n2,sigma,size,PSFT, Lens_op2, lensed, lvl):
ns1,ns2 = n1*size, n2*size
ns1,ns2 = int(n1*size), int(n2*size)
ones = np.ones((n1,n2))
lensed[lensed == 0] = 1
noise = ones*sigma
@@ -1081,7 +1084,7 @@ def level_source(n1,n2,sigma,size,PSFT, Lens_op2, lensed, lvl):

dirac[int(ns1/2),int(ns2/2)] = 1

wave_dirac = tools.wave_transform(dirac, lvl)
wave_dirac, _ = tools.wave_transform(dirac, lvl)
levels = np.zeros(wave_dirac.shape)
for i in range(lvl):
if np.size(noise.shape) > 2:
@@ -1170,7 +1173,7 @@ def mk_bound(Fkappa, n1,n2,size):


def mk_simu(n1,n2,lvl,size, sigma, I_op, transform, n):
storage = np.zeros((lvl,n1*size, n2*size, n))
storage = np.zeros((lvl, int(n1*size), int(n2*size), n))
for i in range(n):
noise = np.random.randn(n1,n2)*sigma
noise_lens = I_op(noise)
@@ -1201,7 +1204,7 @@ def simulate_noise(n1,n2, sigma, size, I_op, transform, lvl, Npar = np.int(mtp.c
n = 500
if Npar>mtp.cpu_count():
Npar = mtp.cpu_count()
ns1,ns2 = n1*size, n2*size
ns1,ns2 = int(n1*size), int(n2*size)
# lvl = np.int(np.log2(ns1))
w_levels = np.zeros((lvl,ns1,ns2))

File renamed without changes.
Loading