Commit aa425018 authored by alexander.dornheim's avatar alexander.dornheim
Browse files

Updated AP algortihm (partially)

Added test for tasse with AP, RAAR with beta_0=beta_max=0.5
parent 043d6ca5
......@@ -23,9 +23,9 @@ class AP(Algorithm):
Dictionary containing the problem configuration.
It must contain the following mappings:
proj1: ProxOperator
prox1: ProxOperator
First ProxOperator (the class, no instance)
proj2: ProxOperator
prox2: ProxOperator
Second ProxOperator (the class, no instance)
beta0: number
Starting relaxation parmater
......@@ -44,17 +44,34 @@ class AP(Algorithm):
dim: int
Size of the product space
"""
self.proj1 = config['proj1'](config); self.proj2 = config['proj2'](config);
self.prox1 = config['proxoperators'][0](config); self.prox2 = config['proxoperators'][1](config);
self.normM = config['normM'];
self.Nx = config['Nx']; self.Ny = config['Ny']; self.Nz = config['Nz'];
self.dim = config['dim'];
self.product_space_dimension = config['product_space_dimension'];
self.iters = 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
"""
proj1 = self.proj1; proj2 = self.proj2;
prox1 = self.prox1; prox2 = self.prox2;
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]
normM = self.normM;
iters = self.iters
......@@ -62,47 +79,60 @@ class AP(Algorithm):
change[0] = 999;
gap = change.copy();
tmp1 = proj2.work(u);
tmp1 = prox2.work(u);
while iters < maxiter and change[iters] >= tol:
iters += 1;
tmp_u = proj1.work(tmp1);
tmp1 = proj2.work(tmp_u);
tmp_u = prox1.work(tmp1);
tmp1 = prox2.work(tmp_u);
tmp_change = 0; tmp_gap = 0;
if self.Ny == 1 or self.Nx == 1:
tmp_change = (norm(u-tmp_u,'fro')/normM)**2;
tmp_gap = (norm(tmp1-tmp_u,'fro')/normM)**2;
elif self.Nz == 1:
for j in range(self.dim):
tmp_change += (norm(u[:,:,j]-tmp_u[:,:,j],'fro')/normM)**2;
tmp_gap += (norm(tmp1[:,:,j]-tmp_u[:,:,j])/normM,'fro')**2;
if p==1 and q==1:
tmp_change= (norm(u-tmp_u, 'fro')/normM)**2
tmp_gap = (norm(tmp1-tmp_u,'fro')/normM)**2
if hasattr(self, 'truth'):
if self.truth_dim[0] == 1:
z=tmp_u[0,:]
elif self.truth_dim[1]==1:
z=tmp_u[:,0]
else:
z=tmp_u
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(tmp1[:,:,j]-tmp_u[:,:,j],'fro')/normM)**2
else:
for j in range(self.dim):
for k in range(self.Nz):
tmp_change += (norm(u[:,:,k,j]-tmp_u[:,:,k,j],'fro')/normM)**2;
tmp_gap += (norm(tmp1[:,:,k,j]-tmp_u[:,:,k,j],'fro')/normM)**2;
for j in range(self.product_space_dimension):
for k in range(self.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(tmp1[:,:,k,j]-tmp_u[:,:,k,j],'fro')/(normM))**2
change[iters] = sqrt(tmp_change);
gap[iters] = sqrt(tmp_gap);
u = tmp_u.copy();
tmp = proj1.work(u);
tmp2 = proj2.work(u);
if self.Ny == 1:
u1 = tmp[:,1];
u2 = tmp2[:,1];
elif self.Nx == 1:
u1 = tmp[1,:];
u2 = tmp2[1,:];
elif self.Nz == 1:
u1 = tmp[:,:,1];
u2 = tmp2[:,:,1];
tmp = prox1.work(u);
tmp2 = 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;
u1 = tmp
u2 = tmp2
change = change[1:iters+1];
gap = gap[1:iters+1];
......
......@@ -166,8 +166,8 @@ class RAAR(Algorithm):
u1 = tmp[0,:];
u2 = tmp2[0,:];
elif self.Nz == 1 and tmp.ndim > 2:
u1 = tmp[:,:,1]
u2 = tmp2[:,:,1]
u1 = tmp[:,:,0]
u2 = tmp2[:,:,0]
else:
u1 = tmp;
u2 = tmp2;
......
......@@ -173,17 +173,6 @@ class Phase(Problem):
proj2 = self.config['proxoperators'][1](self.config)
u_2 = proj2.work(u_1);
print(norm(self.config['u_0']))
print(norm(u_1))
print(norm(u_2))
print(norm(u_1-u_2))
print(self.config['norm_rt_data'])
#print(nonzero(u_1))
print(u_1[51,55]);
#print(size(self.config['support_idx'][1]))
print(u_2[51,55])
print(proj2)
# estimate the gap in the relevant metric
if self.config['Nx'] ==1 or self.config['Ny']==1 :
tmp_gap = square(norm(u_1-u_2)/self.config['norm_rt_data']);
......@@ -196,7 +185,6 @@ class Phase(Problem):
tmp_gap = tmp_gap+(norm(u_1[:,:,j]-u_2[:,:,j])/self.config['norm_rt_data'])**2
gap_0=sqrt(tmp_gap);
print(gap_0)
# sets the set fattening to be a percentage of the
# initial gap to the unfattened set with
......@@ -207,7 +195,7 @@ class Phase(Problem):
# the second tolerance relative to the oder of
# magnitude of the metric
self.config['TOL2'] = self.config['data_ball']*1e-15;
self.config['proxoperators']
self.algorithm = getattr(Algorithms, self.config['algorithm'])(self.config);
......@@ -263,15 +251,30 @@ class Phase(Problem):
return
u1 = f['u1'].value.view(np.complex)
else:
if self.config['MAXIT'] == 1000 :
f = h5py.File('Phase_test_data/tasse_u1_1000.mat')
elif self.config['MAXIT'] == 20:
f = h5py.File('Phase_test_data/tasse_u1_20.mat')
elif self.config['MAXIT'] == 1:
f = h5py.File('Phase_test_data/tasse_u1_1.mat')
else:
print("No file available to compare to.")
return
if self.config['algorithm'] == 'RAAR':
if self.config['beta_0'] == 0.95:
if self.config['MAXIT'] == 1000 :
f = h5py.File('Phase_test_data/tasse_u1_1000.mat')
elif self.config['MAXIT'] == 20:
f = h5py.File('Phase_test_data/tasse_u1_20.mat')
elif self.config['MAXIT'] == 1:
f = h5py.File('Phase_test_data/tasse_u1_1.mat')
else:
print("No file available to compare to.")
return
elif self.config['beta_0'] == 0.50:
if self.config['MAXIT'] == 20:
f = h5py.File('Phase_test_data/tasse_u1_20_beta_0_5.mat')
else:
print("No file available to compare to.")
return
else:
print("No file available to compare to.")
return
elif self.config['algorithm'] == 'AP':
if self.config['MAXIT'] == 20:
f = h5py.File('Phase_test_data/tasse_u1_ap_20.mat')
u1 = f['u1'].value.view(np.float)
u1 =np.array(u1)
......
......@@ -85,11 +85,11 @@ new_config = {
#if(strcmp('problem_family,'Phase'))
## maximum number of iterations and tolerances
'MAXIT' : 1,
'MAXIT' : 20,
'TOL' : 1e-12,
## relaxaton parameters in RAAR, HPR and HAAR
'beta_0' : 0.95, # starting relaxation prameter (only used with
'beta_0' : 0.95, #0.95 # starting relaxation prameter (only used with
# HAAR, HPR and RAAR)
'beta_max' :0.50, # maximum relaxation prameter (only used with
# HAAR, RAAR, and HPR)
......
......@@ -262,7 +262,6 @@ class Approx_P_JWST_Poisson(ProxOperator):
self.exp_illu_minus = exp(-1j*self.illumination_phase)*tile(self.indicator_ampl[...,None],(1,1,self.product_space_dimension-1))
#@profile
def work(self,u):
TOL2 = self.TOL2;
epsilon = self.epsilon;
......
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