SimpleAlgortihm.py 8.96 KB
 jansen31 committed May 10, 2019 1 2 ``````# Simple Algortihm # Analog to AlgortihmWrapper in ProxMatlab `````` 3 4 5 `````` from .algorithms import Algorithm `````` jansen31 committed Jun 12, 2019 6 ``````from numpy import zeros, angle, trace, exp, sqrt, sum `````` 7 8 ``````from numpy.linalg import norm `````` jansen31 committed May 10, 2019 9 `````` `````` jansen31 committed Jun 12, 2019 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 committed May 10, 2019 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 committed Jun 12, 2019 48 49 `````` self.Nx = config['Nx'] self.Ny = config['Ny'] `````` jansen31 committed May 10, 2019 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: `````` markus.meier01 committed Dec 15, 2019 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: `````` Matthijs committed Apr 21, 2020 61 `````` self.diagnostic = config['diagnostic'] `````` 62 63 64 `````` else: self.diagnostic = False `````` jansen31 committed Jun 12, 2019 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 committed May 10, 2019 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 committed May 10, 2019 105 `````` if hasattr(self, 'truth'): `````` 106 107 `````` Relerrs = change.copy() `````` Matthijs committed Apr 21, 2020 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 committed May 10, 2019 116 117 118 `````` config['iter'] = iter # makes call to algorithm: config['u'] = u `````` 119 120 `````` u_new = self.evaluate(u) `````` Matthijs committed Apr 21, 2020 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 committed May 10, 2019 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 committed Jun 12, 2019 136 `````` tmp_change = phase_offset_compensated_norm(u, u_new, norm_type='fro', norm_factor=norm_data) ** 2 `````` Matthijs committed Apr 21, 2020 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 committed Jun 12, 2019 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'): `````` markus.meier01 committed Dec 15, 2019 144 `````` if self.truth_dim[0] == 1: `````` jansen31 committed May 10, 2019 145 `````` z = u1[0, :] `````` markus.meier01 committed Dec 15, 2019 146 `````` elif self.truth_dim[1] == 1: `````` jansen31 committed May 10, 2019 147 148 149 `````` z = u1[:, 0] else: z = u1 `````` markus.meier01 committed Dec 15, 2019 150 `````` Relerrs[iter] = phase_offset_compensated_norm(self.truth, z, norm_factor=self.norm_truth, `````` jansen31 committed Jun 12, 2019 151 `````` norm_type='fro') `````` jansen31 committed May 10, 2019 152 `````` `````` 153 `````` elif q == 1: # more complex: loop over product space dimension `````` 154 `````` for j in range(self.product_space_dimension): `````` jansen31 committed May 10, 2019 155 `````` tmp_change = tmp_change + (norm(u[:, :, j] - u_new[:, :, j], 'fro') / norm_data) ** 2 `````` Matthijs committed Apr 21, 2020 156 `````` if self.diagnostic: `````` 157 `````` # compute (||P_SP_Mx-P_Mx||/norm_data)^2: `````` jansen31 committed May 10, 2019 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 `````` Matthijs committed Apr 21, 2020 160 `````` if self.diagnostic: `````` 161 `````` if hasattr(self, 'truth'): `````` jansen31 committed May 10, 2019 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 `````` markus.meier01 committed Dec 15, 2019 165 `````` `````` 166 `````` else: # loop over product space dimension as well as additional z dimension `````` Matthijs committed Apr 21, 2020 167 `````` if self.diagnostic: `````` alexander.dornheim committed Dec 13, 2017 168 169 `````` if hasattr(self, 'truth'): Relerrs[iter] = 0 `````` 170 `````` for j in range(self.product_space_dimension): `````` jansen31 committed May 10, 2019 171 172 `````` for k in range(self.Nz): tmp_change = tmp_change + (norm(u[:, :, k, j] - u_new[:, :, k, j], 'fro') / norm_data) ** 2 `````` Matthijs committed Apr 21, 2020 173 `````` if self.diagnostic: `````` jansen31 committed May 10, 2019 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 committed Jun 12, 2019 177 `````` norm(u2[:, :, k, j] - shadow[:, :, k, j], 'fro') / (norm_data)) ** 2 `````` jansen31 committed May 10, 2019 178 179 `````` if hasattr(self, 'truth') and (j == 0): Relerrs[iter] = Relerrs[iter] + norm( `````` markus.meier01 committed Dec 15, 2019 180 `````` self.truth - exp(-1j * angle(trace(self.truth.T * u1[:, :, k, 1]))) * u1[:, :, k, 1], `````` markus.meier01 committed Dec 10, 2019 181 `````` 'fro') / self.norm_Truth `````` 182 `````` `````` 183 `````` # Add values to the list of change, gap and shadow_change `````` 184 `````` change[iter] = sqrt(tmp_change) `````` Matthijs committed Apr 21, 2020 185 `````` if self.diagnostic: `````` 186 `````` gap[iter] = sqrt(tmp_gap) `````` jansen31 committed May 10, 2019 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 committed May 10, 2019 193 `````` u = u_new `````` Matthijs committed Apr 21, 2020 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 committed May 10, 2019 198 `````` shadow = u2 `````` 199 200 `````` ##### POSTPROCESSING `````` jansen31 committed May 10, 2019 201 202 `````` u2 = self.prox2.work(u) u1 = self.prox1.work(u2) `````` 203 `````` if self.Nx == 1: `````` jansen31 committed May 10, 2019 204 205 `````` u1 = u1[:, 0] u2 = u2[:, 0] `````` 206 `````` elif self.Ny == 1: `````` jansen31 committed May 10, 2019 207 208 `````` u1 = u1[0, :] u2 = u2[0, :] `````` 209 `````` elif self.Nz == 1 and u1.ndim > 2: `````` jansen31 committed May 10, 2019 210 211 `````` u1 = u1[:, :, 0] u2 = u2[:, :, 0] `````` 212 `````` `````` jansen31 committed May 10, 2019 213 `````` change = change[1:iter + 1] `````` 214 `````` `````` jansen31 committed May 10, 2019 215 `````` output = {'u': u, 'u1': u1, 'u2': u2, 'iter': iter, 'change': change} `````` Matthijs committed Apr 21, 2020 216 `````` if self.diagnostic: `````` jansen31 committed May 10, 2019 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 committed May 10, 2019 222 `````` output['Relerrs'] = Relerrs `````` 223 `` return output``