Dear Gitlab users, due to maintenance reasons, Gitlab will not be available on Thursday 30.09.2021 from 5:00 pm to approximately 5:30 pm.

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

Merge branch 'master' of gitlab.gwdg.de:nam/ProxPython

parents 91fb2337 4eaa8342
......@@ -6,3 +6,4 @@ build
dist
.vscode/launch.json
.vscode/settings.json
tests/orb_tomog_regression_tests.py
\ No newline at end of file
......@@ -115,10 +115,10 @@ to run the demo called "demo_name".
For example,
``` bash
python3 JWST_RAAR.py
python3 JWST_DRl.py
```
will run the RAAR algortihm on the JWST experiment.
will run the DRl algortihm on the JWST experiment.
## Structure
......
import SetProxPythonPath
from proxtoolbox.experiments.orbitaltomography.degenerate_orbits import DegenerateOrbital
exp_pm = DegenerateOrbital(algorithm='CP',
constraint='sparse real',
sparsity_parameter=180)
# exp_pm.plotInputData()
exp_pm.run()
exp_pm.show()
import SetProxPythonPath
from proxtoolbox.experiments.orbitaltomography.orthogonal_orbits import OrthogonalOrbitals
exp_pm = OrthogonalOrbitals(algorithm='CP',
constraint='sparse real')
# exp_pm.plotInputData()
exp_pm.run()
exp_pm.show()
import SetProxPythonPath
from proxtoolbox.experiments.orbitaltomography.planar_molecule import PlanarMolecule
exp_pm = PlanarMolecule(algorithm='DRl',
constraint='sparse real',
sparsity_parameter=40)
# exp_pm.plotInputData()
exp_pm.run()
exp_pm.show()
from proxtoolbox.algorithms.algorithm import Algorithm
from proxtoolbox.algorithms.algorithm import Algorithm
class AP(Algorithm):
"""
......
......@@ -34,4 +34,4 @@ from .PHeBIE_ptychography_objective import *
from .PHeBIE_phase_objective import *
from .Wirtinger import *
from .Wirt_IterateMonitor import *
from .NSLS_sourceloc_objective import *
from .NSLS_sourceloc_objective import *
\ No newline at end of file
from numpy import exp
try:
from tqdm import tqdm, tqdm_notebook
load_tqdm = True
except ModuleNotFoundError:
load_tqdm = False
pass
from numpy import zeros, angle, trace, exp, sqrt, sum
from numpy.linalg import norm
class Algorithm:
"""
......@@ -9,7 +15,7 @@ class Algorithm:
The algorithm object is created when the Experiment object is initialized.
"""
def __init__(self, experiment, iterateMonitor, accelerator = None):
def __init__(self, experiment, iterateMonitor, accelerator=None):
# instantiate prox operators
self.proxOperators = []
......@@ -17,7 +23,7 @@ class Algorithm:
if proxClass != '':
prox = proxClass(experiment)
self.proxOperators.append(prox)
#for convenience
# for convenience
self.prox1 = self.proxOperators[0]
self.prox2 = self.proxOperators[1]
......@@ -40,6 +46,11 @@ class Algorithm:
self.iterateMonitor = iterateMonitor
self.accelerator = accelerator
if hasattr(experiment, 'progressbar'):
self.progressbar = experiment.progressbar
else:
self.progressbar = None
# temporaries for current iteration
self.iter = 0
self.u = None
......@@ -48,14 +59,13 @@ class Algorithm:
# for debugging
self.debug = experiment.debug
def preprocess(self):
"""
The default implementation calls the iterate monitor's
preprocess() method which initializes the arrays that
will be used to collect statistics.
"""
self.iterateMonitor.preprocess(self)
self.iterateMonitor.preprocess(self)
def stoppingCondition(self):
"""
......@@ -149,15 +159,26 @@ class Algorithm:
self.u = u
self.preprocess()
self.displayProgress()
while (self.stoppingCondition()):
if self.progressbar == 'tqdm_notebook':
pbar = tqdm_notebook(total=self.maxIter)
elif self.progressbar == 'tqdm':
pbar = tqdm(total=self.maxIter)
else:
pbar = None
while self.stoppingCondition():
self.iter += 1
if pbar is not None:
pbar.update()
self.u_new = self.evaluate(self.u)
self.updateStatistics()
self.displayProgress()
self.u = self.u_new
if pbar is not None:
pbar.close()
return self.postprocess()
def getDescription(self):
"""
Function returning a string giving the algorithm name
......@@ -176,7 +197,7 @@ class Algorithm:
desc = self.getDescriptionHelper()
return desc
def getDescriptionHelper(self, param_name = None, param_0 = None, param_max = None):
def getDescriptionHelper(self, param_name=None, param_0=None, param_max=None):
"""
Generate a string giving a short description of
the algorithm (the name of its class name and optionally,
......@@ -191,9 +212,9 @@ class Algorithm:
"""
desc = type(self).__name__
if param_name is not None and param_0 is not None \
and param_max is not None:
and param_max is not None:
desc += ", $" + param_name + "$ = " + str(param_0) \
+ " to " + str(param_max)
+ " to " + str(param_max)
return desc
def computeRelaxation(self):
......@@ -206,15 +227,10 @@ class Algorithm:
lmbda: Float64
Relaxation parameter.
"""
iter = self.iter + 1 # add 1 to conform with matlab version
iteration = self.iter + 1 # add 1 to conform with matlab version
# gradually changes lambda from lambda_0 into lambda_max
a = -iter/self.lambda_switch
a *= a*a # compute a^3
a = -iteration / self.lambda_switch
a *= a * a # compute a^3
s = exp(a)
lmbda = s*self.lambda_0 + (1-s)*self.lambda_max
lmbda = s * self.lambda_0 + (1 - s) * self.lambda_max
return lmbda
@staticmethod
def phase_offset_compensated_norm(arr1, arr2, norm_factor=1, norm_type='fro'):
phase_offset = angle(sum(arr1 * arr2.conj()))
return norm(arr1 - arr2 * exp(1j * phase_offset), norm_type) / norm_factor
from numpy import angle, trace, exp, sqrt, matmul, reshape
from numpy.linalg import norm
from proxtoolbox.algorithms.iterateMonitor import IterateMonitor
from proxtoolbox.utils.cell import Cell, isCell
from proxtoolbox.utils.size import size_matlab
import numpy as np
from numpy import zeros, angle, trace, exp, sqrt, sum, matmul, array, reshape
from numpy.linalg import norm
class FeasibilityIterateMonitor(IterateMonitor):
"""
......@@ -38,7 +39,7 @@ class FeasibilityIterateMonitor(IterateMonitor):
super(FeasibilityIterateMonitor, self).preprocess(alg)
nProx = len(self.proxOperators)
# Even if more detailed diagnostics are not run, the mathematical
# exit criteria will be on step-sizes, for which one needs to monitor the
# norm of the difference between successive iterates of the blocks, in
......@@ -53,10 +54,10 @@ class FeasibilityIterateMonitor(IterateMonitor):
# assumes that u0 is an array
self.u_monitor[0] = self.u0
for j in range(1, nProx):
self.u_monitor[j] = self.proxOperators[j].eval(self.u_monitor[j-1], j)
else: # always creates a u_monitor, issue is whether to monitor the intermediate
self.u_monitor[j] = self.proxOperators[j].eval(self.u_monitor[j - 1], j)
else: # always creates a u_monitor, issue is whether to monitor the intermediate
# steps for the gap, as the above prepares for in the case that
# diagnostic ist true.
# diagnostic is true.
self.u_monitor = self.u0
# set up diagnostic arrays
......@@ -104,22 +105,22 @@ class FeasibilityIterateMonitor(IterateMonitor):
# reveals the RELEVANT convergence of the (shadow sequence of the) algorithm:
if isCell(u):
for j in range(len(u)):
self.u_monitor[0][j] = exp(-1j*angle(trace(matmul(prev_u_mon[0][j].T.conj(), u[j]))))*u[j]
self.u_monitor[0][j] = exp(-1j * angle(trace(matmul(prev_u_mon[0][j].T.conj(), u[j])))) * u[j]
else:
# assumes that u is a ndarray
self.u_monitor[0] = exp(-1j*angle(trace(matmul(prev_u_mon[0].T.conj(), u))))*u
self.u_monitor[0] = exp(-1j * angle(trace(matmul(prev_u_mon[0].T.conj(), u)))) * u
u1 = self.proxOperators[nProx-1].eval(self.u_monitor[0]) # Note: to compute u1 which prox_index should we use?
# Matlab takes whatever is defined in method_input, which corresponds
u1 = self.proxOperators[nProx - 1].eval(self.u_monitor[0]) # Note: to compute u1 which prox_index should we use?
# Matlab takes whatever is defined in method_input, which corresponds
# to the last prox evaluation if it was set (not always the case)
# for Douglas Rachford and DRlambda the iterates themselves do not
# converge to the intersection, even if this is nonempty. The Shadows, however,
# do go to the intersection, if this is nonempty.
if self.diagnostic:
self.u_monitor[1] = self.proxOperators[0].eval(u1, 0)
for j in range(1, nProx-1):
self.u_monitor[j+1] = self.proxOperators[j].eval(self.u_monitor[j], j)
for j in range(1, nProx - 1):
self.u_monitor[j + 1] = self.proxOperators[j].eval(self.u_monitor[j], j)
# convert u1 to an array if necessary
if isCell(u1):
_m, _n, p, q = size_matlab(u1[len(u1)-1])
......@@ -133,9 +134,9 @@ class FeasibilityIterateMonitor(IterateMonitor):
if p == 1 and q == 1:
if isCell(self.u_monitor[0]):
for j in range(len(self.u_monitor[0])):
tmp_change = tmp_change + (norm(self.u_monitor[0][j] - prev_u_mon[0][j])/normM)**2
tmp_change = tmp_change + (norm(self.u_monitor[0][j] - prev_u_mon[0][j]) / normM) ** 2
else:
tmp_change = (norm(self.u_monitor[0] - prev_u_mon[0])/normM)**2
tmp_change = (norm(self.u_monitor[0] - prev_u_mon[0]) / normM) ** 2
if self.diagnostic:
# compute shadow and gap
tmp_shadow = tmp_change
......@@ -145,11 +146,11 @@ class FeasibilityIterateMonitor(IterateMonitor):
# (see Bauschke-Combettes-Luke, J. Approx. Theory, 2004)
if isCell(u1):
for jj in range(len(u1)):
tmp_shadow = tmp_shadow + (norm(self.u_monitor[1][jj] - prev_u_mon[1][jj])/normM)**2
tmp_gap = tmp_gap + (norm(self.u_monitor[1][jj] - u1[jj])/normM)**2
tmp_shadow = tmp_shadow + (norm(self.u_monitor[1][jj] - prev_u_mon[1][jj]) / normM) ** 2
tmp_gap = tmp_gap + (norm(self.u_monitor[1][jj] - u1[jj]) / normM) ** 2
else:
tmp_shadow = tmp_shadow + (norm(self.u_monitor[1] - prev_u_mon[1])/normM)**2
tmp_gap = tmp_gap + (norm(self.u_monitor[1] - u1)/normM)**2
tmp_shadow = tmp_shadow + (norm(self.u_monitor[1] - prev_u_mon[1]) / normM) ** 2
tmp_gap = tmp_gap + (norm(self.u_monitor[1] - u1) / normM) ** 2
for j in range(2, nProx):
# For Douglas-Rachford,in general it is appropriate to monitor the
# SHADOWS of the iterates, since in the convex case these converge
......@@ -157,21 +158,21 @@ class FeasibilityIterateMonitor(IterateMonitor):
# (see Bauschke-Combettes-Luke, J. Approx. Theory, 2004)
if isCell(self.u_monitor[j]):
for jj in range(len(self.u_monitor[j])):
tmp_shadow = tmp_shadow + (norm(self.u_monitor[j][jj] - prev_u_mon[j][jj])/normM)**2
tmp_gap = tmp_gap + (norm(self.u_monitor[j][jj] - self.u_monitor[j-1][jj])/normM)**2
tmp_shadow = tmp_shadow + (norm(self.u_monitor[j][jj] - prev_u_mon[j][jj]) / normM) ** 2
tmp_gap = tmp_gap + (norm(self.u_monitor[j][jj] - self.u_monitor[j - 1][jj]) / normM) ** 2
else:
tmp_shadow = tmp_shadow + (norm(self.u_monitor[j] - prev_u_mon[j])/normM)**2
tmp_gap = tmp_gap + (norm(self.u_monitor[j] - self.u_monitor[j-1])/normM)**2
tmp_gap = tmp_gap + (norm(self.u_monitor[j] - self.u_monitor[j-1])/normM)**2
# compute relative error
if self.truth is not None:
# convert u1 to an array if necessary
# TODO need to simplify this...
# TODO need to simplify this...
if isCell(u1):
z = u1[len(u1)-1]
z = u1[len(u1) - 1]
if self.truth_dim[0] == 1:
z = reshape(z, (1,len(z))) # we want a true matrix not just a vector
z = reshape(z, (1, len(z))) # we want a true matrix not just a vector
elif self.truth_dim[1] == 1:
z = reshape(z, (len(z),1)) # we want a true matrix not just a vector
z = reshape(z, (len(z), 1)) # we want a true matrix not just a vector
else:
if u1.ndim == 1 or self.formulation == 'cyclic':
z = u1
......@@ -182,14 +183,14 @@ class FeasibilityIterateMonitor(IterateMonitor):
z = reshape(z, (len(z),1)) # we want a true matrix not just a vector
else: # product space
if self.truth_dim[0] == 1:
z = u1[0,:]
z = reshape(z, (1,len(z))) # we want a true matrix not just a vector
z = u1[0, :]
z = reshape(z, (1, len(z))) # we want a true matrix not just a vector
elif self.truth_dim[1] == 1:
z = u1[:,0]
z = reshape(z, (len(z),1)) # we want a true matrix not just a vector
z = u1[:, 0]
z = reshape(z, (len(z), 1)) # we want a true matrix not just a vector
else:
if u1.ndim == 3:
z = u1[:,:,0]
z = u1[:, :, 0]
else:
z = u1
"""
......@@ -223,35 +224,41 @@ class FeasibilityIterateMonitor(IterateMonitor):
if isCell(self.u_monitor[0]):
for j in range(len(self.u_monitor[0])):
for jj in range(p):
tmp_change = tmp_change + (norm(self.u_monitor[0][j][:,:,jj] - prev_u_mon[0][j][:,:,jj])/normM)**2
tmp_change = tmp_change + (
norm(self.u_monitor[0][j][:, :, jj] - prev_u_mon[0][j][:, :, jj]) / normM) ** 2
if self.diagnostic:
tmp_shadow = tmp_change
for k in range(1, nProx):
# compute (||P_SP_Mx-P_Mx||/normM)^2:
for j in range(len(self.u_monitor[0])):
for jj in range(p):
tmp_gap = tmp_gap + (norm(self.u_monitor[k][j][:,:,jj] - self.u_monitor[k-1][:,:,jj])/normM)**2
tmp_shadow = tmp_shadow + (norm(self.u_monitor[k][j][:,:,jj] - prev_u_mon[k][j][:,:,jj])/normM)**2
tmp_gap = tmp_gap + (norm(
self.u_monitor[k][j][:, :, jj] - self.u_monitor[k - 1][:, :, jj]) / normM) ** 2
tmp_shadow = tmp_shadow + (norm(
self.u_monitor[k][j][:, :, jj] - prev_u_mon[k][j][:, :, jj]) / normM) ** 2
else:
for j in range(p):
tmp_change = tmp_change + (norm(self.u_monitor[0][:,:,j] - prev_u_mon[0][:,:,j])/normM)**2
tmp_change = tmp_change + (norm(self.u_monitor[0][:, :, j] - prev_u_mon[0][:, :, j]) / normM) ** 2
if self.diagnostic:
tmp_shadow = tmp_change
for k in range(1, nProx):
# compute (||P_SP_Mx-P_Mx||/normM)^2:
tmp_gap = tmp_gap + (norm(self.u_monitor[k][:,:,j] - self.u_monitor[k-1][:,:,j])/normM)**2
tmp_shadow = tmp_shadow + (norm(self.u_monitor[k][:,:,j] - prev_u_mon[k][:,:,j])/normM)**2
tmp_gap = tmp_gap + (
norm(self.u_monitor[k][:, :, j] - self.u_monitor[k - 1][:, :, j]) / normM) ** 2
tmp_shadow = tmp_shadow + (
norm(self.u_monitor[k][:, :, j] - prev_u_mon[k][:, :, j]) / normM) ** 2
if self.diagnostic and self.truth is not None:
z = u1[:,:,0]
rel_error = norm(self.truth - exp(-1j*angle(trace(matmul(self.truth.T.conj(), z)))) * z) / self.norm_truth
else: # cells of 4D arrays???
z = u1[:, :, 0]
rel_error = norm(
self.truth - exp(-1j * angle(trace(matmul(self.truth.T.conj(), z)))) * z) / self.norm_truth
else: # cells of 4D arrays???
pass
self.changes[alg.iter] = sqrt(tmp_change)
if self.diagnostic:
self.gaps[alg.iter] = sqrt(tmp_gap)
self.shadow_changes[alg.iter] = sqrt(tmp_shadow) # this is the Euclidean norm of the gap to
self.shadow_changes[alg.iter] = sqrt(tmp_shadow) # this is the Euclidean norm of the gap to
# the unregularized set. To monitor the Euclidean norm of the gap to the
# regularized set is expensive to calculate, so we use this surrogate.
# Since the stopping criteria is on the change in the iterates, this
......@@ -288,8 +295,8 @@ class FeasibilityIterateMonitor(IterateMonitor):
output = super(FeasibilityIterateMonitor, self).postprocess(alg, output)
if self.diagnostic:
stats = output['stats']
stats['gaps'] = self.gaps[1:alg.iter+1]
stats['shadow_changes'] = self.shadow_changes[1:alg.iter+1]
stats['gaps'] = self.gaps[1:alg.iter + 1]
stats['shadow_changes'] = self.shadow_changes[1:alg.iter + 1]
if self.truth is not None:
stats['rel_errors'] = self.rel_errors[1:alg.iter+1]
return output
\ No newline at end of file
stats['rel_errors'] = self.rel_errors[1:alg.iter + 1]
return output
from numpy import zeros, sqrt
from numpy.linalg import norm
from proxtoolbox import algorithms
from proxtoolbox.utils.cell import Cell, isCell
from proxtoolbox.utils.cell import isCell
from proxtoolbox.utils.size import size_matlab
from numpy import zeros, angle, trace, exp, sqrt, sum, mod
from numpy.linalg import norm
class IterateMonitor:
......@@ -14,8 +15,8 @@ class IterateMonitor:
def __init__(self, experiment):
self.u0 = experiment.u0
self.u_monitor = self.u0 # the part of the sequence that is being monitored
# put u0 as default
assert self.u0 is not None, 'No valid initial guess given'
self.u_monitor = self.u0 # the part of the sequence that is being monitored put u0 as default
self.isCell = isCell(self.u0)
self.norm_data = experiment.norm_data
self.truth = experiment.truth
......@@ -122,10 +123,10 @@ class IterateMonitor:
u_m = self.u_monitor
# Python code for [m,n,p] = size(u_m)
if isCell(u_m):
# in matlab this corresponded to a cell
# in matlab this corresponded to a cell
# here we assume that such a cell is m by n where m = 1
# so far this seems to be always the case
# the following tests attemp to match the ndarray case below
# the following tests attempt to match the ndarray case below
m = 1
n = len(u_m)
if n == self.n_product_Prox:
......@@ -161,11 +162,6 @@ class IterateMonitor:
output['stats']['changes'] = self.changes[1:alg.iter + 1]
return output
@staticmethod
def phase_offset_compensated_norm(arr1, arr2, norm_factor=1, norm_type='fro'):
phase_offset = angle(sum(arr1 * arr2.conj()))
return norm(arr1 - arr2 * exp(1j * phase_offset), norm_type) / norm_factor
def getIterateSize(self, u):
"""
Given an iterate, determine p and q parameters
......@@ -207,7 +203,7 @@ class IterateMonitor:
def calculateObjective(self, alg):
"""
Calculate objective value. The dafault implementation
Calculate objective value. The default implementation
uses the optimality monitor if it exists
Parameters
......
from proxtoolbox import algorithms
from proxtoolbox import proxoperators
from proxtoolbox.proxoperators.proxoperator import ProxOperator
from proxtoolbox.utils.cell import Cell, isCell
import json
import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib.pyplot import subplots, show, figure
import sys
import numpy as np
import os
import importlib
import scipy
import time
from datetime import datetime
import json
import scipy
import numpy as np
from math import sqrt
from numpy import square, mod
from numpy.linalg import norm
from math import sqrt
from proxtoolbox import algorithms
from proxtoolbox import proxoperators
class ExperimentMetaClass(type):
......@@ -105,7 +101,7 @@ class Experiment(metaclass=ExperimentMetaClass):
# (including what would have been generated with verbose mode)
graphics=0,
anim=False,
anim_step = 2, # Tell how frequently the animation is updated, i.e., every `anim_step` step(s)
anim_step=2, # Tell how frequently the animation is updated, i.e., every `anim_step` step(s)
anim_colorbar=False, # if True the animated image will include a colorbar (updated at each iteration)
graphics_display=None,
figure_width=9, # figure width in inches
......@@ -129,6 +125,8 @@ class Experiment(metaclass=ExperimentMetaClass):
self.object = object
self.formulation = formulation
self.output = {}
# Those data members will define the final dimensions
# of the dataset to be used in this experiment.
# At this stage we set them to the desired values
......@@ -139,7 +137,7 @@ class Experiment(metaclass=ExperimentMetaClass):
self.Ny = Ny
if Nz is not None:
self.Nz = Nz
else:
else: # TODO: Bad for experiments with 3D data, where Nz may be None to simply follow the data
self.Nz = 1
# we store here the original values specified by the user
self.desired_product_space_dimension = product_space_dimension
......@@ -180,11 +178,11 @@ class Experiment(metaclass=ExperimentMetaClass):
self.accelerator_name = accelerator_name
self.rotate = rotate # tells the iterate monitor to rotate the iterates
# to a fixed global reference
# to a fixed global reference
self.rnd_seed = rnd_seed
# additional attributes
self.data = None
self.data_sq = None
self.norm_data = None
......@@ -195,9 +193,9 @@ class Experiment(metaclass=ExperimentMetaClass):
self.truth = None
self.truth_dim = None
self.norm_truth = None
self.optimality_monitor = None
self.proj_iter = None # not sure what it does
self.proxOperators = []
......@@ -205,8 +203,8 @@ class Experiment(metaclass=ExperimentMetaClass):
self.productProxOperators = []
self.n_product_Prox = None # size of the product space
# TODO: which default: None, 0 or 1 ?
# see how it is used
# TODO: which default: None, 0 or 1 ?
# see how it is used
self.propagator = None
self.inverse_propagator = None
......@@ -217,7 +215,6 @@ class Experiment(metaclass=ExperimentMetaClass):
# for animation:
self.anim_figure = None
self.debug = debug
def initialize(self):
......@@ -235,8 +232,8 @@ class Experiment(metaclass=ExperimentMetaClass):
# load data
self.loadData()
self.reshapeData(self.Nx, self.Ny, self.Nz) # TODO need to revisit this
# - arguments not needed in particular
self.reshapeData(self.Nx, self.Ny, self.Nz) # TODO need to revisit this
# - arguments not needed in particular
if self.TOL2 is None:
self.TOL2 = 1e-20
......@@ -277,7 +274,7 @@ class Experiment(metaclass=ExperimentMetaClass):
t = time.time()
self.output = self.runAlgorithm(self.u0)
self.output = self.runAlgorithm(self.u0) # TODO: append to dict? (in overwrite mode)
self.output['stats']['time'] = time.time() - t
......@@ -286,7 +283,7 @@ class Experiment(metaclass=ExperimentMetaClass):
self.output['stats']['time'], "seconds.")
self.postprocess()
if self.save_output:
self.saveOutput()
......@@ -313,8 +310,8 @@ class Experiment(metaclass=ExperimentMetaClass):
dictionary containing the output of the algorithm. In
particular, it contains the last iterate and various
statistics.