Skip to content

GitLab

Commit f719d5c8 authored by markus.meier01's avatar markus.meier01
Browse files

Added graphical output for 1D, fixed some things in QNAP

parent d7885374
Hide whitespace changes
Inline Side-by-side
proxtoolbox/Algorithms/QNAP.py
View file @ f719d5c8
......@@ -72,14 +72,14 @@ class QNAP(Algorithm):
else:
self.diagnostic = False
self.samsara = Samsara()
self.Samsara = Samsara()
def run(self, u, tol, maxiter):
"""
Runs the algorithm for the specified input data
"""
Samsara = self.Samsara
##### PREPROCESSING
iter = self.iter
......@@ -116,7 +116,6 @@ class QNAP(Algorithm):
while iter < maxiter and change[iter] >= tol:
iter += 1
self.iter =iter
tmp_u = prox1.work(tmp1)
# make call to SAMSARA for acceleration step
# since for the product space Prox1 is the
......@@ -125,33 +124,33 @@ class QNAP(Algorithm):
# array to samsara. Have to reshape the complex
# array as a real-valued vector. Use the gap as the
# objective function, and the negative gradient tmp_u- u
if (p==1) and (q==1):
if (np.all(isreal(u))):
if iter > 4:
uold_vec=u[:,0]
unew_vec=tmp_u[:,0]
gradfold_vec = gradfnew_vec
gradfnew_vec = (u[:,0]- tmp_u[:,0])
unew_vec, uold_vec, gap[iter-1], gradfold_vec, change[iter] = \
self.samsara.run(uold_vec, unew_vec,
gap[iter-2]*self.Nx*self.Ny,
gap[iter-1]*self.Nx*self.Ny,
gradfold_vec, gradfnew_vec)
for j in range(self.product_space_dimension):
tmp_u[:,j]=unew_vec
gap[iter-1]=gap[iter-1]/(self.Nx*self.Ny)
if (p==1) and (q==1):
if (np.all(isreal(u))):
if iter > 4:
uold_vec=u[:,0]
unew_vec=tmp_u[:,0]
gradfold_vec = gradfnew_vec
gradfnew_vec = (u[:,0]- tmp_u[:,0])
unew_vec, uold_vec, gap[iter-1], gradfold_vec, change[iter] = \
Samsara.run(uold_vec, unew_vec,
gap[iter-2]*self.Nx*self.Ny,
gap[iter-1]*self.Nx*self.Ny,
gradfold_vec, gradfnew_vec)
for j in range(self.product_space_dimension):
tmp_u[:,j]=unew_vec
gap[iter-1]=gap[iter-1]/(self.Nx*self.Ny)
else:
unew_vec=tmp_u[:,0]
uold_vec=u[:,0]
gradfnew_vec = (u[:,0]- tmp_u[:,0])
else:
unew_vec=tmp_u[:,0]
uold_vec=u[:,0]
gradfnew_vec = (u[:,0]- tmp_u[:,0])
else:
if iter>3:
uold_vec= np.concatenate(np.real(u[:,0]), np.imag(u[:,0])).reshape(self.Ny*self.Nx*2, 0)
if iter>3:
uold_vec= np.concatenate((np.real(u[:,0]), np.imag(u[:,0]))).reshape(self.Ny*self.Nx*2, 1)
unew_vec=np.concatenate(np.real(tmp_u[:,0]), np.imag(tmp_u[:,0])).reshape(self.Ny*self.Nx*2, 0)
unew_vec=np.concatenate((np.real(tmp_u[:,0]), np.imag(tmp_u[:,0]))).reshape(self.Ny*self.Nx*2, 1)
gradfold_vec = gradfnew_vec
gradfnew_vec = uold_vec-unew_vec
gradfold_vec = gradfnew_vec
gradfnew_vec = uold_vec-unew_vec
# unew_vec,uold_vec,gap[iter-0],gradfold_vec, change[iter]= \
# feval('samsara', uold_vec, unew_vec,
......@@ -159,115 +158,115 @@ class QNAP(Algorithm):
# gap(iter-1)*self.Nx*self.Ny,
# gradfold_vec, gradfnew_vec)
unew_vec, uold_vec, gap[iter], gradfold_vec, change[iter] = Samsara.run(uold_vec, unew_vec, gap(iter-2)*self.Nx*self.Ny, gap(iter-1)*self.Nx*self.Ny, gradfold_vec, gradfnew_vec)
unew_vec, uold_vec, gap[iter], gradfold_vec, change[iter] = Samsara.run(uold_vec, unew_vec, gap[iter-2]*self.Nx*self.Ny, gap[iter-1]*self.Nx*self.Ny, gradfold_vec, gradfnew_vec)
# now reshape unew_vec
tmp_u_vec = unew_vec[0:self.Ny*self.Nx-1]+1j*unew_vec[self.Ny*self.Nx:self.Ny*self.Nx*2-1]
for j in range(self.product_space_dimension):
tmp_u[:,j]=tmp_u_vec
gap[iter-1]=gap[iter-1]/(self.Nx*self.Ny)
else:
uold_vec= np.array(np.real(u[:,0]), np.imag(u[:,0])).reshape(self.Ny*self.Nx*2, 1)
unew_vec= np.array( np.real(tmp_u[:,0]), np.imag(tmp_u[:,0])).reshape(self.Ny*self.Nx*2, 1)
gradfnew_vec = uold_vec-unew_vec
elif (p!=1) and (q==1):
if (np.all(isreal(u))):
tmp2 = u[:,:,0]
uold_vec= tmp2.reshape(self.Nx*self.Ny,1)
tmp2 = tmp_u[:,:,0]
unew_vec = tmp2.reshape(self.Nx*self.Ny,1)
if iter<=3:
gradfnew_vec = uold_vec-unew_vec
else:
gradfold_vec = gradfnew_vec
gradfnew_vec = uold_vec-unew_vec
unew_vec,uold_vec,gap[iter-1],gradfold_vec, change[iter]= \
self.samsara.run(uold_vec, unew_vec,
tmp_u_vec = unew_vec[0:self.Ny*self.Nx-1]+1j*unew_vec[self.Ny*self.Nx:self.Ny*self.Nx*2-1]
for j in range(self.product_space_dimension):
tmp_u[:,j]=tmp_u_vec
gap[iter-1]=gap[iter-1]/(self.Nx*self.Ny)
else:
uold_vec= np.array((np.real(u[:,0]), np.imag(u[:,0]))).reshape(self.Ny*self.Nx*2, 1)
unew_vec= np.array((np.real(tmp_u[:,0]), np.imag(tmp_u[:,0]))).reshape(self.Ny*self.Nx*2, 1)
gradfnew_vec = uold_vec-unew_vec
elif (p!=1) and (q==1):
if (np.all(isreal(u))):
tmp2 = u[:,:,0]
uold_vec= tmp2.reshape(self.Nx*self.Ny,1)
tmp2 = tmp_u[:,:,0]
unew_vec = tmp2.reshape(self.Nx*self.Ny,1)
if iter<=3:
gradfnew_vec = uold_vec-unew_vec
else:
gradfold_vec = gradfnew_vec
gradfnew_vec = uold_vec-unew_vec
unew_vec,uold_vec,gap[iter-1],gradfold_vec, change[iter]= \
Samsara.run(uold_vec, unew_vec,
gap(iter-2)*self.Nx*self.Ny,
gap(iter-1)*self.Nx*self.Ny,
gradfold_vec, gradfnew_vec)
# now reshape and replace u
gap[iter-1]=gap[iter-1]/(self.Nx*self.Ny);
gap[iter-1]=gap[iter-1]/(self.Nx*self.Ny);
tmp2=unew_vec.reshape(self.Ny,self.Nx);
for j in range(self.product_space_dimension):
tmp_u[:,:,j]=tmp2
tmp2=unew_vec.reshape(self.Ny,self.Nx);
for j in range(self.product_space_dimension):
tmp_u[:,:,j]=tmp2
else:
tmp2=np.concatenate((np.real(u[:,:,0]),np.imag(u[:,:,0])), axis=1)
uold_vec=tmp2.reshape(self.Nx*self.Ny*2,1)
tmp2= np.concatenate((np.real(tmp_u[:,:,0]), np.imag(tmp_u[:,:,0])), axis=1)
unew_vec = tmp2.reshape(self.Nx*self.Ny*2,1)
if iter<=3:
gradfnew_vec = uold_vec-unew_vec
else:
gradfold_vec = gradfnew_vec
gradfnew_vec = uold_vec-unew_vec
unew_vec,uold_vec,gap[iter-1],gradfold_vec, change[iter]= \
self.samsara.run( uold_vec, unew_vec,
gap(iter-2)*self.Nx*self.Ny,
gap(iter-1)*self.Nx*self.Ny,
tmp2=np.concatenate((np.real(u[:,:,0]),np.imag(u[:,:,0])), axis=1)
uold_vec=tmp2.reshape(self.Nx*self.Ny*2,1)
tmp2= np.concatenate((np.real(tmp_u[:,:,0]), np.imag(tmp_u[:,:,0])), axis=1)
unew_vec = tmp2.reshape(self.Nx*self.Ny*2,1)
if iter<=3:
gradfnew_vec = uold_vec-unew_vec
else:
gradfold_vec = gradfnew_vec
gradfnew_vec = uold_vec-unew_vec
unew_vec,uold_vec,gap[iter-1],gradfold_vec, change[iter]= \
Samsara.run( uold_vec, unew_vec,
gap[iter-2]*self.Nx*self.Ny,
gap[iter-1]*self.Nx*self.Ny,
gradfold_vec, gradfnew_vec)
# now reshape and replace u
gap[iter-1]=gap[iter-1]/(self.Nx*self.Ny)
tmp2=unew_vec.reshape(self.Ny,self.Nx*2)
unew=tmp2[:,1:self.Nx] +1j*tmp2[:,self.Nx+1:2*self.Nx]
for j in range(self.product_space_dimension):
tmp_u[:,:,j]=unew
else: # product space of 3D arrays
print('Cannot handle 3D arrays on the product space yet')
tmp1 = self.prox2.work(tmp_u)
tmp_change=0; tmp_gap=0
if(p==1)and(q==1):
tmp_change = np.linalg.norm(u-tmp_u, 'fro')/norm_data**2
tmp_gap = np.linalg.norm(tmp1-tmp_u,'fro')/norm_data**2
if hasattr(self, 'diagnostic'): #('diagnostic' in config)
if hasattr(self, "truth"): #('truth' in config)
if self.truth_dim[0]==1:
z=tmp_u[0,:]
elif self.truth_dim[0]==1:
z=tmp_u[:,0]
else:
z=tmp_u;
gap[iter-1]=gap[iter-1]/(self.Nx*self.Ny)
tmp2=unew_vec.reshape(self.Ny,self.Nx*2)
unew=tmp2[:,1:self.Nx] +1j*tmp2[:,self.Nx+1:2*self.Nx]
for j in range(self.product_space_dimension):
tmp_u[:,:,j]=unew
else: # product space of 3D arrays
print('Cannot handle 3D arrays on the product space yet')
tmp1 = self.prox2.work(tmp_u)
tmp_change=0; tmp_gap=0
if(p==1)and(q==1):
tmp_change = np.linalg.norm(u-tmp_u, 'fro')/norm_data**2
tmp_gap = np.linalg.norm(tmp1-tmp_u,'fro')/norm_data**2
if hasattr(self, 'diagnostic'): #('diagnostic' in config)
if hasattr(self, "truth"): #('truth' in config)
if self.truth_dim[0]==1:
z=tmp_u[0,:]
elif self.truth_dim[0]==1:
z=tmp_u[:,0]
else:
z=tmp_u;
Relerrs[iter] = np.linalg.norm(self.truth - np.exp(-1j*np.angle(np.trace(self.truth*z))) * z, 'fro')/self.norm_truth;
Relerrs[iter] = np.linalg.norm(self.truth - np.exp(-1j*np.angle(np.trace(self.truth*z))) * z, 'fro')/self.norm_truth;
elif(q==1):
for j in range(0,self.product_space_dimension):
tmp_change= tmp_change+ np.linalg.norm(u[:,:,j]-tmp_u[:,:,j], 'fro')/norm_data**2;
tmp_gap = tmp_gap+ np.linalg.norm(tmp1[:,:,j]-tmp_u[:,:,j],'fro')/norm_data**2
elif(q==1):
for j in range(0,self.product_space_dimension):
tmp_change= tmp_change+ np.linalg.norm(u[:,:,j]-tmp_u[:,:,j], 'fro')/norm_data**2;
tmp_gap = tmp_gap+ np.linalg.norm(tmp1[:,:,j]-tmp_u[:,:,j],'fro')/norm_data**2
if hasattr(self, "truth"): #'truth' in config
Relerrs[iter] = np.linalg.norm(self.truth - np.exp(-1j*np.angle(np.trace(self.truth*tmp_u[:,:,1]))) * tmp_u[:,:,1], 'fro')/self.norm_truth;
if hasattr(self, "truth"): #'truth' in config
Relerrs[iter] = np.linalg.norm(self.truth - np.exp(-1j*np.angle(np.trace(self.truth*tmp_u[:,:,1]))) * tmp_u[:,:,1], 'fro')/self.norm_truth;
else:
Relerrs[iter]=0;
for j in range(0,self.product_space_dimension):
for k in range(r1,self.Nz):
tmp_change= tmp_change+(np.linalg.norm(u[:,:,k,j]-tmp_u[:,:,k,j], 'fro')/norm_data)**2;
# compute (||P_Sx-P_Mx||/norm_data)^2:
tmp_gap = tmp_gap+(np.linalg.norm(tmp1[:,:,k,j]-tmp_u[:,:,k,j],'fro')/(norm_data))**2;
else:
Relerrs[iter]=0;
for j in range(0,self.product_space_dimension):
for k in range(r1,self.Nz):
tmp_change= tmp_change+(np.linalg.norm(u[:,:,k,j]-tmp_u[:,:,k,j], 'fro')/norm_data)**2;
# compute (||P_Sx-P_Mx||/norm_data)^2:
tmp_gap = tmp_gap+(np.linalg.norm(tmp1[:,:,k,j]-tmp_u[:,:,k,j],'fro')/(norm_data))**2;
if hasattr(self, "truth") and j == 1: # (any(strcmp('truth',fieldnames(method_input))))&&(j==1)
Relerrs[iter] = Relerrs[iter]+np.linalg.norm(self.truth - np.exp(-1j*np.angle(np.trace(self.truth*tmp_u[:,:,:,1]))) * tmp_u[:,:,k,1], 'fro')/self.norm_truth;
if hasattr(self, "truth") and j == 1: # (any(strcmp('truth',fieldnames(method_input))))&&(j==1)
Relerrs[iter] = Relerrs[iter]+np.linalg.norm(self.truth - np.exp(-1j*np.angle(np.trace(self.truth*tmp_u[:,:,:,1]))) * tmp_u[:,:,k,1], 'fro')/self.norm_truth;
change[iter]=np.sqrt(tmp_change);
gap[iter] = np.sqrt(tmp_gap);
u=tmp_u;
change[iter]=np.sqrt(tmp_change);
gap[iter] = np.sqrt(tmp_gap);
u=tmp_u;
if hasattr(self, "diagnostic"): #(any(strcmp('diagnostic', fieldnames(method_input))))
if hasattr(self, "diagnostic"): #(any(strcmp('diagnostic', fieldnames(method_input))))
# graphics
#print(config["anim"])
if (self.anim>=1)and (iter%2==0):
self.u=tmp1;
if (self.anim>=1)and (iter%2==0):
self.u=tmp1;
#self=self.animation(method_output);
#self.animation()
......
proxtoolbox/Algorithms/Wirtinger.py
View file @ f719d5c8
......@@ -143,7 +143,6 @@ class Wirtinger(Algorithm):
gap[iter] = sqrt(tmp_gap)
change[iter] = sqrt(tmp_change)
print(14)
tmp = Prox1.work(u)
tmp2 = Prox2.work(u)
......
proxtoolbox/Problems/Phase/Graphics/Phase_graphics.py
View file @ f719d5c8
......@@ -39,12 +39,13 @@ def Phase_graphics(config, output):
#beta_max = config['beta_max']
u_0 = config['u_0']
if output['u1'].ndim == 2:
if output['u1'].ndim == 1 or output['u1'].ndim == 2:
u = output['u1']
u2 = output['u2']
else:
u = output['u1'][:,:,0]
u2 = output['u2'][:,:,0]
iter = output['iter']
change = output['change']
......@@ -55,38 +56,50 @@ def Phase_graphics(config, output):
f, ((ax1, ax2), (ax3, ax4)) = subplots(2, 2)
if output['u1'].ndim >= 2:
f, ((ax1, ax2), (ax3, ax4)) = subplots(2, 2)
im=ax1.imshow(np.abs(u),cmap='gray')
f.colorbar(im, ax=ax1)
ax1.set_title('best approximation amplitude - physical constraint satisfied')
im=ax1.imshow(np.abs(u),cmap='gray')
f.colorbar(im, ax=ax1)
ax1.set_title('best approximation amplitude - physical constraint satisfied')
im=ax2.imshow(np.real(u),cmap='gray')
f.colorbar(im, ax=ax2)
ax2.set_title('best approximation phase - physical constraint satisfied')
im=ax2.imshow(np.real(u),cmap='gray')
f.colorbar(im, ax=ax2)
ax2.set_title('best approximation phase - physical constraint satisfied')
im=ax3.imshow(np.abs(u2),cmap='gray')
f.colorbar(im, ax=ax3)
ax3.set_title('best approximation amplitude - Fourier constraint satisfied')
im=ax3.imshow(np.abs(u2),cmap='gray')
f.colorbar(im, ax=ax3)
ax3.set_title('best approximation amplitude - Fourier constraint satisfied')
im=ax4.imshow(np.real(u2),cmap='gray')
f.colorbar(im, ax=ax4)
ax4.set_title('best approximation amplitude - Fourier constraint satisfied')
im=ax4.imshow(np.real(u2),cmap='gray')
f.colorbar(im, ax=ax4)
ax4.set_title('best approximation amplitude - Fourier constraint satisfied')
g, ((bx1, bx2), (bx3, bx4)) = subplots(2, 2)
im=bx1.imshow(np.abs(u),cmap='gray')
f.colorbar(im, ax=bx1)
bx1.set_title('best approximation amplitude - physical constraint satisfied')
im = bx2.imshow(np.real(u), cmap='gray')
f.colorbar(im, ax=bx2)
bx2.set_title('best approximation phase - physical constraint satisfied')
bx3.semilogy(change)
bx3.set_xlabel('iteration')
bx3.set_ylabel('$||x^{2k+2}-x^{2k}||$')
if 'diagnostic' in config:
bx4.semilogy(output['gap'])
bx4.set_xlabel('iteration')
bx4.set_ylabel('$||x^{2k+1}-x^{2k}||$')
g, ((bx1, bx2), (bx3, bx4)) = subplots(2, 2)
im=bx1.imshow(np.abs(u),cmap='gray')
f.colorbar(im, ax=bx1)
bx1.set_title('best approximation amplitude - physical constraint satisfied')
im = bx2.imshow(np.real(u), cmap='gray')
f.colorbar(im, ax=bx2)
bx2.set_title('best approximation phase - physical constraint satisfied')
bx3.semilogy(change)
bx3.set_xlabel('iteration')
bx3.set_ylabel('$||x^{2k+2}-x^{2k}||$')
if 'diagnostic' in config:
bx4.semilogy(output['gap'])
bx4.set_xlabel('iteration')
bx4.set_ylabel('$||x^{2k+1}-x^{2k}||$')
else:
f, (ax1,ax2) = subplots(2,1)
ax1.semilogy(change)
ax1.set_xlabel('iteration')
ax1.set_ylabel('$||x^{2k+2}-x^{2k}||$')
if 'diagnostic' in config:
ax2.semilogy(output['gap'])
ax2.set_xlabel('iteration')
ax2.set_ylabel('$||x^{2k+1}-x^{2k}||$')
show()
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