Tried fft from pyfftw (is a bit fast than fft from numpy.fft, but commented out pyfftw)

Tried naive implementation of magproj, should work now, but commented out
Replaced u*conj(u) by u.real**2+u.imag**2 in magproj. Computation time for 500 JWST iterations is down from 1.3s to 0.5 s (much faster!!!). Computation before was complex!!!
Did the same in Approx_P_JWST_Poisson (also somewhat faster, since using real muliplication, log now).
Implemented P_S

proxtoolbox/Problems/Phase/CDI_data_processor.py

@@ -37,8 +37,9 @@
 
 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, ones
+from numpy import log10, log2, ceil, floor, argmax, unravel_index, zeros, where, nonzero, pi, sqrt, ones, float64
 from numpy.fft import fftshift, ifft2
+#from pyfftw.interfaces.scipy_fftpack import fftshift, ifft2
 from numpy.linalg import norm
 from numpy.random import rand
 
@@ -110,7 +111,7 @@ def CDI_data_processor(config):
     config['supp_phase'] = []
 
     # initial guess
-    config['u_0'] = 0.5*S*ones((N,N))#S*rand(N,N)
+    config['u_0'] = 0.5*S*ones((N,N),float64)#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/JWST_in.py

@@ -70,7 +70,7 @@ new_config = {
 # able to control (without too much damage)'''
 
 # Algorithm:
-    'algorithm' : 'AP', #'Accelerated_AP_product_space';  
+    'algorithm' : 'RAAR', #'Accelerated_AP_product_space';  
     'numruns' : 1, # the only time this parameter will
 # be different than 1 is when we are
 # benchmarking...not something a normal user
@@ -88,7 +88,7 @@ new_config = {
 
 
     #maximum number of iterations and tolerances
-    'MAXIT' : 1,
+    'MAXIT' : 500,
     'TOL' : 1e-8,
 
     # relaxaton parameters in RAAR, HPR and HAAR

proxtoolbox/Problems/Phase/phase.py

@@ -213,6 +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'])
 
     
     def _postsolve(self):
@@ -270,10 +271,12 @@ class Phase(Problem):
                 else:
                     print("No file available to compare to.")
                     return
-            elif self.config['algorithm'] == 'AP':
+            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')
+            elif self.config['algorithm'] == 'AP' and self.config['constraint'] == 'nonnegative and support':
                         f = h5py.File('Phase_test_data/tasse_u1_ap_' + str(self.config['MAXIT']) + '.mat')
 
-            u1 = f['u1'].value.view(np.float)
+            u1 = f['u1'].value.view(np.float64)
 
         u1 =np.array(u1)
         u1 = u1.T

proxtoolbox/ProxOperators/__init__.py

@@ -9,4 +9,4 @@ The "ProxOperators"-module contains all proximal operators that are already impl
 
 from .proxoperators import *
 
-__all__ = ["P_diag","P_parallel","magproj", "Approx_P_JWST_Poisson", "P_amp", "P_SP", "Approx_PM_Gaussian", "Approx_PM_Poisson"]
+__all__ = ["P_diag","P_parallel","magproj", "Approx_P_JWST_Poisson", "P_amp", "P_SP", "Approx_PM_Gaussian", "Approx_PM_Poisson", "P_S"]

proxtoolbox/ProxOperators/proxoperators.py

@@ -9,9 +9,10 @@ The "ProxOperators"-module contains various specific operators that do the actua
 
 import numpy as np
 from numpy import conj, dot, empty, ones, sqrt, sum, zeros, exp, nonzero, log, tile, shape, real
+#from pyfftw.interfaces.scipy_fftpack import fft2, ifft2
 from numpy.fft import fft2, ifft2
 
-__all__ = ["P_diag","P_parallel","magproj", "Approx_P_JWST_Poisson", "P_amp", "P_SP", "Approx_PM_Gaussian", "Approx_PM_Poisson"]
+__all__ = ["P_diag","P_parallel","magproj", "Approx_P_JWST_Poisson", "P_amp", "P_SP", "Approx_PM_Gaussian", "Approx_PM_Poisson", "P_S"]
 
 
 
@@ -45,44 +46,48 @@ class ProxOperator:
         raise NotImplementedError("This is just an abstract interface")
 
 
-
+@profile
 def magproj(u, constr):
-	"""
-	Projection operator onto a magnitude constraint
-
-	Parameters
-	----------
-	u : array_like - The function to be projected onto constr (can be complex)
-	constr : array_like - A nonnegative array that is the magnitude constraint
-
-	Returns
-	-------
-	array_like - The projection
-	"""
-
-	"""
-	# naive implementation:
-	mod_u = sqrt(conj(u) * u)
-	return constr*(u/mod_u)
-
-	"""
+    """
+    Projection operator onto a magnitude constraint
 
-	"""  
-	Inexact, but stable implementation of magnitude projection.  
-	See LukeBurkeLyon, SIREV 2002
-	"""
+    Parameters
+    ----------
+    u : array_like - The function to be projected onto constr (can be complex)
+    constr : array_like - A nonnegative array that is the magnitude constraint
 
-	eps = 1e-30
+    Returns
+    -------
+    array_like - The projection
+    """
+    """
+    # naive implementation: should work now, roughly as fast as below
+    mod_u = sqrt(u.real**2+u.imag**2)
+    with np.errstate(divide='ignore', invalid='ignore'):
+        proj = constr/mod_u
+        proj[np.isnan(proj)] = 0 #then mod_u=0 and constr=0
+    proj = proj*u
+    index_inf = np.isinf(proj)
+    proj[index_inf] = constr[index_inf] #then mod_u=0 and constr!=0
+    return  proj
+    """
 
-	modsq_u = conj(u) * u
-	denom = modsq_u+eps
-	denom2 = sqrt(denom)
-	r_eps = (modsq_u/denom2) - constr
-	dr_eps = (denom+eps)/(denom*denom2)
+    """  
+    Inexact, but stable implementation of magnitude projection.  
+    See LukeBurkeLyon, SIREV 2002
+    """
+    
+    eps = 3e-30
 
-	return (1 - (dr_eps*r_eps)) * u
+    modsq_u = u.real**2+u.imag**2 #beaware: for U * conj(U) subsequent calculations are much slower since complex (more than double computation time)
+    denom = modsq_u+eps
+    denom2 = sqrt(denom)
+    r_eps = (modsq_u/denom2) - constr
+    dr_eps = (denom+eps)/(denom*denom2)
 
+    return (1 - (dr_eps*r_eps)) * u
 
+    
 class P_diag(ProxOperator):
     """
     Projection onto the diagonal of a product space
@@ -256,12 +261,13 @@ class Approx_P_JWST_Poisson(ProxOperator):
         self.data = config['data'];
         self.product_space_dimension = config['product_space_dimension'];
         self.abs_illumination = config['abs_illumination'];
+    
 
         #calculate the following once (not every iteration) for better performance (speedup ~20# for 500 Iterations JWST)
         self.exp_illu = exp(1j*self.illumination_phase)*tile(self.indicator_ampl[...,None],(1,1,self.product_space_dimension-1))
         self.exp_illu_minus = exp(-1j*self.illumination_phase)*tile(self.indicator_ampl[...,None],(1,1,self.product_space_dimension-1))
 
-    
+    @profile    
     def work(self,u):
         TOL2 = self.TOL2;
         epsilon = self.epsilon;
@@ -280,7 +286,7 @@ class Approx_P_JWST_Poisson(ProxOperator):
 
         for j in range(product_space_dimension-1):
             U = fft2( exp_illu[:,:,j] * u[:,:,j]);
-            U_sq = U * conj(U);
+            U_sq = U.real**2+U.imag**2 #slightly faster real(U * conj(U)) and be aware: for U * conj(U) subsequent calculations much slower since complex
             tmp = U_sq/data_sq[:,:,j];
             mask = ones(tmp.shape, np.bool);
             mask[data_zeros[j]] = 0;
@@ -288,8 +294,8 @@ class Approx_P_JWST_Poisson(ProxOperator):
             U_sq[mask]=0;
             IU= tmp==0;
             tmp[IU]=1;
-            tmp=log(tmp);
-            hU = sum(sum(U_sq *tmp + data_sq[:,:,j] - U_sq));
+            tmp=log(tmp)
+            hU = sum(U_sq *tmp + data_sq[:,:,j] - U_sq)
             if hU>=epsilon+TOL2:
                 U0 = magproj(U,data[:,:,j]); #argument order changed compared to matlab implementation!!!
                 u_new[:,:,j] = exp_illu_minus[:,:,j] *ifft2(U0);
@@ -358,7 +364,7 @@ class P_SP(ProxOperator):
         self.support_idx = config['support_idx']
 
     def work(self,u):
-        p_SP=zeros(shape(u))
+        p_SP=zeros(shape(u), np.float64)
         p_SP[self.support_idx] = np.maximum(real(u[self.support_idx]),0)
         return p_SP
 
@@ -438,6 +444,7 @@ class Approx_PM_Poisson(ProxOperator):
         self.Ib = config['data_zeros']
         self.epsilon = config['data_ball']
 
+    
     def work(self,u):
         TOL2 = self.TOL2
         M = self.M
@@ -459,6 +466,15 @@ class Approx_PM_Poisson(ProxOperator):
         else:
             return u
 
+class P_S(ProxOperator):
+
+    def __init__(self,config):
+        self.support_idx = config['support_idx']
+
+    def work(self,u):
+        p_S=zeros(shape(u),dtype = np.complex128)
+        p_S[self.support_idx] = u[self.support_idx]
+        return p_S