Commit 1c2e7054 authored by s.gretchko's avatar s.gretchko
Browse files

Added DyRePr algorithm, Gen. accel., JWST demos + animation

parent 0819a3fd
......@@ -2,6 +2,6 @@
import SetProxPythonPath
from proxtoolbox.experiments.CT.ART_Experiment import ART_Experiment
ART = ART_Experiment(algorithm='AP')
ART = ART_Experiment(algorithm='AP', anim=True, anim_step=1)
ART.run()
ART.show()
......@@ -2,6 +2,6 @@
import SetProxPythonPath
from proxtoolbox.experiments.CT.ART_Experiment import ART_Experiment
ART = ART_Experiment(algorithm='CP', formulation ='cyclic')
ART = ART_Experiment(algorithm='CP', formulation ='cyclic', anim=True, anim_step=1)
ART.run()
ART.show()
import SetProxPythonPath
from proxtoolbox.experiments.phase.JWST_Experiment import JWST_Experiment
JWST = JWST_Experiment(algorithm='CDRl', formulation='cyclic',
lambda_0=1.0, lambda_max=1.0, MAXIT=600)
JWST.run()
JWST.show()
import SetProxPythonPath
from proxtoolbox.experiments.phase.JWST_Experiment import JWST_Experiment
JWST = JWST_Experiment(algorithm='CDRl', formulation='cyclic',
lambda_0=1.0, lambda_max=1.0, MAXIT=300, noise=True)
JWST.run()
JWST.show()
# old demo
import SetProxPythonPath
from proxtoolbox.experiments.phase.JWST_Experiment import JWST_Experiment
JWST = JWST_Experiment(algorithm='DRl', lambda_0=1.0, lambda_max=1.0,
MAXIT=300, noise=True)
JWST.run()
JWST.show()
......@@ -2,6 +2,6 @@
import SetProxPythonPath
from proxtoolbox.experiments.phase.JWST_Experiment import JWST_Experiment
JWST = JWST_Experiment(algorithm='DRl', lambda_0=0.55, lambda_max=0.55)
JWST = JWST_Experiment(algorithm='DRl', lambda_0=0.55, lambda_max=0.55, anim=True)
JWST.run()
JWST.show()
# old demo
import SetProxPythonPath
from proxtoolbox.experiments.phase.JWST_Experiment import JWST_Experiment
JWST = JWST_Experiment(algorithm = 'RAAR', lambda_0 = 1.0, lambda_max = 1.0, MAXIT = 1000)
JWST = JWST_Experiment(algorithm='DyRePr', TOL=5e-9, anim=True, MAXIT=20)
JWST.run()
JWST.show()
import SetProxPythonPath
from proxtoolbox.experiments.phase.JWST_Experiment import JWST_Experiment
JWST = JWST_Experiment(algorithm='DyRePr', TOL=5e-9, noise=True)
JWST.run()
JWST.show()
import SetProxPythonPath
from proxtoolbox.experiments.phase.JWST_Experiment import JWST_Experiment
JWST = JWST_Experiment(algorithm='AvP', accelerator_name='GenericAccelerator')
JWST.run()
JWST.show()
import SetProxPythonPath
from proxtoolbox.experiments.phase.JWST_Experiment import JWST_Experiment
JWST = JWST_Experiment(algorithm='AvP', accelerator_name='GenericAccelerator',
noise=True)
JWST.run()
JWST.show()
......@@ -3,6 +3,7 @@
import SetProxPythonPath
from proxtoolbox.experiments.phase.JWST_Experiment import JWST_Experiment
JWST = JWST_Experiment(algorithm = 'HPR', lambda_0 = 0.85, lambda_max = 0.85, MAXIT = 500)
JWST = JWST_Experiment(algorithm = 'HPR', lambda_0 = 0.85, lambda_max = 0.85,
MAXIT = 5000, anim=True)
JWST.run()
JWST.show()
......@@ -15,7 +15,7 @@ warmupParams = {
'iterate_monitor_name': 'IterateMonitor',
'verbose': 1,
'graphics': 0,
'anim': 0
'anim': False
}
Pty = PtychographyExperiment(warmup=True, MAXIT=10, warmup_params=warmupParams)
......
......@@ -121,7 +121,7 @@ warmupParams = {
'iterate_monitor_name': 'IterateMonitor',
'verbose': 1,
'graphics': 0,
'anim': 0
'anim': False
}
Pty = PtychographyExperiment(debug = False, warmup = True, MAXIT=10, warmup_params=warmupParams, verbose=0)
......
......@@ -202,7 +202,7 @@ by defining the static method :meth:`getDefaultParameters`
'rotate': False,
'verbose': 0,
'graphics': 1,
'anim': 2,
'anim': False,
'debug': True
}
return defaultParams
......
......@@ -91,7 +91,7 @@ defines default parameters via its static method :meth:`getDefaultParameters`
'rotate': False,
'verbose': 0,
'graphics': 1,
'anim': 2,
'anim': False,
'debug': True
}
return defaultParams
......
......@@ -27,7 +27,7 @@ class ADMM(Algorithm):
lagrange_mult_update_class = getattr(algorithms, lagrange_mult_update_name)
self.lagrange_mult_update = lagrange_mult_update_class(experiment)
else:
raise('Lagrange multiplier update not defined')
raise AttributeError('Lagrange multiplier update not defined')
def evaluate(self, u):
......
......@@ -47,6 +47,8 @@ class AvP(Algorithm):
proxClass = getattr(proxoperators, 'Prox_product_space')
self.p_product_space = proxClass(exp_AvP)
self.n_product_Prox = len(experiment.proxOperators)
else:
self.n_product_Prox = experiment.n_product_Prox
def evaluate(self, u):
......
......@@ -4,10 +4,6 @@ from proxtoolbox.algorithms.feasibilityIterateMonitor import FeasibilityIterateM
from proxtoolbox.utils.cell import Cell, isCell
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from matplotlib.pyplot import subplots, show, figure
class CT_IterateMonitor(IterateMonitor):
#class CT_IterateMonitor(FeasibilityIterateMonitor):
......@@ -19,51 +15,14 @@ class CT_IterateMonitor(IterateMonitor):
It would be more appropriate to derive from
FeasibilityIterateMonitor, but this would
result in too many calculations
Was used for plotting. Now does nothing, but this
could change...
"""
def __init__(self, experiment):
super(CT_IterateMonitor, self).__init__(experiment)
self.figure_width = experiment.figure_width
self.figure_height = experiment.figure_height
self.figure_dpi = experiment.figure_dpi
self.figure = None
self.ax = None
self.colorbar = None
def displayProgress(self, alg):
"""
Display progress. This method is called after
UpdateStatistics().
Parameters
----------
alg : instance of Algorithm class
Algorithm being monitored.
"""
# the output is a picture in a column vector
u = alg.u_new
if isCell(u):
u = u[0]
elif u.ndim == 2:
u = u[:,0]
N = int(np.sqrt(u.shape[0]))
image = np.reshape(u, (N,N), order='F')
if self.figure == None:
# create figure
plt.ion()
self.figure = plt.figure(figsize = (self.figure_width, self.figure_height),
dpi = self.figure_dpi)
self.ax = self.figure.add_subplot(111)
pass
im = self.ax.imshow(image, cmap='gray')
if self.colorbar is not None:
self.colorbar.remove()
self.colorbar = plt.colorbar(im, ax=self.ax, fraction=0.046, pad=0.04)
title = "Iteration " + str(alg.iter)
self.ax.set_title(title)
self.figure.canvas.draw()
self.figure.canvas.flush_events()
from proxtoolbox.algorithms.AvP import AvP
from proxtoolbox import proxoperators
from proxtoolbox.proxoperators.P_diag import P_diag
from proxtoolbox.proxoperators.Prox_product_space import Prox_product_space
from proxtoolbox.utils.cell import Cell, isCell
from proxtoolbox.utils.size import size_matlab
from numpy import mod, zeros
from numpy.linalg import norm
import copy
from numpy import exp
class DyRePr(AvP):
"""
Dynamically Reweighted Prox mappings as presented in
Luke,Burke,Lyon "OPTICAL WAVEFRONT RECONSTRUCTION: THEORY AND
NUMERICAL METHODS", SIAM Review, 2002, where it is
described as extended least squares: a gradient
descent method for minimizing the log likelihood
objective of the sum of squared distances to sets.
Based on Matlab code written by Russell Luke (Inst. Fuer
Numerische und Angewandte Mathematik, Universitaet
Gottingen) on Aug.18, 2017.
"""
def __init__(self, experiment, iterateMonitor, accelerator):
super(DyRePr, self).__init__(experiment, iterateMonitor, accelerator)
def evaluate(self, u):
"""
Update of the Averaged Projections algorithm
Parameters
----------
u : ndarray or a list of ndarray objects
The current iterate.
Returns
-------
u_new : ndarray or a list of ndarray objects
The new iterate (same type as input parameter `u`).
Notes
-----
We try to figure out how the user is calling this. Either
(1) the user knows that averaged projections is alternating
projectons on the product space with one of the sets being
the diagonal, or (2) the user hasn't preprocessed this
and has just passed the list of prox mappings and left it to us
to rearrange and average. If it is case (1) one of the
prox mappings will be P_Diag, otherwise we will assume
that we are in case (2)
"""
if self.has_p_diag is False:
# Case (2): the user hasn't preprocessed this
# and has just passed the list of prox mappings
# and left it to us to rearrange and average.
u_intermed = Cell(self.n_product_Prox)
for j in range(self.n_product_Prox):
u_intermed[j] = u
u_intermed = self.p_product_space.eval(u_intermed)
if self.accelerator is not None:
u_intermed = self.accelerator.evaluate(u_intermed, self)
for j in range(self.n_product_Prox):
weight = norm(u_intermed[j], ord=2)**2
u_intermed[j] /= weight + 1
u_new = self.p_avg.eval(u_intermed)
else:
# Case (1) the user knows that averaged projections is alternating
# projectons on the product space with one of the sets being
# the diagonal.
# First, evaluate the other prox op
u_tmp = self.proxOperators[self.other_prox_index].eval(u)
for j in range(self.n_product_Prox):
if isCell(u):
weight = norm(u[j]-u_tmp[j], ord=2)**2
u_tmp[j] /= weight + 1
else:
m, n, _p, _q = size_matlab(u)
if m == self.n_product_Prox:
weight = norm(u[j,:]-u_tmp[j,:], ord=2)**2
u_tmp[j,:] /= weight + 1
elif n == self.n_product_Prox:
weight = norm(u[:,j]-u_tmp[:,j], ord=2)**2
u_tmp[:,j] /= weight + 1
else:
weight = norm(u[:,:,j]-u_tmp[:,:,j], ord=2)**2
u_tmp[:,:,j] /= weight + 1
# Finally, evaluate diag prox
u_new = self.proxOperators[self.p_diag_index].eval(u_tmp)
return u_new
import numpy as np
from numpy import zeros, isreal, reshape, real, imag
from proxtoolbox.utils.size import size_matlab
from proxtoolbox.utils.cell import Cell, isCell
import importlib
import copy
class GenericAccelerator:
"""
Generic Nesterov-type accelerator for first-order methods
Based on Matlab code written by Russell Luke (Inst. Fuer
Numerische und Angewandte Mathematik, Universitaet
Gottingen) on April 5, 2018.
"""
def __init__(self, experiment):
self.u_old = None
def evaluate(self, u, alg):
"""
Evaluation of generic Nesterov-type accelerator
for first-order methods
Parameters
----------
u : ndarray or a cell of ndarray objects
The point to be accelerated.
alg : Algorithm object
The algorithm that uses this accelerator
Returns
-------
u_new : ndarray or a cell of ndarray objects
(same type as input parameter u)
The step generated by this acceleration method
(same type as input parameter `u`).
"""
if self.u_old is None:
self.u_old = copy.deepcopy(u) # use deep copy because we may have a cell of arrays
s_k = (alg.iter)/(alg.iter+3)
if isCell(u):
len_u = len(u)
u_tmp = Cell(len_u)
for j in range(len_u):
u_tmp[j] = u[j].copy()
u[j] = (s_k+1)*u_tmp[j] - s_k*self.u_old[j]
else:
u_tmp = u.copy()
u = (s_k+1)*u_tmp - s_k*self.u_old
self.u_old = u_tmp
return u
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