1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
|
# -*- coding: utf-8 -*-
# Copyright 2019 Science Technology Facilities Council
# Copyright 2019 University of Manchester
#
# This work is part of the Core Imaging Library developed by Science Technology
# Facilities Council and University of Manchester
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0.txt
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
from ccpi.optimisation.algorithms import Algorithm
import numpy
class CGLS(Algorithm):
r'''Conjugate Gradient Least Squares algorithm
Problem:
.. math::
\min || A x - b ||^2_2
|
Parameters :
:parameter operator : Linear operator for the inverse problem
:parameter x_init : Initial guess ( Default x_init = 0)
:parameter data : Acquired data to reconstruct
:parameter tolerance: Tolerance/ Stopping Criterion to end CGLS algorithm
Reference:
https://web.stanford.edu/group/SOL/software/cgls/
'''
def __init__(self, **kwargs):
super(CGLS, self).__init__()
x_init = kwargs.get('x_init', None)
operator = kwargs.get('operator', None)
data = kwargs.get('data', None)
tolerance = kwargs.get('tolerance', 1e-6)
if x_init is not None and operator is not None and data is not None:
print(self.__class__.__name__ , "set_up called from creator")
self.set_up(x_init=x_init, operator=operator, data=data, tolerance=tolerance)
def set_up(self, x_init, operator, data, tolerance=1e-6):
self.x = x_init * 0.
self.operator = operator
self.tolerance = tolerance
self.r = data - self.operator.direct(self.x)
self.s = self.operator.adjoint(self.r)
self.p = self.s
self.norms0 = self.s.norm()
self.norms = self.s.norm()
self.gamma = self.norms0**2
self.normx = self.x.norm()
self.xmax = self.normx
self.loss.append(self.r.squared_norm())
self.configured = True
def update(self):
self.q = self.operator.direct(self.p)
delta = self.q.squared_norm()
alpha = self.gamma/delta
self.x += alpha * self.p
self.r -= alpha * self.q
self.s = self.operator.adjoint(self.r)
self.norms = self.s.norm()
self.gamma1 = self.gamma
self.gamma = self.norms**2
self.beta = self.gamma/self.gamma1
self.p = self.s + self.beta * self.p
self.normx = self.x.norm()
self.xmax = numpy.maximum(self.xmax, self.normx)
def update_objective(self):
a = self.r.squared_norm()
if a is numpy.nan:
raise StopIteration()
self.loss.append(a)
def should_stop(self):
return self.flag() or self.max_iteration_stop_cryterion()
def flag(self):
flag = (self.norms <= self.norms0 * self.tolerance) or (self.normx * self.tolerance >= 1)
if flag:
self.update_objective()
if self.iteration > self._iteration[-1]:
print (self.verbose_output())
print('Tolerance is reached: {}'.format(self.tolerance))
return flag
|
