phase.py 15.2 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 import Graphics
8
from numpy.linalg import norm
9
10
import numpy as np
import h5py
11
from numpy import square, sqrt, nonzero, size
12
13


14
class Phase(Problem):
15
    """
16
    Phase Problem
17
18
19
20
21
22
    """
    config = {
    }
    
    def __init__(self, new_config={}):
        """
23
        The initialization of a Phase instance takes the default configuration
24
25
26
27
28
29
        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.
        """
30
        self.config.update(new_config)     
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47

        #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']));

48
        self.config['normM']=self.config['norm_rt_data']; #previously (in matlab) norm_data
49
        if 'Nz' not in self.config:
50
            self.config['Nz'] = 1;
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


        #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':
85
            if self.config['distance'] == 'far field':
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
                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:
121
122
            if prox != '':
                self.config['proxoperators'].append(getattr(ProxOperators, prox))
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142

        # 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

143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
        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)
162
            u_2 = proj2.work(u_1);
163
164
165

        # estimate the gap in the relevant metric
        if self.config['Nx'] ==1 or self.config['Ny']==1 :
166
            tmp_gap = square(norm(u_1-u_2)/self.config['norm_rt_data']);
167
        elif self.config['product_space_dimension'] == 1:
168
            tmp_gap = (norm(u_1-u_2)/self.config['norm_rt_data'])**2
169
170
171
172
        else:
            tmp_gap=0;
            for j in range(self.config['product_space_dimension']):
                # compute (||P_Sx-P_Mx||/normM)^2:
173
                tmp_gap = tmp_gap+(norm(u_1[:,:,j]-u_2[:,:,j])/self.config['norm_rt_data'])**2
174
        
175
176
177
178
179
180
181
        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.
182
        self.config['data_ball']=self.config['data_ball']*gap_0;
183
184
        # the second tolerance relative to the oder of 
        # magnitude of the metric
185
        self.config['TOL2'] = self.config['data_ball']*1e-15; 
186
        self.config['proxoperators']
187
        self.algorithm = getattr(Algorithms, self.config['algorithm'])(self.config);
188
189
190
191
192
193
194
        
    
    
    def _presolve(self):
        """
        Prepares argument for actual solving routine
        """
195
        
196
197
198
199
200
201
202
    
    def _solve(self):
        """
        Runs the algorithm to solve the given sudoku problem
        """
#        algorithm = self.config['algorithm'](self.config)
        
203
        self.output = self.algorithm.run(self.config['u_0'],self.config['TOL'],self.config['MAXIT'])
204
        print('Iterations:' + str(self.output['iter']))
205

206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
    
    def _postsolve(self):
        """
        Processes the solution and generates the output
        """
        
        
        
    def show(self):
        """
        Generates graphical output from the solution
        """
        
        print("Calculation time:")
        print(self.elapsed_time)
221
222
        graphics = getattr(Graphics,self.config['graphics_display'])
        graphics(self.config,self.output)
223
224
225
226
227

    def compare_to_matlab(self):
        """
        Routine to test and verify results by comparing to matlab
        Note that this is only for development and should not be used by a normal user
alexander.dornheim's avatar
alexander.dornheim committed
228
        For result to match u_0 should be chosen as np.multiply(config['abs_illumination'],exp(1j*2*pi*0.5*np.ones(newres)))']
229
        """
230

alexander.dornheim's avatar
alexander.dornheim committed
231
232
        print(self.config['proxoperators'])

233
        if self.config['experiment'] == 'JWST':
234
235
236
237
238
239
240
241
242
243
244
            if self.config['algorithm'] == 'RAAR':
                if self.config['MAXIT'] == 1:
                    f = h5py.File('Phase_test_data/u1_1.mat')
                elif self.config['MAXIT'] == 500 :
                    f = h5py.File('Phase_test_data/u1_500.mat')
                else:
                    print("No file available to compare to.")
                    return
            elif self.config['algorithm'] == 'AP':
                f = h5py.File('Phase_test_data/JWST_u1_ap_' + str(self.config['MAXIT']) + '.mat')

245
            u1 = f['u1'].value.view(np.complex)
246
        elif self.config['data_filename'] == 'Siemens_processor' and self.config['constraint'] == 'amplitude':
247
            f = h5py.File('Phase_test_data/siemens_amplitude_u1_' + str(self.config['MAXIT']) + '.mat')
248
            u1 = f['u1'].value.view(np.complex)
249
250
251
        elif self.config['data_filename'] == 'Siemens_processor' and self.config['constraint'] == 'nonnegative and support':
            f = h5py.File('Phase_test_data/siemens_nonneg_u1_' + str(self.config['MAXIT']) + '.mat')
            u1 = f['u1'].value.view(np.float64)
252
253
254
        elif self.config['data_filename'] == 'Siemens_processor' and self.config['constraint'] == 'real and support':
            f = h5py.File('Phase_test_data/siemens_real_u1_' + str(self.config['MAXIT']) + '.mat')
            u1 = f['u1'].value.view(np.float64)
255
        else:
256
257
258
259
260
261
262
263
264
265
266
267
            if self.config['algorithm'] == 'RAAR':
                if self.config['beta_0'] == 0.95:
                    if self.config['MAXIT'] == 1000 :
                        f = h5py.File('Phase_test_data/tasse_u1_1000.mat')
                    elif self.config['MAXIT'] == 20:
                        f = h5py.File('Phase_test_data/tasse_u1_20.mat')
                    elif self.config['MAXIT'] == 1:
                        f = h5py.File('Phase_test_data/tasse_u1_1.mat')
                    else:
                        print("No file available to compare to.")
                        return
                elif self.config['beta_0'] == 0.50:
alexander.dornheim's avatar
alexander.dornheim committed
268
                    f = h5py.File('Phase_test_data/tasse_u1_'+ str(self.config['MAXIT']) + '_beta_0_5.mat')
269
                else:
alexander.dornheim's avatar
alexander.dornheim committed
270
271
                    print("No file available to compare to.")
                    return
272
273
            elif self.config['algorithm'] == 'AP' and self.config['constraint'] == 'support only':
                        f = h5py.File('Phase_test_data/tasse_supp_u1_ap_' + str(self.config['MAXIT']) + '.mat')
274
            elif ( self.config['algorithm'] == 'AP' or self.config['algorithm'] == 'AP_expert') and self.config['constraint'] == 'nonnegative and support':
alexander.dornheim's avatar
alexander.dornheim committed
275
                        f = h5py.File('Phase_test_data/tasse_u1_ap_' + str(self.config['MAXIT']) + '.mat')
276

277
            u1 = f['u1'].value.view(np.float64)
278

279
280
281
        u1 =np.array(u1)
        u1 = u1.T
        print("Compare u1:")
282
283
284
285
        #print("Nonzero indices matlab:")
        #print(nonzero(u1))
        #print("Nonzero indices python:")
        #print(nonzero(self.output['u1']))
286
287
        print("Nonzero indices equal:")
        print(np.array_equal(nonzero(u1),nonzero(self.output['u1'])))
288
289
290
291
292
293
        #print("Nonzero values matlab:")
        #print(u1[nonzero(u1)])
        #print("Nonzero values python:")
        #print(self.output['u1'][nonzero(self.output['u1'])])
        #print("Difference at nonzero values:")
        #nonz = nonzero(u1)
294
        diff = u1 - self.output['u1']
295
        #print(diff[nonz])
296
        print("Maximum norm of difference:")
297
        print(np.amax(abs(diff)));
298
        print("Frobenius norm of difference:")
299
300
301
302
303
        print(norm(diff))
        print("Frobenius norm of matlab u1:")
        print(norm(u1))
        print("Frobenius norm of python u1:")
        print(norm(self.output['u1']))
304
305
306
307
308
309
310
311
312
313

    def compare_data_to_matlab(self):
        """
        Routine to test and verify results by comparing to matlab
        Note that this is only for development and should not be used by a normal user
        For result to match u_0 should be chosen as np.multiply(config['abs_illumination'],exp(1j*2*pi*0.5*np.ones(newres)))'] =
        """

        if self.config['data_filename'] == 'CDI_data_processor':
            f = h5py.File('Phase_test_data/CDI_data_processor_rt_data.mat')
314
315
316
317
318
319
320
321
322
323
            data_mat = f['rt_data'].value.view(np.float)
            data_py = self.config['rt_data']
        elif self.config['data_filename'] == 'Siemens_processor':
            f = h5py.File('Phase_test_data/Siemens_processor_truth.mat')
            data_mat = f['truth'].value.view(np.float)
            data_py = self.config['truth']

        data_mat =np.array(data_mat)
        data_mat = data_mat.T
        print("Compare:")
324
        print("Nonzero indices equal:")
325
326
        print(np.array_equal(nonzero(data_mat),nonzero(data_py)))
        diff = data_mat - data_py
327
328
329
330
        print("Maximum norm of difference:")
        print(np.amax(abs(diff)));
        print("Frobenius norm of difference:")
        print(norm(diff))
331
332
333
334
        print("Frobenius norm of matlab data:")
        print(norm(data_mat))
        print("Frobenius norm of python data:")
        print(norm(data_py))
335
336
337
338
339

        if self.config['data_filename'] == 'CDI_data_processor':
            f = h5py.File('Phase_test_data/CDI_data_processor_S.mat')
            S = f['S'].value.view(np.float)

340
341
342
343
344
345
346
347
348
349
350
351
352
353
            S =np.array(S)
            S = S.T
            print("Compare S:")
            print("Nonzero indices equal:")
            print(np.array_equal(nonzero(S),nonzero(self.config['abs_illumination'])))
            diff = S - self.config['abs_illumination']
            print("Maximum norm of difference:")
            print(np.amax(abs(diff)));
            print("Frobenius norm of difference:")
            print(norm(diff))
            print("Frobenius norm of matlab S:")
            print(norm(S))
            print("Frobenius norm of python S:")
            print(norm(self.config['abs_illumination']))
354