Commit 12a22454 authored by alexander.dornheim's avatar alexander.dornheim
Browse files

Added Simpple RAAR Algortihm

Adjusted file path in JWST_data_processor
parent 597bbe05
# -*- coding: utf-8 -*-
"""
Created on Mon Dec 14 12:48:26 2015
# RAAR.m
# written on Aug.18 , 2017 by
# Russell Luke
# Inst. Fuer Numerische und Angewandte Mathematik
#
# DESCRIPTION: User-friendly version of the Relaxed Averaged Alternating Reflection algorithm.
# For background see:
# D.R.Luke, Inverse Problems 21:37-50(2005)
# D.R. Luke, SIAM J. Opt. 19(2): 714--739 (2008).
from .SimpleAlgortihm import SimpleAlgorithm
from numpy import exp
class RAAR(SimpleAlgorithm):
def evaluate(self, u):
iter = self.config['iter']
beta0 = self.config['beta_0']
beta_max = self.config['beta_max']
beta_switch = self.config['beta_switch']
tmp1 = 2*self.prox2.work(u)-u
tmp2=self.prox1.work(tmp1)
# update
beta = exp((-iter/beta_switch)**3)*beta0+(1-exp((-iter/beta_switch)**3))*beta_max # unrelaxes as the
unew = (beta*(2*tmp2-tmp1) + (1-beta)*tmp1 + u)/2
return unew
@author: rebecca
"""
from numpy import zeros, angle, trace, exp, sqrt
from numpy.linalg import norm
from .algorithms import Algorithm
class RAAR(Algorithm):
"""
Relaxed Averaged Alternating Reflection algorithm
"""
def __init__(self,config):
"""
Parameters
----------
config : dict
Dictionary containing the problem configuration.
It must contain the following mappings:
proxoperators: 2 ProxOperators
Tuple of ProxOperators (the class, no instance)
beta_0: number
Starting relaxation parmater
beta_max: number
Maximum relaxation parameter
beta_switch: int
Iteration at which beta moves from beta_0 -> 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.prox1 = config['proxoperators'][0](config)
self.prox2 = config['proxoperators'][1](config)
self.normM = config['normM']
self.beta_0 = config['beta_0']
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.product_space_dimension = config['product_space_dimension']
self.iter = 0
if 'truth' in config:
self.truth = config['truth']
self.truth_dim = config['truth_dim']
self.norm_truth = config['norm_truth']
def run(self, u, tol, maxiter):
"""
Runs the algorithm for the specified input data
"""
##### PREPROCESSING
normM = self.normM
beta = self.beta_0
iter = self.iter
if u.ndim < 3:
p = 1
q = 1
elif u.ndim == 3:
p = u.shape[2]
q = 1
else:
p = u.shape[2]
q = u.shape[3]
change = zeros(maxiter+1,dtype=u.dtype)
change[0] = 999
gap = change.copy()
shadow_change = change.copy()
if hasattr(self, 'truth'):
Relerrs = change.copy()
tmp1 = 2*self.prox2.work(u) - u
shadow = self.prox1.work(u)
##### LOOP
while iter < maxiter and change[iter] >= tol:
tmp = exp((-(iter+1)/self.beta_switch)**3);
beta = (tmp*self.beta_0) + ((1-tmp)*self.beta_max)
iter += 1;
tmp3 = self.prox1.work(tmp1)
tmp_u = ((beta*(2*tmp3-tmp1)) + ((1-beta)*tmp1) + u)/2
tmp2 = self.prox2.work(tmp_u)
tmp3 = self.prox1.work(tmp2)
tmp_change = 0; tmp_gap = 0; tmp_shadow=0;
#print(norm(tmp_u))
if p==1 and q==1:
tmp_change= (norm(u-tmp_u, 'fro')/normM)**2
tmp_shadow = (norm(tmp3-shadow,'fro')/normM)**2
tmp_gap = (norm(tmp3-tmp2,'fro')/normM)**2
if hasattr(self, 'truth'):
if self.truth_dim[0] == 1:
z=tmp3[0,:]
elif self.truth_dim[1] == 1:
z=tmp3[:,0]
else:
z=tmp3;
Relerrs[iter] = norm(self.truth - exp(-1j*angle(trace(self.truth.T*z))) * z, 'fro')/self.norm_truth
elif q==1:
for j in range(self.product_space_dimension):
tmp_change= tmp_change+ (norm(u[:,:,j]-tmp_u[:,:,j], 'fro')/normM)**2
# compute (||P_SP_Mx-P_Mx||/normM)^2:
tmp_gap = tmp_gap+(norm(tmp3[:,:,j]-tmp2[:,:,j],'fro')/normM)**2
tmp_shadow = tmp_shadow+(norm(tmp3[:,:,j]-shadow[:,:,j],'fro')/normM)**2
if hasattr(self, 'truth'):
z=tmp3[:,:,0]
Relerrs[iter] = norm(self.truth - exp(-1j*angle(trace(self.truth.T*z))) * z, 'fro')/self.norm_truth
else:
for j in range(self.product_space_dimension):
for k in range(Nz):
tmp_change= tmp_change+ (norm(u[:,:,k,j]-tmp_u[:,:,k,j], 'fro')/normM)**2
# compute (||P_Sx-P_Mx||/normM)^2:
tmp_gap = tmp_gap+(norm(tmp3[:,:,k,j]-tmp2[:,:,k,j],'fro')/(normM))**2
tmp_shadow = tmp_shadow+(norm(tmp3[:,:,k,j]-shadow[:,:,k,j],'fro')/(normM))**2
change[iter] = sqrt(tmp_change)
gap[iter] = sqrt(tmp_gap)
shadow_change[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
# does not matter.
# graphics
u = tmp_u
tmp1 = (2*tmp2) - tmp_u
##### POSTPROCESSING
u = tmp2;
tmp = self.prox1.work(u);
tmp2 = self.prox2.work(u);
if self.Nx == 1:
u1 = tmp[:,0];
u2 = tmp2[:,0];
elif self.Ny == 1:
u1 = tmp[0,:];
u2 = tmp2[0,:];
elif self.Nz == 1 and tmp.ndim > 2:
u1 = tmp[:,:,0]
u2 = tmp2[:,:,0]
else:
u1 = tmp;
u2 = tmp2;
change = change[1:iter+1];
gap = gap[1:iter+1];
shadow_change = shadow_change[1:iter+1]
output = {'u1': u1, 'u2': u2, 'iter': iter, 'change': change, 'gap': gap, 'shadow_change' : shadow_change}
if hasattr(self, 'truth'):
output['Relerrs']=Relerrs
return output
# -*- coding: utf-8 -*-
"""
Created on Mon Dec 14 12:48:26 2015
@author: rebecca
"""
from numpy import zeros, angle, trace, exp, sqrt
from numpy.linalg import norm
from .algorithms import Algorithm
class RAAR_expert(Algorithm):
"""
Relaxed Averaged Alternating Reflection algorithm
"""
def __init__(self,config):
"""
Parameters
----------
config : dict
Dictionary containing the problem configuration.
It must contain the following mappings:
proxoperators: 2 ProxOperators
Tuple of ProxOperators (the class, no instance)
beta_0: number
Starting relaxation parmater
beta_max: number
Maximum relaxation parameter
beta_switch: int
Iteration at which beta moves from beta_0 -> 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.prox1 = config['proxoperators'][0](config)
self.prox2 = config['proxoperators'][1](config)
self.normM = config['normM']
self.beta_0 = config['beta_0']
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.product_space_dimension = config['product_space_dimension']
self.iter = 0
if 'truth' in config:
self.truth = config['truth']
self.truth_dim = config['truth_dim']
self.norm_truth = config['norm_truth']
def run(self, u, tol, maxiter):
"""
Runs the algorithm for the specified input data
"""
##### PREPROCESSING
normM = self.normM
beta = self.beta_0
iter = self.iter
if u.ndim < 3:
p = 1
q = 1
elif u.ndim == 3:
p = u.shape[2]
q = 1
else:
p = u.shape[2]
q = u.shape[3]
change = zeros(maxiter+1,dtype=u.dtype)
change[0] = 999
gap = change.copy()
shadow_change = change.copy()
if hasattr(self, 'truth'):
Relerrs = change.copy()
tmp1 = 2*self.prox2.work(u) - u
shadow = self.prox1.work(u)
##### LOOP
while iter < maxiter and change[iter] >= tol:
tmp = exp((-(iter+1)/self.beta_switch)**3);
beta = (tmp*self.beta_0) + ((1-tmp)*self.beta_max)
iter += 1;
tmp3 = self.prox1.work(tmp1)
tmp_u = ((beta*(2*tmp3-tmp1)) + ((1-beta)*tmp1) + u)/2
tmp2 = self.prox2.work(tmp_u)
tmp3 = self.prox1.work(tmp2)
tmp_change = 0; tmp_gap = 0; tmp_shadow=0;
#print(norm(tmp_u))
if p==1 and q==1:
tmp_change= (norm(u-tmp_u, 'fro')/normM)**2
tmp_shadow = (norm(tmp3-shadow,'fro')/normM)**2
tmp_gap = (norm(tmp3-tmp2,'fro')/normM)**2
if hasattr(self, 'truth'):
if self.truth_dim[0] == 1:
z=tmp3[0,:]
elif self.truth_dim[1] == 1:
z=tmp3[:,0]
else:
z=tmp3;
Relerrs[iter] = norm(self.truth - exp(-1j*angle(trace(self.truth.T*z))) * z, 'fro')/self.norm_truth
elif q==1:
for j in range(self.product_space_dimension):
tmp_change= tmp_change+ (norm(u[:,:,j]-tmp_u[:,:,j], 'fro')/normM)**2
# compute (||P_SP_Mx-P_Mx||/normM)^2:
tmp_gap = tmp_gap+(norm(tmp3[:,:,j]-tmp2[:,:,j],'fro')/normM)**2
tmp_shadow = tmp_shadow+(norm(tmp3[:,:,j]-shadow[:,:,j],'fro')/normM)**2
if hasattr(self, 'truth'):
z=tmp3[:,:,0]
Relerrs[iter] = norm(self.truth - exp(-1j*angle(trace(self.truth.T*z))) * z, 'fro')/self.norm_truth
else:
for j in range(self.product_space_dimension):
for k in range(Nz):
tmp_change= tmp_change+ (norm(u[:,:,k,j]-tmp_u[:,:,k,j], 'fro')/normM)**2
# compute (||P_Sx-P_Mx||/normM)^2:
tmp_gap = tmp_gap+(norm(tmp3[:,:,k,j]-tmp2[:,:,k,j],'fro')/(normM))**2
tmp_shadow = tmp_shadow+(norm(tmp3[:,:,k,j]-shadow[:,:,k,j],'fro')/(normM))**2
change[iter] = sqrt(tmp_change)
gap[iter] = sqrt(tmp_gap)
shadow_change[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
# does not matter.
# graphics
u = tmp_u
tmp1 = (2*tmp2) - tmp_u
##### POSTPROCESSING
u = tmp2;
tmp = self.prox1.work(u);
tmp2 = self.prox2.work(u);
if self.Nx == 1:
u1 = tmp[:,0];
u2 = tmp2[:,0];
elif self.Ny == 1:
u1 = tmp[0,:];
u2 = tmp2[0,:];
elif self.Nz == 1 and tmp.ndim > 2:
u1 = tmp[:,:,0]
u2 = tmp2[:,:,0]
else:
u1 = tmp;
u2 = tmp2;
change = change[1:iter+1];
gap = gap[1:iter+1];
shadow_change = shadow_change[1:iter+1]
output = {'u1': u1, 'u2': u2, 'iter': iter, 'change': change, 'gap': gap, 'shadow_change' : shadow_change}
if hasattr(self, 'truth'):
output['Relerrs']=Relerrs
return output
......@@ -40,9 +40,6 @@ class SimpleAlgorithm:
self.prox1 = config['proxoperators'][0](config)
self.prox2 = config['proxoperators'][1](config)
self.normM = config['normM']
self.beta_0 = config['beta_0']
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.product_space_dimension = config['product_space_dimension']
self.iter = 0
......
......@@ -44,13 +44,13 @@ def JWST_data_processor(config):
diversity=3;
with open('../InputData/phase_retrieval/phase_p37.pmod','r') as fid:
with open('../InputData/Phase/phase_p37.pmod','r') as fid:
# read lower endian float <f
temp1 = np.fromfile(fid, dtype='<f');
temp1 = temp1.astype(np.float64)
temp1 = temp1.reshape(512,512).T;
with open('../InputData/phase_retrieval/phase_m37.pmod','r') as fid:
with open('../InputData/Phase/phase_m37.pmod','r') as fid:
# read lower endian float <f
temp2 = np.fromfile(fid, dtype='<f');
temp2 = temp2.astype(np.float64)
......
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