P_orthonorm.py

"""
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.
"""
from numpy import sum, array, sqrt, zeros_like

from proxtoolbox.proxoperators.proxoperator import ProxOperator

__all__ = ['P_orthonorm', "P_norm"]


class P_orthonorm(ProxOperator):
    def __init__(self, experiment):
        """
        Normalize two iterates and rotate them such that they are orthogonal (sum(u.v) == 0).
        Applies a real-valuedness projection first.

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

        self.p_norm = P_norm(experiment)

    def eval(self, u, prox_idx=None):
        # normalize u[0] and u[1]
        u_norm = self.p_norm.eval(u.real)
        norms = [sqrt(sum(abs(u_i) ** 2)) for u_i in u_norm]
        # determine angle _a_ between u[0] and u[1]
        a = sum(u_norm[0] * u_norm[1]) / (norms[0] * norms[1])
        if a == 0:  # for already orthogonal iterates, apply no change
            return u_norm
        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)
        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])
        if a < 0:  # Not sure if or why this is needed...
            u_new *= -1
        return u_new


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

    def eval(self, u, prox_idx=None):
        # Normalize components of u
        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)])


# Basic test functionality
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)

    checksum = []
    thlist = np.arange(0.11, 6.01, 0.1)
    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")