Source Code for Module minimise.levenberg_marquardt

  1  from LinearAlgebra import solve_linear_equations 
  2  from Numeric import Float64, copy, zeros 
  3   
  4  from generic_minimise import generic_minimise 
  5   
  6   
  7 -class levenberg_marquardt(generic_minimise): 
  8 -        def __init__(self, chi2_func=None, dchi2_func=None, dfunc=None, errors=None, args=(), x0=None, func_tol=1e-5, maxiter=1000, full_output=0, print_flag=0): 
  9                  """Class for Levenberg-Marquardt minimisation specific functions. 
 10   
 11                  Function options 
 12                  ~~~~~~~~~~~~~~~~ 
 13   
 14                  chi2_func:  User supplied chi-squared function which is run with the function parameters and args as options. 
 15                  dchi2_func:  User supplied chi-squared gradient function which is run with the function parameters and args as options. 
 16                  dfunc:  User supplied function which should return a vector of partial derivatives of the function which back calculates 
 17                          values for the chi-squared function. 
 18                  params:  The initial function parameter values. 
 19                  errors:  The experimental errors. 
 20                  args:  A tuple containing the arguments to send to chi2_func and dchi2_func. 
 21                  maxiter:  The maximum number of iterations. 
 22                  full_output:  A flag specifying what should be returned. 
 23   
 24   
 25                  Output 
 26                  ~~~~~~ 
 27   
 28                  If full_output = 0, the parameter values and chi-squared value are returned as a tuple. 
 29                  If full_output = 1, the parameter values, chi-squared value, number of iterations, and the warning flag are returned as a tuple. 
 30   
 31                  """ 
 32   
 33                  self.chi2_func = chi2_func 
 34                  self.dchi2_func = dchi2_func 
 35                  self.dfunc = dfunc 
 36                  self.errors = errors 
 37                  self.args = args 
 38                  self.xk = x0 
 39                  self.func_tol = func_tol 
 40                  self.maxiter = maxiter 
 41                  self.full_output = full_output 
 42                  self.print_flag = print_flag 
 43   
 44                  # Initialise the function, gradient, and hessian evaluation counters. 
 45                  self.f_count = 0 
 46                  self.g_count = 0 
 47                  self.h_count = 0 
 48   
 49                  # Initialise the warning string. 
 50                  self.warning = None 
 51   
 52                  # The initial chi-squared value, chi-squared gradient vector, and derivative function matrix. 
 53                  self.fk, self.f_count = apply(self.chi2_func, (self.xk,)+self.args), self.f_count + 1 
 54                  self.dfk, self.g_count = -0.5 * apply(self.dchi2_func, (self.xk,)+self.args), self.g_count + 1 
 55                  self.df = self.dfunc() 
 56   
 57                  # Initial value of lambda (the Levenberg-Marquardt fudge factor). 
 58                  self.l = 1.0 
 59                  self.n = len(self.xk) 
 60   
 61                  # Minimisation. 
 62                  self.minimise = self.generic_minimise 
 63   
 64   
 65 -        def backup_current_data(self): 
 66                  "Function to backup the current data into fk_last." 
 67   
 68                  self.fk_last = self.fk 
 69   
 70   
 71 -        def create_lm_matrix(self): 
 72                  """Function to create the Levenberg-Marquardt matrix. 
 73   
 74                  The matrix is: 
 75   
 76                                         _n_ 
 77                                         \   /     1        d y(xi)   d y(xi)                \  
 78                          LM_matrix_jk =  >  | ---------- . ------- . ------- . (1 + lambda) | 
 79                                         /__ \ sigma_i**2     dj        dk                   / 
 80                                         i=1 
 81   
 82                          where j == k is one of the function parameters. 
 83   
 84                                         _n_ 
 85                                         \   /     1        d y(xi)   d y(xi) \  
 86                          LM_matrix_jk =  >  | ---------- . ------- . ------- | 
 87                                         /__ \ sigma_i**2     dj        dk    / 
 88                                         i=1 
 89   
 90                          where j != k are function parameters. 
 91                  """ 
 92   
 93                  # Create the Levenberg-Marquardt matrix with elements equal to zero. 
 94                  self.lm_matrix = zeros((self.n, self.n), Float64) 
 95   
 96                  # Loop over the error points from i=1 to n. 
 97                  for i in range(len(self.errors)): 
 98                          # Calculate the inverse of the variance to minimise calculations. 
 99                          i_variance = 1.0 / self.errors[i]**2 
100   
101                          # Loop over all function parameters. 
102                          for param_j in range(self.n): 
103                                  # Loop over the function parameters from the first to 'param_j' to create the Levenberg-Marquardt matrix. 
104                                  for param_k in range(param_j + 1): 
105                                          if param_j == param_k: 
106                                                  matrix_jk = i_variance * self.df[param_j, i] * self.df[param_k, i] * (1.0 + self.l) 
107                                          else: 
108                                                  matrix_jk = i_variance * self.df[param_j, i] * self.df[param_k, i] 
109   
110                                          self.lm_matrix[param_j, param_k] = self.lm_matrix[param_j, param_k] + matrix_jk 
111                                          self.lm_matrix[param_k, param_j] = self.lm_matrix[param_k, param_j] + matrix_jk 
112   
113   
114 -        def new_param_func(self): 
115                  "Find the new parameter vector self.xk_new" 
116   
117                  # Create the Levenberg-Marquardt matrix. 
118                  self.create_lm_matrix() 
119   
120                  # Solve the Levenberg-Marquardt equation to get the vector of function parameter changes. 
121                  self.param_change = solve_linear_equations(self.lm_matrix, self.dfk) 
122   
123                  # Add the current parameter vector to the parameter change vector to find the new parameter vector. 
124                  self.xk_new = zeros(self.n, Float64) 
125                  self.xk_new = self.xk + self.param_change 
126   
127   
128 -        def tests(self): 
129                  "Levenberg-Marquardt tests." 
130   
131                  if self.fk > self.fk_last: 
132                          if self.l <= 1e99: 
133                                  self.l = self.l * 10.0 
134                  else: 
135                          # Finish minimising when the chi-squared difference is insignificant. 
136                          if abs(self.fk_last - self.fk) < self.func_tol: 
137                                  return 1 
138                          if self.l >= 1e-99: 
139                                  self.l = self.l / 10.0 
140                  return 0 
141   
142   
143 -        def update_data(self): 
144                  "Function to update the chi-squared value, chi-squared gradient vector, and derivative function matrix." 
145   
146                  self.xk = copy.deepcopy(self.xk_new) 
147                  self.fk, self.f_count = apply(self.chi2_func, (self.xk,)+self.args), self.f_count + 1 
148                  self.dfk, self.g_count = -0.5 * apply(self.dchi2_func, (self.xk,)+self.args), self.g_count + 1 
149                  self.df = self.dfunc() 
150