small changes

proxtoolbox/Algorithms/SimpleAlgortihm.py

@@ -3,11 +3,15 @@
 
 from .algorithms import Algorithm
 
-from numpy import zeros, angle, trace, exp, sqrt
-
+from numpy import zeros, angle, trace, exp, sqrt, sum
 from numpy.linalg import norm
 
 
+def phase_offset_compensated_norm(arr1, arr2, norm_factor=1, norm_type='fro'):
+    phase_offset = angle(sum(arr1 * arr2.conj()))
+    return norm(arr1 - arr2 * exp(1j * phase_offset), norm_type) / norm_factor
+
+
 class SimpleAlgorithm:
 
     def __init__(self, config):
@@ -41,8 +45,8 @@ class SimpleAlgorithm:
         self.prox1 = config['proxoperators'][0](config)
         self.prox2 = config['proxoperators'][1](config)
         self.norm_data = config['norm_data']
-        self.Nx = config['Nx'];
-        self.Ny = config['Ny'];
+        self.Nx = config['Nx']
+        self.Ny = config['Ny']
         self.Nz = config['Nz']
         self.product_space_dimension = config['product_space_dimension']
         self.iter = 0
@@ -58,6 +62,14 @@ class SimpleAlgorithm:
         else:
             self.diagnostic = False
 
+    def evaluate(self, u):
+        """
+        Method to override in specific algorithm subclass
+        :param u: old iterate
+        :return: new iterate
+        """
+        return u
+
     def run(self, u, tol, maxiter):
         """
         Runs the algorithm for the specified input data
@@ -117,14 +129,18 @@ class SimpleAlgorithm:
             tmp_shadow = 0
 
             if p == 1 and q == 1:
-                tmp_change = (norm(u - u_new, 'fro') / norm_data) ** 2
+                # tmp_change = (norm(u - u_new, 'fro') / norm_data) ** 2
+                tmp_change = phase_offset_compensated_norm(u, u_new, norm_type='fro', norm_factor=norm_data) ** 2
+
                 if 'diagnostic' in self.config:
                     # For Douglas-Rachford,in general it is appropriate to monitor the
                     # SHADOWS of the iterates, since in the convex case these converge
                     # even for beta=1.
                     # (see Bauschke-Combettes-Luke, J. Approx. Theory, 2004)
-                    tmp_shadow = (norm(u2 - shadow, 'fro') / norm_data) ** 2
-                    tmp_gap = (norm(u1 - u2, 'fro') / norm_data) ** 2
+                    # tmp_shadow = (norm(u2 - shadow, 'fro') / norm_data) ** 2
+                    # tmp_gap = (norm(u1 - u2, 'fro') / norm_data) ** 2
+                    tmp_shadow = phase_offset_compensated_norm(u2, shadow, norm_factor=norm_data, norm_type='fro') ** 2
+                    tmp_gap = phase_offset_compensated_norm(u1, u2, norm_factor=norm_data, norm_type='fro') ** 2
 
                     if hasattr(self, 'truth'):
                         if self.truth_dim[0] == 1:
@@ -133,8 +149,10 @@ class SimpleAlgorithm:
                             z = u1[:, 0]
                         else:
                             z = u1
-                        Relerrs[iter] = norm(self.truth - exp(-1j * angle(trace(self.truth.T * z))) * z,
-                                             'fro') / self.norm_truth
+                        # Relerrs[iter] = norm(self.truth - exp(-1j * angle(trace(self.truth.T * z))) * z,
+                        #                      'fro') / self.norm_truth
+                        Relerrs[iter] = phase_offset_compensated_norm(self.truth, z, norm_factor=self.norm_truth,
+                                                                      norm_type='fro')
 
             elif q == 1:
                 for j in range(self.product_space_dimension):
@@ -160,7 +178,7 @@ class SimpleAlgorithm:
                             # compute (||P_Sx-P_Mx||/norm_data)^2:
                             tmp_gap = tmp_gap + (norm(u1[:, :, k, j] - u2[:, :, k, j], 'fro') / (norm_data)) ** 2
                             tmp_shadow = tmp_shadow + (
-                                        norm(u2[:, :, k, j] - shadow[:, :, k, j], 'fro') / (norm_data)) ** 2
+                                    norm(u2[:, :, k, j] - shadow[:, :, k, j], 'fro') / (norm_data)) ** 2
                 if hasattr(self, 'truth') and (j == 0):
                     Relerrs[iter] = Relerrs[iter] + norm(
                         self.truth - exp(-1j * angle(trace(self.truth.T * u1[:, :, k, 1]))) * u1[:, :, k, 1],

proxtoolbox/Problems/OrbitalTomog/coronene_test.py → proxtoolbox/Problems/OrbitalTomog/coronene_config.py

-## This is the input file that the user sees/modifies.  It should be simple, 
-## avoid jargon or acronyms, and should be a model for a menu-driven GUI
+# This is the input file that the user sees/modifies.  It should be simple,
+# avoid jargon or acronyms, and should be a model for a menu-driven GUI
 
 new_config = {
 
@@ -17,6 +17,7 @@ new_config = {
     # What type of constraints do we have?
     # Options are: 'support only', 'real and support', 'nonnegative and support',
     #              'amplitude only', 'sparse real', 'sparse complex', and 'hybrid'
+    #              'symmetric sparse real', 'symmetric sparse complex'
     'constraint': 'real and support',
 
     # What type of measurements are we working with?
@@ -62,14 +63,14 @@ new_config = {
     # benchmarking...not something a normal user
     # would be doing.
 
-    ## The following are parameters specific to RAAR, HPR, and HAAR that the
-    ## user should be able to set/modify.  Surely
-    ## there will be other algorithm specific parameters that a user might
-    ## want to play with.  Don't know how best
-    ## to do this.  Thinking of a GUI interface, we could hard code all the
-    ##  parameters the user might encounter and have the menu options change
-    ## depending on the value of the 'method field.
-    ## do different things depending on the chosen algorithm:
+    # # The following are parameters specific to RAAR, HPR, and HAAR that the
+    # # user should be able to set/modify.  Surely
+    # # there will be other algorithm specific parameters that a user might
+    # # want to play with.  Don't know how best
+    # # to do this.  Thinking of a GUI interface, we could hard code all the
+    # #  parameters the user might encounter and have the menu options change
+    # # depending on the value of the 'method field.
+    # # do different things depending on the chosen algorithm:
     # if(strcmp('method,'RAAR')||strcmp('method,'AP')||...
     #        strcmp('method,'HPR')||strcmp('method,'HAAR'))
     # the following just points this driver to a patch that communicates 
@@ -118,9 +119,9 @@ new_config = {
     #    'N_CG' : 70,
     # end
 
-    # ##==========================================
-    # ## parameters for plotting and diagnostics
-    # ##==========================================
+    # # ==========================================
+    # # parameters for plotting and diagnostics
+    # # ==========================================
     # 'plotWhat.n1=2,
     # 'plotWhat.n2=3,
     # 'plotWhat.plots' : 'PYWpyw',
@@ -135,9 +136,9 @@ new_config = {
     #                             # the graphics.  Otherwise, the user
     #                             # can write their own plotting subroutine
     #
-    #==========================================
-    # parameters for plotting and diagnostics
-    #==========================================
+    # # ==========================================
+    # # parameters for plotting and diagnostics
+    # # ==========================================
     'diagnostic': True,
     'verbose': 1,  # options are 0 or 1
     'graphics': 1,  # whether or not to display figures, options are 0 or 1.
@@ -153,9 +154,9 @@ new_config = {
 
 }
 
-#======================================================================
+# =====================================================================
 #  Technical/software specific parameters
-#======================================================================
+# =====================================================================
 # Given the parameter values above, the following technical/algorithmic
 # parameters are automatically set.  The user does not need to know
 # about these details, and so probably these parameters should be set in

proxtoolbox/Problems/OrbitalTomog/orbitaltomog_data_processor.py

@@ -2,11 +2,15 @@ import numpy as np
 import numpy.fft as fft
 from scipy.io import loadmat
 import matplotlib.pyplot as plt
+import imageio
 
 
 def data_processor(config):
-    inp_data = loadmat('../../../InputData/OrbitalTomog/' + config['data_filename'])
-    inp = inp_data["I2"]
+    if config['data_filename'][-4:] == '.mat':
+        inp_data = loadmat('../../../InputData/OrbitalTomog/' + config['data_filename'])
+        inp = inp_data["I2"]
+    else:
+        inp = imageio.imread('../../../InputData/OrbitalTomog/' + config['data_filename'])
     ny, nx = inp.shape
       
     config['data'] = abs(inp)
@@ -15,16 +19,14 @@ def data_processor(config):
     # Keep the same resolution?
     config['Ny'], config['Nx'] = ny, nx
 
-    # Autokorrelation
+    # Autocorrelation
     config['threshold for support'] = 0.2
     autocorrelation = abs(fft.fftshift(fft.ifft2(fft.ifftshift(inp))))
     config['support'] = (autocorrelation > config['threshold for support'] * np.amax(autocorrelation)).astype(np.uint) # [0,0,1,1,0,..]
     config['abs_illumination'] = np.ones_like(config['support'])
-    # Anfangsbedingungen
+    # Initial guess
     ph_init = 2 * np.pi * np.random.rand(ny, nx)
-    # ph_init = np.angle(np.fft.fft2(ph_init));
     config['u_0'] = inp * np.exp(1j * ph_init)
-#    previous = fft.fftshift(fft.ifft2(fft.ifftshift(config['u_0'])))
 
     if config['dataprocessor_plotting']:
         plt.figure(figsize=(15, 4))

setup.py

@@ -59,6 +59,7 @@ setup(
         'Programming Language :: Python :: 3.2',
         'Programming Language :: Python :: 3.3',
         'Programming Language :: Python :: 3.4',
+        'Programming Language :: Python :: 3.7',
     ],
 
     # What does your project relate to?
@@ -77,6 +78,7 @@ setup(
     # requirements files see:
     # https://packaging.python.org/en/latest/requirements.html
     #install_requires=['peppercorn'],
+    # # June 2019 dependencies: numpy, scipy, matplotlib, imageio
 
     # List additional groups of dependencies here (e.g. development
     # dependencies). You can install these using the following syntax,