RAAR.py 6.08 KB
Newer Older
1
2
3
4
5
6
7
# -*- coding: utf-8 -*-
"""
Created on Mon Dec 14 12:48:26 2015

@author: rebecca
"""

8
from numpy import zeros, angle, trace, exp, sqrt
9

10
from numpy.linalg import norm
11
12
13
14
15
16
17
18
19
20

from .algorithms import Algorithm

class RAAR(Algorithm):
    """
    Relaxed Averaged Alternating Reflection algorithm
    """
    
    def __init__(self,config):
        """
rnahme's avatar
rnahme committed
21
22
23
24
25
        Parameters
        ----------
        config : dict        
                 Dictionary containing the problem configuration.
                 It must contain the following mappings:
26
            
27
28
                proxoperators: 2 ProxOperators
                    Tuple of ProxOperators (the class, no instance)
29
                beta_0: number
30
                    Starting relaxation parmater
rnahme's avatar
rnahme committed
31
                beta_max: number
32
                    Maximum relaxation parameter
rnahme's avatar
rnahme committed
33
                beta_switch: int
34
                    Iteration at which beta moves from beta_0 -> beta_max
rnahme's avatar
rnahme committed
35
                normM: number
36
                    ?
rnahme's avatar
rnahme committed
37
                Nx: int
38
                    Row-dim of the product space elements
rnahme's avatar
rnahme committed
39
                Ny: int
40
                    Column-dim of the product space elements
rnahme's avatar
rnahme committed
41
                Nz: int
42
                    Depth-dim of the product space elements
rnahme's avatar
rnahme committed
43
                dim: int
44
45
                    Size of the product space
        """
46
47
48

        self.prox1 = config['proxoperators'][0](config)
        self.prox2 = config['proxoperators'][1](config)
49
        self.normM = config['normM']
50
        self.beta_0 = config['beta_0']
51
52
        self.beta_max = config['beta_max']
        self.beta_switch = config['beta_switch']
53
        self.Nx = config['Nx']; self.Ny = config['Ny']; self.Nz = config['Nz'];
54
55
56
57
58
59
60
        self.product_space_dimension = config['product_space_dimension']
        self.iter = 0

        if 'truth' in config:
            self.truth = config['truth']
            self.truth_dim = config['truth_dim']
            self.norm_truth = config['norm_truth']
61
62

    def run(self, u, tol, maxiter):
63
64
65
        """
        Runs the algorithm for the specified input data
        """
66
67
        
        ##### PREPROCESSING
68
        normM = self.normM
69
        
70
71
72
73
74
75
76
77
78
79
80
81
        beta = self.beta_0
        iter = self.iter
        if u.ndim < 3:
            p = 1
            q = 1
        elif u.ndim == 3:
            p = u.shape[2]
            q = 1
        else:
            p = u.shape[2]
            q = u.shape[3]

82
83
84
        change = zeros(maxiter+1,dtype=u.dtype)
        change[0] = 999
        gap = change.copy()
85
86
87
88
89
90
        shadow_change = change.copy()
        if  hasattr(self, 'truth'):
            Relerrs = change.copy()

        tmp1 = 2*self.prox2.work(u) - u
        shadow = self.prox1.work(u)
91

92
        
93
94
        
        ##### LOOP
95
96
        while iter < maxiter and shadow_change[iter] >= tol:

97
            tmp = exp((-(iter+1)/self.beta_switch)**3);
98
            beta = (tmp*self.beta_0) + ((1-tmp)*self.beta_max)
99
            
100
            iter += 1;
101
            
102
103
104
            tmp3 = self.prox1.work(tmp1)
            tmp_u = ((beta*(2*tmp3-tmp1)) + ((1-beta)*tmp1) + u)/2
            tmp2 = self.prox2.work(tmp_u)
105
            
106
            tmp3 = self.prox1.work(tmp2)
107
            
108
109
            tmp_change = 0; tmp_gap = 0; tmp_shadow=0;

110
111
            #print(norm(tmp_u))

112
113
114
115
116
117
118
119
            if p==1 and q==1:
              tmp_change= (norm(u-tmp_u, 'fro')/normM)**2
              tmp_shadow = (norm(tmp3-shadow,'fro')/normM)**2
              tmp_gap = (norm(tmp3-tmp2,'fro')/normM)**2

              if hasattr(self, 'truth'):
                  if self.truth_dim[0] == 1:
                      z=tmp3[0,:]
120
                  elif self.truth_dim[1] == 1:
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
                      z=tmp3[:,0]
                  else:
                      z=tmp3;
                  Relerrs[iter] = norm(self.truth - exp(-1j*angle(trace(self.truth.T*z))) * z, 'fro')/self.norm_truth

            elif q==1:
              for j in range(self.product_space_dimension):
                  tmp_change= tmp_change+ (norm(u[:,:,j]-tmp_u[:,:,j], 'fro')/normM)**2
                  # compute (||P_SP_Mx-P_Mx||/normM)^2:
                  tmp_gap = tmp_gap+(norm(tmp3[:,:,j]-tmp2[:,:,j],'fro')/normM)**2
                  tmp_shadow = tmp_shadow+(norm(tmp3[:,:,j]-shadow[:,:,j],'fro')/normM)**2

              if hasattr(self, 'truth'):
                 z=tmp3[:,:,0]
                 Relerrs[iter] = norm(self.truth - exp(-1j*angle(trace(self.truth.T*z))) * z, 'fro')/self.norm_truth

137
            else:
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
              for j in range(self.product_space_dimension):
                  for k in range(Nz):
                      tmp_change= tmp_change+ (norm(u[:,:,k,j]-tmp_u[:,:,k,j], 'fro')/normM)**2
                      # compute (||P_Sx-P_Mx||/normM)^2:
                      tmp_gap = tmp_gap+(norm(tmp3[:,:,k,j]-tmp2[:,:,k,j],'fro')/(normM))**2
                      tmp_shadow = tmp_shadow+(norm(tmp3[:,:,k,j]-shadow[:,:,k,j],'fro')/(normM))**2
   
            change[iter] = sqrt(tmp_change)
            gap[iter] = sqrt(tmp_gap)
            shadow_change[iter] = sqrt(tmp_shadow) # this is the Euclidean norm of the gap to
            # 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.
            # graphics
153
            
154
155
            u = tmp_u
            tmp1 = (2*tmp2) - tmp_u
156
            
157
158
        
        ##### POSTPROCESSING
159
        u = tmp2;
160
161
162
163
164
165
166
167
        tmp = self.prox1.work(u);
        tmp2 = self.prox2.work(u);
        if self.Nx == 1:
            u1 = tmp[:,0];
            u2 = tmp2[:,0];
        elif self.Ny == 1:
            u1 = tmp[0,:];
            u2 = tmp2[0,:];
168
        elif self.Nz == 1 and tmp.ndim > 2:
169
170
            u1 = tmp[:,:,0]
            u2 = tmp2[:,:,0]
171
172
173
        else:
            u1 = tmp;
            u2 = tmp2;
174

175
176
177
178
        change = change[1:iter+1];
        gap = gap[1:iter+1];
        shadow_change = shadow_change[1:iter+1]

179
180
181
182
        output = {'u1': u1, 'u2': u2, 'iter': iter, 'change': change, 'gap': gap, 'shadow_change' : shadow_change}
        if hasattr(self, 'truth'):
            output['Relerrs']=Relerrs
        return output