Commit 284d8efb authored by s.gretchko's avatar s.gretchko
Browse files

Added all source localization demos

parent 6d6c5a08
......@@ -28,7 +28,9 @@ class AvP2(Algorithm):
def __init__(self, experiment, iterateMonitor, accelerator):
super(AvP2, self).__init__(experiment, iterateMonitor, accelerator)
if not hasattr(experiment, 'shift_data'):
if hasattr(experiment, 'shift_data'):
self.shift_data = experiment.shift_data
else:
L = self.product_space_dimension
self.shift_data = Cell(L)
for j in range(L):
......
......@@ -93,9 +93,12 @@ class AvP2_IterateMonitor(IterateMonitor):
#
# method_input.gap(iter) = 999; % just a dummy value.
# for now use dummy value
objValue = 999
return objValue
if self.optimality_monitor is not None:
return self.optimality_monitor.calculateObjective(alg)
else:
# for now use dummy value
objValue = 999
return objValue
def postprocess(self, alg, output):
"""
......
from proxtoolbox.utils.cell import Cell, isCell
from proxtoolbox.utils.size import size_matlab
import numpy as np
from numpy.linalg import norm
class NSLS_sourceloc_objective:
"""
objective function for the NSLS (NonSmooth Least Squares) formulation of the
source localization experiment
Based on Matlab code written by Russell Luke (Inst. Fuer
Numerische und Angewandte Mathematik, Universitaet
Gottingen) on Sept. 12, 2017.
"""
def __init__(self, experiment):
self.sensors = experiment.sensors
self.norm_data = experiment.norm_data
self.data = experiment.data
self.shift_data = experiment.shift_data
def calculateObjective(self, alg):
"""
Evaluate the objective function for the NSLS
Parameters
----------
alg : algorithm instance
The algorithm that is running
Returns
-------
objValue : real
The value of the objective function
"""
objValue = 0
u = alg.u_new[2]
if isCell(u):
for j in range (len(u)):
norm_vec = norm(self.shift_data[j] - u[j])
objValue += (norm_vec - self.data[j])**2
else:
errMsg = "Only cells are supported for now"
raise NotImplementedError(errMsg)
return objValue
......@@ -34,3 +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 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
......@@ -108,7 +109,10 @@ class FeasibilityIterateMonitor(IterateMonitor):
# assumes that u is a ndarray
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])
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.
......@@ -118,11 +122,9 @@ class FeasibilityIterateMonitor(IterateMonitor):
self.u_monitor[j+1] = self.proxOperators[j].eval(self.u_monitor[j], j)
# convert u1 to an array if necessary
if isCell(u1):
p, q = self.getIterateSize(u1[len(u1)-1])
#[~,~,p,q]=size(u1{length(u1)});
_m, _n, p, q = size_matlab(u1[len(u1)-1])
else:
p, q = self.getIterateSize(u1)
#[~,~,p,q]=size(u1);
_m, _n, p, q = size_matlab(u1)
# compute the normalized change in successive iterates:
# change(iter) = sum(sum((feval('P_M',M,u)-tmp).^2))/normM;
tmp_change = 0
......@@ -135,6 +137,7 @@ class FeasibilityIterateMonitor(IterateMonitor):
else:
tmp_change = (norm(self.u_monitor[0] - prev_u_mon[0])/normM)**2
if self.diagnostic:
# compute shadow and gap
tmp_shadow = tmp_change
# For Douglas-Rachford,in general it is appropriate to monitor the
# SHADOWS of the iterates, since in the convex case these converge
......@@ -159,6 +162,7 @@ class FeasibilityIterateMonitor(IterateMonitor):
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
# compute relative error
if self.truth is not None:
# convert u1 to an array if necessary
# TODO need to simplify this...
......@@ -169,10 +173,13 @@ class FeasibilityIterateMonitor(IterateMonitor):
elif self.truth_dim[1] == 1:
z = reshape(z, (len(z),1)) # we want a true matrix not just a vector
else:
if self.formulation == 'cyclic':
if u1.ndim == 1 or self.formulation == 'cyclic':
z = u1
if z.ndim == 1:
z = reshape(z, (len(z),1)) # we want a true matrix not just a vector
if self.truth_dim[0] == 1:
z = reshape(z, (1,len(z))) # we want a true matrix not just a vector
else:
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,:]
......@@ -211,7 +218,7 @@ class FeasibilityIterateMonitor(IterateMonitor):
zcov=fftshift(real(ifft2(fft2(z).*ifft2(z)))./method_input.Nx);
Relerrs=norm(zcov-method_input.cov_true_support);
"""
rel_error = norm(self.truth - exp(-1j*angle(trace(matmul(self.truth.T.conj(), z)))) * z) / self.norm_truth
rel_error = norm(self.truth - exp(-1j*angle(trace(matmul(self.truth.T.conj(), z)))) * z) / self.norm_truth
elif q == 1:
if isCell(self.u_monitor[0]):
for j in range(len(self.u_monitor[0])):
......
......@@ -136,19 +136,23 @@ class IterateMonitor:
u1 = u_m
u2 = u_m
else: # ndarray
m, n, p, _q = size_matlab(u_m)
if n == self.n_product_Prox:
u1 = u_m[:, 0]
u2 = u_m[:, n - 1]
elif m == self.n_product_Prox:
u1 = u_m[0, :]
u2 = u_m[m - 1, :]
elif p == self.n_product_Prox:
u1 = u_m[:, :, 0]
u2 = u_m[:, :, p - 1]
else:
if u_m.ndim == 1:
u1 = u_m
u2 = u_m
else:
m, n, p, _q = size_matlab(u_m)
if n == self.n_product_Prox:
u1 = u_m[:, 0]
u2 = u_m[:, n - 1]
elif m == self.n_product_Prox:
u1 = u_m[0, :]
u2 = u_m[m - 1, :]
elif p == self.n_product_Prox:
u1 = u_m[:, :, 0]
u2 = u_m[:, :, p - 1]
else:
u1 = u_m
u2 = u_m
output['u1'] = u1
output['u2'] = u2
......
......@@ -24,7 +24,7 @@ class SourceLocExperiment(Experiment):
'instance': 'simple random',
'object': 'real',
'constraint': 'diagonal_and_balls',
'formulation': 'feasibility',
'formulation': 'product space',
'product_space_dimension': 1,
'sensors': 3,
'phys_dimension': 3,
......@@ -91,7 +91,7 @@ class SourceLocExperiment(Experiment):
# random numbers as Matlab)
sensorMatrix = 2*self.phys_boundary*(rand(sensors*dimension) - 0.5)
sensorMatrix = np.reshape(sensorMatrix, (sensors, dimension), order='F')
# connvert to cells
# convert to cells
self.shift_data = Cell(sensors)
for i in range(sensors):
self.shift_data[i] = sensorMatrix[i]
......@@ -110,6 +110,13 @@ class SourceLocExperiment(Experiment):
+ " is not yet implemented"
raise NotImplementedError(errMsg)
if self.formulation == "product space":
u0 = self.u0
self.u0 = Cell(self.sensors)
for j in range(self.sensors):
self.u0[j] = u0.copy()
def setupProxOperators(self):
"""
......@@ -117,45 +124,71 @@ class SourceLocExperiment(Experiment):
"""
super(SourceLocExperiment, self).setupProxOperators() # call parent's method
# TODO should not use feasibility - to be consistent with other experiments
# we should use "cyclic" or "product space"
if self.formulation == 'feasibility':
# new code:
#
if self.formulation == 'cyclic':
self.iterate_monitor_name = 'FeasibilityIterateMonitor'
self.product_space_dimension = 1
self.proxOperators = []
self.productProxOperators = []
if self.algorithm_name == 'CP' or self.algorithm_name == 'CDRl':
if self.product_space_dimension == 1:
self.sets = self.sensors
self.nProx = self.sensors
for _j in range(self.nProx):
self.proxOperators.append('Project_on_sphere')
self.sets = self.sensors
self.nProx = self.sensors
for _j in range(self.nProx):
self.proxOperators.append('Project_on_sphere')
elif self.formulation == 'product space':
self.iterate_monitor_name = 'FeasibilityIterateMonitor'
self.product_space_dimension = self.sensors
self.sets = 2
self.nProx = 2
self.proxOperators.append('P_diag')
self.proxOperators.append('Prox_product_space')
# add product prox operators
self.n_product_Prox = self.product_space_dimension
for _j in range(self.n_product_Prox):
self.productProxOperators.append('Project_on_sphere')
elif self.formulation == 'NSLS':
self.optimality_monitor='NSLS_sourceloc_objective'
if self.algorithm_name == 'AvP2':
self.iterate_monitor_name = 'AvP2_IterateMonitor'
self.Prox = 2
self.product_space_dimension = self.sensors
self.n_product_Prox = self.sensors
if self.instance == 'JWST':
self.proxOperators.append('Project_on_JWST_diagonal')
self.proxOperators.append('Project_on_JWST_spheres')
elif self.instance == 'CDP':
self.proxOperators.append('Project_on_CDP_diagonal')
self.proxOperators.append('Project_on_CDP_spheres')
elif self.instance == 'Goettingen':
# not fully implemented
self.proxOperators.append('Project_on_shifted_diagonal')
self.proxOperators.append('Project_on_phase_spheres')
else:
# more sophisticated blocking strategies possible.
# We start with just one block and two prox mappings:
# the product space prox mapping and a projection
# onto the shifted diagonal. This is just averaged projections
# in thw case of CP and a kind of rotating RAAR/DRlambda in the
# case of CAARl.
self.sets = 2
self.product_space_dimension = self.sensors
self.nProx = 2
self.proxOperators.append('P_diag')
self.proxOperators.append('Prox_product_space')
# add product prox operators
self.n_product_Prox = self.product_space_dimension
for _j in range(self.n_product_Prox):
self.productProxOperators.append('Project_on_sphere')
self.proxOperators.append('Project_on_spheres')
u0 = Cell(self.sensors)
for j in range(self.sensors):
u0[j] = self.u0
self.u0 = Cell(3)
self.u0[0] = u0[0]
prox2_class = getattr(proxoperators, self.proxOperators[1])
prox2 = prox2_class(self)
self.u0[2] = prox2.eval(u0 - self.shift_data)
prox1_class = getattr(proxoperators, self.proxOperators[0])
prox1 = prox1_class(self)
u0 = prox1.eval(self.u0[2] + self.shift_data)
self.u0[1] = u0[0]
else:
# not sure there is a difference with code just above...
self.sets = 2
self.product_space_dimension = self.sets
self.nProx = 2
self.proxOperators.append('P_diag')
self.proxOperators.append('Prox_product_space')
# add product prox operators
self.n_product_Prox = self.product_space_dimension
for _j in range(self.n_product_Prox):
self.productProxOperators.append('Project_on_sphere')
errMsg = "formulation " + self.formulation \
+ " for algorithm ' + + ' is not supported"
raise NotImplementedError(errMsg)
else:
errMsg = "formulation " + self.formulation \
+ " is not supported"
raise NotImplementedError(errMsg)
def show(self):
......
......@@ -6,7 +6,21 @@ from proxtoolbox.utils.accessors import accessors
import numpy as np
from numpy.linalg import norm
class Project_on_sphere(ProxOperator):
class SourceLocProx(ProxOperator):
"""
Projection onto the given entries in a Sudoku problem
"""
def __init__(self, experiment):
self.sensors = experiment.sensors
self.data = experiment.data
self.shift_data = experiment.shift_data
self.Nx = experiment.Nx
self.Ny = experiment.Ny
self.product_space_dimension = experiment.product_space_dimension
class Project_on_sphere(SourceLocProx):
"""
Prox operator that projects a vector on the product of
balls with different radius
......@@ -17,12 +31,7 @@ class Project_on_sphere(ProxOperator):
"""
def __init__(self, experiment):
self.sensors = experiment.sensors
self.data = experiment.data
self.shift_data = experiment.shift_data
self.Nx = experiment.Nx
self.Ny = experiment.Ny
self.product_space_dimension = experiment.product_space_dimension
super(Project_on_sphere, self).__init__(experiment)
def eval(self, u, prox_index = None):
"""
......@@ -30,14 +39,14 @@ class Project_on_sphere(ProxOperator):
Parameters
----------
u : Cell or ndarray
u : ndarray
Input data to be projected
prox_idx : int, optional
Index of this prox operator
Returns
-------
u_new : Cell or ndarray
u_new : ndarray
The projection
"""
u_new = u.copy()
......@@ -46,17 +55,49 @@ class Project_on_sphere(ProxOperator):
else:
j = 0
#L = self.product_space_dimension
#get, _set, _data_shape = accessors(self.shift_data, self.Nx, self.Ny, L)
center = self.shift_data[j]
norm_u = norm(u - center)
u_new = (u - center) * self.data[j]/norm_u + center
return u_new
# if isCell(self.shift_data):
# center = self.shift_data[j]
# else:
# m, n, p, _q = size_matlab(self.shift_data)
# if m == self.sensors:
# center = self.
class Project_on_spheres(SourceLocProx):
"""
Prox operator that projects a vector on the product of
balls with different radius
Based on Matlab code written by Russell Luke (Inst. Fuer
Numerische und Angewandte Mathematik, Universitaet
Gottingen) on Sept. 8, 2016.
"""
def __init__(self, experiment):
super(Project_on_spheres, self).__init__(experiment)
def eval(self, u, prox_index = None):
"""
Projection of a vector on the product of balls with different radius
Parameters
----------
u : Cell
Input data to be projected
prox_idx : int, optional
Index of this prox operator
Returns
-------
u_new : Cell
The projection
"""
if isCell(u):
u_new = Cell(len(u))
for j in range (len(u)):
norm_u = norm(u[j])
u_new[j] = u[j] * self.data[j] / norm_u
else:
errMsg = "Only cells are supported for the input data to be projected"
raise NotImplementedError(errMsg)
return u_new
Supports Markdown
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