Uploaded CDP_processor_ADMM, slight changes in CDP_source_loc_processor

proxtoolbox/Problems/Phase/CDP_processor_ADMM.py

+from numpy.random import randn, random_sample
+from numpy.linalg import norm
+from numpy.fft import fft, ifft, ifft2, fft2
+from numpy import sqrt, conj, tile, mean
+import random
+import numpy as np
+
+def CDP_processor_ADMM(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
+    # integrated into the ProxToolbox by 
+    # Russell Luke, September 2016.
+
+    # 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)
+    config['truth'] = x
+    config['norm_truth'] = norm(x, 'fro')
+    config['truth_dim'] = x.shape
+    
+    ## Make masks and linear sampling operators
+    
+    L = config['product_space_dimension']  # 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 = 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
+    
+    # support constraint: none
+    config['indicator_ampl'] = 1
+    config['abs_illumination'] = 1
+    
+    # Not sure if the following code is correct or if it is needed
+    # It is supposed to be the equivalent to input.support_idx=find(x~=0) 
+    # in matlab which searches x for all nonzero elements and puts the
+    # corresponding indices in a vector. The indices are calculated as if x was a column vector
+    nonzerovector = x.nonzero()
+    # Puts the indices of non zero elements in an array
+    
+    nonzerovector = nonzerovector.transpose
+    #transposes the array
+    
+    for i in range(nonzerovector.shape[0]):
+        config['support_idx'][i] = nonzerovector[i,1]*x.shape[0]+ nonzerovector[i,0]
+        # Is supposed to change the array output into the equivalent column vector indices
+        
+    # 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
+    
+    ## Initialization
+    
+    npower_iter = config['warmup_iter']         # Number of power iterations
+    z0 = randn(n1,n2)
+    z0 = z0/norm(n0, 'fro')                     # Initial guess
+    start_time = time.time()                    # Power iterations 
+    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])
+    elif n1==1:
+        Relerrs = norm(x-exp(-1j*angle(trace(z*x)))*z,'fro')/norm(x,'fro')
+        config['u_0'] = tile(z, [L,1])
+    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].size,z[1].size, L])
+        
+    print('Run time of initialization: %.2f seconds ', Times)
+    print('Relative error after initialization: %.2f ', Relerrs)

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
+from numpy import sqrt, conj, reshape, tile, mean, angle, trace, nonzero, transpose
 import numpy as np
 import random
 
@@ -10,6 +10,8 @@ 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
+# integrated into the ProxToolbox by 
+# Russell Luke, September 2016.
 
 # The input data are coded diffraction patterns about a random complex
 # valued image.
@@ -58,9 +60,23 @@ def CDP_source_loc_processor(config):
     # support constraint: none
     config['indicator_ampl'] = 1
     config['abs_illumination'] = 1
-    config ['Illumination'] = 1
-    ##suport_idx is missing###
+    config['Illumination'] = 1
     
+    # Not sure if the following code is correct or if it is needed
+    # It is supposed to be the equivalent to input.support_idx=find(x~=0) 
+    # in matlab which searches x for all nonzero elements and puts the
+    # corresponding indices in a vector. The indices are calculated as if x was a column vector
+    nonzerovector = x.nonzero()
+    # Puts the indices of non zero elements in an array
+    
+    nonzerovector = nonzerovector.transpose
+    #transposes the array
+    
+    for i in range(nonzerovector.shape[0]):
+        config['support_idx'][i] = nonzerovector[i,1]*x.shape[0]+ nonzerovector[i,0]
+        # Is supposed to change the array output into the equivalent column vector indices
+     
+     
     # Data
     Y=abs(A(x))
     config['rt_data'] = Y
@@ -87,9 +103,16 @@ def CDP_source_loc_processor(config):
         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['truth'] = truth[:,0]
         config['norm_truth'] = norm(truth[:,0], 'fro')
         config['truth_dim'] = truth[1,:].size
+    elif n1==1:
+        Relerrs = norm(x-exp(-1j*angle(trace(z*x)))*z, 'fro')/norm(x,'fro')
+        config['u_0'] = tile(z, [L,1])
+        truth = A(x)
+        config['truth'] = truth[:,0]
+        config['norm_truth'] = norm(truth[0,:],'fro')
+        config['truth_dim'] = truth[0,:].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])