Commit fa3112db authored by s.gretchko's avatar s.gretchko
Browse files

Added Kruger Experiment, DataProcessor, AvP demo, and verbose mode

parent 69aaab3b
import SetProxPythonPath
from proxtoolbox.experiments.phase.Krueger_Experiment import Krueger_Experiment
CDI = Krueger_Experiment(algorithm = 'AvP', data_ball = 24)
CDI.run()
CDI.show()
......@@ -71,7 +71,7 @@ class DataProcessor:
self.norm_data = self.norm_rt_data
# The following lines of code correspond to the Cone_and_Sphere branch. However,
# the case where the data is a list (or cell in Matlab) is not treated
# the case where the data is a cell is not treated (in Matlab implementation)
# TODO: Need to check is this is OK.
if not isCell(self.data):
tmp = self.data.shape
......
......@@ -225,7 +225,6 @@ class Experiment(metaclass = ExperimentMetaClass):
self.inverse_propagator = getattr(proxoperators, self.inverse_propagator)
self.TOL2 = 1e-20
self.data_ball = 1e-30
# TODO in Cone_and_Sphere estimation of the parameters above is done differently (set_fattener)
# New matlab code:
......@@ -258,20 +257,25 @@ class Experiment(metaclass = ExperimentMetaClass):
The initial iterate is obtained from the data processor.
TODO: implement multiple runs
"""
print ("Running " + self.algorithm_name + " on " + self.experiment_name + "...")
print("Running " + self.algorithm_name + " on " + self.experiment_name + "...")
if self.verbose:
self.printInput()
t = time.time()
self.output = self.algorithm.run(self.data_processor.u0) # run algorithm
if 'stats' in self.output:
self.output['stats']['time'] = time.time() - t
self.output['stats']['time'] = time.time() - t
print("Took", self.output['stats']['iter'], "iterations and",
self.output['stats']['time'], "seconds.")
if self.save_output:
self.saveOutput()
if self.verbose:
self.printBenchmark()
def setupProxOperators(self):
"""
......@@ -354,6 +358,54 @@ class Experiment(metaclass = ExperimentMetaClass):
except TypeError as e:
errMsg = e.__str__()
print("Unable to serialize object for json output file: ", errMsg)
def printInput(self):
print("Input:")
vars_dict = vars(self)
for var in vars_dict:
obj = vars_dict[var]
# Only print variables that are not an instance of a class
if hasattr(obj, '__dict__') is False:
print('\t', var, ':', obj)
def printBenchmark(self):
benchmark = self.createBenchmark()
print("Benchmark:")
for var in benchmark:
obj = benchmark[var]
print('\t', var, ':', obj)
def createBenchmark(self):
# create a dictionary with relevant results
benchmark = {}
iterCount = self.output['stats']['iter']
benchmark['iter'] = iterCount
benchmark['time'] = self.output['stats']['time']
changes = self.output['stats']['changes']
benchmark['end_change'] = changes[iterCount-1]
if iterCount > 8:
# if algorithm exits in fewer than 8 iterations, this will
# throw an error because it is taking an average over
# the last 10 iterations. Adjust accordingly.
rates = changes[1:iterCount] / changes[0:iterCount-1]
rate = sum(rates[iterCount-5:iterCount-1])/3
benchmark['rate'] = rate
a_posteriori_error = rate/(1-rate)*changes[iterCount-2]
benchmark['a_posteriori_error'] = a_posteriori_error
if self.diagnostic:
self.truth = self.data_processor.truth
if self.truth is not None and 'rel_errors' in self.output['stats']:
rel_errors = self.output['stats']['rel_errors']
if iterCount == 1:
true_error = rel_errors[iterCount]
else:
true_error = rel_errors[iterCount-1]
benchmark['true_error'] = true_error
return benchmark
def estimateGap(self):
"""
......
......@@ -51,7 +51,7 @@ class JWST_Experiment(PhaseExperiment):
self.nProx = self.data_processor.sets
self.product_space_dimension = 1
# add projectors onto amplitude data
for j in range(self.nProx-1):
for _j in range(self.nProx-1):
self.proxOperators.append('Approx_Pphase_FreFra_Poisson') #'Approx_P_cyclic_JWST_Poisson'
self.proxOperators.append('P_amp_support') #projector onto pupil amplitude constraint
else:
......@@ -61,7 +61,7 @@ class JWST_Experiment(PhaseExperiment):
self.proxOperators.append('Prox_product_space')
self.n_product_Prox = self.product_space_dimension
self.productProxOperators = []
for j in range(self.n_product_Prox-1):
for _j in range(self.n_product_Prox-1):
self.productProxOperators.append('Approx_Pphase_FreFra_Poisson')
self.productProxOperators.append('P_amp_support') #projector onto pupil amplitude constraint
#input.Prox0 = input.Prox{1}
......@@ -97,10 +97,10 @@ class JWST_Experiment(PhaseExperiment):
title = "Algorithm " + algo_desc
# figure 904
titles = ["best approximation amplitude - physical constraint satisfied",
"best approximation phase - physical constraint satisfied",
"best approximation amplitude - physical constraint satisfied",
"best approximation phase - physical constraint satisfied"]
titles = ["Best approximation amplitude - physical constraint satisfied",
"Best approximation phase - physical constraint satisfied",
"Best approximation amplitude - Fourier constraint satisfied",
"Best approximation phase - Fourier constraint satisfied"]
f = self.createFourImageFigure(u, u2, titles)
f.suptitle(title)
plt.subplots_adjust(hspace = 0.3) # adjust vertical space (height) between subplots (default = 0.2)
......@@ -121,7 +121,7 @@ class JWST_Experiment(PhaseExperiment):
ax3.set_xlabel(label)
ax3.set_ylabel('Log of change in iterates')
if 'gaps' in self.output:
if 'gaps' in self.output['stats']:
gaps = self.output['stats']['gaps']
ax4.semilogy(gaps)
ax4.set_xlabel(label)
......
from proxtoolbox.experiments.dataProcessor import DataProcessor
import numpy as np
from numpy import exp, sqrt, log2, log10, ceil, floor, unravel_index, argmax, zeros
from scipy.io import loadmat
from numpy.fft import fftshift, ifft2, fft2
from proxtoolbox.Utilities.mypoissonrnd import mypoissonrnd
from proxtoolbox.Utilities.gaussian import gaussian
from proxtoolbox.Utilities.cell import Cell, isCell
class Krueger_DataProcessor(DataProcessor):
"""
Data processor for the Krueger experiment
"""
def __init__(self):
"""
"""
super(Krueger_DataProcessor, self).__init__()
self.farfield = None
self.fresnel_nr = None
self.magn = None
self.FT_conv_kernel = Cell(1)
self.beam = Cell(1)
self.data_zeros = Cell(1)
self.abs_illumination = None
self.support_idx = None
self.product_space_dimension = 1
self.supp_phase = None
def loadData(self, experiment):
"""
Load Krueger dataset. Create the initial iterate.
"""
np.random.seed(1234) # for tests
# Sample: NTT-AT, ATN/XREESCO-50HC
# 500 nm thick Ta structure: amplitude transmission 0.9644,
# phase shift 0.4rad at E = 17.5 keV
# parameters see below
# empty waveguide (WG) beam
WG = loadmat('../InputData/Phase/WG_beam.mat')
WG = WG['WG']
# hologram
print('Loading data hologram_not-normalized.mat')
I_exp = loadmat('../InputData/Phase/hologram_not-normalized.mat')
I_exp = I_exp['I_exp']
I_exp[np.isnan(I_exp)] = 1
# The following is NOT used because it is not clear how the ``normalized" hologram
# is obtained. I suspect it is obtained by simply dividing out the empty beam
# data in the obervation plane. But this is incorrect. Despite the
# efforts of the Hohage team to get around this error by ad hoc regularization,
# we take the following exact approach: back propagate the empty beam, and
# divide this from the unknown object in the object plane. This approach is true
# to the imaging model and does not introduce any approximations. It does, however,
# expose the fact that we do not know the phase of the beam.
# clear I_exp
# load 'data/hologram.mat'
## single image reconstruction
# number of pixels
Ny = I_exp.shape[0]
Nx = I_exp.shape[1]
experiment.Ny = Ny
experiment.Nx = Nx
##########################
# Experimental parameters
##########################
# energy in keV
E = 17.5
# wavelength [m]
lambd = 12.398 / E * 1E-10
k = 2 * np.pi / lambd
# distance source-sample [m]
z1 = 7.48e-3
# distance sample-detector [m]
z2 = 3.09
# effective distance of detector
z_eff = z1 * z2 / (z1 + z2)
# effective magnification factor
M = (z1 + z2) / z1
# pixel size in detector plane
dx_det = 55e-6
dy_det = 55e-6
# magnified coordinate system in detector plane
# [X_det,Y_det] = meshgrid(dx_det*[0:1:Nx-1],dy_det*[0:1:Ny-1]);
# demagnified pixel size in sample plane
dx = dx_det / M
self.dx = dx
dy = dy_det / M
self.dy = dy
# coordinate system in sample plane
# [X,Y] = meshgrid(dx*((1:1:Nx)-floor(Nx/2)-1),dy*((1:1:Ny)-floor(Ny/2)-1));
# magnified coordinate system in detector plane
# [X_det,Y_det] = meshgrid(dx_det*((1:1:Nx)-floor(Nx/2)-1),dy_det*((1:1:Ny)-floor(Ny/2)-1));
# grid conversion in q-space
dqx = 2 * np.pi / (Nx * dx)
dqy = 2 * np.pi / (Ny * dy)
[Qx, Qy] = np.meshgrid(dqx * (range(1, Nx + 1) - np.floor(Nx / 2) - 1),
dqy * (range(1, Ny + 1) - np.floor(Ny / 2) - 1))
# Prepare data for ProxToolbox:
# Fresnel propagator:
self.farfield = False
#config['Nx'] = Nx
#config['Ny'] = Ny
self.fresnel_nr = 1 * 2 * np.pi * Nx # Fresnel number
self.magn = 1 # magnification factor (should this be M above, or
# is it okay since the demagnified pixel size is used?)
kappa = np.sqrt(k**2 - (np.square(Qx) + np.square(Qy)))
self.FT_conv_kernel[0] = fftshift(exp(-1j * kappa * z_eff))
self.beam[0] = abs(ifft2(fft2(sqrt(WG * 7.210116465530644e-04)) / self.FT_conv_kernel[0]))
# the scaling factor of the waveguide array above was
# reverse engineered from the scaling that was apparently
# used in the preparation of the data in the file hologram.mat
# I don't want to use this hologram as the data, preferring instead
# to keep the experimental data as close to the actual observation
# as possible. Also, I am disagreeing with what my colleagues in
# physics think is the correct way to compensate for a structured
# empty beam. According to my reverse engineering, it appears that
# they have divided the EMPTY BEAM MEASUREMENT WG from the object
# IN THE OBJECT PLANE. It is true that you need to do the empty beam
# correction in the object plane, though not with the empty beam, but
# rather the reverse propagated empty beam, as I have done above. The
# reason being that the empty beam is measured in the measurement plane,
# so it does not make sense to divide the object&beam in the object
# plane by the empty beam in the measurement plane - you should divide
# by the empty beam propagated back to the object plane. The difficulty in
# this is that we do not know the phase of the empty beam. The above
# formulation sets the phase of the beam to zero at the object plane.
# The beam is saved as the square root of the empty beam to conform with
# the projeciton operations working on the amplitude instead of the
# magnitude.
# problem data
self.data = Cell(1)
self.data[0] = I_exp
self.rt_data = Cell(1)
self.rt_data[0] = sqrt(I_exp)
#self.data_zeros = self.data == 0
self.data_zeros[0] = np.where(self.data[0] == 0)
self.norm_rt_data = sqrt(sum(sum(self.data[0])))
# use the abs_illumination field to represent the
# support constraint.
self.abs_illumination = np.ones(I_exp.shape)
self.support_idx = np.nonzero(self.abs_illumination)
self.sets = 2
self.product_space_dimension = 1
self.supp_phase = []
# start values
self.u0 = ifft2(fft2(self.rt_data[0] * self.magn) / self.FT_conv_kernel[0]) / self.beam[0]
from proxtoolbox.experiments.phase.phaseExperiment import PhaseExperiment
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.pyplot import subplots, show, figure
class Krueger_Experiment(PhaseExperiment):
'''
Krueger experiment class
'''
def __init__(self,
**kwargs):
"""
"""
# call parent's __init__ method
super(Krueger_Experiment, self).__init__(**kwargs)
# do here any data member initialization
self.TOL2 = 1e-14*self.data_ball
@staticmethod
def getDefaultParameters():
defaultParams = {
'experiment_name' : 'Krueger',
'data_processor_name': 'Krueger_DataProcessor',
'data_processor_package': 'proxtoolbox.experiments.phase',
'object': 'phase',
'constraint': 'phase on support',
'noise': None,
'snr': 10e2,
'farfield': False,
'MAXIT': 1000,
'TOL': 5e-7,
'lambda_0': 0.05,
'lambda_max': 0.05,
'lambda_switch': 4,
'data_ball': 1e-0,
'diagnostic': True,
'iterate_monitor_name': 'FeasibilityIterateMonitor',
'verbose': 1,
'graphics': 1,
'anim': 0
}
return defaultParams
def setupProxOperators(self):
"""
Determine the prox operators to be used for this experiment
"""
super(Krueger_Experiment, self).setupProxOperators() # call parent's method
self.propagator = 'Propagator_FreFra'
self.inverse_propagator = 'InvPropagator_FreFra'
self.proxOperators = []
self.proxOperators.append('P_Amod')
self.proxOperators.append('Approx_Pphase_FreFra_Poisson')
self.sets = self.data_processor.sets
self.nProx = self.sets
def show(self):
# display data
f, ((ax1, ax2), (ax3, ax4)) = subplots(2, 2, \
figsize = (self.figure_width, self.figure_height),
dpi = self.figure_dpi)
self.createImageSubFigure(f, ax1, self.data_processor.beam[0],
'Empty beam - Detector plane')
self.createImageSubFigure(f, ax2, self.data_processor.data[0],
'Uncorrected near field observation')
self.createImageSubFigure(f, ax3, np.abs(self.data_processor.u0),
'Initial guess amplitude')
self.createImageSubFigure(f, ax4, np.angle(self.data_processor.u0),
'Initial guess phase')
plt.subplots_adjust(wspace = 0.3) # adjust horizontal space (width)
# between subplots (default = 0.2)
plt.subplots_adjust(hspace = 0.3) # adjust vertical space (height)
# between subplots (default = 0.2)
f.suptitle('Krueger Data')
# call parent to display the other plots
super(Krueger_Experiment, self).show()
def createImageSubFigure(self, f, ax, u, title = None):
# Note: this method overides the method defined in
# class PhaseExperiment
x_max = self.Nx*self.data_processor.dx*1E6
y_max = self.Ny*self.data_processor.dy*1E6
im = ax.imshow(u, extent=[0, x_max, y_max, 0], cmap='gray')
# The "magic" values for fraction and pad adjust the
# size of the color bar so that its height is comparable
# to the plot:
f.colorbar(im, ax=ax, fraction=0.046, pad=0.04)
ax.set_xlabel(r'$x$ [$\mathrm{\mu m}$]')
ax.set_ylabel(r'$y$ [$\mathrm{\mu m}$]')
if title is not None:
ax.set_title(title)
......@@ -53,10 +53,10 @@ class PhaseExperiment(Experiment):
title = "Algorithm: " + algo_desc
# figure 904
titles = ["best approximation amplitude - physical constraint satisfied",
"best approximation phase - physical constraint satisfied",
"best approximation amplitude - physical constraint satisfied",
"best approximation phase - physical constraint satisfied"]
titles = ["Best approximation amplitude - Physical constraint satisfied",
"Best approximation phase - Physical constraint satisfied",
"Best approximation amplitude - Fourier constraint satisfied",
"Best approximation phase - Fourier constraint satisfied"]
if m == self.product_space_dimension or m == 1:
f = self.createFourGraphFigure(u[0:], u2[0,:], titles)
......
from proxtoolbox.proxoperators.proxoperator import ProxOperator
import numpy as np
from numpy import abs
class P_Amod(ProxOperator):
"""
Projection onto support constraints
Based on Matlab code written by Russell Luke (Inst. Fuer
Numerische und Angewandte Mathematik, Universitaet
Gottingen) on June 1, 2017.
"""
def __init__(self, experiment):
self.support_idx = experiment.data_processor.support_idx
def eval(self, u, prox_idx = None):
"""
Projects the input data onto support constraints.
Parameters
----------
u : ndarray
Function in the physical domain to be projected
prox_idx : int, optional
Index of this prox operator
Returns
-------
p_A : ndarray
The projection in the physical (time) domain.
Leaves the phase alone on the support, resets
amplitude to 1 and sets amplitude to 1 with 0
phase outside support
"""
p_A = np.ones_like(u)
p_A[self.support_idx] = u[self.support_idx] / abs(u[self.support_idx])
return p_A
......@@ -15,6 +15,7 @@ from .P_avg import *
from .P_diag import *
from .P_amp_support import *
from .P_SP import *
from .P_Amod import *
'''
from .P_Sparsity import *
......
......@@ -98,9 +98,9 @@ class Propagator_FreFra(Propagator_FreFra_Base):
u_hat = FFT(u)/sqrt(self.Nx*self.Ny)
else: # near field
if self.beam is not None:
u_hat = IFFT(self.FT_conv_kernel[j]*FFT(u*self.beam[j]))/self.magn[j]
u_hat = IFFT(self.FT_conv_kernel[j]*FFT(u*self.beam[j]))/self.magn
else:
u_hat = IFFT(self.FT_conv_kernel[j]*FFT(u))/self.magn[j]
u_hat = IFFT(self.FT_conv_kernel[j]*FFT(u))/self.magn
return u_hat
......@@ -158,9 +158,9 @@ class InvPropagator_FreFra(Propagator_FreFra_Base):
u_new = IFFT(p_Mhat)*sqrt(self.Nx*self.Ny)
else: # near field
if self.beam is not None:
u_new = IFFT(conj(self.FT_conv_kernel[j])*FFT(p_Mhat*self.magn[j])) / self.beam[j]
u_new = IFFT(conj(self.FT_conv_kernel[j])*FFT(p_Mhat*self.magn)) / self.beam[j]
else:
u_new = IFFT(conj(self.FT_conv_kernel[j])*FFT(p_Mhat*self.magn[j]))
u_new = IFFT(conj(self.FT_conv_kernel[j])*FFT(p_Mhat*self.magn))
return u_new
......
......@@ -575,37 +575,6 @@ class Q_Heau(ProxOperator):
return y + (nu / rho) * (xsi * (x - y) + nu * (z - y))
# P_Amod.m
# written on June 1, 2017 by
# Russell Luke
# Inst. Fuer Numerische und Angewandte Mathematik
# Universitaet Gottingen
#
# DESCRIPTION: Projection subroutine for projecting onto
# support constraints
#
# INPUT: input.supp_ampl = OBJECT (TIME) DOMAIN CONSTRAINT: indeces of the
# nonzero elements in the object domain.
# u = function in the physical domain to be projected
#
# OUTPUT: p_A = the projection IN THE PHYSICAL (time) DOMAIN,
# leaves the phase alone on the support, resets
# amplitude to 1 and sets amplitude to 1 with 0
# phase outside support.
#
# USAGE: p_A = P_Amod(input,u)
#
# Nonstandard Matlab function calls:
class P_Amod(ProxOperator):
def __init__(self, experiment):
self.support_idx = config['support_idx']
def eval(self, u, prox_idx = None):
p_A = np.ones_like(u)
p_A[self.support_idx] = u[self.support_idx] / abs(u[self.support_idx])
return (p_A)
Markdown is supported
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