Commit 7f014975 authored by Matthijs's avatar Matthijs
Browse files

higher-dimensional norm function

parent 7de32940
......@@ -3,17 +3,34 @@
from .algorithms import Algorithm
from numpy import zeros, angle, trace, exp, sqrt, sum
from numpy import zeros, angle, trace, exp, sqrt, sum, ndarray
from numpy.linalg import norm
def phase_offset_compensated_norm(arr1, arr2, norm_factor=1, norm_type='fro'):
def nd_norm_difference(arr1: ndarray, arr2: ndarray, norm_factor: float = 1.):
"""
Root-sum-square norm of the difference of 2 numpy arrays, equal to Frobenius norm for arr.ndim == 2
"""
return sqrt(sum((arr1 - arr2) ** 2)) / norm_factor
def phase_offset_compensated_norm(arr1: ndarray, arr2: ndarray, norm_factor: float = 1.,
norm_type: str = 'fro') -> float:
"""
Norm of the difference between 2 numpy arrays.
:param arr1: numpy array
:param arr2: numpy array like arr1
:param norm_factor: normalization factor
:param norm_type: for 1d and 2d arrays, pass norm_type to numpy.linalg.norm
:return: float
"""
if arr1.ndim < 3:
phase_offset = angle(sum(arr1 * arr2.conj()))
return norm(arr1 - arr2 * exp(1j * phase_offset), norm_type) / norm_factor
else:
phase_offset = angle(sum(arr1 * arr2.conj()))
return sqrt(sum((arr1 - arr2 * exp(1j * phase_offset))**2)) / norm_factor
return nd_norm_difference(arr1, arr2 * exp(1j * phase_offset), norm_factor=norm_factor)
class SimpleAlgorithm:
......@@ -155,7 +172,7 @@ class SimpleAlgorithm:
norm_type='fro')
elif q == 1: # more complex: loop over product space dimension
for j in range(self.product_space_dimension):
for j in range(self.product_space_dimension): # this loop might be avoided using the 3d norm function
tmp_change = tmp_change + (norm(u[:, :, j] - u_new[:, :, j], 'fro') / norm_data) ** 2
if self.diagnostic:
# compute (||P_SP_Mx-P_Mx||/norm_data)^2:
......@@ -171,8 +188,8 @@ class SimpleAlgorithm:
if self.diagnostic:
if hasattr(self, 'truth'):
Relerrs[iter] = 0
for j in range(self.product_space_dimension):
for k in range(self.Nz):
for j in range(self.product_space_dimension): # this loop might be avoided using a 4d norm function
for k in range(self.Nz): # this loop might be avoided using the 3d norm function
tmp_change = tmp_change + (norm(u[:, :, k, j] - u_new[:, :, k, j], 'fro') / norm_data) ** 2
if self.diagnostic:
# compute (||P_Sx-P_Mx||/norm_data)^2:
......
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