Updated RAAR algortihm (should now work identically to matlab version), small...

Updated RAAR algortihm (should now work identically to matlab version), small changes in other files

proxtoolbox/Algorithms/RAAR.py

@@ -5,10 +5,9 @@ Created on Mon Dec 14 12:48:26 2015
 @author: rebecca
 """
 
-from math import exp, sqrt
-from numpy import zeros
+from numpy import zeros, angle, trace, exp, sqrt
 
-from scipy.linalg import norm
+from numpy.linalg import norm
 
 from .algorithms import Algorithm
 
@@ -27,12 +26,12 @@ class RAAR(Algorithm):
             
                 proxoperators: 2 ProxOperators
                     Tuple of ProxOperators (the class, no instance)
-                beta0: number
+                beta_0: number
                     Starting relaxation parmater
                 beta_max: number
                     Maximum relaxation parameter
                 beta_switch: int
-                    Iteration at which beta moves from beta0 -> beta_max
+                    Iteration at which beta moves from beta_0 -> beta_max
                 normM: number
                     ?
                 Nx: int
@@ -44,15 +43,21 @@ class RAAR(Algorithm):
                 dim: int
                     Size of the product space
         """
-        self.proj1 = config['proxoperators'][0](config)
-        self.proj2 = config['proxoperators'][1](config)
+
+        self.prox1 = config['proxoperators'][0](config)
+        self.prox2 = config['proxoperators'][1](config)
         self.normM = config['normM']
-        self.beta0 = config['beta0']
+        self.beta_0 = config['beta_0']
         self.beta_max = config['beta_max']
         self.beta_switch = config['beta_switch']
         self.Nx = config['Nx']; self.Ny = config['Ny']; self.Nz = config['Nz'];
-        self.dim = config['dim']
-        self.iters = 0
+        self.product_space_dimension = config['product_space_dimension']
+        self.iter = 0
+
+        if 'truth' in config:
+            self.truth = config['truth']
+            self.truth_dim = config['truth_dim']
+            self.norm_truth = config['norm_truth']
 
     def run(self, u, tol, maxiter):
         """
@@ -62,65 +67,110 @@ class RAAR(Algorithm):
         ##### PREPROCESSING
         normM = self.normM
         
-        beta = self.beta0
-        iters = self.iters
+        beta = self.beta_0
+        iter = self.iter
+        if u.ndim < 3:
+            p = 1
+            q = 1
+        elif u.ndim == 3:
+            p = u.shape[2]
+            q = 1
+        else:
+            p = u.shape[2]
+            q = u.shape[3]
+
         change = zeros(maxiter+1,dtype=u.dtype)
         change[0] = 999
         gap = change.copy()
-        
-        tmp1 = 2*self.proj2.work(u) - u
+        shadow_change = change.copy()
+        if  hasattr(self, 'truth'):
+            Relerrs = change.copy()
+
+        tmp1 = 2*self.prox2.work(u) - u
+        shadow = self.prox1.work(u)
         
         
         ##### LOOP
-        while iters < maxiter and change[iters] >= tol:
-            tmp = exp((-iters/self.beta_switch)**3);
-            beta = (tmp*self.beta0) + ((1-tmp)*self.beta_max);
-            iters += 1;
+        while iter < maxiter and shadow_change[iter] >= tol:
+
+            tmp = exp((-iter/self.beta_switch)**3);
+            beta = (tmp*self.beta_0) + ((1-tmp)*self.beta_max)
+            iter += 1;
             
-            tmp3 = self.proj1.work(tmp1);
-            tmp_u = ((beta*(2*tmp3-tmp1)) + ((1-beta)*tmp1) + u)/2;
-            tmp2 = self.proj2.work(tmp_u);
+            tmp3 = self.prox1.work(tmp1)
+            tmp_u = ((beta*(2*tmp3-tmp1)) + ((1-beta)*tmp1) + u)/2
+            tmp2 = self.prox2.work(tmp_u)
             
-            tmp3 = self.proj1.work(tmp2);
+            tmp3 = self.prox1.work(tmp2)
             
-            tmp_change = 0; tmp_gap = 0;
-            if self.Ny == 1 or self.Nx == 1:
-                tmp_change = (norm(u-tmp_u,'fro')/normM)**2;
-                tmp_gap = (norm(tmp3-tmp2,'fro')/normM)**2;
-            elif self.Nz == 1:
-                for j in range(self.dim):
-                    tmp_change += (norm(u[:,:,j]-tmp_u[:,:,j],'fro')/normM)**2;
-                    tmp_gap += (norm(tmp3[:,:,j]-tmp2[:,:,j],'fro')/normM)**2;
+            tmp_change = 0; tmp_gap = 0; tmp_shadow=0;
+
+            if p==1 and q==1:
+              tmp_change= (norm(u-tmp_u, 'fro')/normM)**2
+              tmp_shadow = (norm(tmp3-shadow,'fro')/normM)**2
+              tmp_gap = (norm(tmp3-tmp2,'fro')/normM)**2
+
+              if hasattr(self, 'truth'):
+                  if self.truth_dim[0] == 1:
+                      z=tmp3[0,:]
+                  elif slef.truth_dim[1] == 1:
+                      z=tmp3[:,0]
+                  else:
+                      z=tmp3;
+                  Relerrs[iter] = norm(self.truth - exp(-1j*angle(trace(self.truth.T*z))) * z, 'fro')/self.norm_truth
+
+            elif q==1:
+              for j in range(self.product_space_dimension):
+                  tmp_change= tmp_change+ (norm(u[:,:,j]-tmp_u[:,:,j], 'fro')/normM)**2
+                  # compute (||P_SP_Mx-P_Mx||/normM)^2:
+                  tmp_gap = tmp_gap+(norm(tmp3[:,:,j]-tmp2[:,:,j],'fro')/normM)**2
+                  tmp_shadow = tmp_shadow+(norm(tmp3[:,:,j]-shadow[:,:,j],'fro')/normM)**2
+
+              if hasattr(self, 'truth'):
+                 z=tmp3[:,:,0]
+                 Relerrs[iter] = norm(self.truth - exp(-1j*angle(trace(self.truth.T*z))) * z, 'fro')/self.norm_truth
+
             else:
-                for j in range(self.dim):
-                    for k in range(self.Nz):
-                        tmp_change += (norm(u[:,:,k,j]-tmp_u[:,:,k,j],'fro')/normM)**2;
-                        tmp_gap += (norm(tmp3[:,:,k,j]-tmp2[:,:,k,j],'fro')/normM)**2;
-            
-            change[iters] = sqrt(tmp_change);
-            gap[iters] = sqrt(tmp_gap);
+              for j in range(self.product_space_dimension):
+                  for k in range(Nz):
+                      tmp_change= tmp_change+ (norm(u[:,:,k,j]-tmp_u[:,:,k,j], 'fro')/normM)**2
+                      # compute (||P_Sx-P_Mx||/normM)^2:
+                      tmp_gap = tmp_gap+(norm(tmp3[:,:,k,j]-tmp2[:,:,k,j],'fro')/(normM))**2
+                      tmp_shadow = tmp_shadow+(norm(tmp3[:,:,k,j]-shadow[:,:,k,j],'fro')/(normM))**2
+   
+            change[iter] = sqrt(tmp_change)
+            gap[iter] = sqrt(tmp_gap)
+            shadow_change[iter] = sqrt(tmp_shadow) # this is the Euclidean norm of the gap to
+            # the unregularized set.  To monitor the Euclidean norm of the gap to the
+            # regularized set is expensive to calculate, so we use this surrogate.
+            # Since the stopping criteria is on the change in the iterates, this
+            # does not matter.
+            # graphics
             
-            u = tmp_u;
-            tmp1 = (2*tmp2) - tmp_u;
+            u = tmp_u
+            tmp1 = (2*tmp2) - tmp_u
             
         
         ##### POSTPROCESSING
         u = tmp2;
-        tmp = self.proj1.work(u);
-        tmp2 = self.proj2.work(u);
-        if self.Ny == 1:
-            u1 = tmp[:,1];
-            u2 = tmp2[:,1];
-        elif self.Nx == 1:
-            u1 = tmp[1,:];
-            u2 = tmp2[1,:];
+        tmp = self.prox1.work(u);
+        tmp2 = self.prox2.work(u);
+        if self.Nx == 1:
+            u1 = tmp[:,0];
+            u2 = tmp2[:,0];
+        elif self.Ny == 1:
+            u1 = tmp[0,:];
+            u2 = tmp2[0,:];
         elif self.Nz == 1:
-            u1 = tmp[:,:,1];
-            u2 = tmp2[:,:,1];
+            u1 = tmp[:,:,0];
+            u2 = tmp2[:,:,0];
         else:
             u1 = tmp;
             u2 = tmp2;
-        change = change[1:iters+1];
-        gap = gap[1:iters+1];
+        change = change[1:iter+1];
+        gap = gap[1:iter+1];
+        shadow_change = shadow_change[1:iter+1]
+        print(iter)
+
         
-        return u1, u2, iters, change, gap
+        return u1, u2, iter, change, gap#, shadow_change

proxtoolbox/Problems/Phase/CDI_data_processor.py

@@ -108,7 +108,7 @@ def CDI_data_processor(config):
     config['supp_phase'] = []
 
     # initial guess
-    config['u_0'] = S*rand(N)
+    config['u_0'] = 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

@@ -92,7 +92,7 @@ new_config = {
     'TOL' : 1e-8,
 
     # relaxaton parameters in RAAR, HPR and HAAR
-    'beta0' : 1.0,                # starting relaxation prameter (only used with
+    'beta_0' : 1.0,                # starting relaxation prameter (only used with
     # HAAR, HPR and RAAR)
     'beta_max' : 1.0,             # maximum relaxation prameter (only used with
     # HAAR, RAAR, and HPR)

proxtoolbox/Problems/Phase/phase.py

@@ -130,7 +130,8 @@ class Phase(Problem):
         self.config['proxoperators'] = [];
 
         for prox in proxoperators:
-            self.config['proxoperators'].append(getattr(ProxOperators, prox))
+            if prox != '':
+                self.config['proxoperators'].append(getattr(ProxOperators, prox))
 
         # input.Proj1_input.F=F;  % is it any more expensive to pass everything
         # into the projectors rather than just a selection of data and
@@ -176,6 +177,8 @@ class Phase(Problem):
         # 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']);
+        elif self.config['product_space_dimension'] == 1:
+            tmp_gap = (norm(u_1-u_2)/square(self.config['norm_rt_data']));
         else:
             tmp_gap=0;
             for j in range(self.config['product_space_dimension']):

proxtoolbox/Problems/Phase/tasse_raar_in.py

@@ -54,7 +54,7 @@ new_config = {
 ## able to control (without too much damage)
 
 ## Algorithm:
-'method' : 'RAAR',  # used to be 'Projection', 
+'algorithm' : 'RAAR',  # used to be 'Projection', 
 'numruns' : 1, # the only time this parameter will
 # be different than 1 is when we are
 # benchmarking...not something a normal user

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"]
+__all__ = ["P_diag","P_parallel","magproj", "Approx_P_JWST_Poisson", "P_amp", "P_SP", "Approx_PM_Gaussian", "Approx_PM_Poisson"]

proxtoolbox/ProxOperators/proxoperators.py

@@ -8,10 +8,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
+from numpy import conj, dot, empty, ones, sqrt, sum, zeros, exp, nonzero, log, tile, shape, real
 from numpy.fft import fft2, ifft2
 
-__all__ = ["P_diag","P_parallel","magproj", "Approx_P_JWST_Poisson", "P_amp", "P_SP"]
+__all__ = ["P_diag","P_parallel","magproj", "Approx_P_JWST_Poisson", "P_amp", "P_SP", "Approx_PM_Gaussian", "Approx_PM_Poisson"]
 
 
 
@@ -356,8 +356,111 @@ class P_amp(ProxOperator):
 class P_SP(ProxOperator):
 
     def __init__(self,config):
-        self.support_idx = config['support_idx'];
+        self.support_idx = config['support_idx']
 
     def work(self,u):
-        p_SP=zeros(shape(u));
-        p_SP[self.support_idx] = max(real(u[self.support_idx]),0);
+        p_SP=zeros(shape(u))
+        p_SP[self.support_idx] = np.max(real(u[self.support_idx]),0)
+        return p_SP
+
+#                      Approx_PM_Poisson.m
+#             written on Feb. 18, 2011 by 
+#                   Russell Luke
+#   Inst. Fuer Numerische und Angewandte Mathematik
+#                Universitaet Gottingen
+#
+# DESCRIPTION:  Projection subroutine for projecting onto Fourier
+#               magnitude constraints
+#
+# INPUT:        Func_params = a data structure with .data = nonegative real FOURIER DOMAIN CONSTRAINT
+#                             .data_ball is the regularization parameter described in 
+#                             D. R. Luke, Nonlinear Analysis 75 (2012) 1531–1546.
+#                            .TOL2 is an extra tolerance. 
+#               u = function in the physical domain to be projected
+#
+# OUTPUT:       p_M    = the projection IN THE PHYSICAL (time) DOMAIN
+#               phat_M = projection IN THE FOURIER DOMAIN
+#               
+# USAGE: u_epsilon = Approx_PM_Gaussian(Func_params,u)
+#
+# Nonstandard Matlab function calls:  MagProj
+
+class Approx_PM_Gaussian(ProxOperator):
+
+    def __init__(self,config):
+        self.TOL2 = config['TOL2']
+        self.b = config['data'] * config['data']
+        self.epsilon = config['data_ball']
+        self.data = config['data']
+
+    def work(self,u):
+        TOL2 = self.TOL2;
+        b = self.b
+        epsilon = self.epsilon;
+        U = fft2(u)
+        U0 = magproj(U, self.data)
+        U0_sq = U0 * conj(U0)
+        tmp = U0_sq - b
+        h=sum(sum(tmp * conj(tmp)))
+        if h>=epsilon+TOL2:
+            return ifft(U0)
+        else:
+            return u
+
+#                      Approx_PM_Poisson.m
+#             written on Feb. 18, 2011 by 
+#                   Russell Luke
+#   Inst. Fuer Numerische und Angewandte Mathematik
+#                Universitaet Gottingen
+#
+# DESCRIPTION:  Projection subroutine for projecting onto Fourier
+#               magnitude constraints
+#
+# INPUT:       Func_params = a data structure with 
+#                             .data = nonegative real FOURIER DOMAIN CONSTRAINT
+#                             .data_sq = the elementwise square of .data
+#                             .data_ball is the regularization parameter described in 
+#                             D. R. Luke, Nonlinear Analysis 75 (2012) 1531–1546.
+#                            .TOL2 is an extra tolerance. 
+#               u = function in the physical domain to be projected
+#
+# OUTPUT:       
+#               
+# USAGE: u_epsilon = P_M(Func_params,u)
+#
+# Nonstandard Matlab function calls:  MagProj
+
+class Approx_PM_Poisson(ProxOperator):
+
+    def __init__(self,config):
+        self.TOL2 = config['TOL2']
+        self.M = config['data']
+        self.b = config['data_sq']
+        self.Ib = config['data_zeros']
+        self.epsilon = config['data_ball']
+
+    def work(self,u):
+        TOL2 = self.TOL2
+        M = self.M
+        b = self.b
+        Ib = self.Ib
+        epsilon = self.epsilon
+
+        U = ifft2(u)
+        U_sq = U * conj(U)
+        tmp = U_sq/b; tmp[Ib] = 1
+        U_sq[Ib]=0
+        IU = tmp==0
+        tmp[IU]=1
+        tmp=log(tmp)
+        hU = sum(sum(U_sq * tmp + b - U_sq))
+        if hU>=epsilon+TOL2:
+            U0 = magproj(U,M)
+            return ifft2(U0)
+        else:
+            return u
+
+        
+
+
+