phase.py 16.3 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

        #call data processor, read data

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

alexander.dornheim's avatar
alexander.dornheim committed
38
39
40
41
42
43
        # reshape and rename the data
        if('data_sq' in self.config):
            pass    # already put data into required format, skip
        else:
            self.config['data_sq'] = self.config['data']
            self.config['data'] = self.config['rt_data']
44
45
            if('norm_data' in self.config):
                self.config['norm_data_sq']= self.config['norm_data']
alexander.dornheim's avatar
alexander.dornheim committed
46
47
            self.config['norm_data']=self.config['norm_rt_data']
            
48
            tmp = self.config['data'].shape
alexander.dornheim's avatar
alexander.dornheim committed
49
50
51
52
53
54
55
            if(tmp[0]==1 or tmp[1]==1):
                    self.config['data_sq'] = self.config['data_sq'].reshape((self.config['Nx'],self.config['Ny']))
                    #the prox algorithms work with the square root of the measurement:
                    self.config['data'] = self.config['data'].reshape((self.config['Nx'],self.config['Ny']))
            
            if 'Nz' not in self.config:
                self.config['Nz'] = 1
56
57
58
59
60


        #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:
61
            formulation = self.config['formulation']
62
        else:
63
            formulation = 'product space'
64
65
66

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

        if self.config['constraint'] == 'hybrid':
70
            proxoperators[0] = 'P_cP' # This will be problem specific
71
        elif self.config['constraint'] == 'support only':
72
            proxoperators[0] = 'P_S'
73
        elif self.config['constraint'] == 'real and support':
74
            proxoperators[0] ='P_S_real'
75
        elif self.config['constraint'] =='nonnegative and support':
76
            proxoperators[0] ='P_SP'
77
        elif self.config['constraint'] =='amplitude only':
78
            proxoperators[0] ='P_amp'
alexander.dornheim's avatar
alexander.dornheim committed
79
        elif self.config['constraint'] == 'phase on support':
alexander.dornheim's avatar
alexander.dornheim committed
80
            proxoperators[0] ='P_Amod'
81
        elif self.config['constraint'] =='minimum amplitude':
82
            proxoperators[0] = 'P_min_amp'
83
        elif self.config['constraint'] =='sparse':
84
            proxoperators[0] = 'not in yet' 
85
        elif self.config['constraint'] =='phaselift':
86
            proxoperators[0] = 'P_mean_SP'
87
        elif self.config['constraint'] =='phaselift2':
88
89
            proxoperators[0] ='P_liftM'
            proxoperators[2] ='Approx_PM_Poisson' # Patrick: This is just to monitor the change of phases!  
90

alexander.dornheim's avatar
alexander.dornheim committed
91
        if self.config['experiment'] == 'single diffraction' or self.config['experiment'] == 'CDI':
92
            if self.config['distance'] == 'far field':
93
                if self.config['constraint'] == 'phaselift':
94
                    proxoperators[1] = 'P_Rank1'
95
                elif self.config['constraint'] == 'phaselift2':
96
                    proxoperators[1] = 'P_rank1_SR'
97
98
                else:
                    if self.config['noise'] == 'Poisson':
99
                        proxoperators[1] ='Approx_PM_Poisson'
100
                    else:
101
                        proxoperators[1] ='Approx_PM_Gaussian'
102
            else:
103
104
105
106
107
108
109
110
111
112
113
114
115
                if self.config['noise'] == 'Poisson':
                    proxoperators[1]='Approx_P_FreFra_Poisson'
        # possibly not necessary, but we set the prox operators here for specific
        # named data sets
        elif (self.config['experiment'] == 'Krueger') or (self.config['experiment'] == 'Near_field_cell_syn'):
             proxoperators[0]='P_Amod'
             proxoperators[1]='Approx_P_FreFra_Poisson'
        # The following selects the prox mappings for diversity diffraction not
        # performed in the product space. So far only used for RCAAR.
        elif self.config['experiment'] in ('dict', 'dictyM103_stx6_600frames','xenopus','living_worm'):
            proxoperators[0]='P_Amod' # Not sure this is the appropriate prox operator for these
                                   # experiments...
            proxoperators[1]='Approx_P_FreFra_Poisson'
116
117
118
119

        # 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':
120
121
            proxoperators[1] = 'Approx_P_RCAAR_JWST_Poisson'
            proxoperators[0] = proxoperators[1]
122
        elif self.config['experiment'] == 'JWST': 
123
124
125
            proxoperators[1] = 'Approx_P_JWST_Poisson'  
            proxoperators[2] = proxoperators[0]
            proxoperators[0] = 'P_diag'
126
        elif self.config['experiment'] == 'CDP':
127
128
129
            proxoperators[1] = 'P_CDP'
            proxoperators[2] = proxoperators[0]
            proxoperators[0] = 'P_diag'
130
        elif self.config['experiment'] == 'ptychography':
131
            proxoperators[1] = 'not in yet'
132
        elif self.config['experiment'] == 'complex':
133
            proxoperators[1] = 'not in yet'
134
        elif self.config['constraint'] == 'phaselift':
135
            proxoperators[1] ='P_PL_lowrank'
136

137
        self.config['proxoperators'] = []
138
139

        for prox in proxoperators:
140
141
            if prox != '':
                self.config['proxoperators'].append(getattr(ProxOperators, prox))
142
143
144
145
146
147
148

        # 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:
149
            self.config['product_space_dimension'] = 1
150
151

        # set the animation program:
152
        self.config['animation']='Phase_animation'
153
154
155
156
157
158
159
160
161
        #
        # 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

162
        self.config['TOL2'] = 1e-15
163
164
165
166
167
168
169

        #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']):
170
                self.config['proj_iter'] =j
171
                proj1 = self.config['proxoperators'][0](self.config)
172
173
                u_1[:,:,j]= proj1.work(self.config['u_0'])
                self.config['proj_iter'] = mod(j,config['product_space_dimension'])+1
174
                proj1 = self.config['proxoperators'][0](self.config)
175
176
                u_1[:,:,j]= proj1.work(self.config['u_0'])
            end
177
178
        else: #i.e. formulation=='product space'
            proj1 = self.config['proxoperators'][0](self.config)
179
            u_1 = proj1.work(self.config['u_0'])
180
            proj2 = self.config['proxoperators'][1](self.config)
181
            u_2 = proj2.work(u_1)
182
183
184

        # estimate the gap in the relevant metric
        if self.config['Nx'] ==1 or self.config['Ny']==1 :
185
            tmp_gap = square(norm(u_1-u_2)/self.config['norm_rt_data'])
186
        elif self.config['product_space_dimension'] == 1:
187
            tmp_gap = (norm(u_1-u_2)/self.config['norm_rt_data'])**2
188
        else:
189
            tmp_gap=0
190
            for j in range(self.config['product_space_dimension']):
191
                # compute (||P_Sx-P_Mx||/norm_data)^2:
192
                tmp_gap = tmp_gap+(norm(u_1[:,:,j]-u_2[:,:,j])/self.config['norm_rt_data'])**2
193
        
194
        gap_0=sqrt(tmp_gap)
195
196
197
198
199
200

        # 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.
201
        self.config['data_ball']=self.config['data_ball']*gap_0
202
203
        # the second tolerance relative to the oder of 
        # magnitude of the metric
204
        self.config['TOL2'] = self.config['data_ball']*1e-15
205
        self.config['proxoperators']
206
        self.algorithm = getattr(Algorithms, self.config['algorithm'])(self.config)
207
208
209
210
211
212
213
        
    
    
    def _presolve(self):
        """
        Prepares argument for actual solving routine
        """
214
        
215
216
217
218
219
220
221
    
    def _solve(self):
        """
        Runs the algorithm to solve the given sudoku problem
        """
#        algorithm = self.config['algorithm'](self.config)
        
222
        self.output = self.algorithm.run(self.config['u_0'],self.config['TOL'],self.config['MAXIT'])
223
        print('Iterations:' + str(self.output['iter']))
224

225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
    
    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)
240
241
        graphics = getattr(Graphics,self.config['graphics_display'])
        graphics(self.config,self.output)
242
243
244
245
246

    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
247
        For result to match u_0 should be chosen as np.multiply(config['abs_illumination'],exp(1j*2*pi*0.5*np.ones(newres)))']
248
        """
249

alexander.dornheim's avatar
alexander.dornheim committed
250
251
        print(self.config['proxoperators'])

252
        if self.config['experiment'] == 'JWST':
253
254
255
256
257
258
259
260
261
262
263
            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')

264
            u1 = f['u1'].value.view(np.complex)
265
        elif self.config['data_filename'] == 'Siemens_processor' and self.config['constraint'] == 'amplitude':
266
            f = h5py.File('Phase_test_data/siemens_amplitude_u1_' + str(self.config['MAXIT']) + '.mat')
267
            u1 = f['u1'].value.view(np.complex)
268
269
270
        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)
271
272
273
        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)
274
        else:
275
276
277
278
279
280
281
282
283
284
285
286
            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
287
                    f = h5py.File('Phase_test_data/tasse_u1_'+ str(self.config['MAXIT']) + '_beta_0_5.mat')
288
                else:
alexander.dornheim's avatar
alexander.dornheim committed
289
290
                    print("No file available to compare to.")
                    return
291
292
            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')
293
            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
294
                        f = h5py.File('Phase_test_data/tasse_u1_ap_' + str(self.config['MAXIT']) + '.mat')
295

296
            u1 = f['u1'].value.view(np.float64)
297

298
299
300
        u1 =np.array(u1)
        u1 = u1.T
        print("Compare u1:")
301
302
303
304
        #print("Nonzero indices matlab:")
        #print(nonzero(u1))
        #print("Nonzero indices python:")
        #print(nonzero(self.output['u1']))
305
306
        print("Nonzero indices equal:")
        print(np.array_equal(nonzero(u1),nonzero(self.output['u1'])))
307
308
309
310
311
312
        #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)
313
        diff = u1 - self.output['u1']
314
        #print(diff[nonz])
315
        print("Maximum norm of difference:")
316
        print(np.amax(abs(diff)))
317
        print("Frobenius norm of difference:")
318
319
320
321
322
        print(norm(diff))
        print("Frobenius norm of matlab u1:")
        print(norm(u1))
        print("Frobenius norm of python u1:")
        print(norm(self.output['u1']))
323
324
325
326
327
328
329
330
331
332

    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')
333
334
335
336
337
338
339
340
341
342
            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:")
343
        print("Nonzero indices equal:")
344
345
        print(np.array_equal(nonzero(data_mat),nonzero(data_py)))
        diff = data_mat - data_py
346
        print("Maximum norm of difference:")
347
        print(np.amax(abs(diff)))
348
349
        print("Frobenius norm of difference:")
        print(norm(diff))
350
351
352
353
        print("Frobenius norm of matlab data:")
        print(norm(data_mat))
        print("Frobenius norm of python data:")
        print(norm(data_py))
354
355
356
357
358

        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)

359
360
361
362
363
364
365
            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:")
366
            print(np.amax(abs(diff)))
367
368
369
370
371
372
            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']))
373