Commit 0b9d1e00 authored by alexander.dornheim's avatar alexander.dornheim
Browse files

Added first version of Q_Heau

parent 3cb7f67f
......@@ -145,19 +145,19 @@ class HAAR(Algorithm):
T_output = self.T.run(self.T_config['u_0'],self.T_config['TOL'],self.T_config['MAXIT'])
Ty = T_output['u']
if(p==1) and (q==1) and (n!=1) and (m!=1):
Ty = Ty.reshape((n*m,1))
tmp_y= y.reshape((n*m,1))
elif (p>1) and (q==1):
tmp_y = zeros((n*m*p,1,1))
if dim_number == 2:
Ty = Ty.reshape((s[0]*s[1],1))
tmp_y= y.reshape((s[0]*s[1],1))
elif dim_number == 3:
tmp_y = zeros((s[0]*s[1]*s[2],1,1))
tmp_Ty=tmp_y
for j in range(p):
ytmp = y[:,:,j].reshape((n*m,1))
tmp_y[n*m*j:m*n*(j+1)] = ytmp
ytmp = Ty[:,:,j].reshape((n*m,1))
tmp_Ty[n*m*j:m*n*(j+1)] = ytmp
ytmp = y[:,:,j].reshape((s[0]*s[1],1))
tmp_y[s[0]*s[1]*j:s[0]*s[1]*(j+1)] = ytmp
ytmp = Ty[:,:,j].reshape((s[0]*s[1],1))
tmp_Ty[s[0]*s[1]*j:s[0]*s[1]*(j+1)] = ytmp
Ty=tmp_Ty
elif(q>1):
elif dim_number == 4:
print('not ready for 4D arrays')
y_new=feval('Q_Heau',y0,tmp_y,(1-mu)*tmp_y+mu*Ty)
......
......@@ -654,3 +654,36 @@ class P_block_sequential_hyperplane(ProxOperator):
# but one step is close enough if the blocks consist of orthogonal
# elemenets...
return v0
# Q_Heau.m
# written on MAY 31, 2006 by
# Russell Luke & Daniel Gempesaw
# University of Delaware
#
# DESCRIPTION: Q operator for the Heaugaseau-like algorithm as
# described in Bauschke Combettes&Luke Jour. Approx. Thry., 2006
#
# INPUT: x,y,z - all column vectors, not nec. real
#
# OUTPUT: y_new = Q(x,y,z)
#
# USAGE: y_new = Q_Heau(x,y,z)
#
# Nonstandard Matlab function calls: none
class Q_Heau(ProxOperator):
def run(x,y,z):
xsi = real((y-z).T @ (x-y))+imag((y-z).T @ (x-y))
eta = (x-y).T @ (x-y)
nu = (y-z).T @ (y-z)
rho = eta*nu-xsi**2
if((rho<=0) and (xsi>=0)):
return z
elif((rho>0) and (xsi*nu>=rho)):
return x+(1+xsi/nu)*(z-y)
elif((rho>0) and (xsi*nu<rho)):
return y+(nu/rho)*(xsi*(x-y)+nu*(z-y))
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