from proxtoolbox.algorithms.algorithm import Algorithm
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 numpy import mod
import copy
from numpy import exp
[docs]class AvP(Algorithm):
"""
Averaged Projections algorithm
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(AvP, self).__init__(experiment, iterateMonitor, accelerator)
# see if P_Diag is in the list of prox operators
# prox operators have already been instanciated by parent (Algorithm class)
p_diag_index = 0
self.has_p_diag = False
for prox in self.proxOperators:
if isinstance(prox, P_diag):
self.has_p_diag = True
self.p_diag_index = p_diag_index
self.other_prox_index = (p_diag_index+1) % 2
break
p_diag_index += 1
if self.has_p_diag is False:
# instantiate P_avg prox
proxClass = getattr(proxoperators, 'P_avg')
self.p_avg = proxClass(experiment)
# Instantiate Prox_product_space operator:
# The constructor of this class expects the experiment object
# to have a nonempty data member productProxOperators which is
# not the case here. Hence, we need to buid a copy of 'experiment'
# which will fit the requirements.
exp_AvP = copy.copy(experiment)
exp_AvP.productProxOperators = experiment.proxOperators
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
[docs] 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)
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 = self.proxOperators[self.other_prox_index].eval(u)
if self.accelerator is not None:
u = self.accelerator.evaluate(u, self)
# Finally, evaluate diag prox
u_new = self.proxOperators[self.p_diag_index].eval(u)
return u_new