Commit fa4b60ad authored by Russell Luke's avatar Russell Luke
Browse files

running...now to physics!

parent f91e1f67
......@@ -33,26 +33,23 @@ import numpy as np
def Phase_graphics(config, output):
algortihm=config['algorithm']
#algortihm=config['algorithm']
#beta0 = config['beta_0']
#beta_max = config['beta_max']
u_0 = config['u_0']
#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']
#iter = output['iter']
change = output['change']
if 'time' in output:
time = output['time']
else:
time=999
# if 'time' in output:
# time = output['time']
# else:
# time=999
f, ((ax1, ax2), (ax3, ax4)) = subplots(2, 2)
......@@ -87,5 +84,4 @@ def Phase_graphics(config, output):
bx4.set_xlabel('iteration')
bx4.set_ylabel('$||x^{2k+1}-x^{2k}||$')
show()
from .JWST_graphics import JWST_graphics
from .Phase_graphics import Phase_graphics
__all__ = ["Phase_graphics", "JWST_graphics"]
__all__ = ["Phase_graphics"]
......@@ -17,7 +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'
'constraint': 'support only',
'constraint': 'real and support',
# What type of measurements are we working with?
# Options are: 'single diffraction', 'diversity diffraction',
......@@ -56,7 +56,7 @@ new_config = {
# able to control (without too much damage)
# Algorithm:
'algorithm': 'AP', # used to be 'Projection',
'algorithm': 'DRl', # 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
......@@ -85,15 +85,18 @@ new_config = {
# if(strcmp('problem_family,'Phase'))
# maximum number of iterations and tolerances
'MAXIT': 5000,
'MAXIT': 500,
'TOL': 1e-10,
# relaxaton parameters in RAAR, HPR and HAAR
'beta_0': 0.85, # 0.95 # starting relaxation prameter (only used with
'lambda_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)
'lambda_max': 0.50,
'beta_switch': 30, # iteration at which beta moves from beta_0 -> beta_max
'lambda_switch': 100,
# parameter for the data regularization
# need to discuss how/whether the user should
......@@ -142,7 +145,8 @@ new_config = {
'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
'graphics_display': 'Phase_graphics', # unless specified, a default
'dataprocessor_plotting':True
# plotting subroutine will generate
# the graphics. Otherwise, the user
# can write their own plotting subroutine
......
......@@ -8,41 +8,46 @@ def data_processor(config):
inp_data = loadmat('../../../InputData/OrbitalTomog/' + config['data_filename'])
inp = inp_data["I2"]
ny, nx = inp.shape
config['data'] = inp
config['data'] = abs(inp)
config['norm_data'] = np.sqrt(np.sum(config['data'] ** 2))
# Keep the same resolution?
config['Ny'], config['Nx'] = ny, nx
# Autokorrelation
config['threshold for support'] = 0.1
autocorrelation = fft.fftshift(fft.ifft2(fft.ifftshift(inp)))
config['support'] = (autocorrelation < config['threshold for support'] * np.amax(autocorrelation)).astype(np.uint)
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
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'])))
# previous = fft.fftshift(fft.ifft2(fft.ifftshift(config['u_0'])))
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()
if config['dataprocessor_plotting']:
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("iFFT(data)")
plt.subplot(133)
plt.imshow(config['support'])
plt.colorbar()
plt.title("Initial support")
plt.show()
# Other settings
config['fresnel_nr'] = 0
config['FT_conv_kernel'] = 1
config['use_farfield_formula'] = True
config['magn'] = 1
config['data_sq'] = abs(config['data'])**2
config['data_zeros'] = 1
return config
......@@ -4,7 +4,7 @@ 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 proxtoolbox.Problems.OrbitalTomog import Graphics
# from .proxoperators import ProxOperator
#from phase_retrieval.back_end_utilities import norm
from numpy.linalg import norm
......@@ -124,7 +124,7 @@ class Phase(Problem):
# the second tolerance relative to the order 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)#, self.back_end)
self.algorithm = getattr(Algorithms, self.config['algorithm'])(self.config)#, self.back_end)
def _presolve(self):
"""
......@@ -155,5 +155,6 @@ class Phase(Problem):
print("Calculation time:")
print(self.elapsed_time)
_graphics = getattr(graphics, self.config['graphics_display'])
_graphics = getattr(Graphics, self.config['graphics_display'])
print("Retrieved graphics function")
_graphics(self.config, self.output)
......@@ -821,15 +821,23 @@ class Approx_P_FreFra_Poisson(ProxOperator):
fhat = ifft2(self.FT_conv_kernel * fft2(f * self.beam)) / self.magn
else:
fhat = ifft2(self.FT_conv_kernel * fft2(f)) / self.magn
#
# fhat_sq = fhat.real ** 2 + fhat.imag ** 2
# tmp = fhat_sq / self.data_sq
# tmp[self.data_zeros] = 1
# fhat_sq[self.data_zeros] = 0
# Ifhat = tmp == 0
# tmp[Ifhat] = 1
# tmp = log(tmp)
# hfhat = sum(fhat_sq * tmp + self.data_sq - fhat_sq)
fhat_sq = fhat.real ** 2 + fhat.imag ** 2
tmp = fhat_sq / self.data_sq
tmp[self.data_zeros] = 1
fhat_sq[self.data_zeros] = 0
tmp = self.data_sq / fhat_sq
Ifhat = tmp == 0
tmp[Ifhat] = 1
tmp = log(tmp)
hfhat = sum(fhat_sq * tmp + self.data_sq - fhat_sq)
hfhat = sum(self.data_sq * tmp + fhat_sq - self.data_sq)
if hfhat >= self.data_ball + self.TOL2:
p_Mhat = magproj(fhat, self.data)
if self.use_farfield_formula:
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment