phase.py 9 KB
Newer Older
1
2
# -*- coding: utf-8 -*-

3
from proxtoolbox.Problems.problems import Problem
4
5
6
from proxtoolbox import Algorithms
from proxtoolbox import ProxOperators
from proxtoolbox.ProxOperators.proxoperators import ProxOperator
7
from proxtoolbox.Problems.Phase.JWST_graphics import JWST_graphics
8
from numpy.linalg import norm
9
from numpy import square, sqrt, nonzero
10
11


12
class Phase(Problem):
13
    """
14
    Phase Problem
15
16
17
18
19
20
    """
    config = {
    }
    
    def __init__(self, new_config={}):
        """
21
        The initialization of a Phase instance takes the default configuration
22
23
24
25
26
27
28
        and updates the parameters with the arguments in new_config.
        
        Parameters
        ----------
        new_config : dict, optional - Parameters to initialize the problem. If unspecified, the default config is used.
        """
        self.config.update(new_config)
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46


        #call data processor, read data

        module = __import__(self.config['data_filename'])
        data_processor = getattr(module, self.config['data_filename'])
        data_processor(self.config)
        

        #reshape and rename the data 
        self.config['data_sq'] = self.config['data'];
        self.config['data'] = self.config['rt_data'];
        tmp = self.config['data'].shape;
        if(tmp[0]==1 or tmp[1]==1):
            self.config['data_sq'] = self.config['data_sq'].reshape((self.config['Nx'],self.config['Ny']));
            #the projection algorithms work with the square root of the measurement:
            self.config['data'] = self.config['data'].reshape((self.config['Nx'],self.config['Ny']));

47
        self.config['normM']=self.config['norm_rt_data']; #previously (in matlab) norm_data
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
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
        if 'Nz' not in self.config:
            self.config[Nz] = 1;


        #If method_config[formulation is does not exist, i.e. not specified in 
        #the *_in.m file, use the product space as the default.
        if 'formulation' in self.config:
            formulation = self.config['formulation'];
        else:
            formulation = 'product space';

        # Set the projectors and inputs based on the types of constraints and 
        # experiments
        proxoperators = ['','',''];

        if self.config['constraint'] == 'hybrid':
            proxoperators[0] = 'P_cP'; # This will be problem specific
        elif self.config['constraint'] == 'support only':
            proxoperators[0] = 'P_S';
        elif self.config['constraint'] == 'real and support':
            proxoperators[0] ='P_S_real';
        elif self.config['constraint'] =='nonnegative and support':
            proxoperators[0] ='P_SP';
        elif self.config['constraint'] =='amplitude only':
            proxoperators[0] ='P_amp';
        elif self.config['constraint'] =='minimum amplitude':
            proxoperators[0] = 'P_min_amp';
        elif self.config['constraint'] =='sparse':
            proxoperators[0] = 'not in yet';  
        elif self.config['constraint'] =='phaselift':
            proxoperators[0] = 'P_mean_SP';
        elif self.config['constraint'] =='phaselift2':
            proxoperators[0] ='P_liftM';
            proxoperators[2] ='Approx_PM_Poisson'; # Patrick: This is just to monitor the change of phases!  

        if self.config['experiment'] == 'single diffraction':
            if self.config[distance] == 'far field':
                if self.config['constraint'] == 'phaselift':
                    proxoperators[1] = 'P_Rank1';
                elif self.config['constraint'] == 'phaselift2':
                    proxoperators[1] = 'P_rank1_SR';
                else:
                    if self.config['noise'] == 'Poisson':
                        proxoperators[1] ='Approx_PM_Poisson';
                    else:
                        proxoperators[1] ='Approx_PM_Gaussian';
            else:
                proxoperators[1]='P_Fresnel';

        # The following selects the projectors for diversity diffraction not
        # performed in the product space. So far only used for RCAAR.
        elif self.config['experiment'] == 'diversity diffraction' and formulation == 'sequential':
            proxoperators[1] = 'Approx_P_RCAAR_JWST_Poisson';
            proxoperators[0] = proxoperators[1];
        elif self.config['experiment'] == 'JWST': 
            proxoperators[1] = 'Approx_P_JWST_Poisson';  
            proxoperators[2] = proxoperators[0];
            proxoperators[0] = 'P_diag';
        elif self.config['experiment'] == 'CDP':
            proxoperators[1] = 'P_CDP';  
            proxoperators[2] = proxoperators[0];
            proxoperators[0] = 'P_diag';
        elif self.config['experiment'] == 'ptychography':
            proxoperators[1] = 'not in yet';
        elif self.config['experiment'] == 'complex':
            proxoperators[1] = 'not in yet';
        elif self.config['constraint'] == 'phaselift':
            proxoperators[1] ='P_PL_lowrank';

        self.config['proxoperators'] = [];

        for prox in proxoperators:
            self.config['proxoperators'].append(getattr(ProxOperators, prox))

        # input.Proj1_input.F=F;  % is it any more expensive to pass everything
        # into the projectors rather than just a selection of data and
        # parameters?  If not, and we pass everything anyway, there is no need
        # to create a new structure element.

        if 'product_space_dimension' not in self.config:
            self.config['product_space_dimension'] = 1;

        # set the animation program:
        self.config['animation']='Phase_animation';
        #
        # if you are only working with two sets but
        # want to do averaged projections
        # (= alternating projections on the product space)
        # or RAAR on the product space (=swarming), then
        # you will want to change product_space_dimension=2
        # and adjust your input files and projectors accordingly. 
        # you could also do this within the data processor

141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
        self.config['TOL2'] = 1e-15;

        #To estimate the gap in the sequential formulation, we build the
        # appropriate point in the product space. This allows for code reuse.
        # Note for sequential diversity diffraction, input.Proj1 is the "RCAAR"
        # version of the function.
        if formulation == 'sequential':
            for j in range(self.config['product_space_dimension']):
                self.config['proj_iter'] =j;
                proj1 = self.config['proxoperators'][0](self.config)
                u_1[:,:,j]= proj1.work(self.config['u_0']);
                self.config['proj_iter'] = mod(j,config['product_space_dimension'])+1;
                proj1 = self.config['proxoperators'][0](self.config)
                u_1[:,:,j]= proj1.work(self.config['u_0']);
            end;
        else: #i.e. formulation=='product space'
            proj1 = self.config['proxoperators'][0](self.config)
            u_1 = proj1.work(self.config['u_0']);
            proj2 = self.config['proxoperators'][1](self.config)
160
            u_2 = proj2.work(u_1);
161

162

163
164
        # estimate the gap in the relevant metric
        if self.config['Nx'] ==1 or self.config['Ny']==1 :
165
            tmp_gap = square(norm(u_1-u_2)/self.config['norm_rt_data']);
166
167
168
169
        else:
            tmp_gap=0;
            for j in range(self.config['product_space_dimension']):
                # compute (||P_Sx-P_Mx||/normM)^2:
170
                tmp_gap = tmp_gap+(norm(u_1[:,:,j]-u_2[:,:,j])/square(self.config['norm_rt_data']));
171
172
173
174
175
176
177
178

        gap_0=sqrt(tmp_gap);

        # sets the set fattening to be a percentage of the
        # initial gap to the unfattened set with 
        # respect to the relevant metric (KL or L2), 
        # that percentage given by
        # input.data_ball input by the user.
179
        self.config['data_ball']=self.config['data_ball']*gap_0;
180
181
        # the second tolerance relative to the oder of 
        # magnitude of the metric
182
        self.config['TOL2'] = self.config['data_ball']*1e-15; 
183
184

        self.algorithm = getattr(Algorithms, self.config['algorithm'])(self.config);
185
186
187
188
189
190
191
        
    
    
    def _presolve(self):
        """
        Prepares argument for actual solving routine
        """
192
        
193
194
195
196
197
198
199
    
    def _solve(self):
        """
        Runs the algorithm to solve the given sudoku problem
        """
#        algorithm = self.config['algorithm'](self.config)
        
200
201
202
        self.output = dict();
        
        self.output['u1'],self.output['u2'],self.output['iter'],self.output['change'],self.output['gap'] = \
203
            self.algorithm.run(self.config['u_0'],self.config['TOL'],self.config['MAXIT'])
204

205
        #print(self.iters)
206
207
        #print(nonzero(self.u1))
        #print(self.u1[self.u1.nonzero()]);
208
209
210
211
        #print(norm(self.u1));
        #print(self.u1[64,16]);
        #print(norm(self.u2));
        #print(self.u1[71,46]);
212
213
        #print(self.gap);
        #print(self.change);
214
215
216
217
218
219
220
221
222
223
224
225
226
    
    
    def _postsolve(self):
        """
        Processes the solution and generates the output
        """
        
        
        
    def show(self):
        """
        Generates graphical output from the solution
        """
227
228
229
230

        print(self.elapsed_time)
        JWST_graphics(self.config,self.output)
        
231