Fixed a bug in CDI_data_prcessor (index adjustment), now rt_data matches matlab result.

Added a test function to check CDI_data_processors to phase

proxtoolbox/Algorithms/RAAR.py

@@ -161,12 +161,10 @@ class RAAR(Algorithm):
         elif self.Ny == 1:
             u1 = tmp[0,:];
             u2 = tmp2[0,:];
-        elif self.Nz == 1:
-            u1 = tmp[:,:,0];
-            u2 = tmp2[:,:,0];
         else:
             u1 = tmp;
             u2 = tmp2;
+
         change = change[1:iter+1];
         gap = gap[1:iter+1];
         shadow_change = shadow_change[1:iter+1]

proxtoolbox/Problems/Phase/CDI_data_processor.py

@@ -37,7 +37,7 @@
 
 from scipy.io import loadmat
 from matplotlib.pyplot import colorbar, show, imshow
-from numpy import log10, log2, ceil, floor, argmax, unravel_index, zeros, where, nonzero, pi, sqrt
+from numpy import log10, log2, ceil, floor, argmax, unravel_index, zeros, where, nonzero, pi, sqrt, ones
 from numpy.fft import fftshift, ifft2
 from numpy.linalg import norm
 from numpy.random import rand
@@ -64,14 +64,14 @@ def CDI_data_processor(config):
     workres = 2**(step_up)*2**(floor(log2(orig_res)))#desired array size
 
     N=int(workres)
-    print(N)
+
     ## center data and use region of interest around center
     #central pixel
     #find max data value
     argmx = unravel_index(argmax(dp), dp.shape)
-    print(argmx)
-    Xi = argmx[0]
-    Xj = argmx[1]
+
+    Xi = argmx[0]+1
+    Xj = argmx[1]+1
 
     #get desired roi:
     #necessary conversion
@@ -83,13 +83,15 @@ def CDI_data_processor(config):
     Nj = 2*(Xj+Dj*(Dj<0)-1)
 
     tmp = zeros((N,N))
-    tmp[int(Di*(Di>0)+1):int(Di*(Di>0)+Ni),int(Dj*(Dj>0)+1):int(Dj*(Dj>0)+Nj)] = dp[int(Xi-Ni/2):int(Xi+Ni/2-1),int(Xj-Nj/2):int(Xj+Nj/2-1)]
+
+    tmp[int(Di*(Di>0)):int(Di*(Di>0)+Ni),int(Dj*(Dj>0)):int(Dj*(Dj>0)+Nj)] = dp[int(Xi-Ni/2)-1:int(Xi+Ni/2-1),int(Xj-Nj/2)-1:int(Xj+Nj/2-1)]
+
     M=(fftshift(tmp))**(.5)
     # M=tmp.^(.5)
     ## define support
     DX = ceil(N/7)
     S = zeros((N,N))
-    S[int(N/4+1+DX):int(3/4*N-1-DX),int(N/4+1+DX):int(3/4*N-1-DX)] = 1
+    S[int(N/4+1+DX)-1:int(3/4*N-1-DX)-1,int(N/4+1+DX):int(3/4*N-1-DX)] = 1
 
     # config['data_ball=config['Nx*config['Ny*data_ball
     config['rt_data']=M
@@ -108,7 +110,7 @@ def CDI_data_processor(config):
     config['supp_phase'] = []
 
     # initial guess
-    config['u_0'] = S*rand(N,N)
+    config['u_0'] = 0.5*S*ones((N,N))#S*rand(N,N)
     config['u_0'] = config['u_0']/norm(config['u_0'],'fro')*config['norm_rt_data'] 
 
     if config['fresnel_nr']*config['magn'] <=  2*pi*sqrt(config['Nx']*config['Ny']):

proxtoolbox/Problems/Phase/phase.py

@@ -241,31 +241,45 @@ class Phase(Problem):
         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['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')
+
+        if self.config['experiment'] == 'JWST':
+            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
+            u1 = f['u1'].value.view(np.complex)
         else:
-            print("No file available to compare to.")
-            return
-        u1 = f['u1'].value.view(np.complex)
+            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
+            u1 = f['u1'].value.view(np.float)
+
         u1 =np.array(u1)
         u1 = u1.T
         print("Compare u1:")
-        print("Nonzero indices matlab:")
-        print(nonzero(u1))
-        print("Nonzero indices python:")
-        print(nonzero(self.output['u1']))
+        #print("Nonzero indices matlab:")
+        #print(nonzero(u1))
+        #print("Nonzero indices python:")
+        #print(nonzero(self.output['u1']))
         print("Nonzero indices equal:")
         print(np.array_equal(nonzero(u1),nonzero(self.output['u1'])))
-        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)
+        #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)
         diff = u1 - self.output['u1']
-        print(diff[nonz])
+        #print(diff[nonz])
         print("Maximum norm of difference:")
         print(np.amax(abs(diff)));
         print("Frobenius norm of difference:")
@@ -274,4 +288,30 @@ class Phase(Problem):
         print(norm(u1))
         print("Frobenius norm of python u1:")
         print(norm(self.output['u1']))
+
+    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')
+            rt_data = f['rt_data'].value.view(np.float)
+
+        rt_data =np.array(rt_data)
+        rt_data = rt_data.T
+        print("Compare rt_data:")
+        print("Nonzero indices equal:")
+        print(np.array_equal(nonzero(rt_data),nonzero(self.config['rt_data'])))
+        diff = rt_data - self.config['rt_data']
+        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']))
         

proxtoolbox/Problems/Phase/tasse_raar_in.py

@@ -85,7 +85,7 @@ new_config = {
 
 #if(strcmp('problem_family,'Phase'))
     ## maximum number of iterations and tolerances
-    'MAXIT' : 1000,
+    'MAXIT' : 1,
     'TOL' : 1e-12,
 
     ## relaxaton parameters in RAAR, HPR and HAAR