SimpleAlgortihm.py 8.96 KB
Newer Older
jansen31's avatar
jansen31 committed
1
2
# Simple Algortihm
# Analog to AlgortihmWrapper in ProxMatlab
3
4
5

from .algorithms import Algorithm

jansen31's avatar
jansen31 committed
6
from numpy import zeros, angle, trace, exp, sqrt, sum
7
8
from numpy.linalg import norm

jansen31's avatar
jansen31 committed
9

jansen31's avatar
jansen31 committed
10
11
12
13
14
def phase_offset_compensated_norm(arr1, arr2, norm_factor=1, norm_type='fro'):
    phase_offset = angle(sum(arr1 * arr2.conj()))
    return norm(arr1 - arr2 * exp(1j * phase_offset), norm_type) / norm_factor


15
16
class SimpleAlgorithm:

jansen31's avatar
jansen31 committed
17
    def __init__(self, config):
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
        """
        Parameters
        ----------
        config : dict        
         Dictionary containing the problem configuration.
         It must contain the following mappings:

        proxoperators: 2 ProxOperators
            Tuple of ProxOperators (the class, no instance)
        beta_0: number
            Starting relaxation parmater
        beta_max: number
            Maximum relaxation parameter
        beta_switch: int
            Iteration at which beta moves from beta_0 -> beta_max
33
        norm_data: number
34
35
36
37
38
39
40
41
42
43
44
45
46
            ?
        Nx: int
            Row-dim of the product space elements
        Ny: int
            Column-dim of the product space elements
        Nz: int
            Depth-dim of the product space elements
        dim: int
            Size of the product space
        """

        self.prox1 = config['proxoperators'][0](config)
        self.prox2 = config['proxoperators'][1](config)
47
        self.norm_data = config['norm_data']
jansen31's avatar
jansen31 committed
48
49
        self.Nx = config['Nx']
        self.Ny = config['Ny']
jansen31's avatar
jansen31 committed
50
        self.Nz = config['Nz']
51
52
53
54
55
        self.product_space_dimension = config['product_space_dimension']
        self.iter = 0
        self.config = config

        if 'truth' in config:
56
57
58
            self.truth = config['truth']
            self.truth_dim = config['truth_dim']
            self.norm_truth = config['norm_truth']
59
60

        if 'diagnostic' in config:
61
            self.diagnostic = config['diagnostic']
62
63
64
        else:
            self.diagnostic = False

jansen31's avatar
jansen31 committed
65
66
67
68
69
70
71
72
    def evaluate(self, u):
        """
        Method to override in specific algorithm subclass
        :param u: old iterate
        :return: new iterate
        """
        return u

73
74
75
76
77
78
79
    def run(self, u, tol, maxiter):
        """
        Runs the algorithm for the specified input data
        """

        ##### PREPROCESSING
        config = self.config
80
        norm_data = self.norm_data
81
82

        iter = self.iter
83
84
85
86

        # TODO: select the right case also for a 3d problem.
        #  This could involve using the parameter self.product_space_dimension
        if u.ndim < 3:  # temp note: as used in (e.g.) 2d CDI, or orbital tomography
87
88
89
90
91
92
93
94
95
            p = 1
            q = 1
        elif u.ndim == 3:
            p = u.shape[2]
            q = 1
        else:
            p = u.shape[2]
            q = u.shape[3]

jansen31's avatar
jansen31 committed
96
        change = zeros(maxiter + 1, dtype=u.dtype)
97
98
99
100
101
102
103
104
        change[0] = 999

        ######################################
        #  set up diagnostic arrays
        ######################################
        if self.diagnostic:
            gap = change.copy()
            shadow_change = change.copy()
jansen31's avatar
jansen31 committed
105
        if hasattr(self, 'truth'):
106
107
            Relerrs = change.copy()

108
        # Different algorithms will have different qualities that one
109
110
111
112
113
114
115
        # should monitor.  Since we take a fixed point perspective, the main
        # thing to monitor is the change in the iterates.
        shadow = self.prox1.work(u)

        ##### LOOP
        while iter < maxiter and change[iter] >= tol:
            iter += 1
jansen31's avatar
jansen31 committed
116
117
118
            config['iter'] = iter
            # makes call to algorithm:
            config['u'] = u
119
120
            u_new = self.evaluate(u)

121
            if self.diagnostic:
122
123
124
125
126
127
128
                # the next prox operation only gets used in the computation of
                # the size of the gap.  This extra step is not
                # required in alternating projections, which makes RAAR
                u2 = self.prox2.work(u_new)
                u1 = self.prox1.work(u2)

            # compute the normalized change in successive iterates: 
129
            # change(iter) = sum(sum((feval('P_M',M,u)-tmp).^2))/norm_data;
jansen31's avatar
jansen31 committed
130
131
132
            tmp_change = 0
            tmp_gap = 0
            tmp_shadow = 0
133

134
135
            # the simple case, where everything can be calculated in one difference
            if self.product_space_dimension == 1 or (p == 1 and q == 1):
jansen31's avatar
jansen31 committed
136
                tmp_change = phase_offset_compensated_norm(u, u_new, norm_type='fro', norm_factor=norm_data) ** 2
137
                if self.diagnostic:
138
139
140
                    # For Douglas-Rachford,in general it is appropriate to monitor the SHADOWS of the iterates,
                    # since in the convex case these converge even for beta=1. (see Bauschke-Combettes-Luke,
                    # J. Approx. Theory, 2004)
jansen31's avatar
jansen31 committed
141
142
                    tmp_shadow = phase_offset_compensated_norm(u2, shadow, norm_factor=norm_data, norm_type='fro') ** 2
                    tmp_gap = phase_offset_compensated_norm(u1, u2, norm_factor=norm_data, norm_type='fro') ** 2
143
                    if hasattr(self, 'truth'):
144
                        if self.truth_dim[0] == 1:
jansen31's avatar
jansen31 committed
145
                            z = u1[0, :]
146
                        elif self.truth_dim[1] == 1:
jansen31's avatar
jansen31 committed
147
148
149
                            z = u1[:, 0]
                        else:
                            z = u1
150
                        Relerrs[iter] = phase_offset_compensated_norm(self.truth, z, norm_factor=self.norm_truth,
jansen31's avatar
jansen31 committed
151
                                                                      norm_type='fro')
jansen31's avatar
jansen31 committed
152

153
            elif q == 1:  # more complex: loop over product space dimension
154
                for j in range(self.product_space_dimension):
jansen31's avatar
jansen31 committed
155
                    tmp_change = tmp_change + (norm(u[:, :, j] - u_new[:, :, j], 'fro') / norm_data) ** 2
156
                    if self.diagnostic:
157
                        # compute (||P_SP_Mx-P_Mx||/norm_data)^2:
jansen31's avatar
jansen31 committed
158
159
                        tmp_gap = tmp_gap + (norm(u1[:, :, j] - u2[:, :, j], 'fro') / norm_data) ** 2
                        tmp_shadow = tmp_shadow + (norm(u2[:, :, j] - shadow[:, :, j], 'fro') / norm_data) ** 2
160
                if self.diagnostic:
161
                    if hasattr(self, 'truth'):
jansen31's avatar
jansen31 committed
162
                        z = u1[:, :, 0]
163
164
                        Relerrs[iter] = norm((self.truth - exp(-1j * angle(trace(self.truth.T.transpose() * z))) * z),
                                             'fro') / self.norm_truth
165

166
            else:  # loop over product space dimension as well as additional z dimension
167
                if self.diagnostic:
168
169
                    if hasattr(self, 'truth'):
                        Relerrs[iter] = 0
170
                for j in range(self.product_space_dimension):
jansen31's avatar
jansen31 committed
171
172
                    for k in range(self.Nz):
                        tmp_change = tmp_change + (norm(u[:, :, k, j] - u_new[:, :, k, j], 'fro') / norm_data) ** 2
173
                        if self.diagnostic:
jansen31's avatar
jansen31 committed
174
175
176
                            # compute (||P_Sx-P_Mx||/norm_data)^2:
                            tmp_gap = tmp_gap + (norm(u1[:, :, k, j] - u2[:, :, k, j], 'fro') / (norm_data)) ** 2
                            tmp_shadow = tmp_shadow + (
jansen31's avatar
jansen31 committed
177
                                    norm(u2[:, :, k, j] - shadow[:, :, k, j], 'fro') / (norm_data)) ** 2
jansen31's avatar
jansen31 committed
178
179
                if hasattr(self, 'truth') and (j == 0):
                    Relerrs[iter] = Relerrs[iter] + norm(
180
                        self.truth - exp(-1j * angle(trace(self.truth.T * u1[:, :, k, 1]))) * u1[:, :, k, 1],
181
                        'fro') / self.norm_Truth
182

183
            # Add values to the list of change, gap and shadow_change
184
            change[iter] = sqrt(tmp_change)
185
            if self.diagnostic:
186
                gap[iter] = sqrt(tmp_gap)
jansen31's avatar
jansen31 committed
187
                shadow_change[iter] = sqrt(tmp_shadow)  # this is the Euclidean norm of the gap to
188
189
190
                # the unregularized set. To monitor the Euclidean norm of the gap to the regularized set is
                # expensive to calculate, so we use this surrogate. Since the stopping criteria is on the change
                # in the iterates, this does not matter.
191

192
            # update iterate
jansen31's avatar
jansen31 committed
193
            u = u_new
194
            if self.diagnostic:
195
196
                # For Douglas-Rachford,in general it is appropriate to monitor the SHADOWS of the iterates, since in
                # the convex case these converge even for beta=1.
197
                # (see Bauschke-Combettes-Luke, J. Approx. Theory, 2004)
jansen31's avatar
jansen31 committed
198
                shadow = u2
199
200

        ##### POSTPROCESSING
jansen31's avatar
jansen31 committed
201
202
        u2 = self.prox2.work(u)
        u1 = self.prox1.work(u2)
203
        if self.Nx == 1:
jansen31's avatar
jansen31 committed
204
205
            u1 = u1[:, 0]
            u2 = u2[:, 0]
206
        elif self.Ny == 1:
jansen31's avatar
jansen31 committed
207
208
            u1 = u1[0, :]
            u2 = u2[0, :]
209
        elif self.Nz == 1 and u1.ndim > 2:
jansen31's avatar
jansen31 committed
210
211
            u1 = u1[:, :, 0]
            u2 = u2[:, :, 0]
212

jansen31's avatar
jansen31 committed
213
        change = change[1:iter + 1]
214

jansen31's avatar
jansen31 committed
215
        output = {'u': u, 'u1': u1, 'u2': u2, 'iter': iter, 'change': change}
216
        if self.diagnostic:
jansen31's avatar
jansen31 committed
217
218
219
220
            gap = gap[1:iter + 1]
            shadow_change = shadow_change[1:iter + 1]
            output['gap'] = gap
            output['shadow_change'] = shadow_change
221
            if hasattr(self, 'truth'):
jansen31's avatar
jansen31 committed
222
                output['Relerrs'] = Relerrs
223
        return output