model_free

11 - def __init__(self, mf):

12 "Class used to create and process input and output for the program Modelfree 4." 13 14 self.mf = mf 15 16 self.update_data() 17 self.ask_stage() 18 if match('1', self.mf.data.stage): 19 self.minimisation() 20 elif match('^2', self.mf.data.stage): 21 self.model_selection()

22 23

24 - def minimisation(self):

25 self.extract_relax_data() 26 self.mf.data.calc_frq() 27 self.mf.data.calc_constants() 28 29 self.mf.file_ops.mkdir(self.mf.data.model) 30 self.mf.results = open(self.mf.data.model + '/results', 'w') 31 32 print "\n\n[ Model-free fitting ]\n" 33 diff_type = 'iso' 34 diff_params = (float(10.0 * 1e-9)) 35 mf_model = self.mf.data.model 36 37 # Flags. 38 ######## 39 40 # Debugging flag 41 self.mf.min_debug = 1 42 43 # Flag for the diagonal scaling of the model-free parameters. 44 self.scaling_flag = 1 45 46 # Minimisation options. 47 chi2_tol = 1e-15 48 max_iterations = 5000 49 50 # Set the function, gradient, and hessian functions. 51 if self.scaling_flag: 52 self.functions = mf_trans_functions(self.mf) 53 func = self.functions.chi2 54 dfunc = self.functions.dchi2 55 d2func = self.functions.d2chi2 56 else: 57 self.functions = mf_functions(self.mf) 58 func = self.functions.chi2 59 dfunc = self.functions.dchi2 60 d2func = self.functions.d2chi2 61 62 # The value of the constant. 63 self.c = 1e10 64 65 # Initialise the grid options. 66 if match('[Gg]rid', self.mf.usr_param.init_params[0]): 67 init_params = [self.mf.usr_param.init_params[0], self.init_grid_ops()] 68 69 # Initialise the fixed model-free parameter options. 70 elif match('[Ff]ixed', self.mf.usr_param.init_params[0]): 71 init_params = [self.mf.usr_param.init_params[0], self.init_fixed_params()] 72 73 # Initialise the Levenberg-Marquardt minimisation options. 74 if match('[Ll][Mm]$', self.mf.usr_param.minimiser[0]) or match('[Ll]evenburg-[Mm]arquardt$', self.mf.usr_param.minimiser[0]): 75 if self.scaling_flag: 76 self.mf.usr_param.minimiser.append(mf_trans_functions.lm_dri) 77 else: 78 self.mf.usr_param.minimiser.append(mf_functions.lm_dri) 79 self.mf.usr_param.minimiser.append([]) 80 81 # Print a header into the results file. 82 self.print_header() 83 84 # Loop over all data sets. 85 for self.res in range(len(self.mf.data.relax_data[0])): 86 if self.mf.min_debug >= 1: 87 print "\n\n<<< Fitting to residue: " + self.mf.data.relax_data[0][self.res][0] + " " + self.mf.data.relax_data[0][self.res][1] + " >>>" 88 else: 89 print "Residue: " + self.mf.data.relax_data[0][self.res][0] + " " + self.mf.data.relax_data[0][self.res][1] 90 91 # Initialise the iteration counter and function, gradient, and hessian call counters. 92 self.iter_count = 0 93 self.f_count = 0 94 self.g_count = 0 95 self.h_count = 0 96 97 # Initialise the relaxation data and error data structures. 98 relax_data = [] 99 errors = [] 100 for data in range(len(self.mf.data.relax_data)): 101 relax_data.append(self.mf.data.relax_data[data][self.res][2]) 102 errors.append(self.mf.data.relax_data[data][self.res][3]) 103 function_ops = (diff_type, diff_params, mf_model, relax_data, errors) 104 if match('^[Ll][Mm]$', self.mf.usr_param.minimiser[0]) or match('^[Ll]evenburg-[Mm]arquardt$', self.mf.usr_param.minimiser[0]): 105 self.mf.usr_param.minimiser[2] = errors 106 107 # Initialisation of model-free parameter values. 108 results = self.mf.minimise(func, dfunc=dfunc, d2func=d2func, args=function_ops, x0=None, minimiser=init_params, full_output=1, print_flag=self.mf.min_debug) 109 self.params, self.chi2, iter, fc, gc, hc, self.warning = results 110 111 # Minimisation. 112 results = self.mf.minimise(func, dfunc=dfunc, d2func=d2func, args=function_ops, x0=self.params, minimiser=self.mf.usr_param.minimiser, func_tol=chi2_tol, maxiter=max_iterations, full_output=1, print_flag=self.mf.min_debug) 113 self.params, self.chi2, iter, fc, gc, hc, self.warning = results 114 self.iter_count = self.iter_count + iter 115 self.f_count = self.f_count + fc 116 self.g_count = self.g_count + gc 117 self.h_count = self.h_count + hc 118 119 if self.mf.min_debug: 120 print "\n\n<<< Finished minimiser >>>" 121 122 # Write the results to file. 123 self.print_results() 124 125 print "\n[ Done ]\n\n" 126 self.mf.results.close()

127 128

129 - def model_selection(self):

130 "Model selection stage."

131

132 - def ask_stage(self):

133 "User input of stage number." 134 135 print "\n[ Select the stage for model-free analysis ]\n" 136 print "The stages are:" 137 print " Stage 1 (1): Minimisation." 138 print " Stage 2 (2): Model selection." 139 140 while 1: 141 input = raw_input('> ') 142 valid_stages = ['1', '2'] 143 if input in valid_stages: 144 self.mf.data.stage = input 145 break 146 else: 147 print "Invalid stage number. Choose either 1 or 2." 148 if match('1', self.mf.data.stage): 149 while 1: 150 print "Please select which model-free model to fit to." 151 print " (1): Model-free model 1, {S2}." 152 print " (2): Model-free model 2, {S2, te}." 153 print " (3): Model-free model 3, {S2, Rex}." 154 print " (4): Model-free model 4, {S2, te, Rex}." 155 print " (5): Model-free model 5, {S2f, S2s, ts}." 156 input = raw_input('> ') 157 valid_stages = ['1', '2', '3', '4', '5'] 158 if input in valid_stages: 159 self.mf.data.model = 'm' + input 160 break 161 else: 162 print "Invalid model-free model." 163 elif match('2', self.mf.data.stage): 164 while 1: 165 print "Stage 2 has the following two options for the final run:" 166 print " (a): No optimization of the diffusion tensor." 167 print " (b): Optimization of the diffusion tensor." 168 input = raw_input('> ') 169 valid_stages = ['a', 'b'] 170 if input in valid_stages: 171 self.mf.data.stage = self.mf.data.stage + input 172 break 173 else: 174 print "Invalid option, choose either a or b." 175 176 print "The stage chosen is " + self.mf.data.stage + "\n" 177 if match('1', self.mf.data.stage): 178 print "The model-free model chosen is " + self.mf.data.model 179 print "\n"

180 181

182 - def init_fixed_params(self):

183 """Fix the inital model-free parameter values. 184 185 The function needs to be modified so the user can specify the fixed values. 186 """ 187 188 if match('m1', self.mf.data.model): 189 params = array([0.5], Float64) 190 elif match('m2', self.mf.data.model): 191 if self.scaling_flag: 192 params = array([0.5, 100.0 * 1e-12 * self.c], Float64) 193 else: 194 params = array([0.5, 100.0 * 1e-12], Float64) 195 elif match('m3', self.mf.data.model): 196 params = array([0.5, 0.0], Float64) 197 elif match('m4', self.mf.data.model): 198 if self.scaling_flag: 199 params = array([0.5, 0.0, 100.0 * 1e-12 * self.c], Float64) 200 else: 201 params = array([0.5, 0.0, 100.0 * 1e-12], Float64) 202 elif match('m5', self.mf.data.model): 203 if self.scaling_flag: 204 params = array([0.5, 0.5, 100.0 * 1e-12 * self.c], Float64) 205 else: 206 params = array([0.5, 0.5, 100.0 * 1e-12], Float64) 207 208 return params

209 210

211 - def init_grid_ops(self):

212 """Generate the data structure of grid options for the grid search. 213 214 The function needs to be modified so the user can specify the number of iterations for each parameter. 215 """ 216 217 inc = self.mf.usr_param.init_params[1] 218 grid_ops = [] 219 if match('m1', self.mf.data.model): 220 grid_ops.append([inc, 0.0, 1.0]) 221 elif match('m2', self.mf.data.model): 222 grid_ops.append([inc, 0.0, 1.0]) 223 if self.scaling_flag: 224 grid_ops.append([inc, 0.0, 10000.0 * 1e-12 * self.c]) 225 else: 226 grid_ops.append([inc, 0.0, 10000.0 * 1e-12]) 227 elif match('m3', self.mf.data.model): 228 grid_ops.append([inc, 0.0, 1.0]) 229 grid_ops.append([inc, 0.0, 30.0 / (1e-8 * self.mf.data.frq[0])**2]) 230 elif match('m4', self.mf.data.model): 231 grid_ops.append([inc, 0.0, 1.0]) 232 if self.scaling_flag: 233 grid_ops.append([inc, 0.0, 10000.0 * 1e-12 * self.c]) 234 else: 235 grid_ops.append([inc, 0.0, 10000.0 * 1e-12]) 236 grid_ops.append([inc, 0.0, 30.0 / (1e-8 * self.mf.data.frq[0])**2]) 237 elif match('m5', self.mf.data.model): 238 grid_ops.append([inc, 0.0, 1.0]) 239 grid_ops.append([inc, 0.0, 1.0]) 240 if self.scaling_flag: 241 grid_ops.append([inc, 0.0, 10000.0 * 1e-12 * self.c]) 242 else: 243 grid_ops.append([inc, 0.0, 10000.0 * 1e-12]) 244 245 return grid_ops

246 247

248 - def print_header(self):

249 self.mf.results.write("%-5s" % "Num") 250 self.mf.results.write("%-6s" % "Name") 251 if match('m1', self.mf.data.model): 252 self.mf.results.write("%-30s" % "S2") 253 elif match('m2', self.mf.data.model): 254 self.mf.results.write("%-30s" % "S2") 255 self.mf.results.write("%-30s" % "te (ps)") 256 elif match('m3', self.mf.data.model): 257 self.mf.results.write("%-30s" % "S2") 258 self.mf.results.write("%-30s" % ("Rex (" + self.mf.data.frq_label[0] + " MHz)")) 259 elif match('m4', self.mf.data.model): 260 self.mf.results.write("%-30s" % "S2") 261 self.mf.results.write("%-30s" % "te (ps)") 262 self.mf.results.write("%-30s" % ("Rex (" + self.mf.data.frq_label[0] + " MHz)")) 263 elif match('m5', self.mf.data.model): 264 self.mf.results.write("%-30s" % "S2f") 265 self.mf.results.write("%-30s" % "S2s") 266 self.mf.results.write("%-30s" % "ts (ps)") 267 self.mf.results.write("%-30s" % "Chi-squared statistic") 268 self.mf.results.write("%-15s" % "Iterations") 269 self.mf.results.write("%-15s" % "Func calls") 270 self.mf.results.write("%-15s" % "Grad calls") 271 self.mf.results.write("%-15s" % "Hessian calls") 272 self.mf.results.write("%-30s\n" % "Warning")

273 274

275 - def print_results(self):

276 self.mf.results.write("%-5s" % self.mf.data.relax_data[0][self.res][0]) 277 self.mf.results.write("%-6s" % self.mf.data.relax_data[0][self.res][1]) 278 if match('m1', self.mf.data.model): 279 s2 = self.params[0] 280 print "S2: " + `s2` + " Chi2: " + `self.chi2` 281 self.mf.results.write("%-30s" % `s2`) 282 elif match('m2', self.mf.data.model): 283 s2 = self.params[0] 284 if self.scaling_flag: 285 te = self.params[1] * 1e12 / self.c 286 else: 287 te = self.params[1] * 1e12 288 print "S2: " + `s2` + " te: " + `te` + " Chi2: " + `self.chi2` 289 self.mf.results.write("%-30s" % `s2`) 290 self.mf.results.write("%-30s" % `te`) 291 elif match('m3', self.mf.data.model): 292 s2 = self.params[0] 293 rex = self.params[1] * (1e-8 * self.mf.data.frq[0])**2 294 print "S2: " + `s2` + " Rex: " + `rex` + " Chi2: " + `self.chi2` 295 self.mf.results.write("%-30s" % `s2`) 296 self.mf.results.write("%-30s" % `rex`) 297 elif match('m4', self.mf.data.model): 298 s2 = self.params[0] 299 if self.scaling_flag: 300 te = self.params[1] * 1e12 / self.c 301 else: 302 te = self.params[1] * 1e12 303 rex = self.params[2] * (1e-8 * self.mf.data.frq[0])**2 304 print "S2: " + `s2` + " te: " + `te` + " Rex: " + `rex` + " Chi2: " + `self.chi2` 305 self.mf.results.write("%-30s" % `s2`) 306 self.mf.results.write("%-30s" % `te`) 307 self.mf.results.write("%-30s" % `rex`) 308 elif match('m5', self.mf.data.model): 309 s2f = self.params[0] 310 s2s = self.params[1] 311 if self.scaling_flag: 312 ts = self.params[2] * 1e12 / self.c 313 else: 314 ts = self.params[2] * 1e12 315 print "S2f: " + `s2f` + " S2s: " + `s2s` + " ts: " + `ts` + " Chi2: " + `self.chi2` 316 self.mf.results.write("%-30s" % `s2f`) 317 self.mf.results.write("%-30s" % `s2s`) 318 self.mf.results.write("%-30s" % `ts`) 319 self.mf.results.write("%-30s" % `self.chi2`) 320 self.mf.results.write("%-15s" % `self.iter_count`) 321 self.mf.results.write("%-15s" % `self.f_count`) 322 self.mf.results.write("%-15s" % `self.g_count`) 323 self.mf.results.write("%-15s" % `self.h_count`) 324 if self.warning: 325 self.mf.results.write("%-30s\n" % `self.warning`) 326 else: 327 self.mf.results.write("\n")

Source Code for Module model_free