Commit e111433a authored by s.gretchko's avatar s.gretchko
Browse files

Added output saving (.json and .mat format)

parent 365d6dea
......@@ -3,6 +3,6 @@
import SetProxPythonPath
from proxtoolbox.experiments.phase.JWST_Experiment import JWST_Experiment
JWST = JWST_Experiment(algorithm = 'RAAR', beta_0 = 1.0, beta_max = 1.0, MAXIT = 1000)
JWST = JWST_Experiment(algorithm = 'RAAR', lambda_0 = 1.0, lambda_max = 1.0, MAXIT = 1000)
JWST.run()
JWST.show()
......@@ -2,6 +2,6 @@
import SetProxPythonPath
from proxtoolbox.experiments.phase.JWST_Experiment import JWST_Experiment
JWST = JWST_Experiment(algorithm = 'DRAP', beta_0 = 0.65, beta_max = 0.65)
JWST = JWST_Experiment(algorithm = 'DRAP', lambda_0 = 0.65, lambda_max = 0.65)
JWST.run()
JWST.show()
......@@ -3,6 +3,6 @@
import SetProxPythonPath
from proxtoolbox.experiments.phase.JWST_Experiment import JWST_Experiment
JWST = JWST_Experiment(algorithm = 'HPR', beta_0 = 0.85, beta_max = 0.85, MAXIT = 500)
JWST = JWST_Experiment(algorithm = 'HPR', lambda_0 = 0.85, lambda_max = 0.85, MAXIT = 500)
JWST.run()
JWST.show()
......@@ -3,6 +3,6 @@
import SetProxPythonPath
from proxtoolbox.experiments.phase.JWST_Experiment import JWST_Experiment
JWST = JWST_Experiment(algorithm = 'RAAR', beta_0 = 0.85, beta_max = 0.85, MAXIT = 1000)
JWST = JWST_Experiment(algorithm = 'RAAR', lambda_0 = 0.85, lambda_max = 0.85, MAXIT = 1000)
JWST.run()
JWST.show()
......@@ -31,10 +31,10 @@ class DRAP(Algorithm):
The new iterate.
"""
beta = self.computeBeta()
lmbd = self.computeRelaxation()
if not isinstance(u, list):
tmp1 = self.prox2.eval(u)
tmp2 = beta*(tmp1 - u)
tmp2 = lmbd*(tmp1 - u)
tmp3 = self.prox1.eval(tmp1 + tmp2)
u_new = tmp3 - tmp2
else:
......@@ -43,7 +43,7 @@ class DRAP(Algorithm):
tmp3 = []
u_new = []
for tmp1_j, u_j in zip(tmp1, u):
tmp2_j = beta*(tmp1_j - u_j)
tmp2_j = lmbd*(tmp1_j - u_j)
tmp3_j = tmp1_j + tmp2_j
tmp2.append(tmp2_j)
tmp3.append(tmp3_j)
......@@ -54,4 +54,4 @@ class DRAP(Algorithm):
return u_new
def getDescription(self):
return self.getDescriptionHelper("\\beta", self.beta_0, self.beta_max)
\ No newline at end of file
return self.getDescriptionHelper("\\lambda", self.lambda_0, self.lambda_max)
\ No newline at end of file
......@@ -31,23 +31,23 @@ class HPR(Algorithm):
The new iterate.
"""
beta = self.computeBeta()
lmbd = self.computeRelaxation()
if not isinstance(u, list):
tmp = self.prox2.eval(u)
tmp2 = (2+beta-1)*tmp - u
u_new = (2*self.prox1.eval(tmp2) - tmp2 + u + (1-beta)*tmp)/2
tmp2 = (2+lmbd-1)*tmp - u
u_new = (2*self.prox1.eval(tmp2) - tmp2 + u + (1-lmbd)*tmp)/2
else:
tmp = self.prox2.eval(u)
tmp2 = []
u_new = []
for j in range(len(tmp)):
tmp2_j = (2+beta-1)*tmp[j]- u[j]
tmp2_j = (2+lmbd-1)*tmp[j]- u[j]
tmp2.append(tmp2_j)
tmp3 = self.prox1.eval(tmp2)
for j in range(len(tmp)):
u_new_j = (2*tmp3[j]-tmp2[j] + u[j] +(1-beta)*tmp[j])/2
u_new_j = (2*tmp3[j]-tmp2[j] + u[j] +(1-lmbd)*tmp[j])/2
u_new.append(u_new_j)
return u_new
def getDescription(self):
return self.getDescriptionHelper("\\beta", self.beta_0, self.beta_max)
\ No newline at end of file
return self.getDescriptionHelper("\\lambda", self.lambda_0, self.lambda_max)
\ No newline at end of file
......@@ -33,12 +33,12 @@ class RAAR(Algorithm):
(same type as input parameter u)
The new iterate.
"""
beta = self.computeBeta()
lmbd = self.computeRelaxation()
if not isinstance(u, list):
tmp1 = 2*self.prox2.eval(u) - u # Reflection using prox2 on u
tmp2 = self.prox1.eval(tmp1) # Projection of the reflection
# update
u_new = (beta*(2*tmp2 - tmp1) + (1 - beta)*tmp1 + u)/2
u_new = (lmbd*(2*tmp2 - tmp1) + (1 - lmbd)*tmp1 + u)/2
else:
u_new = self.prox2.eval(u)
tmp1 = []
......@@ -48,8 +48,8 @@ class RAAR(Algorithm):
tmp2 = self.prox1.eval(tmp1)
# update
for j in range(len(u)):
u_new[j] = (beta*(2*tmp2[j]-tmp1[j]) + (1-beta)*tmp1[j] + u[j])/2
u_new[j] = (lmbd*(2*tmp2[j]-tmp1[j]) + (1-lmbd)*tmp1[j] + u[j])/2
return u_new
def getDescription(self):
return self.getDescriptionHelper("\\beta", self.beta_0, self.beta_max)
\ No newline at end of file
return self.getDescriptionHelper("\\lambda", self.lambda_0, self.lambda_max)
\ No newline at end of file
......@@ -34,9 +34,9 @@ class Algorithm:
self.maxIter = experiment.MAXIT
self.tol = experiment.TOL
self.beta_0 = experiment.beta_0
self.beta_max = experiment.beta_max
self.beta_switch = experiment.beta_switch
self.lambda_0 = experiment.lambda_0
self.lambda_max = experiment.lambda_max
self.lambda_switch = experiment.lambda_switch
self.iterateMonitor = iterateMonitor
self.accelerator = accelerator
......@@ -108,7 +108,7 @@ class Algorithm:
The number of iterations the algorithm performed
additional entries given by the iterate monitor
"""
output = {'u': self.u, 'iter': self.iter}
output = {'stats': {'iter': self.iter}, 'u': self.u}
self.iterateMonitor.postprocess(self, output)
return output
......@@ -146,7 +146,7 @@ class Algorithm:
"""
Function returning a string giving the algorithm name
and optionally a list of parameters that characterize
this algorithm (i.e., beta_0, and beta_max constants).
this algorithm (i.e., lambda_0, and lambda_max constants).
This string is used for graphical output.
The default implementation only returns the name of the
class implementing this algorithm. Derived classes may
......@@ -168,23 +168,23 @@ class Algorithm:
+ " to " + str(param_max)
return desc
def computeBeta(self):
def computeRelaxation(self):
"""
Compute the relaxation parameters beta based on the
Compute the relaxation parameter based on the
current iteration count.
Returns
-------
beta: Float64
lmbda: Float64
Relaxation parameter.
"""
iter = self.iter + 1 # add 1 to conform with matlab version
# gradually changes beta from beta0 into beta_max
a = -iter/self.beta_switch
# gradually changes lambda from lambda_0 into lambda_max
a = -iter/self.lambda_switch
a *= a*a # compute a^3
s = exp(a)
beta = s*self.beta_0 + (1-s)*self.beta_max
return beta
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'):
......
......@@ -131,7 +131,7 @@ class FeasibilityIterateMonitor(IterateMonitor):
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
# even for beta=1.
# even for lambda=1.
# (see Bauschke-Combettes-Luke, J. Approx. Theory, 2004)
if isinstance(u1, list):
for jj in range(len(u1)):
......@@ -260,6 +260,8 @@ class FeasibilityIterateMonitor(IterateMonitor):
"""
output = super(FeasibilityIterateMonitor, self).postprocess(alg, output)
if self.diagnostic:
output['gaps'] = self.gaps[1:alg.iter+1]
output['shadow_changes'] = self.shadow_changes[1:alg.iter+1]
stats = output['stats']
stats['gaps'] = self.gaps[1:alg.iter+1]
stats['shadow_changes'] = self.shadow_changes[1:alg.iter+1]
stats['rel_errors'] = self.rel_errors[1:alg.iter+1]
return output
\ No newline at end of file
......@@ -136,7 +136,8 @@ class IterateMonitor:
output['u1'] = u1
output['u2'] = u2
output['changes'] = self.changes[1:alg.iter+1]
if 'stats' in output:
output['stats']['changes'] = self.changes[1:alg.iter+1]
return output
@staticmethod
......
#SG
from proxtoolbox.experiments.dataProcessor import DataProcessor
from proxtoolbox import algorithms
......@@ -9,7 +8,11 @@ import sys
import os
import importlib
import time
from datetime import datetime
import json
import scipy
import numpy as np
from numpy import square, mod
from numpy.linalg import norm
from math import sqrt
......@@ -28,6 +31,23 @@ class ExperimentMetaClass(type):
obj.initialize() # initialize experiment object
return obj
class JSONNumpyEncoder(json.JSONEncoder):
""" json encoder for numpy types """
def default(self, obj): # pylint: disable=E0202
if isinstance(obj, (np.int_, np.intc, np.intp, np.int8,
np.int16, np.int32, np.int64, np.uint8,
np.uint16, np.uint32, np.uint64)):
return int(obj)
elif isinstance(obj, (np.float_, np.float16, np.float32,
np.float64)):
return float(obj)
elif isinstance(obj, complex):
return [obj.real, obj.imag]
elif isinstance(obj,(np.ndarray,)):
return obj.tolist()
return json.JSONEncoder.default(self, obj)
class Experiment(metaclass = ExperimentMetaClass):
'''
......@@ -61,9 +81,9 @@ class Experiment(metaclass = ExperimentMetaClass):
numruns = 1,
MAXIT = 1000,
TOL = 1e-8,
beta_0 = 0.85, # Initial beta value for first iterations
beta_max = 0.85, # Final beta value for later iterations
beta_switch = 20, # Timescale over which to change beta
lambda_0 = 0.85, # Initial relaxation value for first iterations
lambda_max = 0.85, # Final relaxation value for later iterations
lambda_switch = 20, # Timescale over which to change relaxation parameter
data_ball = 1e-3,
diagnostic = True,
iterate_monitor_name = 'IterateMonitor',
......@@ -77,6 +97,9 @@ class Experiment(metaclass = ExperimentMetaClass):
figure_height = 7, # figure height in inches
figure_dpi = 100, # figure dpi for display
debug = False,
save_output = True, # Set to true for tests, but default should be false
json_output = True,
matlab_output = True,
**kwargs):
"""
......@@ -115,9 +138,14 @@ class Experiment(metaclass = ExperimentMetaClass):
self.numruns = numruns
self.MAXIT = MAXIT
self.TOL = TOL
self.beta_0 = beta_0
self.beta_max = beta_max
self.beta_switch = beta_switch
self.lambda_0 = lambda_0
self.lambda_max = lambda_max
self.lambda_switch = lambda_switch
self.save_output = save_output
self.json_output = json_output
self.matlab_output = matlab_output
self.iterate_monitor_name = iterate_monitor_name
self.accelerator_name = accelerator_name
......@@ -233,14 +261,17 @@ class Experiment(metaclass = ExperimentMetaClass):
print ("Running " + self.algorithm_name + " on " + self.experiment_name + "...")
t = time.time()
self.output = self.algorithm.run(self.data_processor.u0)
self.output['time'] = time.time() - t
self.output = self.algorithm.run(self.data_processor.u0) # run algorithm
if 'stats' in self.output:
self.output['stats']['time'] = time.time() - t
print("Took", self.output['iter'], "iterations and",
self.output['time'], "seconds.")
print("Took", self.output['stats']['iter'], "iterations and",
self.output['stats']['time'], "seconds.")
# TODO save results
if self.save_output:
self.saveOutput()
def setupProxOperators(self):
"""
......@@ -274,7 +305,56 @@ class Experiment(metaclass = ExperimentMetaClass):
self.proxOperators.append(None)
self.proxOperators.append('Approx_PM_Poisson') # Patrick: This is just to monitor the change of phases!
def saveOutput(self):
# Create directory (if needed)
# does not raise an exception
# if the directory already exists
date_time = datetime.now() # current date and time
date_str = date_time.strftime("%Y_%m_%d")
directory = '../OutputData/' + date_str + '/'
os.makedirs(directory, exist_ok=True)
filename_root = self.experiment_name + '_' + self.algorithm_name + '.output'
file_path = directory + filename_root
# prepare a special version of the output before saving it
output = {}
for key in self.output:
if key == 'u1' or key == 'u2':
pass
elif key == 'u_monitor':
# replace list with array
length = len(self.output[key])
u_mon = np.zeros(length, dtype=np.object)
for i in range(length):
u_mon[i] = self.output[key][i]
output[key] = u_mon
else:
output[key] = self.output[key]
if self.matlab_output:
# save .mat file
file_path_ext = file_path + '.mat'
# wrap all the variables in a dictionary ('python_output') as
# in ProxMatlab, give it a different name to avoid name collision
# with ProxMatlab
scipy.io.savemat(file_path_ext, {'python_output': output})
if self.json_output:
# save json file
file_path_ext = file_path + '.json'
try:
with open(file_path_ext, 'w') as file:
json.dump(output, file, indent=4, cls=JSONNumpyEncoder)
except IOError as e:
errMsg = e.__str__()
print("Unable to write json output file ", file_path_ext, " : ", errMsg)
except TypeError as e:
errMsg = e.__str__()
print("Unable to serialize object for json output file: ", errMsg)
def estimateGap(self):
"""
Estimate the gap in the sequential formulation.
......@@ -348,9 +428,9 @@ class Experiment(metaclass = ExperimentMetaClass):
'product_space_dimension': 4,
'MAXIT': 6000,
'TOL': 5e-5,
'beta_0': 0.85,
'beta_max': 0.85,
'beta_switch': 20,
'lambda_0': 0.85,
'lambda_max': 0.85,
'lambda_switch': 20,
'data_ball': 1e-30,
'diagnostic': True,
'iterate_monitor_name': 'FeasibilityIterateMonitor',
......
......@@ -36,9 +36,9 @@ class CDI_Experiment(PhaseExperiment):
'farfield': True,
'MAXIT': 6000,
'TOL': 1e-8,
'beta_0': 0.5,
'beta_max': 0.50,
'beta_switch': 30,
'lambda_0': 0.5,
'lambda_max': 0.50,
'lambda_switch': 30,
'data_ball': .999826e-30,
'diagnostic': True,
'iterate_monitor_name': 'FeasibilityIterateMonitor',
......
......@@ -36,9 +36,9 @@ class CDP_Experiment(PhaseExperiment):
'farfield': True,
'MAXIT': 10000,
'TOL': 1e-10,
'beta_0': 0.55,
'beta_max': 0.55,
'beta_switch': 20,
'lambda_0': 0.55,
'lambda_max': 0.55,
'lambda_switch': 20,
'data_ball': 1e-30,
'diagnostic': True,
'iterate_monitor_name': 'FeasibilityIterateMonitor',
......
......@@ -27,9 +27,9 @@ class JWST_Experiment(PhaseExperiment):
'product_space_dimension': 4,
'MAXIT': 6000,
'TOL': 5e-5,
'beta_0': 0.85,
'beta_max': 0.85,
'beta_switch': 20,
'lambda_0': 0.85,
'lambda_max': 0.85,
'lambda_switch': 20,
'data_ball': 1e-30,
'diagnostic': True,
'iterate_monitor_name': 'FeasibilityIterateMonitor',
......
......@@ -40,9 +40,9 @@ class PhaseExperiment(Experiment):
u2 = u2[:,:,0]
changes = self.output['changes']
if 'time' in self.output:
time = self.output['time']
changes = self.output['stats']['changes']
if 'time' in self.output['stats']:
time = self.output['stats']['time']
else:
time = 999
m = u2.shape[0]
......@@ -105,8 +105,8 @@ class PhaseExperiment(Experiment):
f.suptitle(title)
#figure 902
if 'gaps' in self.output:
gaps = self.output['gaps']
if 'gaps' in self.output['stats']:
gaps = self.output['stats']['gaps']
f = plt.figure(figsize = (self.figure_width, self.figure_height))
plt.semilogy(gaps)
plt.xlabel(label)
......
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