Commit 771ed8f3 authored by markus.meier01's avatar markus.meier01
Browse files

Uploaded P_cP and CDP_candes, attempted dim error fix in CDP_processor.

parent 943962eb
from numpy.random import randn, random_sample
from numpy.matlib import repmat
from numpy.linalg import norm
from numpy import sqrt, conj, reshape
from numpy.fft import fft2, ifft2
import matplotlib.pyplot as plt
import numpy as np
import random
import time
# Implementation of the Wirtinger Flow (WF) algorithm presented in the paper
# "Phase Retrieval via Wirtinger Flow: Theory and Algorithms"
# by E. J. Candes, X. Li, and M. Soltanolkotabi
# The input data are coded diffraction patterns about a random complex
# valued image.
def CDP_Candes(config):
## Make image
n1 = 128,
n2 = 256,
x = randn(n1,n2) + 1j*andn(n1,n2)
# Make masks and linear sampling operators
L = 10, # Number of masks
Masks = np.zeros(n1,n2,L); # Storage for L masks, each of dim n1 x n2
# Sample phases: each symbol in alphabet {1, -1, i , -i} has equal prob.
Randomarray = [1j, -1j, 1, -1]
for ll in range(L):
Masks[:,:,11] = np.random.choice(Randomarray,(n1,n2))
# Sample magnitudes and make masks
temp = random_sample(Masks.shape)
Masks = Masks*((temp <= 0.2)*sqrt(3) + (temp > 0.2)/sqrt(2))
# Make linear operators; A is forward map and
# At its scaled adjoint (At(Y)*numel(Y) is the adjoint)
A = lambda I: fft2(conj(Masks)*reshape(repmat(I, 1, L), np.size(I,1), np.size(I,2), L)) # Input is n1 x n2 image, output is n1 x n2 x L array
At = lambda Y: mean(Masks*ifft2,3) # Input is n1 x n2 X L array, output is n1 x n2 image
# Data
Y = abs(A(x))**2
## Initialization
npower_iter = 50 # Number of power iterations
z0 = z0/norm(z0,'fro') # Initial guess
z0 = randn(n1,n2)
start_time = time.time() # Power iterations
for tt in range (0, npower_iter):
z0 = At(Y*A(z0))
z0 = z0/norm(z0,'fro')
end_time = time.time()
Times = end_time - start_time
normest = sqrt(sum(Y.flatten())/Y.size) # Estimate norm to scale eigenvector
z = normest * z0 # Apply scaling
Relerrs = norm(x - exp(-1j*angle(trace(x*z)))*z, 'fro')/norm(x, 'fro')
T = 500 # Max number of iterations
tau0 = 330 # Time constant for step size
mu = lambda t: min(1-exp(-t(tau0), 0.4)) # Schedule for step size
start_timeTwo = time.time()
for t in range(T):
Bz = A(z)
C = (abs(Bz)**2 -Y) *Bz
grad = At(C) # Wirtinger gradient
z = z - mu(t)/normest**2 * grad # Gradient update
Relerrs = [Relerrs, norm(x - xp(-1j*angle(trace(x*z)))*z, 'fro')/norm(x, 'fro')]
end_timeTwo = time.time()
Times = end_timeTwo - start_timeTwo
print('All done!')
## Check results
print('Run time of initialization: %.2f seconds', Times(1))
print('Relative error after initialization: %f', Relerrs(1))
print('Loop run time: %.2f', Times(T+1))
print('Relative error after %d iterations: #f', T, Relerrs(T+1))
X = range(T)
Y = Relerrs
fig = plt.figure()
ax.plot(X, Y)
plt.title('Relative error vs. iteration count')
plt.xlabel('Iteration')
plt.ylabel('Relative error (log10)')
ax.set_yscale('log')
plt.show()
from numpy.random import randn, random_sample
from numpy.linalg import norm
from numpy import sqrt, conj, tile, mean, exp, angle, trace
from numpy.fft import fft, ifft
from numpy import sqrt, conj, tile, mean, exp, angle, trace, zeros, reshape
from numpy.fft import fft, ifft, fft2
import numpy as np
import time
......@@ -61,7 +61,11 @@ def CDP_processor(config):
config['rt_data']=Y
Y=Y**2
config['data']=Y
config['norm_data']=sum(sum(Y))/Y.size
### next line newly added ###
config['norm_data']=(sum(Y.flatten()))/Y.size
### config['norm_data']=sum(sum(Y))/Y.size ###
normest = sqrt(config['norm_data']) # Estimate norm to scale eigenvector
config['norm_rt_data']=normest
......@@ -76,16 +80,18 @@ def CDP_processor(config):
z0 = z0/norm(z0)
toc = time.time()
z = normest * z0 # Apply scaling
if n2==1:
Relerrs = norm(x - exp(-1j*angle(trace(x.T*z))) * z, 'fro')/norm(x,'fro')
Relerrs = norm(x - exp(-1j*angle(trace(x*z))) * z, 'fro')/norm(x,'fro')
config['u_0'] = tile(z,[1,L])
elif n1==1:
Relerrs = norm(x - exp(-1j*angle(trace(z.T*x))) * z, 'fro')/norm(x,'fro')
Relerrs = norm(x - exp(-1j*angle(trace(z*x))) * z, 'fro')/norm(x,'fro')
config['u_0'] = tilet(z,[L,1])
else:
Relerrs = norm(x - exp(-1j*angle(trace(x.T*z))) * z, 'fro')/norm(x,'fro')
Relerrs = norm(x - exp(-1j*angle(trace(x*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 seconds', toc-tic)
......
from .proxoperators import ProxOperator, magproj
import numpy as np
import math
class P_cP(ProxOperator):
"""
# Projection operator onto a magnitude constraint
"""
def __init__(self,config):
self.support_idx = config['support_idx']
self.abs_illumination = config['abs_illumination']
self.illumination = config['illumination']
def work(self,u):
unew = self.illumination*np.ones((u,u)) ###purpose of this line?###
###uabs = math.sqrt(np.conjugate(u)*u) not needed?###
unew = magproj(self.abs_illumination[self.support_idx]*math.exp(-0.0045*np.angle(u[self.support_idx]/self.illumination[self.support_idx])),u[self.support_idx])
return unew
\ No newline at end of file
......@@ -22,6 +22,7 @@ from .Prox_nonneg_l0 import *
from .Prox_nonneg_l1 import *
from .P_sequential_hyperplane import *
from .P_CDP import *
from .P_cP import *
__all__ = ["P_diag",
......
Markdown is supported
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