phase.py 15.6 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
61
62
63
64
65
66


        #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
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
74
75
        elif self.config['constraint'] == 'real and support':
            proxoperators[0] ='P_S_real';
        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
80
81
        elif self.config['constraint'] == 'phase on support':
            pass
            #proxoperators[0] ='P_Amod' currently not working in ProxPython
82
        elif self.config['constraint'] =='minimum amplitude':
83
            proxoperators[0] = 'P_min_amp'
84
        elif self.config['constraint'] =='sparse':
85
            proxoperators[0] = 'not in yet' 
86
        elif self.config['constraint'] =='phaselift':
87
            proxoperators[0] = 'P_mean_SP'
88
        elif self.config['constraint'] =='phaselift2':
89
90
            proxoperators[0] ='P_liftM'
            proxoperators[2] ='Approx_PM_Poisson' # Patrick: This is just to monitor the change of phases!  
91

alexander.dornheim's avatar
alexander.dornheim committed
92
        if self.config['experiment'] == 'single diffraction' or self.config['experiment'] == 'CDI':
93
            if self.config['distance'] == 'far field':
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
                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:
129
130
            if prox != '':
                self.config['proxoperators'].append(getattr(ProxOperators, prox))
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150

        # 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

151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
        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)
170
            u_2 = proj2.work(u_1);
171
172
173

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

214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
    
    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)
229
230
        graphics = getattr(Graphics,self.config['graphics_display'])
        graphics(self.config,self.output)
231
232
233
234
235

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

alexander.dornheim's avatar
alexander.dornheim committed
239
240
        print(self.config['proxoperators'])

241
        if self.config['experiment'] == 'JWST':
242
243
244
245
246
247
248
249
250
251
252
            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')

253
            u1 = f['u1'].value.view(np.complex)
254
        elif self.config['data_filename'] == 'Siemens_processor' and self.config['constraint'] == 'amplitude':
255
            f = h5py.File('Phase_test_data/siemens_amplitude_u1_' + str(self.config['MAXIT']) + '.mat')
256
            u1 = f['u1'].value.view(np.complex)
257
258
259
        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)
260
261
262
        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)
263
        else:
264
265
266
267
268
269
270
271
272
273
274
275
            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
276
                    f = h5py.File('Phase_test_data/tasse_u1_'+ str(self.config['MAXIT']) + '_beta_0_5.mat')
277
                else:
alexander.dornheim's avatar
alexander.dornheim committed
278
279
                    print("No file available to compare to.")
                    return
280
281
            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')
282
            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
283
                        f = h5py.File('Phase_test_data/tasse_u1_ap_' + str(self.config['MAXIT']) + '.mat')
284

285
            u1 = f['u1'].value.view(np.float64)
286

287
288
289
        u1 =np.array(u1)
        u1 = u1.T
        print("Compare u1:")
290
291
292
293
        #print("Nonzero indices matlab:")
        #print(nonzero(u1))
        #print("Nonzero indices python:")
        #print(nonzero(self.output['u1']))
294
295
        print("Nonzero indices equal:")
        print(np.array_equal(nonzero(u1),nonzero(self.output['u1'])))
296
297
298
299
300
301
        #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)
302
        diff = u1 - self.output['u1']
303
        #print(diff[nonz])
304
        print("Maximum norm of difference:")
305
        print(np.amax(abs(diff)));
306
        print("Frobenius norm of difference:")
307
308
309
310
311
        print(norm(diff))
        print("Frobenius norm of matlab u1:")
        print(norm(u1))
        print("Frobenius norm of python u1:")
        print(norm(self.output['u1']))
312
313
314
315
316
317
318
319
320
321

    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')
322
323
324
325
326
327
328
329
330
331
            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:")
332
        print("Nonzero indices equal:")
333
334
        print(np.array_equal(nonzero(data_mat),nonzero(data_py)))
        diff = data_mat - data_py
335
336
337
338
        print("Maximum norm of difference:")
        print(np.amax(abs(diff)));
        print("Frobenius norm of difference:")
        print(norm(diff))
339
340
341
342
        print("Frobenius norm of matlab data:")
        print(norm(data_mat))
        print("Frobenius norm of python data:")
        print(norm(data_py))
343
344
345
346
347

        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)

348
349
350
351
352
353
354
355
356
357
358
359
360
361
            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']))
362