CDP_processor.py 3.68 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
def CDP_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. 
    
    ## 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'] = s.shape
    
    
    ## Make masks and linear sampling operators
    
    L = config['product_space_dimension']                  # Number of masks 
    if n2==1:
        Masks = np.random.choice(np.array([1j, -1j, 1, -1]),(n1,L))
    elif n1==1:
        Masks = np.random.choice(np.array([1j, -1j, 1, -1]),(L,n2))
    else:
        Masks = 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(np.array([1j, -1j, 1, -1]),(n1,n2))
    
    # Sample magnitudes and make masks 
    temp = rand(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 = @(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)
        # 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
    end
    
    # Data 
54
55
56
57
58
59
60
    Y = abs(A(x))
    config['rt_data']=Y;
    Y=Y**2
    config['data']=Y
    config['norm_data']=sum(Y(:))/numel(Y)
    normest = sqrt(config['norm_data']) # Estimate norm to scale eigenvector 
    config['norm_rt_data']=normest
61
62
63
64
     
    
    ## Initialization
    
65
    npower_iter = config['warmup_iter'];                          # Number of power iterations 
66
67
    z0 = randn(n1,n2); z0 = z0/norm(z0,'fro'); # Initial guess 
    tic                                        # Power iterations 
68
    for tt in range(npower_iter): 
69
        z0 = At(Y.*A(z0)); z0 = z0/norm(z0,'fro');
70
71

    Times = toc 
72
73
    
    z = normest * z0;                   # Apply scaling 
74
    if n2==1:
75
        Relerrs = norm(x - exp(-1i*angle(trace(x'*z))) * z, 'fro')/norm(x,'fro');
76
77
        config['u_0'] = repmat(z,1,L);
    elif n1==1:
78
        Relerrs = norm(x - exp(-1i*angle(trace(z'*x))) * z, 'fro')/norm(x,'fro');
79
80
        config['u_0'] = repmat(z,L,1);
    else:
81
        Relerrs = norm(x - exp(-1i*angle(trace(x'*z))) * z, 'fro')/norm(x,'fro');
82
        config['u_0']=reshape(repmat(z,[1 L]), size(z,1), size(z,2), L);
83
    
84
85
86
     print('Run time of initialization: #.2f\n', Times)
     print('Relative error after initialization: #f\n', Relerrs)
     print('\n')