generic_fns.monte

30 - def __init__(self, relax):

31 """Class containing functions for Monte Carlo simulations.""" 32 33 self.relax = relax

34 35

36 - def create_data(self, run=None, method=None):

37 """Function for creating simulation data. 38 39 It is assumed that all data types are residue specific. 40 """ 41 42 # Arguments. 43 self.run = run 44 45 # Test if the run exists. 46 if not self.run in self.relax.data.run_names: 47 raise RelaxNoRunError, self.run 48 49 # Test if simulations have been set up. 50 if not hasattr(self.relax.data, 'sim_state'): 51 raise RelaxError, "Monte Carlo simulations for the run " + `self.run` + " have not been set up." 52 53 # Test if sequence data is loaded. 54 if not self.relax.data.res.has_key(self.run): 55 raise RelaxNoSequenceError, self.run 56 57 # Test the method argument. 58 valid_methods = ['back_calc', 'direct'] 59 if method not in valid_methods: 60 raise RelaxError, "The simulation creation method " + `method` + " is not valid." 61 62 # Function type. 63 function_type = self.relax.data.run_types[self.relax.data.run_names.index(self.run)] 64 65 # Specific Monte Carlo data creation, data return, and error return function setup. 66 create_mc_data = self.relax.specific_setup.setup('create_mc_data', function_type) 67 return_data = self.relax.specific_setup.setup('return_data', function_type) 68 return_error = self.relax.specific_setup.setup('return_error', function_type) 69 pack_sim_data = self.relax.specific_setup.setup('pack_sim_data', function_type) 70 71 # Loop over the sequence. 72 for i in xrange(len(self.relax.data.res[self.run])): 73 # Skip unselected residues. 74 if not self.relax.data.res[self.run][i].select: 75 continue 76 77 # Create the Monte Carlo data. 78 if method == 'back_calc': 79 data = create_mc_data(self.run, i) 80 81 # Get the original data. 82 else: 83 data = return_data(self.run, i) 84 85 # Get the errors. 86 error = return_error(self.run, i) 87 88 # Loop over the Monte Carlo simulations. 89 random = [] 90 for j in xrange(self.relax.data.sim_number[self.run]): 91 # Randomise the data. 92 random.append([]) 93 for k in xrange(len(data)): 94 # No data or errors. 95 if data[k] == None or error[k] == None: 96 random[j].append(None) 97 continue 98 99 # Gaussian randomisation. 100 random[j].append(gauss(data[k], error[k])) 101 102 # Pack the simulation data. 103 pack_sim_data(self.run, i, random)

104 105

106 - def error_analysis(self, run=None, prune=0.0):

107 """Function for calculating errors from the Monte Carlo simulations. 108 109 The standard deviation formula used to calculate the errors is the square root of the 110 bias-corrected variance, given by the formula: 111 112 ____________________________ 113 / 1 114 sd = / ----- * sum({Xi - Xav}^2)] 115 \/ n - 1 116 117 where: 118 n is the total number of simulations. 119 Xi is the parameter value for simulation i. 120 Xav is the mean parameter value for all simulations. 121 """ 122 123 # Arguments. 124 self.run = run 125 126 # Test if the run exists. 127 if not self.run in self.relax.data.run_names: 128 raise RelaxNoRunError, self.run 129 130 # Test if simulations have been set up. 131 if not hasattr(self.relax.data, 'sim_state'): 132 raise RelaxError, "Monte Carlo simulations for the run " + `self.run` + " have not been set up." 133 134 # Function type. 135 function_type = self.relax.data.run_types[self.relax.data.run_names.index(self.run)] 136 137 # Specific number of instances, return simulation chi2 array, return selected simulation array, return simulation parameter array, and set error functions. 138 count_num_instances = self.relax.specific_setup.setup('num_instances', function_type) 139 if prune > 0.0: 140 return_sim_chi2 = self.relax.specific_setup.setup('return_sim_chi2', function_type) 141 return_selected_sim = self.relax.specific_setup.setup('return_selected_sim', function_type) 142 return_sim_param = self.relax.specific_setup.setup('return_sim_param', function_type) 143 set_error = self.relax.specific_setup.setup('set_error', function_type) 144 145 # Count the number of instances. 146 num_instances = count_num_instances(self.run) 147 148 # Loop over the instances. 149 for instance in xrange(num_instances): 150 # Get the selected simulation array. 151 select_sim = return_selected_sim(self.run, instance) 152 153 # Initialise an array of indecies to prune (an empty array means no prunning). 154 indecies_to_skip = [] 155 156 # Pruning. 157 if prune > 0.0: 158 # Get the array of simulation chi-squared values. 159 chi2_array = return_sim_chi2(self.run, instance) 160 161 # The total number of simulations. 162 n = len(chi2_array) 163 164 # Create a sorted array of chi-squared values. 165 chi2_sorted = deepcopy(chi2_array) 166 chi2_sorted.sort() 167 168 # Number of indecies to remove from one side of the chi2 distribution. 169 num = int(float(n) * 0.5 * prune) 170 171 # Remove the lower tail. 172 for i in xrange(num): 173 indecies_to_skip.append(chi2_array.index(chi2_sorted[i])) 174 175 # Remove the upper tail. 176 for i in xrange(n-num, n): 177 indecies_to_skip.append(chi2_array.index(chi2_sorted[i])) 178 179 # Loop over the parameters. 180 index = 0 181 while 1: 182 # Get the array of simulation parameters for the index. 183 param_array = return_sim_param(self.run, instance, index) 184 185 # Break (no more parameters). 186 if param_array == None: 187 break 188 189 # Simulation parameters with values (ie not None). 190 if param_array[0] != None: 191 # The total number of simulations. 192 n = 0 193 for i in xrange(len(param_array)): 194 # Skip unselected simulations. 195 if not select_sim[i]: 196 continue 197 198 # Prune. 199 if i in indecies_to_skip: 200 continue 201 202 # Increment n. 203 n = n + 1 204 205 # Calculate the sum of the parameter value for all simulations. 206 Xsum = 0.0 207 for i in xrange(len(param_array)): 208 # Skip unselected simulations. 209 if not select_sim[i]: 210 continue 211 212 # Prune. 213 if i in indecies_to_skip: 214 continue 215 216 # Sum. 217 Xsum = Xsum + param_array[i] 218 219 # Calculate the mean parameter value for all simulations. 220 if n == 0: 221 Xav = 0.0 222 else: 223 Xav = Xsum / float(n) 224 225 # Calculate the sum part of the standard deviation. 226 sd = 0.0 227 for i in xrange(len(param_array)): 228 # Skip unselected simulations. 229 if not select_sim[i]: 230 continue 231 232 # Prune. 233 if i in indecies_to_skip: 234 continue 235 236 # Sum. 237 sd = sd + (param_array[i] - Xav)**2 238 239 # Calculate the standard deviation. 240 if n <= 1: 241 sd = 0.0 242 else: 243 sd = sqrt(sd / (float(n) - 1.0)) 244 245 # Simulation parameters with the value None. 246 else: 247 sd = None 248 249 # Set the parameter error. 250 set_error(self.run, instance, index, sd) 251 252 # Increment the parameter index. 253 index = index + 1

254 255

256 - def initial_values(self, run=None):

257 """Function for setting the initial simulation parameter values.""" 258 259 # Arguments. 260 self.run = run 261 262 # Test if the run exists. 263 if not self.run in self.relax.data.run_names: 264 raise RelaxNoRunError, self.run 265 266 # Test if simulations have been set up. 267 if not hasattr(self.relax.data, 'sim_state'): 268 raise RelaxError, "Monte Carlo simulations for the run " + `self.run` + " have not been set up." 269 270 # Function type. 271 function_type = self.relax.data.run_types[self.relax.data.run_names.index(self.run)] 272 273 # Specific initial Monte Carlo parameter value function setup. 274 init_sim_values = self.relax.specific_setup.setup('init_sim_values', function_type) 275 276 # Set the initial parameter values. 277 init_sim_values(self.run)

278 279

280 - def off(self, run=None):

281 """Function for turning simulations off.""" 282 283 # Arguments. 284 self.run = run 285 286 # Test if the run exists. 287 if not self.run in self.relax.data.run_names: 288 raise RelaxNoRunError, self.run 289 290 # Test if simulations have been set up. 291 if not hasattr(self.relax.data, 'sim_state'): 292 raise RelaxError, "Monte Carlo simulations for the run " + `self.run` + " have not been set up." 293 294 # Turn simulations off. 295 self.relax.data.sim_state[self.run] = 0

296 297

298 - def on(self, run=None):

299 """Function for turning simulations on.""" 300 301 # Arguments. 302 self.run = run 303 304 # Test if the run exists. 305 if not self.run in self.relax.data.run_names: 306 raise RelaxNoRunError, self.run 307 308 # Test if simulations have been set up. 309 if not hasattr(self.relax.data, 'sim_state'): 310 raise RelaxError, "Monte Carlo simulations for the run " + `self.run` + " have not been set up." 311 312 # Turn simulations on. 313 self.relax.data.sim_state[self.run] = 1

314 315

316 - def select_all_sims(self, number=None, all_select_sim=None):

317 """Function for setting the select flag of all simulations of all instances to one.""" 318 319 # Function type. 320 function_type = self.relax.data.run_types[self.relax.data.run_names.index(self.run)] 321 322 # Specific number of instances and set the selected simulation array functions. 323 count_num_instances = self.relax.specific_setup.setup('num_instances', function_type) 324 set_selected_sim = self.relax.specific_setup.setup('set_selected_sim', function_type) 325 326 # Count the number of instances. 327 num_instances = count_num_instances(self.run) 328 329 # Create the selected simulation array with all simulations selected. 330 if all_select_sim == None: 331 select_sim = ones(number) 332 333 # Loop over the instances. 334 for instance in xrange(num_instances): 335 # Set up the selected simulation array. 336 if all_select_sim != None: 337 select_sim = all_select_sim[instance].tolist() 338 339 # Set the selected simulation array. 340 set_selected_sim(self.run, instance, select_sim)

341 342

343 - def setup(self, run=None, number=None, all_select_sim=None):

344 """Function for setting up Monte Carlo simulations. 345 346 @param run: The name of the run. 347 @type run: str 348 @param number: The number of Monte Carlo simulations to set up. 349 @type number: int 350 @params all_select_sim: The selection status of the Monte Carlo simulations. The first 351 dimension of this matrix corresponds to the simulation and the second corresponds to the 352 instance. 353 @type all_select_sim: Numeric matrix (int) 354 """ 355 356 # Arguments. 357 self.run = run 358 359 # Test if the run exists. 360 if not self.run in self.relax.data.run_names: 361 raise RelaxNoRunError, self.run 362 363 # Test if Monte Carlo simulations have already been set up for the given run. 364 if hasattr(self.relax.data, 'sim_number') and self.relax.data.sim_number.has_key(self.run): 365 raise RelaxError, "Monte Carlo simulations for the run " + `self.run` + " have already been set up." 366 367 # Create the data structure 'sim_number' if it doesn't exist. 368 if not hasattr(self.relax.data, 'sim_number'): 369 self.relax.data.sim_number = {} 370 371 # Add the simulation number. 372 self.relax.data.sim_number[self.run] = number 373 374 # Create the data structure 'sim_state'. 375 if not hasattr(self.relax.data, 'sim_state'): 376 self.relax.data.sim_state = {} 377 378 # Turn simulations on. 379 self.relax.data.sim_state[self.run] = 1 380 381 # Select all simulations. 382 self.select_all_sims(number=number, all_select_sim=all_select_sim)

Source Code for Module generic_fns.monte_carlo