CDP_processor.py 3.69 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
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
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 
    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;
     
    
    ## Initialization
    
    npower_iter = config['warmup_iter;                          # Number of power iterations 
    z0 = randn(n1,n2); z0 = z0/norm(z0,'fro'); # Initial guess 
    tic                                        # Power iterations 
    for tt = 1:npower_iter 
        z0 = At(Y.*A(z0)); z0 = z0/norm(z0,'fro');
    end
    Times = toc; 
    
    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);
    elseif(n1==1)
        Relerrs = norm(x - exp(-1i*angle(trace(z'*x))) * z, 'fro')/norm(x,'fro');
        config['u_0 = repmat(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);
    end
    
     fprintf('Run time of initialization: #.2f\n', Times)
     fprintf('Relative error after initialization: #f\n', Relerrs)
     fprintf('\n')
    end