Theo Steininger committed Apr 13, 2017 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 ``````# NIFTy # Copyright (C) 2017 Theo Steininger # # Author: Theo Steininger # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # You should have received a copy of the GNU General Public License # along with this program. If not, see . `````` Theo Steininger committed Sep 17, 2016 18 `````` `````` Theo Steininger committed Oct 14, 2016 19 20 21 ``````from __future__ import division import numpy as np `````` Theo Steininger committed Oct 23, 2016 22 ``````from keepers import Loggable `````` Theo Steininger committed Oct 14, 2016 23 `````` `````` Theo Steininger committed Oct 23, 2016 24 `````` `````` Theo Steininger committed Nov 02, 2016 25 ``````class ConjugateGradient(Loggable, object): `````` 26 27 `````` def __init__(self, convergence_tolerance=1E-4, convergence_level=3, iteration_limit=None, reset_count=None, `````` Theo Steininger committed Oct 14, 2016 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 `````` preconditioner=None, callback=None): """ Initializes the conjugate_gradient and sets the attributes (except for `x`). Parameters ---------- A : {operator, function} Operator `A` applicable to a field. b : field Resulting field of the operation `A(x)`. W : {operator, function}, *optional* Operator `W` that is a preconditioner on `A` and is applicable to a field (default: None). spam : function, *optional* Callback function which is given the current `x` and iteration counter each iteration (default: None). reset : integer, *optional* Number of iterations after which to restart; i.e., forget previous conjugated directions (default: sqrt(b.dim)). note : bool, *optional* Indicates whether notes are printed or not (default: False). """ self.convergence_tolerance = np.float(convergence_tolerance) self.convergence_level = np.float(convergence_level) `````` 54 55 56 57 58 59 60 61 `````` if iteration_limit is not None: iteration_limit = int(iteration_limit) self.iteration_limit = iteration_limit if reset_count is not None: reset_count = int(reset_count) self.reset_count = reset_count `````` Theo Steininger committed Oct 14, 2016 62 63 64 65 66 67 68 `````` if preconditioner is None: preconditioner = lambda z: z self.preconditioner = preconditioner self.callback = callback `````` 69 `````` def __call__(self, A, b, x0): `````` Theo Steininger committed Oct 14, 2016 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 `````` """ Runs the conjugate gradient minimization. Parameters ---------- x0 : field, *optional* Starting guess for the minimization. tol : scalar, *optional* Tolerance specifying convergence; measured by current relative residual (default: 1E-4). clevel : integer, *optional* Number of times the tolerance should be undershot before exiting (default: 1). limii : integer, *optional* Maximum number of iterations performed (default: 10 * b.dim). Returns ------- x : field Latest `x` of the minimization. convergence : integer Latest convergence level indicating whether the minimization has converged or not. """ r = b - A(x0) d = self.preconditioner(r) previous_gamma = r.dot(d) if previous_gamma == 0: `````` Theo Steininger committed Oct 18, 2016 99 100 `````` self.logger.info("The starting guess is already perfect solution " "for the inverse problem.") `````` Theo Steininger committed Oct 14, 2016 101 `````` return x0, self.convergence_level+1 `````` 102 `````` norm_b = np.sqrt(b.dot(b)) `````` Theo Steininger committed Oct 14, 2016 103 104 105 `````` x = x0 convergence = 0 iteration_number = 1 `````` Theo Steininger committed Oct 18, 2016 106 `````` self.logger.info("Starting conjugate gradient.") `````` Theo Steininger committed Oct 14, 2016 107 `````` `````` Theo Steininger committed Oct 18, 2016 108 `````` while True: `````` Theo Steininger committed Oct 14, 2016 109 110 111 112 113 114 115 `````` if self.callback is not None: self.callback(x, iteration_number) q = A(d) alpha = previous_gamma/d.dot(q) if not np.isfinite(alpha): `````` Theo Steininger committed Oct 18, 2016 116 `````` self.logger.error("Alpha became infinite! Stopping.") `````` Theo Steininger committed Oct 14, 2016 117 118 119 120 121 122 `````` return x0, 0 x += d * alpha reset = False if alpha.real < 0: `````` Theo Steininger committed Oct 18, 2016 123 `````` self.logger.warn("Positive definiteness of A violated!") `````` Theo Steininger committed Oct 14, 2016 124 `````` reset = True `````` 125 126 127 `````` if self.reset_count is not None: reset += (iteration_number % self.reset_count == 0) if reset: `````` Theo Steininger committed Oct 18, 2016 128 `````` self.logger.info("Resetting conjugate directions.") `````` Theo Steininger committed Oct 14, 2016 129 130 131 132 133 134 135 136 `````` r = b - A(x) else: r -= q * alpha s = self.preconditioner(r) gamma = r.dot(s) if gamma.real < 0: `````` Theo Steininger committed Oct 18, 2016 137 138 `````` self.logger.warn("Positive definitness of preconditioner " "violated!") `````` Theo Steininger committed Oct 14, 2016 139 140 141 `````` beta = max(0, gamma/previous_gamma) `````` 142 `````` delta = np.sqrt(gamma)/norm_b `````` Theo Steininger committed Oct 14, 2016 143 `````` `````` Theo Steininger committed Oct 18, 2016 144 145 146 147 148 149 `````` self.logger.debug("Iteration : %08u alpha = %3.1E " "beta = %3.1E delta = %3.1E" % (iteration_number, np.real(alpha), np.real(beta), np.real(delta))) `````` Theo Steininger committed Oct 14, 2016 150 151 152 `````` if gamma == 0: convergence = self.convergence_level+1 `````` Theo Steininger committed Oct 18, 2016 153 `````` self.logger.info("Reached infinite convergence.") `````` Theo Steininger committed Oct 14, 2016 154 155 156 `````` break elif abs(delta) < self.convergence_tolerance: convergence += 1 `````` Theo Steininger committed Oct 18, 2016 157 158 `````` self.logger.info("Updated convergence level to: %u" % convergence) `````` Theo Steininger committed Oct 14, 2016 159 `````` if convergence == self.convergence_level: `````` Theo Steininger committed Oct 18, 2016 160 `````` self.logger.info("Reached target convergence level.") `````` Theo Steininger committed Oct 14, 2016 161 162 163 164 `````` break else: convergence = max(0, convergence-1) `````` 165 166 `````` if self.iteration_limit is not None: if iteration_number == self.iteration_limit: `````` Theo Steininger committed Oct 18, 2016 167 `````` self.logger.warn("Reached iteration limit. Stopping.") `````` 168 169 170 `````` break d = s + d * beta `````` Theo Steininger committed Oct 14, 2016 171 172 173 174 175 `````` iteration_number += 1 previous_gamma = gamma return x, convergence``````