Package specific_analyses :: Package relax_fit :: Module optimisation
[hide private]
[frames] | no frames]

Source Code for Module specific_analyses.relax_fit.optimisation

  1  ############################################################################### 
  2  #                                                                             # 
  3  # Copyright (C) 2004-2014 Edward d'Auvergne                                   # 
  4  #                                                                             # 
  5  # This file is part of the program relax (http://www.nmr-relax.com).          # 
  6  #                                                                             # 
  7  # This program is free software: you can redistribute it and/or modify        # 
  8  # it under the terms of the GNU General Public License as published by        # 
  9  # the Free Software Foundation, either version 3 of the License, or           # 
 10  # (at your option) any later version.                                         # 
 11  #                                                                             # 
 12  # This program is distributed in the hope that it will be useful,             # 
 13  # but WITHOUT ANY WARRANTY; without even the implied warranty of              # 
 14  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               # 
 15  # GNU General Public License for more details.                                # 
 16  #                                                                             # 
 17  # You should have received a copy of the GNU General Public License           # 
 18  # along with this program.  If not, see <http://www.gnu.org/licenses/>.       # 
 19  #                                                                             # 
 20  ############################################################################### 
 21   
 22  # Module docstring. 
 23  """The R1 and R2 exponential relaxation curve fitting optimisation functions.""" 
 24   
 25  # Python module imports. 
 26  from numpy import average 
 27  from re import search 
 28   
 29  # relax module imports. 
 30  from dep_check import C_module_exp_fn 
 31  from lib.errors import RelaxError, RelaxLenError 
 32  from specific_analyses.relax_fit.parameters import assemble_param_vector, assemble_scaling_matrix 
 33   
 34  # C modules. 
 35  if C_module_exp_fn: 
 36      from target_functions.relax_fit import setup, func, back_calc_I 
 37   
 38   
39 -def back_calc(spin=None, relax_time_id=None):
40 """Back-calculation of peak intensity for the given relaxation time. 41 42 @keyword spin: The spin container. 43 @type spin: SpinContainer instance 44 @keyword relax_time_id: The ID string for the desired relaxation time. 45 @type relax_time_id: str 46 @return: The peak intensity for the desired relaxation time. 47 @rtype: float 48 """ 49 50 # Create the initial parameter vector. 51 param_vector = assemble_param_vector(spin=spin) 52 53 # Create a scaling matrix. 54 scaling_matrix = assemble_scaling_matrix(spin=spin, scaling=False) 55 56 # The keys. 57 keys = list(spin.peak_intensity.keys()) 58 59 # The peak intensities and times. 60 values = [] 61 errors = [] 62 times = [] 63 for key in keys: 64 values.append(spin.peak_intensity[key]) 65 errors.append(spin.peak_intensity_err[key]) 66 times.append(cdp.relax_times[key]) 67 68 # The scaling matrix in a diagonalised list form. 69 scaling_list = [] 70 for i in range(len(scaling_matrix)): 71 scaling_list.append(scaling_matrix[i, i]) 72 73 # Initialise the relaxation fit functions. 74 setup(num_params=len(spin.params), num_times=len(cdp.relax_times), values=values, sd=errors, relax_times=times, scaling_matrix=scaling_list) 75 76 # Make a single function call. This will cause back calculation and the data will be stored in the C module. 77 func_wrapper(param_vector) 78 79 # Get the data back. 80 results = back_calc_I() 81 82 # Return the correct peak height. 83 return results[keys.index(relax_time_id)]
84 85
86 -def func_wrapper(params):
87 """Wrapper function for the C module, for converting numpy arrays. 88 89 @param params: The parameter array from the minimisation code. 90 @type params: numpy array 91 @return: The function value generated by the C module. 92 @rtype: float 93 """ 94 95 # Call the C code. 96 chi2 = func(params.tolist()) 97 98 # Return the chi2 value. 99 return chi2
100 101
102 -def dfunc_wrapper(params):
103 """Wrapper function for the C module, for converting numpy arrays. 104 105 The currently does nothing. 106 """
107 108
109 -def d2func_wrapper(params):
110 """Wrapper function for the C module, for converting numpy arrays. 111 112 The currently does nothing. 113 """
114 115
116 -def grid_search_setup(spin=None, param_vector=None, lower=None, upper=None, inc=None, scaling_matrix=None):
117 """The grid search setup function. 118 119 @keyword spin: The spin data container. 120 @type spin: SpinContainer instance 121 @keyword param_vector: The parameter vector. 122 @type param_vector: numpy array 123 @keyword lower: The lower bounds of the grid search which must be equal to the 124 number of parameters in the model. This optional argument is 125 only used when doing a grid search. 126 @type lower: array of numbers 127 @keyword upper: The upper bounds of the grid search which must be equal to the 128 number of parameters in the model. This optional argument is 129 only used when doing a grid search. 130 @type upper: array of numbers 131 @keyword inc: The increments for each dimension of the space for the grid 132 search. The number of elements in the array must equal to the 133 number of parameters in the model. This argument is only used 134 when doing a grid search. 135 @type inc: array of int 136 @keyword scaling_matrix: The scaling matrix. 137 @type scaling_matrix: numpy diagonal matrix 138 @return: A tuple of the grid size and the minimisation options. For the 139 minimisation options, the first dimension corresponds to the 140 model parameter. The second dimension is a list of the number 141 of increments, the lower bound, and upper bound. 142 @rtype: (int, list of lists [int, float, float]) 143 """ 144 145 # The length of the parameter array. 146 n = len(param_vector) 147 148 # Make sure that the length of the parameter array is > 0. 149 if n == 0: 150 raise RelaxError("Cannot run a grid search on a model with zero parameters.") 151 152 # Lower bounds. 153 if lower != None and len(lower) != n: 154 raise RelaxLenError('lower bounds', n) 155 156 # Upper bounds. 157 if upper != None and len(upper) != n: 158 raise RelaxLenError('upper bounds', n) 159 160 # Increments. 161 if isinstance(inc, list) and len(inc) != n: 162 raise RelaxLenError('increment', n) 163 elif isinstance(inc, int): 164 inc = [inc]*n 165 166 # Set up the default bounds. 167 if not lower: 168 # Init. 169 lower = [] 170 upper = [] 171 172 # Loop over the parameters. 173 for i in range(n): 174 # Relaxation rate (from 0 to 20 s^-1). 175 if spin.params[i] == 'rx': 176 lower.append(0.0) 177 upper.append(20.0) 178 179 # Intensity 180 elif search('^i', spin.params[i]): 181 # Find the ID of the first time point. 182 min_time = min(cdp.relax_times.values()) 183 for key in list(cdp.relax_times.keys()): 184 if cdp.relax_times[key] == min_time: 185 id = key 186 break 187 188 # Defaults. 189 lower.append(0.0) 190 upper.append(average(spin.peak_intensity[id])) 191 192 # Diagonal scaling of minimisation options. 193 lower_new = [] 194 upper_new = [] 195 for i in range(n): 196 lower_new.append(lower[i] / scaling_matrix[i, i]) 197 upper_new.append(upper[i] / scaling_matrix[i, i]) 198 199 # Return the minimisation options. 200 return inc, lower_new, upper_new
201