P_orthonorm.py 3.05 KB
 jansen31 committed Nov 04, 2020 1 2 3 4 5 6 7 ``````""" Written by Matthijs Jansen, October 2020 together with Russell Luke Contains a proxoperator which takes a set of two real-valued vectors (iterate shaped like [[i_1,j_1, .., z_1], [i_2, j_2, .., z_2]] and rotates these such that they are orthogonal. Before this is done, the vectors are normalized and projected onto real values. """ `````` jansen31 committed Oct 26, 2020 8 9 ``````from numpy import sum, array, sqrt, zeros_like `````` jansen31 committed Oct 26, 2020 10 11 ``````from proxtoolbox.proxoperators.proxoperator import ProxOperator `````` jansen31 committed Oct 26, 2020 12 13 ``````__all__ = ['P_orthonorm', "P_norm"] `````` jansen31 committed Oct 26, 2020 14 15 16 17 `````` class P_orthonorm(ProxOperator): def __init__(self, experiment): """ `````` jansen31 committed Nov 03, 2020 18 19 `````` Normalize two iterates and rotate them such that they are orthogonal (sum(u.v) == 0). Applies a real-valuedness projection first. `````` jansen31 committed Oct 26, 2020 20 21 22 23 24 25 26 `````` @param experiment: experiment class """ super(P_orthonorm, self).__init__(experiment) self.p_norm = P_norm(experiment) `````` jansen31 committed Nov 02, 2020 27 `````` def eval(self, u, prox_idx=None): `````` jansen31 committed Oct 26, 2020 28 `````` # normalize u[0] and u[1] `````` jansen31 committed Nov 02, 2020 29 30 `````` u_norm = self.p_norm.eval(u.real) norms = [sqrt(sum(abs(u_i) ** 2)) for u_i in u_norm] `````` jansen31 committed Oct 26, 2020 31 `````` # determine angle _a_ between u[0] and u[1] `````` jansen31 committed Oct 26, 2020 32 `````` a = sum(u_norm[0] * u_norm[1]) / (norms[0] * norms[1]) `````` jansen31 committed Oct 26, 2020 33 `````` if a == 0: # for already orthogonal iterates, apply no change `````` jansen31 committed Oct 26, 2020 34 `````` return u_norm `````` jansen31 committed Oct 26, 2020 35 36 37 38 39 40 41 42 43 44 45 46 `````` elif a == -1 or a == 1: # Cannot do anything for parallel vectors return u_norm # determine roots of y^2 - 2/a y + 1 = 0 elif a > 0: y = 1 / a + sqrt(1 / a ** 2 - 1) elif a < 0: y = 1 / a - sqrt(1 / a ** 2 - 1) else: raise Exception("This should never rise, check calculation of a") # apply projection u_new = zeros_like(u_norm) `````` jansen31 committed Nov 04, 2020 47 48 `````` u_new[1] = u_norm[0] - (y / (y ** 2 - 1)) * (y * u_norm[0] - u_norm[1]) u_new[0] = (1 / (y ** 2 - 1)) * (y * u_norm[0] - u_norm[1]) `````` jansen31 committed Oct 26, 2020 49 `````` return u_new `````` jansen31 committed Oct 26, 2020 50 51 52 53 54 `````` class P_norm(ProxOperator): def __init__(self, experiment): """ `````` jansen31 committed Oct 26, 2020 55 `````` Normalize iterates, such the incoherent sum [sum(abs(u)**2, axis=0)] adds up to experiment.norm `````` jansen31 committed Oct 26, 2020 56 57 58 `````` """ super(P_norm, self).__init__(experiment) self.norm = experiment.norm_data `````` jansen31 committed Nov 03, 2020 59 `````` # Each incoherent component is normalized to norm/sqrt(n) where n is the number of componentsss `````` jansen31 committed Oct 26, 2020 60 `````` `````` jansen31 committed Nov 02, 2020 61 `````` def eval(self, u, prox_idx=None): `````` jansen31 committed Oct 26, 2020 62 `````` # Normalize components of u `````` jansen31 committed Nov 02, 2020 63 64 65 `````` norms = array([sqrt(sum(abs(u_n) ** 2)) for u_n in u]) n_comp = len(u) return array([self.norm * u[i] / norms[i] / sqrt(n_comp) for i in range(n_comp)]) `````` jansen31 committed Oct 26, 2020 66 67 `````` `````` jansen31 committed Nov 03, 2020 68 ``````# Basic test functionality `````` jansen31 committed Oct 26, 2020 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 ``````if __name__ == "__main__": import numpy as np class DummyExperiment: def __init__(self, data=1, norm=sqrt(2)): self.data = data self.norm_data = norm exp = DummyExperiment() portho = P_orthonorm(exp) pnorm = P_norm(exp) `````` jansen31 committed Oct 26, 2020 84 `````` checksum = [] `````` jansen31 committed Nov 04, 2020 85 `````` thlist = np.arange(0.11, 6.01, 0.1) `````` jansen31 committed Oct 26, 2020 86 87 88 89 90 91 `````` for th in thlist: inp = np.array([[1, 0], [np.cos(th * np.pi), np.sin(th * np.pi)]]) out = portho.eval(inp) checksum += [np.sum(out[0] * out[1])] assert np.sum(checksum) < 1e-14, "Test failed" print("Test passed")``````