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

Ptychography is working with warmup algorithm

parent a5cb1bf4
......@@ -2,6 +2,22 @@
import SetProxPythonPath
from proxtoolbox.experiments.ptychography.ptychographyExperiment import PtychographyExperiment
Pty = PtychographyExperiment(debug = False, warmup = False)
warmupParams = {
'algorithm': 'DRl',
'MAXIT': 10,
'TOL': 5e-5,
'lambda_0': 0.75,
'lambda_max': 0.75,
'lambda_switch': 20,
'data_ball': 1e-30,
'diagnostic': True,
'rotate': False,
'iterate_monitor_name': 'IterateMonitor',
'verbose': 1,
'graphics': 0,
'anim': 0
}
Pty = PtychographyExperiment(debug = False, warmup = True, MAXIT=10, warmup_params=warmupParams)
Pty.run()
Pty.show()
......@@ -59,7 +59,7 @@ class PHeBIE(Algorithm):
# Loop through the block-wise implicit-explicit update steps.
K = len(self.proxOperators)
for s in range(K):
u_new[s] = self.explicit[s].eval(u, s) # explicit prox evaluation
u_new[s] = self.explicit[s].eval(u_new, s) # explicit prox evaluation
u_new[s] = self.proxOperators[s].eval(u_new[s], s) # implicit prox evauation
# apply first order acceleration if required
......
......@@ -258,7 +258,7 @@ class Experiment(metaclass=ExperimentMetaClass):
t = time.time()
self.output = self.algorithm.run(self.u0) # run algorithm
self.output = self.runAlgorithm(self.u0)
self.output['stats']['time'] = time.time() - t
......@@ -271,6 +271,26 @@ class Experiment(metaclass=ExperimentMetaClass):
if self.verbose:
self.printBenchmark()
def runAlgorithm(self, u):
"""
Run the algorithm associated with this experiment.
Called by the run() method. May be overriden to provide
additional behaviour, such as running a warmup algorithm
before the main algorithm.
Parameters
----------
u: ndarray or cell of ndarray objects
Initial iterate.
Returns
-------
dictionary containing the output of the algorithm. In
particular, it contains the last iterate and various
statistics.
"""
return self.algorithm.run(u)
def loadData(self):
"""
Load or generate the dataset that will be used for
......
......@@ -73,7 +73,7 @@ class PtychographyExperiment(Experiment):
defaultWarmupParams = {
'algorithm': 'DRl',
'MAXIT': 20,
'MAXIT': 2,
'TOL': 5e-5,
'lambda_0': 0.75,
'lambda_max': 0.75,
......@@ -203,7 +203,7 @@ class PtychographyExperiment(Experiment):
else:
setattr(self, var, value)
# another issue is that the integer variables contained
# Another issue is that the integer variables contained
# in the data file are converted to float by .mat file
# reader (the problem lies in the library h5py that provides
# the services to read the .mat file)
......@@ -618,7 +618,7 @@ class PtychographyExperiment(Experiment):
This is a helper method called during the initialization
process. It is overriden here to deal with an additional
category of prox operators.
category of prox operators : explicit prox operators.
'''
super(PtychographyExperiment, self).retrieveProxOperatorClasses() # call parent's method
......@@ -641,11 +641,67 @@ class PtychographyExperiment(Experiment):
if self.warmup:
# initialize warmup algorithm
warmupExp = copy.copy(self) # create a copy of this experiment
# update this copy with warmup parameters
for var, value in self.warmup_params.items():
if var == 'algorithm':
setattr(warmupExp, 'algorithm_name', value)
else:
setattr(warmupExp, var, value)
# TODO
# set-up prox operators
warmupExp.proxOperators = []
warmupExp.productProxOperators = []
warmupExp.FT_conv_kernel = []
if warmupExp.algorithm_name == 'DRl':
warmupExp.nProx = 2
# prox 0
warmupExp.proxOperators.append(getattr(proxoperators,
'Prox_ptychography_ptwise_diag'))
# prox 1
warmupExp.proxOperators.append(getattr(proxoperators,
'Prox_ptychography_product_space'))
# prox operators for product space
warmupExp.n_product_Prox = self.N_pie
warmupExp.FT_conv_kernel = Cell(self.N_pie)
for j in range(warmupExp.n_product_Prox):
warmupExp.productProxOperators.append(getattr(proxoperators,
'Approx_Pphase_FreFra_Poisson'))
warmupExp.FT_conv_kernel[j] = self.u0[0]
# initial iterate
warmupExp.u0 = Cell(2)
warmupExp.u0[0] = self.u0[1]
warmupExp.u0[1] = self.u0[2]
else:
errMsg = "Warmup algorithm " + warmupExp.algorithm_name \
+ " not supported yet"
raise NotImplementedError(errMsg)
# instantiate algorithm
warmupExp.iterate_monitor_name = 'IterateMonitor'
warmupExp.instanciateAlgorithm()
self.warmupExp = warmupExp
def runAlgorithm(self, u):
if self.warmup:
# run warmup algorithm
print("Running warmup algorithm " + self.warmupExp.algorithm_name)
warmup_out = self.warmupExp.algorithm.run(self.warmupExp.u0)
warmup_u = warmup_out['u']
u0 = Cell(3)
u0[0] = self.u0[0]
u0[1] = warmup_u[0]
u0[2] = warmup_u[1]
else:
u0 = self.u0
print("Running main algorithm (" + self.algorithm_name + ")")
output = super(PtychographyExperiment, self).runAlgorithm(u0) # call parent's method
return output
def show(self):
......@@ -766,21 +822,6 @@ class PtychographyExperiment(Experiment):
title = "Objective value, " + algo_label
plt.title(title)
"""
if isfield(output.stats,'gap')
figure(701)
if(strcmp(input.method,'PHeBIE'))
semilogy(gap),xlabel('iteration'),ylabel('PHeBIE Objective Function','Interpreter','LaTex')
label = ['Algorithm: ',method]; %, ', \beta=',num2str(beta0),' to ',num2str(beta_max)];
title(['objective value, ' label])
else
semilogy(gap),xlabel('iteration'),ylabel('gao','Interpreter','LaTex')
label = ['Algorithm: ',method]; %, ', \beta=',num2str(beta0),' to ',num2str(beta_max)];
title(['gap, ' label])
end
end
"""
plt.show()
def replaceZeroByMin(self, u):
......
from proxtoolbox.proxoperators.proxoperator import ProxOperator
from proxtoolbox import proxoperators
from proxtoolbox.Utilities.cell import Cell, isCell
import numpy as np
from numpy import zeros, ones, zeros_like
......@@ -273,16 +274,150 @@ class P_supp_amp_band(ProxOperator):
u_new = u.copy()
# Apply the support constraint
if self.switch_object_support_constraint:
u_new = u_new * self.object_support
# this may be less efficient that previous code;
# however multiplication by zero keeps the sign instead
# of properly setting the value to zero which creates
# issues when taking the phase (angle)
u_new = np.where(self.object_support == 0, 0, u_new)
# oldcode: u_new = u_new * self.object_support
# Enforce that the object has values in [minTrue,maxTrue]
maxTrue = self.trans_max_true
minTrue = self.trans_min_true
absObject = np.abs(u)
absObject = np.abs(u_new)
# TODO can use np.where() to achieve the same result
high = (absObject > maxTrue).astype(absObject.dtype)
low = (absObject < minTrue).astype(absObject.dtype)
u_new *= (1.0 - low)*(1.0 - high) + (low*minTrue + high*maxTrue)/(absObject + 1e-30)
return u_new
class Prox_ptychography_ptwise_diag(ExplicitPtychographyProxOperator):
"""
The object update for the PHeBIE algorithm
Ref.: "Proximal Heterogeneous Block Implicit-Explicit Method
and Application to Blind Ptychographic Diffraction Imaging",
R. Hesse, D. R. Luke, S. Sabach and M. K. Tam,
SIAM J. on Imaging Sciences, 8(1):426--457 (2015).
Based on Matlab code written by Russell Luke (Inst. Fuer
Numerische und Angewandte Mathematik, Universitaet
Gottingen) on 16th February 2019.
"""
def __init__(self, experiment):
super(Prox_ptychography_ptwise_diag, self).__init__(experiment)
self.FT_conv_kernel = experiment.FT_conv_kernel
proxClass = getattr(proxoperators, 'P_supp_amp_band')
self.prox_supp_amp_band = proxClass(experiment)
def eval(self, u, prox_idx=None):
"""
Evaluate this prox operator
Parameters
----------
u : Cell
Cell of arrays in the domain to be projected.
Assumes u[1] is the object, u[2] is the
block of variables (product space) satisfying the intensity
measurements
prox_idx : int, optional
Index of this prox operator
Returns
-------
u_new : Cell
The projection.
"""
xSum = zeros_like(u[0])
phiSum = zeros(u[0].shape, dtype = np.dtype(np.complex128)) # must be complex
#conjProbe = np.conj(u[0])
#modSqProbe = np.real(conjProbe * u[0])
isCell_u1 = isCell(u[1])
for xPos in range(self.N_pie):
indy = (self.positions[xPos,0] + self.rangeNy).astype(int)
indx = (self.positions[xPos,1] + self.rangeNx).astype(int)
if isCell_u1:
u1_values = u[1][xPos]
else:
u1_values = u[1][:,:,xPos]
conv_kernel = self.FT_conv_kernel[xPos]
conjProbe = np.conj(conv_kernel)
modSqProbe = np.real(conjProbe * conv_kernel)
for y in range(self.Ny):
xSum[indy[y],indx] += modSqProbe[y,:]
phiSum[indy[y],indx] += conjProbe[y,:] * u1_values[y,:]
xsum_denom = self.overrelax * np.maximum(xSum, 1e-30)
u_new = Cell(u.shape)
u_new[0] = phiSum/xsum_denom
u_new[0] = self.prox_supp_amp_band.eval(u_new[0])
u_new[1] = u[1].copy()
return u_new
class Prox_ptychography_product_space(ExplicitPtychographyProxOperator):
"""
Projection subroutine onto constraints arranged in a product space
Based on Matlab code written by Russell Luke (Inst. Fuer
Numerische und Angewandte Mathematik, Universitaet
Gottingen) on Oct 26, 2017.
"""
def __init__(self, experiment):
super(Prox_ptychography_product_space, self).__init__(experiment)
self.FT_conv_kernel = experiment.FT_conv_kernel
# instantiate product prox operators
self.proxOps = []
for proxClass in experiment.productProxOperators:
if proxClass != '':
prox = proxClass(experiment)
self.proxOps.append(prox)
def eval(self, u, prox_idx=None):
"""
Evaluate this prox operator
Parameters
----------
u : Cell
Cell of arrays in the domain to be projected.
Assumes u[1] is the object, u[2] is the
block of variables (product space) satisfying the intensity
measurements.
prox_idx : int, optional
Index of this prox operator
Returns
-------
u_new : Cell
The projection, a cell of arrays in the product space
"""
isCell_u1 = isCell(u[1])
u_new = u.copy() #TODO may have to check if a deep copy is needed
for xPos in range(self.N_pie):
indy = (self.positions[xPos,0] + self.rangeNy).astype(int)
indx = (self.positions[xPos,1] + self.rangeNx).astype(int)
tmp = self.FT_conv_kernel[xPos] * u[0][indy,:][:,indx]
if isCell_u1:
u_new[1][xPos] = tmp
else:
u_new[1][:,:,xPos] = tmp
K = len(self.proxOps)
for k in range(K):
if isCell_u1:
u_new[1][k] = self.proxOps[k].eval(u[1][k], k)
else:
u_new[1][:,:,k] = self.proxOps[k].eval(u[1][:,:,k], k)
return u_new
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