Source Code for Module functions.d2chi2

  1  import sys 
  2  from Numeric import Float64, zeros 
  3   
  4 -class d2chi2: 
  5 -        def __init__(self): 
  6                  "Function to create the chi-squared hessian." 
  7   
  8   
  9 -        def d2chi2(self, params, diff_type, diff_params, model, relax_data, errors, print_flag=0): 
 10                  """Function to create the chi-squared hessian. 
 11   
 12                  Function arguments 
 13                  ~~~~~~~~~~~~~~~~~~ 
 14   
 15                  1:  params - a list containing the parameter values specific for the given model. 
 16                  The order of parameters must be as follows: 
 17                          m1 - {S2} 
 18                          m2 - {S2, te} 
 19                          m3 - {S2, Rex} 
 20                          m4 - {S2, te, Rex} 
 21                          m5 - {S2f, S2s, ts} 
 22                  2:  diff_type - string.  The diffusion tensor, ie 'iso', 'axial', 'aniso' 
 23                  3:  diff_params - array.  An array with the diffusion parameters 
 24                  4:  model - string.  The model 
 25                  5:  relax_data - array.  An array containing the experimental relaxation values. 
 26                  6:  errors - array.  An array containing the experimental errors. 
 27   
 28   
 29                  The chi-sqared hessian 
 30                  ~~~~~~~~~~~~~~~~~~~~~~ 
 31   
 32                  Data structure:  self.data.d2chi2 
 33                  Dimension:  2D, (parameters, parameters) 
 34                  Type:  Numeric array, Float64 
 35                  Dependencies:  self.data.ri, self.data.dri, self.data.d2ri 
 36                  Required by:  None 
 37   
 38   
 39                  Formula 
 40                  ~~~~~~~ 
 41                                                _n_ 
 42                               d2chi2           \       1      /  dRi()     dRi()                         d2Ri()     \  
 43                          ---------------  =  2  >  ---------- | ------- . -------  -  (Ri - Ri()) . --------------- | 
 44                          dthetaj.dthetak       /__ sigma_i**2 \ dthetaj   dthetak                   dthetaj.dthetak /  
 45                                                i=1 
 46                  """ 
 47   
 48                  # debug. 
 49                  if print_flag == 2: 
 50                          print "\n< d2chi2 >" 
 51                          print "Params: " + `params` 
 52   
 53                  self.data.params = params 
 54                  self.data.diff_type = diff_type 
 55                  self.data.diff_params = diff_params 
 56                  self.data.model = model 
 57                  self.data.relax_data = relax_data 
 58                  self.data.errors = errors 
 59   
 60                  # Test to see if relaxation values and spectral density values have previously been calculated for the current parameter values, 
 61                  # If not, calculate the relaxation values which will then calculate the spectral density values, and store the current data and parameter 
 62                  # values in self.data.relax_test and self.data.gradient_test. 
 63                  test = [ self.data.relax_data[:], self.data.errors[:], self.data.params.tolist(), self.data.model ] 
 64                  if test != self.data.relax_test: 
 65                          self.init_param_types() 
 66                          self.Ri() 
 67                          self.data.relax_test = test[:] 
 68                  if test != self.data.gradient_test: 
 69                          self.init_param_types() 
 70                          self.dRi() 
 71                          self.data.gradient_test = test[:] 
 72   
 73                  # Calculate the relaxation hessians. 
 74                  self.d2Ri() 
 75   
 76                  # Initialise the chi-squared hessian. 
 77                  self.data.d2chi2 = zeros((len(self.data.params), len(self.data.params)), Float64) 
 78   
 79                  # Calculate the chi-squared hessian. 
 80                  for i in range(len(self.data.relax_data)): 
 81                          if self.data.errors[i] != 0.0: 
 82                                  # Parameter independent terms. 
 83                                  a = 2.0 / (self.data.errors[i]**2) 
 84                                  b = self.data.relax_data[i] - self.data.ri[i] 
 85   
 86                                  # Loop over the parameters. 
 87                                  for param_j in range(len(self.data.params)): 
 88                                          # Loop over the parameters from 1 to param_k. 
 89                                          for param_k in range(param_j + 1): 
 90                                                  self.data.d2chi2[param_j, param_k] = self.data.d2chi2[param_j, param_k] + a * (self.data.dri[param_j, i] * self.data.dri[param_k, i] - b * self.data.d2ri[param_j, param_k, i]) 
 91                                                  if param_j != param_k: 
 92                                                          self.data.d2chi2[param_k, param_j] = self.data.d2chi2[param_j, param_k] 
 93                          else: 
 94                                  # Loop over the parameters. 
 95                                  for param_k in range(len(self.data.params)): 
 96                                          # Loop over the parameters from 1 to param_k. 
 97                                          for param_j in range(param_k + 1): 
 98                                                  self.data.d2chi2[param_j, param_j] = 1e99 
 99                                                  if param_j != param_k: 
100                                                          self.data.d2chi2[param_k, param_j] = 1e99 
101                                  break 
102   
103                  # debug. 
104                  if print_flag == 2: 
105                          print "J(w):   " + `self.data.jw` 
106                          print "dJ(w):  " + `self.data.djw` 
107                          print "d2J(w): " + `self.data.d2jw` 
108                          print "Ri:     " + `self.data.ri` 
109                          print "dRi:    " + `self.data.dri` 
110                          print "d2Ri:   " + `self.data.d2ri` 
111                          print "d2chi2: " + `self.data.d2chi2` 
112                          print "\n" 
113   
114                  return self.data.d2chi2 
115