P_orthonorm.py 3.05 KB
Newer Older
jansen31's avatar
jansen31 committed
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's avatar
jansen31 committed
8
9
from numpy import sum, array, sqrt, zeros_like

10
11
from proxtoolbox.proxoperators.proxoperator import ProxOperator

jansen31's avatar
jansen31 committed
12
13
__all__ = ['P_orthonorm', "P_norm"]

14
15
16
17

class P_orthonorm(ProxOperator):
    def __init__(self, experiment):
        """
jansen31's avatar
jansen31 committed
18
19
        Normalize two iterates and rotate them such that they are orthogonal (sum(u.v) == 0).
        Applies a real-valuedness projection first.
20
21
22
23
24
25
26

        @param experiment: experiment class
        """
        super(P_orthonorm, self).__init__(experiment)

        self.p_norm = P_norm(experiment)

27
    def eval(self, u, prox_idx=None):
28
        # normalize u[0] and u[1]
29
30
        u_norm = self.p_norm.eval(u.real)
        norms = [sqrt(sum(abs(u_i) ** 2)) for u_i in u_norm]
31
        # determine angle _a_ between u[0] and u[1]
jansen31's avatar
jansen31 committed
32
        a = sum(u_norm[0] * u_norm[1]) / (norms[0] * norms[1])
jansen31's avatar
jansen31 committed
33
        if a == 0:  # for already orthogonal iterates, apply no change
jansen31's avatar
jansen31 committed
34
            return u_norm
jansen31's avatar
jansen31 committed
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's avatar
jansen31 committed
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's avatar
jansen31 committed
49
        return u_new
50
51
52
53
54


class P_norm(ProxOperator):
    def __init__(self, experiment):
        """
jansen31's avatar
jansen31 committed
55
        Normalize iterates,  such the incoherent sum [sum(abs(u)**2, axis=0)] adds up to experiment.norm
56
57
58
        """
        super(P_norm, self).__init__(experiment)
        self.norm = experiment.norm_data
jansen31's avatar
jansen31 committed
59
        # Each incoherent component is normalized to norm/sqrt(n) where n is the number of componentsss
60

61
    def eval(self, u, prox_idx=None):
62
        # Normalize components of u
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's avatar
jansen31 committed
66
67


jansen31's avatar
jansen31 committed
68
# Basic test functionality
jansen31's avatar
jansen31 committed
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's avatar
jansen31 committed
84
    checksum = []
jansen31's avatar
jansen31 committed
85
    thlist = np.arange(0.11, 6.01, 0.1)
jansen31's avatar
jansen31 committed
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")