SimpleAlgortihm.py 10.2 KB
Newer Older
jansen31's avatar
jansen31 committed
1
2
# Simple Algorithm
# Analog to AlgorithmWrapper in ProxMatlab
3

Matthijs's avatar
Matthijs committed
4
from numpy import zeros, angle, trace, exp, sqrt, sum, ndarray
5
6
from numpy.linalg import norm

jansen31's avatar
jansen31 committed
7

Matthijs's avatar
Matthijs committed
8
9
10
11
def nd_norm_difference(arr1: ndarray, arr2: ndarray, norm_factor: float = 1.):
    """
    Root-sum-square norm of the difference of 2 numpy arrays, equal to Frobenius norm for arr.ndim == 2
    """
12
    return sqrt(sum(abs(arr1 - arr2) ** 2)) / norm_factor
Matthijs's avatar
Matthijs committed
13
14
15
16
17
18
19
20
21
22
23
24
25


def phase_offset_compensated_norm(arr1: ndarray, arr2: ndarray, norm_factor: float = 1.,
                                  norm_type: str = 'fro') -> float:
    """
    Norm of the difference between 2 numpy arrays.

    :param arr1: numpy array
    :param arr2: numpy array like arr1
    :param norm_factor: normalization factor
    :param norm_type: for 1d and 2d arrays, pass norm_type to numpy.linalg.norm
    :return: float
    """
26
27
28
29
30
    if arr1.ndim < 3:
        phase_offset = angle(sum(arr1 * arr2.conj()))
        return norm(arr1 - arr2 * exp(1j * phase_offset), norm_type) / norm_factor
    else:
        phase_offset = angle(sum(arr1 * arr2.conj()))
Matthijs's avatar
Matthijs committed
31
        return nd_norm_difference(arr1, arr2 * exp(1j * phase_offset), norm_factor=norm_factor)
jansen31's avatar
jansen31 committed
32
33


34
35
class SimpleAlgorithm:

jansen31's avatar
jansen31 committed
36
    def __init__(self, config):
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
        """
        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
52
        norm_data: number
53
54
55
56
57
58
59
60
61
62
63
64
65
            ?
        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)
66
        self.norm_data = config['norm_data']
jansen31's avatar
jansen31 committed
67
68
        self.Nx = config['Nx']
        self.Ny = config['Ny']
jansen31's avatar
jansen31 committed
69
        self.Nz = config['Nz']
70
71
72
73
74
        self.product_space_dimension = config['product_space_dimension']
        self.iter = 0
        self.config = config

        if 'truth' in config:
75
76
77
            self.truth = config['truth']
            self.truth_dim = config['truth_dim']
            self.norm_truth = config['norm_truth']
78
79

        if 'diagnostic' in config:
80
            self.diagnostic = config['diagnostic']
81
82
83
        else:
            self.diagnostic = False

jansen31's avatar
jansen31 committed
84
85
86
87
88
89
90
91
    def evaluate(self, u):
        """
        Method to override in specific algorithm subclass
        :param u: old iterate
        :return: new iterate
        """
        return u

92
93
94
    def run(self, u, tol, maxiter):
        """
        Runs the algorithm for the specified input data
Matthijs's avatar
Matthijs committed
95
96
97
98
99

        :param u: iterate
        :param tol: tolerance for convergence, if change<tol, stop
        :param maxiter:
        :return:
100
101
102
        """
        ##### PREPROCESSING
        config = self.config
103
        norm_data = self.norm_data
104
105

        iter = self.iter
106
107
108
109

        # 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
110
111
112
113
114
115
116
117
118
            p = 1
            q = 1
        elif u.ndim == 3:
            p = u.shape[2]
            q = 1
        else:
            p = u.shape[2]
            q = u.shape[3]

119
        change = zeros(maxiter + 1)
120
121
122
123
124
125
126
127
        change[0] = 999

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

131
132
        # Different algorithms will have different qualities that one should monitor.
        # Since we take a fixed point perspective, the main thing to monitor is the change in the iterates.
133
134
135
        shadow = self.prox1.work(u)

        ##### LOOP
136
137
138
        if 'tqdm_progressbar' in self.config and self.config['tqdm_progressbar'] is not None:
            pbar = self.config['tqdm_progressbar'](total=maxiter)

139
140
        while iter < maxiter and change[iter] >= tol:
            iter += 1
jansen31's avatar
jansen31 committed
141
            config['iter'] = iter
142
143
            if 'tqdm_progressbar' in self.config and self.config['tqdm_progressbar'] is not None:
                pbar.update()
jansen31's avatar
jansen31 committed
144
145
            # makes call to algorithm:
            config['u'] = u
146
147
            u_new = self.evaluate(u)

148
            if self.diagnostic:
149
150
                # 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 a bit slower
151
152
153
                u2 = self.prox2.work(u_new)
                u1 = self.prox1.work(u2)

154
            # compute the normalized change in successive iterates
jansen31's avatar
jansen31 committed
155
156
157
            tmp_change = 0
            tmp_gap = 0
            tmp_shadow = 0
158
159
            # 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
160
                tmp_change = phase_offset_compensated_norm(u, u_new, norm_type='fro', norm_factor=norm_data) ** 2
161
                if self.diagnostic:
162
163
164
                    # 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
165
166
                    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
167
                    if hasattr(self, 'truth'):
168
                        if self.truth_dim[0] == 1:
jansen31's avatar
jansen31 committed
169
                            z = u1[0, :]
170
                        elif self.truth_dim[1] == 1:
jansen31's avatar
jansen31 committed
171
172
173
                            z = u1[:, 0]
                        else:
                            z = u1
174
                        Relerrs[iter] = phase_offset_compensated_norm(self.truth, z, norm_factor=self.norm_truth,
jansen31's avatar
jansen31 committed
175
                                                                      norm_type='fro')
jansen31's avatar
jansen31 committed
176

177
            elif q == 1:  # more complex: loop over product space dimension
Matthijs's avatar
Matthijs committed
178
                for j in range(self.product_space_dimension):  # this loop might be avoided using the 3d norm function
jansen31's avatar
jansen31 committed
179
                    tmp_change = tmp_change + (norm(u[:, :, j] - u_new[:, :, j], 'fro') / norm_data) ** 2
180
                    if self.diagnostic:
181
                        # compute (||P_SP_Mx-P_Mx||/norm_data)^2:
jansen31's avatar
jansen31 committed
182
183
                        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
184
                if self.diagnostic:
185
                    if hasattr(self, 'truth'):
jansen31's avatar
jansen31 committed
186
                        z = u1[:, :, 0]
187
188
                        Relerrs[iter] = norm((self.truth - exp(-1j * angle(trace(self.truth.T.transpose() * z))) * z),
                                             'fro') / self.norm_truth
189

190
            else:  # loop over product space dimension as well as additional z dimension
191
                if self.diagnostic:
192
193
                    if hasattr(self, 'truth'):
                        Relerrs[iter] = 0
Matthijs's avatar
Matthijs committed
194
195
                for j in range(self.product_space_dimension):  # this loop might be avoided using a 4d norm function
                    for k in range(self.Nz):  # this loop might be avoided using the 3d norm function
jansen31's avatar
jansen31 committed
196
                        tmp_change = tmp_change + (norm(u[:, :, k, j] - u_new[:, :, k, j], 'fro') / norm_data) ** 2
197
                        if self.diagnostic:
jansen31's avatar
jansen31 committed
198
199
200
                            # 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
201
                                    norm(u2[:, :, k, j] - shadow[:, :, k, j], 'fro') / (norm_data)) ** 2
jansen31's avatar
jansen31 committed
202
203
                if hasattr(self, 'truth') and (j == 0):
                    Relerrs[iter] = Relerrs[iter] + norm(
204
                        self.truth - exp(-1j * angle(trace(self.truth.T * u1[:, :, k, 1]))) * u1[:, :, k, 1],
205
                        'fro') / self.norm_Truth
206

207
            # Add values to the list of change, gap and shadow_change
208
            change[iter] = sqrt(tmp_change)
209
            if self.diagnostic:
210
                gap[iter] = sqrt(tmp_gap)
jansen31's avatar
jansen31 committed
211
                shadow_change[iter] = sqrt(tmp_shadow)  # this is the Euclidean norm of the gap to
212
213
214
                # 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.
215

216
            # update iterate
jansen31's avatar
jansen31 committed
217
            u = u_new
218
            if self.diagnostic:
219
220
                # 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
221
                shadow = u2
222
223

        ##### POSTPROCESSING
jansen31's avatar
jansen31 committed
224
225
        u2 = self.prox2.work(u)
        u1 = self.prox1.work(u2)
226
        if self.Nx == 1:
jansen31's avatar
jansen31 committed
227
228
            u1 = u1[:, 0]
            u2 = u2[:, 0]
229
        elif self.Ny == 1:
jansen31's avatar
jansen31 committed
230
231
            u1 = u1[0, :]
            u2 = u2[0, :]
232
        elif self.Nz == 1 and u1.ndim > 2:
jansen31's avatar
jansen31 committed
233
234
            u1 = u1[:, :, 0]
            u2 = u2[:, :, 0]
235

jansen31's avatar
jansen31 committed
236
237
        change = change[1:iter + 1]
        output = {'u': u, 'u1': u1, 'u2': u2, 'iter': iter, 'change': change}
238
        if self.diagnostic:
jansen31's avatar
jansen31 committed
239
240
241
242
            gap = gap[1:iter + 1]
            shadow_change = shadow_change[1:iter + 1]
            output['gap'] = gap
            output['shadow_change'] = shadow_change
243
            if hasattr(self, 'truth'):
jansen31's avatar
jansen31 committed
244
                output['Relerrs'] = Relerrs
245
        return output