generic_fns.monte

29 - def __init__(self, relax):

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

33 34

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

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

103 104

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

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

247 248

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

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

271 272

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

274 """Function for turning simulations off.""" 275 276 # Arguments. 277 self.run = run 278 279 # Test if the run exists. 280 if not self.run in self.relax.data.run_names: 281 raise RelaxNoRunError, self.run 282 283 # Test if simulations have been set up. 284 if not hasattr(self.relax.data, 'sim_state'): 285 raise RelaxError, "Monte Carlo simulations for the run " + `self.run` + " have not been set up." 286 287 # Turn simulations off. 288 self.relax.data.sim_state[self.run] = 0

289 290

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

292 """Function for turning simulations on.""" 293 294 # Arguments. 295 self.run = run 296 297 # Test if the run exists. 298 if not self.run in self.relax.data.run_names: 299 raise RelaxNoRunError, self.run 300 301 # Test if simulations have been set up. 302 if not hasattr(self.relax.data, 'sim_state'): 303 raise RelaxError, "Monte Carlo simulations for the run " + `self.run` + " have not been set up." 304 305 # Turn simulations on. 306 self.relax.data.sim_state[self.run] = 1

307 308

309 - def select_all_sims(self, number=None, select_sim=None):

310 """Function for setting the select flag of all simulations of all instances to one.""" 311 312 # Function type. 313 function_type = self.relax.data.run_types[self.relax.data.run_names.index(self.run)] 314 315 # Specific number of instances and set the selected simulation array functions. 316 count_num_instances = self.relax.specific_setup.setup('num_instances', function_type) 317 set_selected_sim = self.relax.specific_setup.setup('set_selected_sim', function_type) 318 319 # Count the number of instances. 320 num_instances = count_num_instances(self.run) 321 322 # Create the selected simulation array with all simulations selected. 323 if select_sim == None: 324 select_sim = [] 325 for i in xrange(number): 326 select_sim.append(1) 327 328 # Loop over the instances. 329 for instance in xrange(num_instances): 330 # Set the selected simulation array. 331 set_selected_sim(self.run, instance, select_sim)

332 333

334 - def setup(self, run=None, number=None, select_sim=None):

335 """Function for setting up Monte Carlo simulations.""" 336 337 # Arguments. 338 self.run = run 339 340 # Test if the run exists. 341 if not self.run in self.relax.data.run_names: 342 raise RelaxNoRunError, self.run 343 344 # Test if Monte Carlo simulations have already been set up for the given run. 345 if hasattr(self.relax.data, 'sim_number') and self.relax.data.sim_number.has_key(self.run): 346 raise RelaxError, "Monte Carlo simulations for the run " + `self.run` + " have already been set up." 347 348 # Create the data structure 'sim_number' if it doesn't exist. 349 if not hasattr(self.relax.data, 'sim_number'): 350 self.relax.data.sim_number = {} 351 352 # Add the simulation number. 353 if select_sim: 354 self.relax.data.sim_number[self.run] = len(select_sim) 355 else: 356 self.relax.data.sim_number[self.run] = number 357 358 # Create the data structure 'sim_state'. 359 if not hasattr(self.relax.data, 'sim_state'): 360 self.relax.data.sim_state = {} 361 362 # Turn simulations on. 363 self.relax.data.sim_state[self.run] = 1 364 365 # Select all simulations. 366 self.select_all_sims(number=number, select_sim=select_sim)

Source Code for Module generic_fns.monte_carlo