Commit 19274460 authored by rluke1's avatar rluke1
Browse files

git-svn-id:...

git-svn-id: https://svn.num.math.uni-goettingen.de/svn/russell/Regularity/ProxPython_sandbox@1913 f329657e-034e-462b-9bba-5e1f86a05649
parent 764e981b
#
import algorithms, problem, proxoperators
from sudoku import Sudoku
from ptychography import Ptychography
from ptychography import Ptychography_NTT_01_26210
__all__ = ["proxoperators","algorithms","problem","sudoku","ptychography"]
# -*- coding: utf-8 -*-
"""
Created on Wed Oct 1 10:40:18 2014
@author: stefan
"""
import numpy, scipy.linalg, math
__all__ = ["Algorithm","RAAR","HPR","AP"]
class Algorithm:
"""
Generic interface for ProxToolbox algorithms.
"""
def __init__(self):
raise NotImplementedError,"This is just an abstract interface"
def run(self, u, tol, maxiter):
"""
Runs the algorithm one some input data
Parameters
----------
u : array_like
Input data
tol : number
Tolerance
maxiter : int
Maximum number of iterations
Returns
-------
u1 : array_like
Result
u2 : array_like
Result
iters : int
Number of iterations the algorithm performed
change : array_like
The percentage change in the norm
gap : array_like
Squared gap distance normalized by the magnitude constraint
"""
raise NotImplementedError,"This is just an abstract interface"
class RAAR(Algorithm):
"""
Relaxed Averaged Alternating Reflection algorithm
See Also
--------
Algorithm : Generic interface for ProxToolbox algorithms
"""
def __init__(self,config):
"""
Parameters
----------
config : dict
Dictionary containing the problem configuration. It must contain the
following mappings:
proj1:ProxOperator
First ProxOperator (the class, no instance)
proj2:ProxOperator
Second ProxOperator (the class, no instance)
beta0:number
Starting relaxation parmater
beta_max:number
Maximum relaxation parameter
beta_switch:int
Iteration at which beta moves from beta0 -> beta_max
normM:number
?
Nx:int
Row-dim of the product space elements
Ny:int
Column-dim of the product space elements
Nz:int
Depth-dim of the product space elements
dim:int
Size of the product space
"""
self.proj1 = config['proj1'](config);
self.proj2 = config['proj2'](config);
self.normM = config['normM'];
self.beta0 = config['beta0'];
self.beta_max = config['beta_max'];
self.beta_switch = config['beta_switch'];
self.Nx = config['Nx']; self.Ny = config['Ny']; self.Nz = config['Nz'];
self.dim = config['dim'];
def run(self, u, tol, maxiter):
beta0 = self.beta0; beta_max = self.beta_max; beta_switch = self.beta_switch;
proj1 = self.proj1; proj2 = self.proj2;
normM = self.normM;
Nx = self.Nx; Ny = self.Ny; Nz = self.Nz; dim = self.dim;
beta = beta0;
iters = 0;
change = numpy.zeros(maxiter+1,dtype=u.dtype);
change[0] = 999;
gap = change.copy();
tmp1 = 2*proj2.work(u) - u;
while iters < maxiter and change[iters] >= tol:
tmp = math.exp((-iters/beta_switch)**3);
beta = (tmp*beta0) + ((1-tmp)*beta_max);
iters += 1;
tmp3 = proj1.work(tmp1);
tmp_u = ((beta*(2*tmp3-tmp1)) + ((1-beta)*tmp1) + u)/2;
tmp2 = proj2.work(tmp_u);
tmp3 = proj1.work(tmp2);
tmp_change = 0; tmp_gap = 0;
if Ny == 1 or Nx == 1:
tmp_change = (scipy.linalg.norm(u-tmp_u,'fro')/normM)**2;
tmp_gap = (scipy.linalg.norm(tmp3-tmp2,'fro')/normM)**2;
elif Nz == 1:
for j in xrange(dim):
tmp_change += (scipy.linalg.norm(u[:,:,j]-tmp_u[:,:,j],'fro')/normM)**2;
tmp_gap += (scipy.linalg.norm(tmp3[:,:,j]-tmp2[:,:,j])/normM,'fro')**2;
else:
for j in xrange(dim):
for k in xrange(Nz):
tmp_change += (scipy.linalg.norm(u[:,:,k,j]-tmp_u[:,:,k,j],'fro')/normM)**2;
tmp_gap += (scipy.linalg.norm(tmp3[:,:,k,j]-tmp2[:,:,k,j],'fro')/normM)**2;
change[iters] = math.sqrt(tmp_change);
gap[iters] = math.sqrt(tmp_gap);
u = tmp_u;
tmp1 = (2*tmp2) - tmp_u;
u = tmp2;
tmp = proj1.work(u);
tmp2 = proj2.work(u);
if Ny == 1:
u1 = tmp[:,1];
u2 = tmp2[:,1];
elif Nx == 1:
u1 = tmp[1,:];
u2 = tmp2[1,:];
elif Nz == 1:
u1 = tmp[:,:,1];
u2 = tmp2[:,:,1];
else:
u1 = tmp;
u2 = tmp2;
change = change[1:iters+1];
gap = gap[1:iters+1];
return u1, u2, iters, change, gap;
class HPR(Algorithm):
"""
Heinz-Patrick-Russell algorithm
See Also
--------
Algorithm : Generic interface for ProxToolbox algorithms
References
----------
Bauschke, Combettes & Luke, Journal of the Optical Society of America,
20(6):1025-1034, 2003
"""
def __init__(self, config):
"""
Parameters
----------
config : dict
Dictionary containing the problem configuration. It must contain the
following mappings:
proj1:ProxOperator
First ProxOperator (the class, no instance)
proj2:ProxOperator
Second ProxOperator (the class, no instance)
beta0:number
Starting relaxation parmater
beta_max:number
Maximum relaxation parameter
beta_switch:int
Iteration at which beta moves from beta0 -> beta_max
normM:number
?
Nx:int
Row-dim of the product space elements
Ny:int
Column-dim of the product space elements
Nz:int
Depth-dim of the product space elements
dim:int
Size of the product space
"""
self.proj1 = config['proj1'](config);
self.proj2 = config['proj2'](config);
self.normM = config['normM'];
self.beta0 = config['beta0'];
self.beta_max = config['beta_max'];
self.beta_switch = config['beta_switch'];
self.Nx = config['Nx']; self.Ny = config['Ny']; self.Nz = config['Nz'];
self.dim = config['dim'];
def run(self, u, tol, maxiter):
beta0 = self.beta0; beta_max = self.beta_max; beta_switch = self.beta_switch;
proj1 = self.proj1; proj2 = self.proj2;
normM = self.normM;
Nx = self.Nx; Ny = self.Ny; Nz = self.Nz; dim = self.dim;
beta = beta0;
iters = 0;
change = numpy.zeros(maxiter+1,dtype=u.dtype);
change[0] = 999;
gap = change.copy();
tmp = proj2.work(u);
while iters < maxiter and change[iters] >= tol:
a = math.exp((-iters/beta_switch)**3);
beta = (a*beta0) + ((1-a)*beta_max);
iters += 1;
tmp2 = 2*proj2.work(u) - u + (beta-1)*tmp;
tmp_u = (2*proj1.work(tmp2) - tmp2 + u + (1-beta)*tmp)/2;
tmp = proj2.work(tmp_u);
tmp3 = proj1.work(tmp);
tmp_change = 0; tmp_gap = 0;
if Ny == 1 or Nx == 1:
tmp_change = (scipy.linalg.norm(u-tmp_u,'fro')/normM)**2;
tmp_gap = (scipy.linalg.norm(tmp3-tmp2,'fro')/normM)**2;
elif Nz == 1:
for j in xrange(dim):
tmp_change += (scipy.linalg.norm(u[:,:,j]-tmp_u[:,:,j],'fro')/normM)**2;
tmp_gap += (scipy.linalg.norm(tmp3[:,:,j]-tmp2[:,:,j])/normM,'fro')**2;
else:
for j in xrange(dim):
for k in xrange(Nz):
tmp_change += (scipy.linalg.norm(u[:,:,k,j]-tmp_u[:,:,k,j],'fro')/normM)**2;
tmp_gap += (scipy.linalg.norm(tmp3[:,:,k,j]-tmp2[:,:,k,j],'fro')/normM)**2;
change[iters] = math.sqrt(tmp_change);
gap[iters] = math.sqrt(tmp_gap);
u = tmp_u;
tmp = proj1.work(u);
tmp2 = proj2.work(u);
if Ny == 1:
u1 = tmp[:,1];
u2 = tmp2[:,1];
elif Nx == 1:
u1 = tmp[1,:];
u2 = tmp2[1,:];
elif Nz == 1:
u1 = tmp[:,:,1];
u2 = tmp2[:,:,1];
else:
u1 = tmp;
u2 = tmp2;
change = change[1:iters+1];
gap = gap[1:iters+1];
return u1, u2, iters, change, gap;
class AP(Algorithm):
"""
Alternating Projections
See Also
--------
Algorithm : Generic interface for ProxToolbox algorithms
"""
def __init__(self, config):
"""
config : dict
Dictionary containing the problem configuration. It must contain the
following mappings:
proj1:ProxOperator
First ProxOperator (the class, no instance)
proj2:ProxOperator
Second ProxOperator (the class, no instance)
beta0:number
Starting relaxation parmater
beta_max:number
Maximum relaxation parameter
beta_switch:int
Iteration at which beta moves from beta0 -> beta_max
normM:number
?
Nx:int
Row-dim of the product space elements
Ny:int
Column-dim of the product space elements
Nz:int
Depth-dim of the product space elements
dim:int
Size of the product space
"""
self.proj1 = config['proj1'](config); self.proj2 = config['proj2'](config);
self.normM = config['normM'];
self.Nx = config['Nx']; self.Ny = config['Ny']; self.Nz = config['Nz'];
self.dim = config['dim'];
def run(self, u, tol, maxiter):
proj1 = self.proj1; proj2 = self.proj2;
normM = self.normM;
Nx = self.Nx; Ny = self.Ny; Nz = self.Nz; dim = self.dim;
iters = 0;
change = numpy.zeros(maxiter+1);
change[0] = 999;
gap = change.copy();
tmp1 = proj2.work(u);
while iters < maxiter and change[iters] >= tol:
iters += 1;
tmp_u = proj1.work(tmp1);
tmp1 = proj2.work(tmp_u);
tmp_change = 0; tmp_gap = 0;
if Ny == 1 or Nx == 1:
tmp_change = (scipy.linalg.norm(u-tmp_u,'fro')/normM)**2;
tmp_gap = (scipy.linalg.norm(tmp1-tmp_u,'fro')/normM)**2;
elif Nz == 1:
for j in xrange(dim):
tmp_change += (scipy.linalg.norm(u[:,:,j]-tmp_u[:,:,j],'fro')/normM)**2;
tmp_gap += (scipy.linalg.norm(tmp1[:,:,j]-tmp_u[:,:,j])/normM,'fro')**2;
else:
for j in xrange(dim):
for k in xrange(Nz):
tmp_change += (scipy.linalg.norm(u[:,:,k,j]-tmp_u[:,:,k,j],'fro')/normM)**2;
tmp_gap += (scipy.linalg.norm(tmp1[:,:,k,j]-tmp_u[:,:,k,j],'fro')/normM)**2;
change[iters] = math.sqrt(tmp_change);
gap[iters] = math.sqrt(tmp_gap);
u = tmp_u.copy();
tmp = proj1.work(u);
tmp2 = proj2.work(u);
if Ny == 1:
u1 = tmp[:,1];
u2 = tmp2[:,1];
elif Nx == 1:
u1 = tmp[1,:];
u2 = tmp2[1,:];
elif Nz == 1:
u1 = tmp[:,:,1];
u2 = tmp2[:,:,1];
else:
u1 = tmp;
u2 = tmp2;
change = change[1:iters+1];
gap = gap[1:iters+1];
return u1, u2, iters, change, gap;
# -*- coding: utf-8 -*-
"""
Created on Thu Dec 18 09:52:17 2014
@author: stefan
Just a small benchmarking script
"""
from matplotlib import pyplot
from ptychography import Ptychography_NTT_01_26210
import numpy, matplotlib
# Sort out some functions useless for our benchmark
class BenchmarkPtychography(Ptychography_NTT_01_26210):
def _postsolve(self):
return;
def go():
config = BenchmarkPtychography.default_config.copy();
config['blocking_switch'] = True;
config['maxiter'] = 250;
config['block_maxiter'] = 250;
# This might still need some thinking
config['tol'] = 1e-5;
blocking_sizes = (2,4,8,16,32,64);
blocking_schemes = [('sequential','none'),('parallel','none'),
('sequential','sequential'),('parallel','sequential')];
algorithms = ['PALM','PALMregPhiPtwise','Rodenburg','Thibault'];
for i in xrange(6):
config['block_rows'] = blocking_sizes[i];
config['block_cols'] = blocking_sizes[i];
for j in xrange(4):
config['between_blocks_scheme'] = blocking_schemes[j][0];
config['within_blocks_scheme'] = blocking_schemes[j][1];
for k in xrange(4):
config['ptychography_prox'] = algorithms[k];
problem = BenchmarkPtychography(config);
print 'Experiment', k + 4*j + 16*i + 1;
print 'block_rows/block_cols:', blocking_sizes[i];
print 'between_blocks_scheme: %s' % blocking_schemes[j][0];
print 'within_blocks_scheme: %s' % blocking_schemes[j][1];
print 'ptychography_prox : %s' % algorithms[k];
problem.solve();
print '';
pyplot.imsave('%d%d%d_amplitude.png' % (i,j,k), \
numpy.abs(problem.u_final.obj));
pyplot.imsave('%d%d%d_phase.png' % (i,j,k), \
numpy.angle(problem.u_final.obj));
pyplot.imsave('%d%d%d_probe.png' % (i,j,k), \
matplotlib.colors.hsv_to_rgb(problem._im2hsv(problem.u_final.probe,1.0)));
numpy.savez_compressed('%d%d%d_data' % (i,j,k), \
phi=problem.u_final.phi, \
obj=problem.u_final.obj, \
probe=problem.u_final.probe, \
change=problem.change, \
rms_err_obj=problem.custom_errors[0,:], \
rms_err_probe=problem.custom_errors[1,:], \
phi_change=problem.custom_errors[2,:], \
rfactor=problem.custom_errors[3,:], \
objective_value=problem.custom_errors[4,:], \
block_iters=problem.algorithm.block_iters);
\ No newline at end of file
import sys, os
from PyQt4 import QtGui, QtCore, uic
import json
class MainWindow(QtGui.QMainWindow):
def __init__(self, parent=None):
super(MainWindow, self).__init__(parent)
uic.loadUi("blocking_strategies_gui.ui", self)
# Connection to button
self.connect(self.okButton,QtCore.SIGNAL("clicked()"), self.ok)
def ok(self):
# get values from gui and store them into a json file
if str(self.blocking_switch.currentText()) == 'true':
config = {
'blocking_switch':str(self.blocking_switch.currentText()),
'blocking_scheme':str(self.blocking_scheme.currentText()),
'between_blocks_scheme':str(self.between_blocks_scheme.currentText()),
'within_blocks_scheme':str(self.within_blocks_scheme.currentText()),
'block_maxiter':int(self.block_maxiter.toPlainText()),
'block_rows':int(self.block_rows.toPlainText()),
'block_cols':int(self.block_cols.toPlainText()),};
else:
config = {'blocking_switch':str(self.blocking_switch.currentText()),}
with open('data_blocking.json','wb') as outfile:
json.dump(config, outfile)
self.close()
if __name__ == '__main__':
app = QtGui.QApplication(sys.argv)
window = MainWindow()
window.show()
app.exec_()
\ No newline at end of file
<
<?xml version="1.0" encoding="UTF-8"?>
<ui version="4.0">
<class>MainWindow</class>
<widget class="QMainWindow" name="MainWindow">
<property name="geometry">
<rect>
<x>0</x>
<y>0</y>
<width>280</width>
<height>300</height>
</rect>
</property>
<property name="windowTitle">
<string>MainWindow</string>
</property>
<widget class="QWidget" name="centralwidget">
<widget class="QGroupBox" name="groupBox">
<property name="geometry">
<rect>
<x>0</x>
<y>0</y>
<width>280</width>
<height>220</height>
</rect>
</property>
<property name="title">
<string>Parameters relating to blocking strategies</string>
</property>
<layout class="QFormLayout" name="formLayout">
<item row="0" column="0">
<widget class="QLabel" name="label">
<property name="text">
<string>use</string>
</property>
</widget>
</item>
<item row="0" column="1">
<widget class="QComboBox" name="blocking_switch">
<item>
<property name="text">
<string>true</string>
</property>
</item>
<item>
<property name="text">
<string>false</string>
</property>
</item>
</widget>
</item>
<item row="1" column="0">
<widget class="QLabel" name="label_2">
<property name="text">
<string>Blocking scheme</string>
</property>
</widget>
</item>
<item row="1" column="1">
<widget class="QComboBox" name="blocking_scheme">
<item>
<property name="text">
<string>one</string>
</property>
</item>
<item>
<property name="text">
<string>single_view</string>
</property>
</item>
<item>
<property name="text">
<string>greedy</string>
</property>
</item>
<item>
<property name="text">
<string>divide</string>
</property>
</item>
<item>
<property name="text">
<string>split</string>
</property>
</item>
<item>
<property name="text">
<string>orthogonal</string>
</property>
</item>
</widget>
</item>
<item row="2" column="0">
<widget class="QLabel" name="label_3">
<property name="text">
<string>Scheme between blocks</string>
</property>
</widget>
</item>
<item row="2" column="1">