getting to the technical details of OrbitalTomog: phase.py

docs/manual-update-website

+
+echo "Uploading website."
+rsync -av _build/html/ rluke1@webvm.num.math.uni-goettingen.de:public_html/proxtoolbox

proxtoolbox/Problems/OrbitalTomog/Coronene_demo.py

+import sys
+# sys.path.append('../proxtoolbox/Problems/Phase')
+sys.path.append('..')
+import coronene_test
+from phase import Phase
+
+
+orbit = Phase(coronene_test.new_config)
+orbit.solve()
+orbit.show()
\ No newline at end of file

proxtoolbox/Problems/OrbitalTomog/Graphics/Phase_graphics.py

+#                      Phase_graphics.m
+#                  written on May 23, 2012 by 
+#                         Russell Luke
+#       Inst. Fuer Numerische und Angewandte Mathematik
+#                    Universitaet Gottingen
+#
+#
+# DESCRIPTION:  Script driver for viewing results from projection
+#               algorithms on various toy problems
+#
+# INPUT:  
+#              method = character string for the algorithm used.
+#         true_object = the original image
+#                 u_0 = the initial guess
+#                   u = the algorithm "fixed point"
+#              change = the norm square of the change in the
+#                              iterates
+#              error  = squared set distance error at each
+#                              iteration
+#              noneg = the norm square of the nonnegativity/support
+#                              constraint at each iteration
+#              gap  = the norm square of the gap distance, that is
+#                     the distance between the projections of the
+#                     iterates to the sets
+#
+# OUTPUT:       graphics
+# USAGE: Phase_graphics(config,output)
+#
+#############################################################
+
+from matplotlib.pyplot import subplots, show
+import numpy as np
+
+
+def Phase_graphics(config, output):
+              
+    algortihm=config['algorithm']
+    #beta0 = config['beta_0']
+    #beta_max = config['beta_max']
+    u_0 = config['u_0']
+    if output['u1'].ndim == 2:
+        u = output['u1']
+        u2 = output['u2']
+    else:
+        u = output['u1'][:,:,0]
+        u2 = output['u2'][:,:,0]
+    iter = output['iter']
+    change = output['change']
+
+    if 'time' in output:
+        time = output['time']
+    else:
+        time=999
+
+
+
+    f, ((ax1, ax2), (ax3, ax4)) = subplots(2, 2)
+
+    im=ax1.imshow(np.abs(u),cmap='gray')
+    f.colorbar(im, ax=ax1)
+    ax1.set_title('best approximation amplitude - physical constraint satisfied')
+
+    im=ax2.imshow(np.real(u),cmap='gray')
+    f.colorbar(im, ax=ax2)
+    ax2.set_title('best approximation phase - physical constraint satisfied')
+
+    im=ax3.imshow(np.abs(u2),cmap='gray')
+    f.colorbar(im, ax=ax3)
+    ax3.set_title('best approximation amplitude - Fourier constraint satisfied')
+
+    im=ax4.imshow(np.real(u2),cmap='gray')
+    f.colorbar(im, ax=ax4)
+    ax4.set_title('best approximation amplitude - Fourier constraint satisfied')
+
+    g, ((bx1, bx2), (bx3, bx4)) = subplots(2, 2)
+    im=bx1.imshow(np.abs(u),cmap='gray')
+    f.colorbar(im, ax=bx1)
+    bx1.set_title('best approximation amplitude - physical constraint satisfied')
+    im = bx2.imshow(np.real(u), cmap='gray')
+    f.colorbar(im, ax=bx2)
+    bx2.set_title('best approximation phase - physical constraint satisfied')
+    bx3.semilogy(change)
+    bx3.set_xlabel('iteration')
+    bx3.set_ylabel('$||x^{2k+2}-x^{2k}||$')
+    if 'diagnostic' in config:
+        bx4.semilogy(output['gap'])
+        bx4.set_xlabel('iteration')
+        bx4.set_ylabel('$||x^{2k+1}-x^{2k}||$')
+
+
+    show()

proxtoolbox/Problems/OrbitalTomog/Graphics/__init__.py

+from .JWST_graphics import JWST_graphics
+from .Phase_graphics import Phase_graphics
+
+__all__ = ["Phase_graphics", "JWST_graphics"]

proxtoolbox/Problems/OrbitalTomog/coronene_test.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
+
+new_config = {
+
+## We start very general.
+##==========================================
+## Problem parameters
+##==========================================
+## What is the name of the data file?
+'data_filename' : 'coronen_homo1_fourier_noise15.mat',
+
+
+## What type of object are we working with?
+## Options are: 'phase', 'real', 'nonnegative', 'complex'
+'object' : 'complex',
+
+## 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'
+'constraint' : 'support only',
+
+## What type of measurements are we working with?
+## Options are: 'single diffraction', 'diversity diffraction', 
+##              'ptychography', and 'complex'
+## Options are: '2D_ARPES', '3D_ARPES', 
+##              '2D_time', and '3D_time'
+'experiment' : '2D_ARPES',
+
+## Next we move to things that most of our users will know 
+## better than we will.  Some of these may be overwritten in the 
+## data processor file which the user will most likely write. 
+## Are the measurements in the far field or near field?
+## Options are: 'far field' or 'near field', 
+'distance' : 'far field',
+
+## What are the dimensions of the measurements?
+'Nx' : 128,
+'Ny' : 128,
+
+#'fresnel_nr' : 0,
+
+#moved this to phase
+#if(strcmp('distance,'near field'))
+#    'fresnel_nr' : 1*2*pi*'Nx,
+#else
+#    'fresnel_nr' :  0, #1*2*pi*'Nx,
+
+
+#'magn' : 1,
+
+## What are the noise characteristics (Poisson or Gaussian)?
+'noise' : 'Poisson',
+##==========================================
+##  Algorithm parameters
+##==========================================
+## Now set some algorithm parameters that the user should be 
+## able to control (without too much damage)
+
+## Algorithm:
+'algorithm' : 'AP',  # 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
+# 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:
+#if(strcmp('method,'RAAR')||strcmp('method,'AP')||...
+#        strcmp('method,'HPR')||strcmp('method,'HAAR'))
+    # the following just points this driver to a patch that communicates 
+    # the parameters defined at this level to the structures used in the 
+    # algorithms developed by Russell.  
+    #'problem_family' : 'Phase',
+#else # moreregularization parameters for Thorsten's algorithms:
+#   prbl' : complete_itreg_par(prbl),
+    # This is just a patch to Thorsten's tools.  May want to change this
+    # later
+    #'problem_family='Hohage',
+#end
+
+
+#if(strcmp('problem_family,'Phase'))
+    ## maximum number of iterations and tolerances
+    'MAXIT' : 5000,
+    'TOL' : 1e-10,
+
+    ## relaxaton parameters in RAAR, HPR and HAAR
+    'beta_0' : 0.85,    #0.95              # starting relaxation prameter (only used with
+    # HAAR, HPR and RAAR)
+    'beta_max' :0.50,             # maximum relaxation prameter (only used with
+    # HAAR, RAAR, and HPR)
+    'beta_switch' : 30,           # iteration at which beta moves from beta_0 -> beta_max
+
+    ## parameter for the data regularization 
+    ## need to discuss how/whether the user should
+    ## put in information about the noise
+    'data_ball' : .999826,
+    # 'data_ball' : .9998261e-0,
+    # the above is the percentage of the gap
+    # between the measured data and the
+    # initial guess satisfying the
+    # qualitative constraints.  For a number 
+    # very close to one, the gap is not expected 
+    # to improve much.  For a number closer to 0
+    # the gap is expected to improve a lot.  
+    # Ultimately the size of the gap depends
+    # on the inconsistency of the measurement model 
+    # with the qualitative constraints.
+#elseif(strcmp('problem_family,'Hohage'))
+#    'alpha0' : 1e4,
+#    'N_CG' : 70,
+#end
+
+# ##==========================================
+# ## parameters for plotting and diagnostics
+# ##==========================================
+# 'plotWhat.n1=2,
+# 'plotWhat.n2=3,
+# 'plotWhat.plots' : 'PYWpyw',
+# 'verbose' : 1, # options are 0 or 1
+# 'graphics' : 1, # whether or not to display figures, options are 0 or 1.
+#                    # default is 1.
+# 'anim' : 1,  # whether or not to disaply ``real time" reconstructions
+#                 # options are 0=no, 1=yes, 2=make a movie
+#                 # default is 1.
+# 'graphics_display' : [], # unless specified, a default 
+#                             # plotting subroutine will generate 
+#                             # the graphics.  Otherwise, the user
+#                             # can write their own plotting subroutine
+# 
+##==========================================
+## 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.
+                   # default is 1.
+'anim' : 2,  # whether or not to disaply ``real time" reconstructions
+                # options are 0=no, 1=yes, 2=make a movie
+                # default is 1.
+'graphics_display' : 'Phase_graphics' # unless specified, a default 
+                            # plotting subroutine will generate 
+                            # the graphics.  Otherwise, the user
+                            # can write their own plotting subroutine
+
+}
+
+##======================================================================
+##  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 
+## a module one level below this one.  
+

proxtoolbox/Problems/OrbitalTomog/orbitaltomog_data_processor.py

@@ -2,33 +2,39 @@ import numpy as np
 from scipy.io import loadmat
 import matplotlib.pyplot as plt
 
-inp_data=loadmat('../../../InputData/OrbitalTomog/coronen_homo1_fourier_noise15.mat')
-inp=inp_data["I2"]
-ny,nx = inp.shape
+def data_processor(config):
+    inp_data=loadmat('../../../InputData/OrbitalTomog/'+config['data_filename'])
+    inp=inp_data["I2"]
+    ny,nx = inp.shape
+    self.config['data'] = inp
+    self.config['norm_data'] = np.sqrt(np.sum(self.config['data']**2))
+    
+    
+    # Keep the same resolution?
+    config['Ny'],config['Nx'] = ny,nx
 
-# Autokorrelation
-threshold_autocorr=0.1
-autocorrelation=np.fft.fftshift(np.fft.ifft2(np.fft.ifftshift(inp)))
-init_support=(autocorrelation < threshold_autocorr * np.amax(autocorrelation)).astype(np.uint)
+    # Autokorrelation
+    config['threshold for support'] = 0.1
+    autocorrelation=np.fft.fftshift(np.fft.ifft2(np.fft.ifftshift(inp)))
+    config['support'] = (autocorrelation < config['threshold for support'] * np.amax(autocorrelation)).astype(np.uint)
 
-# Anfangsbedingungen
-support=autocorrelation
-ph_init = 2*np.pi*np.random.rand(ny,nx)
-# ph_init = np.angle(np.fft.fft2(ph_init));
-u0 = inp * np.exp(1j*ph_init)
-previous = np.fft.fftshift(np.fft.ifft2(np.fft.ifftshift(u0)))
+    # Anfangsbedingungen
+    ph_init = 2*np.pi*np.random.rand(ny,nx)
+    # ph_init = np.angle(np.fft.fft2(ph_init));
+    config['u0'] = inp * np.exp(1j*ph_init)
+    previous = np.fft.fftshift(np.fft.ifft2(np.fft.ifftshift(u0)))
 
-plt.figure(figsize=(15,4))
-plt.subplot(131)
-plt.imshow(inp)
-plt.colorbar()
-plt.title("Photoelectron spectrum")
-plt.subplot(132)
-plt.imshow(autocorrelation.imag)
-plt.colorbar()
-plt.title("Orbital autocorrelation, real part")
-plt.subplot(133)
-plt.imshow(init_support)
-plt.colorbar()
-plt.title("Initial support")
-plt.show()
\ No newline at end of file
+    plt.figure(figsize=(15,4))
+    plt.subplot(131)
+    plt.imshow(inp)
+    plt.colorbar()
+    plt.title("Photoelectron spectrum")
+    plt.subplot(132)
+    plt.imshow(autocorrelation.real)
+    plt.colorbar()
+    plt.title("Orbital autocorrelation, real part")
+    plt.subplot(133)
+    plt.imshow(config['support'])
+    plt.colorbar()
+    plt.title("Initial support")
+    plt.show()
\ No newline at end of file

proxtoolbox/Problems/OrbitalTomog/phase.py

+# -*- coding: utf-8 -*-
+
+from proxtoolbox.Problems.problems import Problem
+from proxtoolbox import Algorithms
+from proxtoolbox import ProxOperators
+from proxtoolbox.ProxOperators.proxoperators import ProxOperator
+from proxtoolbox.Problems.Phase import Graphics
+from numpy.linalg import norm
+import numpy as np
+import h5py
+from numpy import square, sqrt, nonzero, size
+
+
+class Phase(Problem):
+    """
+    Phase Problem
+    """
+    config = {
+    }
+    
+    def __init__(self, new_config={}):
+        """
+        The initialization of a Phase instance takes the default configuration
+        and updates the parameters with the arguments in new_config.
+        
+        Parameters
+        ----------
+        new_config : dict, optional - Parameters to initialize the problem. If unspecified, the default config is used.
+        """
+        self.config = new_config     
+
+        #call data processor, read data
+
+        module = __import__(self.config['data_filename'])
+        data_processor = getattr(module, self.config['data_filename'])
+        data_processor(self.config)
+
+        if 'Nz' not in self.config:
+            self.config['Nz'] = 1
+
+
+        #If method_config[formulation is does not exist, i.e. not specified in 
+        #the *_in.m file, use the product space as the default.
+        if 'formulation' in self.config:
+            formulation = self.config['formulation']
+        else:
+            formulation = 'product space'
+
+        # Set the projectors and inputs based on the types of constraints and 
+        # experiments
+        proxoperators = ['','','']
+
+        if self.config['constraint'] == 'support only':
+            proxoperators[0] = 'P_S'
+        elif self.config['constraint'] == 'real and support':
+            proxoperators[0] ='P_S_real'
+        elif self.config['constraint'] =='nonnegative and support':
+            proxoperators[0] ='P_SP'
+        elif self.config['constraint'] =='amplitude only':
+            proxoperators[0] ='P_amp'
+        elif self.config['constraint'] == 'phase on support':
+            proxoperators[0] ='P_Amod'
+        elif self.config['constraint'] =='minimum amplitude':
+            proxoperators[0] = 'P_min_amp'
+
+        proxoperators[1]='Approx_P_FreFra_Poisson'
+
+        self.config['proxoperators'] = []
+
+        for prox in proxoperators:
+            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
+        # parameters?  If not, and we pass everything anyway, there is no need
+        # to create a new structure element.
+
+        if 'product_space_dimension' not in self.config:
+            self.config['product_space_dimension'] = 1
+
+        # set the animation program:
+        self.config['animation']='Phase_animation'
+        #
+        # if you are only working with two sets but
+        # want to do averaged projections
+        # (= alternating projections on the product space)
+        # or RAAR on the product space (=swarming), then
+        # you will want to change product_space_dimension=2
+        # and adjust your input files and projectors accordingly. 
+        # you could also do this within the data processor
+
+        self.config['TOL2'] = 1e-15
+
+        #To estimate the gap in the sequential formulation, we build the
+        # appropriate point in the product space. This allows for code reuse.
+        # Note for sequential diversity diffraction, input.Proj1 is the "RCAAR"
+        # version of the function.
+        if formulation == 'sequential':
+            for j in range(self.config['product_space_dimension']):
+                self.config['proj_iter'] =j
+                proj1 = self.config['proxoperators'][0](self.config)
+                u_1[:,:,j]= proj1.work(self.config['u_0'])
+                self.config['proj_iter'] = mod(j,config['product_space_dimension'])+1
+                proj1 = self.config['proxoperators'][0](self.config)
+                u_1[:,:,j]= proj1.work(self.config['u_0'])
+            end
+        else: #i.e. formulation=='product space'
+            proj1 = self.config['proxoperators'][0](self.config)
+            u_1 = proj1.work(self.config['u_0'])
+            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'])
+        elif self.config['product_space_dimension'] == 1:
+            tmp_gap = (norm(u_1-u_2)/self.config['norm_rt_data'])**2
+        else:
+            tmp_gap=0
+            for j in range(self.config['product_space_dimension']):
+                # compute (||P_Sx-P_Mx||/norm_data)^2:
+                tmp_gap = tmp_gap+(norm(u_1[:,:,j]-u_2[:,:,j])/self.config['norm_rt_data'])**2
+        
+        gap_0=sqrt(tmp_gap)
+
+        # sets the set fattening to be a percentage of the
+        # initial gap to the unfattened set with 
+        # respect to the relevant metric (KL or L2), 
+        # that percentage given by
+        # input.data_ball input by the user.
+        self.config['data_ball']=self.config['data_ball']*gap_0
+        # the second tolerance relative to the oder of 
+        # magnitude of the metric
+        self.config['TOL2'] = self.config['data_ball']*1e-15
+        self.config['proxoperators']
+        self.algorithm = getattr(Algorithms, self.config['algorithm'])(self.config)
+        
+    
+    
+    def _presolve(self):
+        """
+        Prepares argument for actual solving routine
+        """
+        
+    
+    def _solve(self):
+        """
+        Runs the algorithm to solve the given sudoku problem
+        """
+#        algorithm = self.config['algorithm'](self.config)
+        
+        self.output = self.algorithm.run(self.config['u_0'],self.config['TOL'],self.config['MAXIT'])
+        print('Iterations:' + str(self.output['iter']))
+
+    
+    def _postsolve(self):
+        """
+        Processes the solution and generates the output
+        """
+        
+        
+        
+    def show(self):
+        """
+        Generates graphical output from the solution
+        """
+        
+        print("Calculation time:")
+        print(self.elapsed_time)
+        graphics = getattr(Graphics,self.config['graphics_display'])
+        graphics(self.config,self.output)
+
+    def compare_to_matlab(self):
+        """
+        Routine to test and verify results by comparing to matlab
+        Note that this is only for development and should not be used by a normal user
+        For result to match u_0 should be chosen as np.multiply(config['abs_illumination'],exp(1j*2*pi*0.5*np.ones(newres)))']
+        """
+
+        print(self.config['proxoperators'])
+
+        if self.config['experiment'] == 'JWST':
+            if self.config['algorithm'] == 'RAAR':
+                if self.config['MAXIT'] == 1:
+                    f = h5py.File('Phase_test_data/u1_1.mat')
+                elif self.config['MAXIT'] == 500 :
+                    f = h5py.File('Phase_test_data/u1_500.mat')
+                else:
+                    print("No file available to compare to.")
+                    return
+            elif self.config['algorithm'] == 'AP':
+                f = h5py.File('Phase_test_data/JWST_u1_ap_' + str(self.config['MAXIT']) + '.mat')
+
+            u1 = f['u1'].value.view(np.complex)
+        elif self.config['data_filename'] == 'Siemens_processor' and self.config['constraint'] == 'amplitude':
+            f = h5py.File('Phase_test_data/siemens_amplitude_u1_' + str(self.config['MAXIT']) + '.mat')
+            u1 = f['u1'].value.view(np.complex)
+        elif self.config['data_filename'] == 'Siemens_processor' and self.config['constraint'] == 'nonnegative and support':
+            f = h5py.File('Phase_test_data/siemens_nonneg_u1_' + str(self.config['MAXIT']) + '.mat')
+            u1 = f['u1'].value.view(np.float64)
+        elif self.config['data_filename'] == 'Siemens_processor' and self.config['constraint'] == 'real and support':
+            f = h5py.File('Phase_test_data/siemens_real_u1_' + str(self.config['MAXIT']) + '.mat')
+            u1 = f['u1'].value.view(np.float64)
+        else:
+            if self.config['algorithm'] == 'RAAR':
+                if self.config['beta_0'] == 0.95:
+                    if self.config['MAXIT'] == 1000 :
+                        f = h5py.File('Phase_test_data/tasse_u1_1000.mat')
+                    elif self.config['MAXIT'] == 20:
+                        f = h5py.File('Phase_test_data/tasse_u1_20.mat')
+                    elif self.config['MAXIT'] == 1:
+                        f = h5py.File('Phase_test_data/tasse_u1_1.mat')
+                    else:
+                        print("No file available to compare to.")
+                        return
+                elif self.config['beta_0'] == 0.50:
+                    f = h5py.File('Phase_test_data/tasse_u1_'+ str(self.config['MAXIT']) + '_beta_0_5.mat')
+                else:
+                    print("No file available to compare to.")
+                    return
+            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' or self.config['algorithm'] == 'AP_expert') 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.float64)
+
+        u1 =np.array(u1)
+        u1 = u1.T
+        print("Compare u1:")
+        #print("Nonzero indices matlab:")
+        #print(nonzero(u1))
+        #print("Nonzero indices python:")
+        #print(nonzero(self.output['u1']))
+        print("Nonzero indices equal:")
+        print(np.array_equal(nonzero(u1),nonzero(self.output['u1'])))
+        #print("Nonzero values matlab:")
+        #print(u1[nonzero(u1)])
+        #print("Nonzero values python:")
+        #print(self.output['u1'][nonzero(self.output['u1'])])
+        #print("Difference at nonzero values:")
+        #nonz = nonzero(u1)
+        diff = u1 - self.output['u1']
+        #print(diff[nonz])
+        print("Maximum norm of difference:")