Uploaded CDP_source_loc_processor, fixed bugs in CDP_candes and CDP_processor

proxtoolbox/Problems/Phase/CDP_candes.py

 from numpy.random import randn, random_sample
 from numpy.linalg import norm
 from numpy.fft import fft2, ifft2
-from numpy import sqrt, conj, reshape, mean
+from numpy import sqrt, conj, reshape, mean, tile
 import matplotlib.pyplot as plt
 import numpy as np
 import random
@@ -19,8 +19,8 @@ import time
 def CDP_Candes(config):
 
 ## Make image
-    n1 = 128,
-    n2 = 256,
+    n1 = 128
+    n2 = 256
     x = randn(n1,n2) + 1j*randn(n1,n2)
     
     # Make masks and linear sampling operators
@@ -53,7 +53,7 @@ def CDP_Candes(config):
 
     start_time = time.time()                         # Power iterations
 
-    for tt in range (0, npower_iter):
+    for tt in range (npower_iter):
         z0 = At(Y*A(z0))
         z0 = z0/norm(z0,'fro')
 
@@ -75,7 +75,7 @@ def CDP_Candes(config):
         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')]
+        Relerrs = [Relerrs, norm(x -exp(-1j*angle(trace(x*z)))*z, 'fro')/norm(x, 'fro')]
     end_timeTwo = time.time()
     Times = end_timeTwo - start_timeTwo
 

proxtoolbox/Problems/Phase/CDP_processor.py

 from numpy.random import randn, random_sample
 from numpy.linalg import norm
-from numpy import sqrt, conj, tile, mean, exp, angle, trace, zeros, reshape
-from numpy.fft import fft, ifft, fft2
+from numpy import sqrt, conj, tile, mean, exp, angle, trace, reshape
+from numpy.fft import fft, ifft, fft2, ifft2
 import numpy as np
 import time
 

proxtoolbox/Problems/Phase/CDP_source_loc_processor.py

+from numpy.random import randn, random_sample
+from numpy.linalg import norm
+from numpy.fft import fft, ifft, fft2, ifft2
+from numpy import sqrt, conj, reshape, tile, mean, angle, trace
+import numpy as np
+import random
+
+
+def CDP_source_loc_processor(config):
+# 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.
+
+    ## Make image
+    n1 = config['Ny']
+    n2 = config['Nx']     # for 1D signals, this will be 1
+    x = randn(n1,n2) + 1j*randn(n1,n2)
+    
+    ## Make masks and linear sampling operators
+    
+    L = config['sets']       # Number of masks 
+    Randomarray = [1j,-1j,1,-1]
+    
+    if n2==1:
+        Masks = np.random.choice(Randomarray,(n1,L))
+    elif n1==1:
+        Masks = np.random.choice(Randomarray,(L,n2))
+    else:
+        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.
+        for ll in range (L):
+            Masks[:,:,ll] = 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))
+    config['Masks'] = conj(Masks)
+    # Saving the conjugate of the mask saves on computing the conjugate
+    # every time the mapping A (below) is applied.
+
+    
+    
+    if n2==1:
+        # Make linear operators; A is forward map and At its scaled adjoint (At(Y)*numel(Y) is the adjoint)
+        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), axis=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 = 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), axis=1)     # 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), axis=2)                                              # Input is n1 x n2 X L array, output is n1 x n2 image
+    
+    # support constraint: none
+    config['indicator_ampl'] = 1
+    config['abs_illumination'] = 1
+    config ['Illumination'] = 1
+    ##suport_idx is missing###
+    
+    # Data
+    Y=abs(A(x))
+    config['rt_data'] = Y
+    Y = Y**2
+    config['data'] = Y
+    config['norm_data'] = (sum(Y.flatten()))/Y.size
+    normest = sqrt(config['norm_data'])    # Estimate norm to scale eigenvector
+    config['norm_rt_data'] = normest
+    
+    npower_iter = config['warmup_iter']    # Number of power iterations
+    z0 = randn(n1,n2)        ## Initial guess
+    z0 = z0/norm(z0,'fro')   ## Power iterations 
+    
+    start_time = time.time()
+    for tt in range(npower_iter):
+        z0 = At(Y*A(z0))
+        z0 = z0/norm(z0,'fro')
+        
+    end_time = time.time()
+    Times = end_time - start_time
+    
+    z = normest * z0                     # Apply scaling 
+    if n2==1:
+        Relerrs = norm(x-exp(-1j*angle(trace(x*z))) *z, 'fro')/norm(x,'fro')
+        config['u_0'] = tile(z,[1,L])
+        truth = A(x)
+        config['truth'] = truth[:,0] ####?####
+        config['norm_truth'] = norm(truth[:,0], 'fro')
+        config['truth_dim'] = truth[1,:].size
+    else:
+        Relerrs = norm(x-exp(-1j*angle(trace(x*z))) *z, 'fro')/norm(x,'fro')
+        config['u_0'] = reshape(tile(z,[1,L]),[z[0].shape,z[1].shape,L])
+        truth = A(x)
+        config['truth'] = truth[:,:,1]
+        config['norm_truth'] = norm(truth[:,:,1],'fro')
+        config['trut_dim'] = truth[:,:,1].shape
+    print('Run time of initialization: %.2f',Times)
+    print('Relative error after initialization: %f', Relerrs)
+    
+    
+    
\ No newline at end of file