P_Sparsity.py 3.29 KB
Newer Older
1
from numpy import zeros_like, unravel_index
robin.requadt's avatar
robin.requadt committed
2
import numpy as np
jansen31's avatar
jansen31 committed
3
from .proxoperators import ProxOperator
jansen31's avatar
jansen31 committed
4
5


robin.requadt's avatar
robin.requadt committed
6
class P_Sparsity(ProxOperator):
7
8
9
    """
    Projection subroutine for projecting onto support constraints
    USAGE: p_Sparsity = P_Sparsity(config)
10
           p_Sparsity.work(u)
11
    """
jansen31's avatar
jansen31 committed
12
13
14
15
    def __init__(self, config):
        """
        Initialization

robin.requadt's avatar
robin.requadt committed
16
17
        Parameters
        ----------
18
        config : dict - Dictionary containing the problem configuration. It must contain the following mappings:
jansen31's avatar
jansen31 committed
19
        'sparsity_parameter' : int
20
        'Ny' : int
jansen31's avatar
jansen31 committed
21
22
        """
        self.sparsity_parameter = config['sparsity_parameter']
23
24
        self.ny = config['Ny']

jansen31's avatar
jansen31 committed
25
26
27
28
29
        if 'sparsity_support' in config:
            self.support = config['sparsity_support'].real.astype(np.uint)
        else:
            self.support = 1

30
31
        if self.sparsity_parameter > 30:
            def value_selection(original, indices, sparsity_parameter):
32
                idx_for_threshold = unravel_index(indices[-sparsity_parameter], original.shape)
33
34
35
36
37
38
                threshold_val = abs(original[idx_for_threshold].get())
                return (abs(original) >= threshold_val) * original
        else:
            def value_selection(original, indices, sparsity_parameter):
                out = zeros_like(original)
                hits = indices[-sparsity_parameter:].get()
39
                hit_idx = [unravel_index(hit, original.shape) for hit in hits]
40
41
42
43
44
                for _idx in hit_idx:
                    out[_idx[0], _idx[1]] = original[_idx[0], _idx[1]]
                return out

        self.value_selection = value_selection
jansen31's avatar
jansen31 committed
45

jansen31's avatar
jansen31 committed
46
47
48
    def work(self, u):
        """
        Parameters
robin.requadt's avatar
robin.requadt committed
49
        ----------
50
        u : array_like - Input, array to be projected
jansen31's avatar
jansen31 committed
51

robin.requadt's avatar
robin.requadt committed
52
53
        Returns
        -------
54
        p_Sparsity : array_like, the projection
jansen31's avatar
jansen31 committed
55
        """
jansen31's avatar
jansen31 committed
56
        u *= self.support # apply support
jansen31's avatar
jansen31 committed
57
        p_Sparsity = 0 * u
58
        sorting = np.argsort(abs(u), axis=None)  # gives indices of sorted array in ascending order
59
        indices = np.asarray([unravel_index(sorting[i], u) for i in range(-1 * self.sparsity_parameter, 0)])
60
        p_Sparsity[indices[:, 0], indices[:, 1]] = u[indices[:, 0], indices[:, 1]]
robin.requadt's avatar
robin.requadt committed
61

jansen31's avatar
jansen31 committed
62
        return p_Sparsity
63
64
65
66
67
68
69
70


class P_Sparsity_real(P_Sparsity):
    def __init__(self, config):
        super(P_Sparsity_real, self).__init__(config)

    def work(self, u):
        return super(P_Sparsity_real, self).work(u.real)
71
72
73
74
75


class P_Sparsity_Symmetric(P_Sparsity):
    def __init__(self, config):
        super(P_Sparsity_Symmetric, self).__init__(config=config)
jansen31's avatar
jansen31 committed
76
77
        self.symmetry = config['symmetry_type'] # antisymmetric = -1, symmetric = 1
        self.symmetry_axis = config['symmetry_axis']  # -1 for last, 0 for first, 'both', 'all' or None for full flip
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94

    def work(self, u):
        if self.symmetry_axis in ['both', 'all']:
            mirror = np.flip(u, axis=None)
        else:
            mirror = np.flip(u, axis=self.symmetry_axis)
        inp = (u + self.symmetry * mirror) / 2
        out = super(P_Sparsity_Symmetric, self).work(inp)
        return out


class P_Sparsity_Symmetric_real(P_Sparsity_Symmetric):
    def __init__(self, config):
        super(P_Sparsity_Symmetric_real, self).__init__(config=config)

    def work(self, u):
        out = super(P_Sparsity_Symmetric_real, self).work(u.real)
95
        return out