Commit 813e667f authored by alexander.dornheim's avatar alexander.dornheim
Browse files

Finished first version of CDP_processor. Still needs testing.

parent 7d39028b
......@@ -39,48 +39,48 @@ def CDP_processor(config):
if n2==1:
# Make linear operators; A is forward map and At its scaled adjoint (At(Y)*numel(Y) is the adjoint)
A = @(I) fft(conj(Masks) .* repmat(I,[1 L])); # Input is n x 1 signal, output is n x L array
At = @(Y) mean(Masks .* ifft(Y), 2); # Input is n x L array, output is n x 1 signal
elseif(n1==1)
A = lambda I: fft(conj(Masks) * tile(I,[1, L])) # Input is n x 1 signal, output is n x L array
At = lambda Y: mean(Masks * ifft(Y), 1) # Input is n x L array, output is n x 1 signal
elif n1==1 :
# Make linear operators; A is forward map and At its scaled adjoint (At(Y)*numel(Y) is the adjoint)
A = @(I) fft(conj(Masks) .* repmat(I,[L 1])); # Input is 1 x n signal, output is L x n array
At = @(Y) mean(Masks .* ifft(Y), 1); # Input is L x n array, output is 1 x n signal
else
A = @(I) fft2(config['Masks .* reshape(repmat(I,[1 L]), size(I,1), size(I,2), L)); # Input is n1 x n2 image, output is n1 x n2 x L array
At = @(Y) mean(Masks .* ifft2(Y), 3); # Input is n1 x n2 X L array, output is n1 x n2 image
A = lambda I: fft(conj(Masks) * tile(I,[L, 1])) # Input is 1 x n signal, output is L x n array
At = lambda Y: mean(Masks * ifft(Y), 0) # Input is L x n array, output is 1 x n signal
else:
A = lambda I: fft2(config['Masks'] * reshape(tile(I,[1 L]), (I.shape[0]I.shape[1], L)) # Input is n1 x n2 image, output is n1 x n2 x L array
At = lambda Y: mean(Masks * ifft2(Y), 2) # Input is n1 x n2 X L array, output is n1 x n2 image
end
# Data
Y = abs(A(x))
config['rt_data']=Y;
config['rt_data']=Y
Y=Y**2
config['data']=Y
config['norm_data']=sum(Y(:))/numel(Y)
config['norm_data']=sum(Y)/Y.size
normest = sqrt(config['norm_data']) # Estimate norm to scale eigenvector
config['norm_rt_data']=normest
## Initialization
npower_iter = config['warmup_iter']; # Number of power iterations
z0 = randn(n1,n2); z0 = z0/norm(z0,'fro'); # Initial guess
tic # Power iterations
npower_iter = config['warmup_iter']; # Number of power iterations
z0 = randn((n1,n2)); z0 = z0/norm(z0,'fro') # Initial guess
tic = time.time() # Power iterations
for tt in range(npower_iter):
z0 = At(Y.*A(z0)); z0 = z0/norm(z0,'fro');
z0 = At(Y*A(z0)); z0 = z0/norm(z0,'fro')
Times = toc
toc = time.time()
z = normest * z0; # Apply scaling
z = normest * z0 # Apply scaling
if n2==1:
Relerrs = norm(x - exp(-1i*angle(trace(x'*z))) * z, 'fro')/norm(x,'fro');
config['u_0'] = repmat(z,1,L);
Relerrs = norm(x - exp(-1j*angle(trace(x.T*z))) * z, 'fro')/norm(x,'fro')
config['u_0'] = tile(z,[1,L])
elif n1==1:
Relerrs = norm(x - exp(-1i*angle(trace(z'*x))) * z, 'fro')/norm(x,'fro');
config['u_0'] = repmat(z,L,1);
Relerrs = norm(x - exp(-1j*angle(trace(z.T*x))) * z, 'fro')/norm(x,'fro')
config['u_0'] = tilet(z,[L,1])
else:
Relerrs = norm(x - exp(-1i*angle(trace(x'*z))) * z, 'fro')/norm(x,'fro');
config['u_0']=reshape(repmat(z,[1 L]), size(z,1), size(z,2), L);
Relerrs = norm(x - exp(-1j*angle(trace(x.T*z))) * z, 'fro')/norm(x,'fro')
config['u_0']=reshape(tile(z,[1, L]), (z.shape[0], z.shape[1], L))
print('Run time of initialization: #.2f\n', Times)
print('Relative error after initialization: #f\n', Relerrs)
print('\n')
print('Run time of initialization: %.2f seconds', toc-tic)
print('Relative error after initialization: %.2f', Relerrs)
print('\n')
\ No newline at end of file
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment