Solved bug in Proxoperator Approx_PM_Poisson. It was using ifft2 instead of fft2.

proxtoolbox/Problems/Phase/CDI_data_processor.py

@@ -91,7 +91,7 @@ def CDI_data_processor(config):
     ## define support
     DX = ceil(N/7)
     S = zeros((N,N))
-    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
+    S[int(N/4+1+DX)-1:int(3/4*N-1-DX),int(N/4+1+DX)-1:int(3/4*N-1-DX)] = 1
 
     # config['data_ball=config['Nx*config['Ny*data_ball
     config['rt_data']=M

proxtoolbox/Problems/Phase/phase.py

@@ -8,7 +8,7 @@ from proxtoolbox.Problems.Phase.JWST_graphics import JWST_graphics
 from numpy.linalg import norm
 import numpy as np
 import h5py
-from numpy import square, sqrt, nonzero
+from numpy import square, sqrt, nonzero, size
 
 
 class Phase(Problem):
@@ -173,7 +173,6 @@ class Phase(Problem):
             proj2 = self.config['proxoperators'][1](self.config)
             u_2 = proj2.work(u_1);
 
-
         # estimate the gap in the relevant metric
         if self.config['Nx'] ==1 or self.config['Ny']==1 :
             tmp_gap = square(norm(u_1-u_2)/self.config['norm_rt_data']);
@@ -184,8 +183,9 @@ class Phase(Problem):
             for j in range(self.config['product_space_dimension']):
                 # compute (||P_Sx-P_Mx||/normM)^2:
                 tmp_gap = tmp_gap+(norm(u_1[:,:,j]-u_2[:,:,j])/square(self.config['norm_rt_data']));
-
+        
         gap_0=sqrt(tmp_gap);
+        print(gap_0)
 
         # sets the set fattening to be a percentage of the
         # initial gap to the unfattened set with 

proxtoolbox/ProxOperators/proxoperators.py

@@ -360,7 +360,7 @@ class P_SP(ProxOperator):
 
     def work(self,u):
         p_SP=zeros(shape(u))
-        p_SP[self.support_idx] = np.max(real(u[self.support_idx]),0)
+        p_SP[self.support_idx] = np.maximum(real(u[self.support_idx]),0)
         return p_SP
 
 #                      Approx_PM_Poisson.m
@@ -446,7 +446,7 @@ class Approx_PM_Poisson(ProxOperator):
         Ib = self.Ib
         epsilon = self.epsilon
 
-        U = ifft2(u)
+        U = fft2(u)
         U_sq = U * conj(U)
         tmp = U_sq/b; tmp[Ib] = 1
         U_sq[Ib]=0