degenerate orbits experiment running + converging, however yielding physically...

degenerate orbits experiment running + converging, however yielding physically incorrect solutions, which are mathematically not unique.

proxtoolbox/algorithms/algorithm.py

-
-from numpy import zeros, angle, trace, exp, sqrt, sum
+from numpy import angle, exp, sum
 from numpy.linalg import norm
+try:
+    from tqdm import tqdm, tqdm_notebook
+    load_tqdm = True
+except ModuleNotFoundError:
+    load_tqdm = False
+    pass
+
 
 class Algorithm:
     """
@@ -17,7 +23,7 @@ class Algorithm:
             if proxClass != '':
                 prox = proxClass(experiment)
                 self.proxOperators.append(prox)
-        #for convenience
+        # for convenience
         self.prox1 = self.proxOperators[0]
         self.prox2 = self.proxOperators[1]
 
@@ -40,6 +46,11 @@ class Algorithm:
         self.iterateMonitor = iterateMonitor
         self.accelerator = accelerator
 
+        if hasattr(experiment, 'progressbar'):
+            self.progressbar = experiment.progressbar
+        else:
+            self.progressbar = None
+
         # temporaries for current iteration
         self.iter = 0
         self.u = None
@@ -48,14 +59,13 @@ class Algorithm:
         # for debugging
         self.debug = experiment.debug
 
-
     def preprocess(self):
         """
         The default implementation calls the iterate monitor's
         preprocess() method which initializes the arrays that 
         will be used to collect statistics.      
         """
-        self.iterateMonitor.preprocess(self)        
+        self.iterateMonitor.preprocess(self)
 
     def stoppingCondition(self):
         """
@@ -149,15 +159,26 @@ class Algorithm:
         self.u = u
         self.preprocess()
         self.displayProgress()
-        while (self.stoppingCondition()):
+
+        if self.progressbar == 'tqdm_notebook':
+            pbar = tqdm_notebook(total=self.maxIter)
+        elif self.progressbar == 'tqdm':
+            pbar = tqdm(total=self.maxIter)
+        else:
+            pbar = None
+
+        while self.stoppingCondition():
             self.iter += 1
+            if pbar is not None:
+                pbar.update()
             self.u_new = self.evaluate(self.u)
             self.updateStatistics()
             self.displayProgress()
             self.u = self.u_new
+        if pbar is not None:
+            pbar.close()
         return self.postprocess()
 
-
     def getDescription(self):
         """
         Function returning a string giving the algorithm name
@@ -176,7 +197,7 @@ class Algorithm:
         desc = self.getDescriptionHelper()
         return desc
 
-    def getDescriptionHelper(self, param_name = None, param_0 = None, param_max = None):
+    def getDescriptionHelper(self, param_name=None, param_0=None, param_max=None):
         """
         Generate a string giving a short description of
         the algorithm (the name of its class name and optionally,
@@ -191,9 +212,9 @@ class Algorithm:
         """
         desc = type(self).__name__
         if param_name is not None and param_0 is not None \
-           and param_max is not None:
+                and param_max is not None:
             desc += ", $" + param_name + "$ = " + str(param_0) \
-                     + " to " + str(param_max)
+                    + " to " + str(param_max)
         return desc
 
     def computeRelaxation(self):
@@ -206,10 +227,10 @@ class Algorithm:
         lmbda: Float64 
             Relaxation parameter.
         """
-        iter = self.iter + 1 # add 1 to conform with matlab version
+        iteration = self.iter + 1  # add 1 to conform with matlab version
         # gradually changes lambda from lambda_0 into lambda_max
-        a = -iter/self.lambda_switch
-        a *= a*a  # compute a^3
+        a = -iteration / self.lambda_switch
+        a *= a * a  # compute a^3
         s = exp(a)
-        lmbda = s*self.lambda_0 + (1-s)*self.lambda_max
+        lmbda = s * self.lambda_0 + (1 - s) * self.lambda_max
         return lmbda

proxtoolbox/experiments/orbitaltomography/degenerate_orbits.py

@@ -17,8 +17,8 @@ class DegenerateOrbital(PlanarMolecule):
             'degeneracy': 2,  # Number of degenerate states to reconstruct
             'constraint': 'sparse real',
             'sparsity_parameter': 40,
-            'use_sparsity_with_support': False,
-            'threshold_for_support': 0.1,
+            'use_sparsity_with_support': True,
+            'threshold_for_support': 0.05,
             'support_filename': None,
             'Nx': None,
             'Ny': None,
@@ -38,6 +38,7 @@ class DegenerateOrbital(PlanarMolecule):
             'graphics': 1,
             'interpolate_and_zoom': True,
             'debug': True,
+            'progressbar': None
         }
         return defaultParams
 
@@ -76,33 +77,16 @@ class DegenerateOrbital(PlanarMolecule):
         - self.proxOperators
         - self.propagator and self.inverse_propagator
         """
-        # Apply orthonormality constraint first
-        self.proxOperators.append('P_orthonorm')
-
         # Select the right real space operator sparsity-based proxoperators
-        if self.constraint == 'sparse real':
-            self.proxOperators.append('P_Sparsity_real')
-        elif self.constraint == 'sparse complex':
-            self.proxOperators.append('P_Sparsity')
-        elif self.constraint in ['symmetric sparse real', 'sparse symmetric real']:
-            self.proxOperators.append('P_Sparsity_Symmetric_real')
-        elif self.constraint in ['symmetric sparse complex', 'symmetric sparse complex']:
-            self.proxOperators.append('P_Sparsity_Symmetric')
+        self.proxOperators.append('P_Sparsity_real_incoherent')
+
+        # Apply orthonormality constraint
+        self.proxOperators.append('P_orthonorm')
 
         # Modulus proxoperator (normally the second operator)
-        if self.experiment_name == '3D ARPES':
-            self.proxOperators.append('P_M_masked')
-            self.propagator = 'PropagatorFFTn'
-            self.inverse_propagator = 'InvPropagatorFFTn'
-        elif self.experiment_name == 'noisy 2D ARPES':
-            # Apply modulus constraint when the difference is large enough (Approx_ prefix)
-            self.proxOperators.append('Approx_P_M')
-            self.propagator = 'PropagatorFFT2'
-            self.inverse_propagator = 'InvPropagatorFFT2'
-        else:  # No noise case: always apply the modulus constraint
-            self.proxOperators.append('P_M')
-            self.propagator = 'PropagatorFFT2'
-            self.inverse_propagator = 'InvPropagatorFFT2'
+        self.proxOperators.append('P_M_incoherent')
+        self.propagator = 'PropagatorFFT2'
+        self.inverse_propagator = 'InvPropagatorFFT2'
 
         self.nProx = len(self.proxOperators)
 
@@ -128,83 +112,45 @@ class DegenerateOrbital(PlanarMolecule):
         """
         Create basic result plots of the phase retrieval procedure
         """
-        self.output['plots'] = []
+        super(PlanarMolecule, self).show()
+
         self.output['u1'] = self.algorithm.prox1.eval(self.algorithm.u)
         self.output['u2'] = self.algorithm.prox2.eval(self.algorithm.u)
-        # if self.interpolate_and_zoom:
-        #     for key in ["u1", "u2"]:
-        #         center = tuple([s // 2 for s in self.output[key].shape])
-        #         self.output[key] = roll_to_pos(self.output[key], pos=center, move_maximum=True)
-        #         self.output[key] = roll_to_pos(self.output[key], pos=center)
-        #         self.output[key] = fourier_interpolate(self.output[key], 2)
-        #         zmy, zmx = tuple([s // 4 for s in self.output[key].shape])
-        #         self.output[key] = self.output[key][zmy:-zmy, zmx:-zmx]
-        u1 = self.output['u1']
-        u2 = self.output['u2']
-        change = self.output['stats']['changes']
-
-        f1 = self.plot_guess(u1, name='Best approximation: physical constraint satisfied', show=False)
-        f2 = self.plot_guess(u2, name='Best approximation: Fourier constraint satisfied', show=False)
-
-        # f, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(9, 6))
-        # im = ax1.imshow(abs(u1), cmap='gray')
-        # f.colorbar(im, ax=ax1)
-        # ax1.set_title('best approximation amplitude: \nphysical constraint satisfied')
-        # ax2.imshow(complex_to_rgb(u1))
-        # ax2.set_title('best approximation phase: \nphysical constraint satisfied')
-        # im = ax3.imshow(abs(u2), cmap='gray')
-        # f.colorbar(im, ax=ax3)
-        # ax3.set_title('best approximation amplitude: \nFourier constraint satisfied')
-        # ax4.imshow(complex_to_rgb(u2))
-        # ax4.set_title('best approximation amplitude: \nFourier constraint satisfied')
-        # f.tight_layout()
-
-        g, ((bx1, bx2), (bx3, bx4)) = plt.subplots(2, 2, figsize=(9, 6))
-        im = bx1.imshow(abs(u1), cmap='gray')
-        g.colorbar(im, ax=bx1)
-        bx1.set_title('best approximation amplitude: \nphysical constraint satisfied')
-        im = bx2.imshow(u1.real, cmap='gray')
-        g.colorbar(im, ax=bx2)
-        bx2.set_title('best approximation phase: \nphysical constraint satisfied')
-        bx3.semilogy(change)
-        bx3.set_xlabel('iteration')
-        bx3.set_ylabel('Change')
-        if 'gaps' in self.output['stats']:
-            gaps = self.output['stats']['gaps']
-            bx4.semilogy(gaps)
-            bx4.set_xlabel('iteration')
-            bx4.set_ylabel('Gap')
-        g.tight_layout()
-
-        h, ax = plt.subplots(1, 3, figsize=(9, 3))
-        ax[0].imshow(self.data)
-        ax[0].set_title("Measured data")
-        prop = self.propagator.eval(self)
-        u_hat = prop.eval(self.algorithm.prox1.eval(self.algorithm.u))
-        ax[1].imshow(abs(u_hat))
-        ax[1].set_title("Predicted measurement intensity")
-        ax[2].imshow(complex_to_rgb(u_hat))
-        ax[2].set_title("Predicted phase (by color)")
-        h.tight_layout()
-        plt.show()
 
-        self.output['plots'].append(f1)
-        self.output['plots'].append(f2)
-        self.output['plots'].append(g)
-        self.output['plots'].append(h)
+        for i, operator in enumerate(self.algorithm.proxOperators):
+            operator_name = self.proxOperators[i].__name__
+            f = self.plot_guess(operator.eval(self.algorithm.u), name="%s satisfied" % operator_name, show=False)
+            self.output['plots'].append(f)
+
+        # f1 = self.plot_guess(self.output['u1'], name='Best approximation: physical constraint satisfied', show=False)
+        # f2 = self.plot_guess(self.output['u2'], name='Best approximation: Fourier constraint satisfied', show=False)
+        # prop = self.propagator(self)
+        # u_hat = prop.eval(self.algorithm.prox1.eval(self.algorithm.u))
+        # h = self.plot_guess(u_hat, show=False, name="Fourier domain measurement projection")
+
+        # self.output['plots'].append(f1)
+        # self.output['plots'].append(f2)
+        # self.output['plots'].append(h)
+        plt.show()
 
     # def saveOutput(self, **kwargs):
     #     super(PlanarMolecule, self).saveOutput(**kwargs)
 
     def plot_guess(self, u, name=None, show=True):
         """"Given a list of fields, plot the individual fields and the combined intensity"""
-        fig, ax = plt.subplots(1, len(u) + 1, figsize=(12, 3.5), num=name)
+        fig, ax = plt.subplots(1, len(u) + 2, figsize=(12, 3), num=name)
         for ii in range(self.degeneracy):
             im = ax[ii].imshow(complex_to_rgb(u[ii]))
-            plt.colorbar(im, ax=ax[ii])
             ax[ii].set_title("Degenerate orbit %d" % ii)
-        ax[-1].imshow(np.sum(abs(u) ** 2, axis=0))
-        ax[-1].set_title("Local density of states")
+        im = ax[-2].imshow(np.sum(abs(u) ** 2, axis=0))
+        ax[-2].set_title("Local density of states")
+        plt.colorbar(im, ax=ax[-2])
+        prop = self.propagator(self)
+        u_hat = prop.eval(u)
+        fourier_intensity = np.sqrt(np.sum(abs(u_hat)**2, axis=0))
+        im = ax[-1].imshow(fourier_intensity)
+        ax[-1].set_title("Fourier domain intensity")
+        plt.colorbar(im, ax=ax[-1])
         plt.tight_layout()
         if show:
             plt.show()

proxtoolbox/experiments/orbitaltomography/orbitExperiment.py

@@ -92,7 +92,11 @@ class OrbitalTomographyExperiment(Experiment):
         """
         self.output['plots'] = []
         convergence_plots, ax = subplots(1, 2, figsize=(7, 3))
-        ax[0].semilogy(self.output['stats']['change'])
+        if 'change' in self.output['stats']:
+            change_key = 'change'
+        else:
+            change_key = 'changes'
+        ax[0].semilogy(self.output['stats'][change_key])
         ax[0].set_title('Change')
         ax[0].set_xlabel('iteration')
         ax[0].set_ylabel('$||x^{2k+2}-x^{2k}||$')

proxtoolbox/experiments/orbitaltomography/planar_molecule.py

@@ -39,6 +39,7 @@ class PlanarMolecule(OrbitalTomographyExperiment):
             'graphics': 1,
             'interpolate_and_zoom': True,
             'debug': True,
+            'progressbar': None  # Valid options: None, 'tqdm' or 'tqdm_notebook'
         }
         return defaultParams
 
@@ -53,6 +54,7 @@ class PlanarMolecule(OrbitalTomographyExperiment):
         self.use_sparsity_with_support = kwargs['use_sparsity_with_support']
         self.threshold_for_support = kwargs['threshold_for_support']  # to determine support from the autocorrelation.
         self.interpolate_and_zoom = kwargs['interpolate_and_zoom']  # For plotting
+        self.progressbar = kwargs['progressbar']
 
         # the following data members are set by loadData(), in addition to those specified in parent classes
         self.support = None  # support

proxtoolbox/proxoperators/P_M.py

 from proxtoolbox.proxoperators.proxoperator import ProxOperator, magproj
-from numpy import where, log, exp, sum
+from numpy import where, log, sum
 
 __all__ = ['P_M', 'P_M_masked', "Approx_P_M", 'Approx_P_M_masked']
 

proxtoolbox/proxoperators/P_incoherent.py

+from numpy import sum
+from numpy.core._multiarray_umath import sqrt, array
+
+from proxtoolbox.proxoperators import P_M, P_Sparsity_real
+
+__all__ = ['P_M_incoherent', 'P_Sparsity_real_incoherent']
+
+
+class P_M_incoherent(P_M):
+    """
+    Apply the Fourier-domain modulus constraint to a set of incoherent fields:
+    scale the amplitudes such that the combined intensity matches the data
+    """
+    def __init__(self, experiment):
+        super(P_M_incoherent, self).__init__(experiment)
+        self.intensity = self.data ** 2
+        self.eps = 1e-15
+
+    def eval(self, u, prox_idx=None):
+        amplitudes = self.prop.eval(u)  # Propagate to measurement domain
+        pred = sum(abs(amplitudes)**2, axis=0)  # Predicted total intensity
+        divisor = sqrt(pred + self.eps**2)  #
+        scaling = 1 - (pred + 2*self.eps**2) / divisor**3 * (pred / divisor - self.data)
+        # scaling = self.data / divisor
+        out = array([scaling * u_i for u_i in amplitudes])
+        return self.invprop.eval(out)  # Propagate back
+
+
+class P_Sparsity_real_incoherent(P_Sparsity_real):
+    def __init__(self, experiment):
+        super(P_Sparsity_real_incoherent, self).__init__(experiment)
+
+    def eval(self, u, prox_idx=None):
+        out = array([super(P_Sparsity_real_incoherent, self).eval(u_i) for u_i in u])
+        return out

proxtoolbox/proxoperators/P_orthonorm.py

 from numpy import sum, array, sqrt, zeros_like
+import numpy as np
 
 from proxtoolbox.proxoperators.proxoperator import ProxOperator
 
@@ -8,7 +9,8 @@ __all__ = ['P_orthonorm', "P_norm"]
 class P_orthonorm(ProxOperator):
     def __init__(self, experiment):
         """
-        Normalize two iterates and rotate them such that they are orthogonal
+        Normalize two iterates and rotate them such that they are orthogonal.
+        Also applies a real-valuedness projection
 
         @param experiment: experiment class
         """
@@ -16,13 +18,12 @@ class P_orthonorm(ProxOperator):
 
         self.p_norm = P_norm(experiment)
 
-    def eval(self, u, **kwargs):
+    def eval(self, u, prox_idx=None):
         # normalize u[0] and u[1]
-        u_norm = self.p_norm.eval(u)
-        norms = [sqrt(sum(abs(u) ** 2)) for u in u_norm]
+        u_norm = self.p_norm.eval(u.real)
+        norms = [sqrt(sum(abs(u_i) ** 2)) for u_i in u_norm]
         # determine angle _a_ between u[0] and u[1]
         a = sum(u_norm[0] * u_norm[1]) / (norms[0] * norms[1])
-
         if a == 0:  # for already orthogonal iterates, apply no change
             return u_norm
         elif a == -1 or a == 1:
@@ -52,9 +53,11 @@ class P_norm(ProxOperator):
         self.norm = experiment.norm_data
         # TODO: define number of iterate components here, divide norm by sqrt(n)
 
-    def eval(self, u, **kwargs):
+    def eval(self, u, prox_idx=None):
         # Normalize components of u
-        return array([self.norm * u_n / sqrt(sum(abs(u_n) ** 2)) / sqrt(len(u)) for u_n in u])
+        norms = array([sqrt(sum(abs(u_n) ** 2)) for u_n in u])
+        n_comp = len(u)
+        return array([self.norm * u[i] / norms[i] / sqrt(n_comp) for i in range(n_comp)])
 
 
 if __name__ == "__main__":

proxtoolbox/proxoperators/__init__.py

@@ -27,4 +27,5 @@ from .P_S import *
 from .P_CDP_cyclic import *
 from .sourceLocProx import *
 from .Pphase_phasepack import *
-from .P_orthonorm import *
\ No newline at end of file
+from .P_orthonorm import *
+from .P_incoherent import *
\ No newline at end of file