Commit 6d6c5a08 authored by s.gretchko's avatar s.gretchko
Browse files

Added source localization experiment (feasibility part)

parent 262f00de
from proxtoolbox.experiments.experiment import Experiment
from proxtoolbox import proxoperators
from proxtoolbox.utils.cell import Cell, isCell
import numpy as np
from numpy.random import rand
from numpy.linalg import norm
import matplotlib
import matplotlib.pyplot as plt
from matplotlib.pyplot import subplots, show, figure
class SourceLocExperiment(Experiment):
'''
Source localization experiment class
'''
@staticmethod
def getDefaultParameters():
defaultParams = {
'experiment_name': 'source location',
'instance': 'simple random',
'object': 'real',
'constraint': 'diagonal_and_balls',
'formulation': 'feasibility',
'product_space_dimension': 1,
'sensors': 3,
'phys_dimension': 3,
'phys_boundary': 100,
'radius_bound': 100,
'noise': False,
'MAXIT': 10000,
'TOL': 1e-11,
'lambda_0': 0.33,
'lambda_max': 0.33,
'lambda_switch': 1,
'diagnostic': True,
'verbose': 0,
'graphics': 1,
'anim': False
}
return defaultParams
def __init__(self,
instance='simple random',
sensors=3,
phys_dimension=3,
phys_boundary=100,
radius_bound=100,
**kwargs):
"""
"""
# call parent's __init__ method
super(SourceLocExperiment, self).__init__(**kwargs)
# do here any data member initialization
self.instance = instance
self.sensors = sensors
self.phys_dimension = phys_dimension
self.phys_boundary = phys_boundary
self.radius_bound = radius_bound
# internal data members
self.shift_data = None # default receiver location
def loadData(self):
"""
Load source localization dataset. Create the initial iterate.
"""
self.Nx = 1
self.Ny = self.phys_dimension
# old code
#self.Nx = self.phys_dimension
#self.Ny = self.sensors
self.Nz = 1
# we load data depending on the instance parameter
if self.instance == 'simple random':
dimension = self.phys_dimension
sensors = self.sensors
# generate the true position of the source
self.truth = 2*self.phys_boundary*(rand(1, dimension) - 0.5)
self.truth_dim = self.truth.shape
self.norm_truth = norm(self.truth)
# generate the sensors
# each row corresponds to a sensor
# (this reshaping procedure is done so that we get the same
# 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
self.shift_data = Cell(sensors)
for i in range(sensors):
self.shift_data[i] = sensorMatrix[i]
self.receiver_data_scale = self.phys_boundary*sensors*dimension
# calculate the distances from the source to the sensors
self.data = np.zeros((sensors))
for i in range(sensors):
self.data[i] = norm(self.shift_data[i] - self.truth)
if self.noise:
self.data[i] += 20*(rand(1)-0.5)
self.norm_data = norm(self.data)
# initial guess
self.u0 = self.truth.flatten() + 0.1*self.phys_boundary*(rand(dimension)-0.5)
else:
errMsg = "Loading data: instance " + self.instance \
+ " is not yet implemented"
raise NotImplementedError(errMsg)
def setupProxOperators(self):
"""
Determine the prox operators to be used for this 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':
self.iterate_monitor_name = 'FeasibilityIterateMonitor'
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')
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')
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')
def show(self):
stats = self.output['stats']
changes = stats['changes']
changesSubplot = 111
rel_errors = None
if 'rel_errors' in stats:
rel_errors = stats['rel_errors']
changesSubplot = 211
time = self.output['stats']['time']
time_str = "{:.{}f}".format(time, 5) # 5 is precision
iterations = np.arange(1, len(changes)+1, 1)
algo_desc = self.algorithm.getDescription()
algo_label = "Algorithm " + algo_desc + " (time = " + time_str + " s)"
title = algo_label
fig = plt.figure(figsize = (self.figure_width, self.figure_height),
dpi = self.figure_dpi)
fig.suptitle(title)
ax1 = fig.add_subplot(changesSubplot)
ax1.semilogy(iterations, changes)
ax1.set_xlabel('Iteration')
yLabel = r'$|\!| x^{k+1}-x^{k} |\!|$'
ax1.set_ylabel(yLabel)
ax1.set_title("Change in iterates")
if rel_errors is not None:
ax2 = fig.add_subplot(212)
ax2.semilogy(iterations, rel_errors)
ax2.set_xlabel('Iteration')
# yLabel = r'$|\!| P_1(x^{k})-P_2(x^{k}) |\!|$' # This is what Matlab does but this not right
yLabel = 'Relative error'
ax2.set_ylabel(yLabel)
ax2.set_title("Relative errors")
plt.subplots_adjust(hspace = 0.3) # adjust vertical space (height) between subplots (default = 0.2)
plt.show()
......@@ -25,3 +25,4 @@ from .Approx_Pphase_JWST_Wirt import *
from .P_CDP_ADMM import *
from .P_S import *
from .P_CDP_cyclic import *
from .sourceLocProx import *
from proxtoolbox.utils.cell import Cell, isCell
from proxtoolbox.proxoperators.proxoperator import ProxOperator
from proxtoolbox.utils.size import size_matlab
from proxtoolbox.utils.accessors import accessors
import numpy as np
from numpy.linalg import norm
class Project_on_sphere(ProxOperator):
"""
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):
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
def eval(self, u, prox_index = None):
"""
Projection of a vector on the product of balls with different radius
Parameters
----------
u : Cell or ndarray
Input data to be projected
prox_idx : int, optional
Index of this prox operator
Returns
-------
u_new : Cell or ndarray
The projection
"""
u_new = u.copy()
if prox_index is not None:
j = prox_index
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.
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