Added tests for Siemens_amplitude

Fixed 2 bugs in data_processor (S, matrix multiplication)
Fixed one bug in Approx_PM_Gaussion (use iift2 instead of ifft)

TestSuite/Phase_Siemens_amplitude_demo.py

@@ -7,5 +7,6 @@ from phase import Phase
 
 Siemens = Phase(Siemens_amplitude_in.new_config)
 Siemens.solve()
-#JWST.compare_to_matlab()
+Siemens.compare_to_matlab()
+Siemens.compare_data_to_matlab()
 #JWST.show()

proxtoolbox/Algorithms/RAAR.py

@@ -150,7 +150,7 @@ class RAAR(Algorithm):
             # Since the stopping criteria is on the change in the iterates, this
             # does not matter.
             # graphics
-            
+            print(shadow_change[iter])
             u = tmp_u
             tmp1 = (2*tmp2) - tmp_u
             

proxtoolbox/Problems/Phase/Siemens_processor.py

@@ -60,8 +60,8 @@ def Siemens_processor(config):
         # support constraint.
         config['amplitude'] = config['norm_rt_data']*S/(norm(S,'fro'))
         config['amplitude'] = config['amplitude']/norm(config['amplitude'],'fro')*config['norm_rt_data'][0]
-        config['supp_phase'] = find(ones(m,n))
-        config['illumination_phase'] = find(ones(m,n))
+        config['supp_phase'] = np.nonzero(ones(m,n))
+        config['illumination_phase'] = np.nonzero(ones(m,n))
     elif config['object'] == 'complex':
         # put some phase across S
         points = config['Nx']
@@ -72,8 +72,10 @@ def Siemens_processor(config):
         config['amplitude'] = config['norm_rt_data']*S/(norm(S,'fro'))
         # input.amplitude = input.amplitude/norm(input.amplitude,'fro')*input.norm_rt_data(1)
         G=np.zeros(S.shape)# Gaussian(points,10,[points/2+1,points/2+1])
-        W=S*np.exp(1j*2*np.pi*G)
+        W=np.matmul(S,np.exp(1j*2*np.pi*G))
+        print(W[86,135])
         M=abs(fft2(W))
+        print(M)
         # config['data_ball']=config['Nx']*config['Ny']*data_ball
         config['rt_data']=M
         # standard for the main program is that
@@ -105,7 +107,7 @@ def Siemens_processor(config):
         config['illumination_phase'] = S    
     
     # initial guess
-    config['u_0'] = ifft2(M*np.exp(2*np.pi*1j*rand(dim[0],dim[1])))
+    config['u_0'] = ifft2(M*np.exp(2*np.pi*1j*0.5*np.ones((dim[0],dim[1]))))# rand(dim[0],dim[1])))
     # config['u_0 = ifft2(M.*exp(2*pi*1i*Gaussian(N,N/2,[N/2+1,N/2+1]))).*S
     # config['u_0 = config['u_0/norm(config['u_0,'fro')*config['abs_illumination 
 

proxtoolbox/Problems/Phase/phase.py

@@ -213,7 +213,7 @@ class Phase(Problem):
 #        algorithm = self.config['algorithm'](self.config)
         
         self.output = self.algorithm.run(self.config['u_0'],self.config['TOL'],self.config['MAXIT'])
-        print(self.output['iter'])
+        print('Iterations:' + str(self.output['iter']))
 
     
     def _postsolve(self):
@@ -253,6 +253,9 @@ class Phase(Problem):
             elif self.config['algorithm'] == 'AP':
                 f = h5py.File('Phase_test_data/JWST_u1_ap_' + str(self.config['MAXIT']) + '.mat')
 
+            u1 = f['u1'].value.view(np.complex)
+        elif self.config['data_filename'] == 'Siemens_processor':
+            f = h5py.File('Phase_test_data/siemens_amplitude_u1_' + str(self.config['MAXIT']) + '.mat')
             u1 = f['u1'].value.view(np.complex)
         else:
             if self.config['algorithm'] == 'RAAR':
@@ -313,39 +316,44 @@ class Phase(Problem):
 
         if self.config['data_filename'] == 'CDI_data_processor':
             f = h5py.File('Phase_test_data/CDI_data_processor_rt_data.mat')
-            rt_data = f['rt_data'].value.view(np.float)
-
-        rt_data =np.array(rt_data)
-        rt_data = rt_data.T
-        print("Compare rt_data:")
+            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:")
         print("Nonzero indices equal:")
-        print(np.array_equal(nonzero(rt_data),nonzero(self.config['rt_data'])))
-        diff = rt_data - self.config['rt_data']
+        print(np.array_equal(nonzero(data_mat),nonzero(data_py)))
+        diff = data_mat - data_py
         print("Maximum norm of difference:")
         print(np.amax(abs(diff)));
         print("Frobenius norm of difference:")
         print(norm(diff))
-        print("Frobenius norm of matlab rt_data:")
-        print(norm(rt_data))
-        print("Frobenius norm of python rt_data:")
-        print(norm(self.config['rt_data']))
+        print("Frobenius norm of matlab data:")
+        print(norm(data_mat))
+        print("Frobenius norm of python data:")
+        print(norm(data_py))
 
         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)
 
-        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']))
+            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']))
         

proxtoolbox/ProxOperators/proxoperators.py

@@ -408,7 +408,7 @@ class Approx_PM_Gaussian(ProxOperator):
         tmp = U0_sq - b
         h=sum(sum(tmp * conj(tmp)))
         if h>=epsilon+TOL2:
-            return ifft(U0)
+            return ifft2(U0)
         else:
             return u