generic_fns.monte

1 ############################################################################### 2 # # 3 # Copyright (C) 2004-2005, 2007-2009 Edward d'Auvergne # 4 # # 5 # This file is part of the program relax. # 6 # # 7 # relax 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 2 of the License, or # 10 # (at your option) any later version. # 11 # # 12 # relax 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 relax; if not, write to the Free Software # 19 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA # 20 # # 21 ############################################################################### 22 23 # Module docstring. 24 """Module for performing Monte Carlo simulations for error analysis.""" 25 26 # Python module imports. 27 from copy import deepcopy 28 from math import sqrt 29 from numpy import ndarray 30 from random import gauss 31 32 # relax module imports. 33 from generic_fns.mol_res_spin import exists_mol_res_spin_data 34 from generic_fns import pipes 35 from relax_errors import RelaxError, RelaxNoSequenceError 36 from specific_fns.setup import get_specific_fn 37 38

39 -def create_data(method=None):

40 """Function for creating simulation data. 41 42 @keyword method: The type of Monte Carlo simulation to perform. 43 @type method: str 44 """ 45 46 # Test if the current data pipe exists. 47 pipes.test() 48 49 # Test if simulations have been set up. 50 if not hasattr(cdp, 'sim_state'): 51 raise RelaxError("Monte Carlo simulations have not been set up.") 52 53 # Test the method argument. 54 valid_methods = ['back_calc', 'direct'] 55 if method not in valid_methods: 56 raise RelaxError("The simulation creation method " + repr(method) + " is not valid.") 57 58 # Specific Monte Carlo data creation, data return, and error return function setup. 59 base_data_loop = get_specific_fn('base_data_loop', cdp.pipe_type) 60 if method == 'back_calc': 61 create_mc_data = get_specific_fn('create_mc_data', cdp.pipe_type) 62 return_data = get_specific_fn('return_data', cdp.pipe_type) 63 return_error = get_specific_fn('return_error', cdp.pipe_type) 64 pack_sim_data = get_specific_fn('pack_sim_data', cdp.pipe_type) 65 66 # Loop over the models. 67 for data_index in base_data_loop(): 68 # Create the Monte Carlo data. 69 if method == 'back_calc': 70 data = create_mc_data(data_index) 71 72 # Get the original data. 73 else: 74 data = return_data(data_index) 75 76 # Get the errors. 77 error = return_error(data_index) 78 79 # List type data. 80 if isinstance(data, list) or isinstance(data, ndarray): 81 # Loop over the Monte Carlo simulations. 82 random = [] 83 for j in xrange(cdp.sim_number): 84 # Randomise the data. 85 random.append([]) 86 for k in xrange(len(data)): 87 # No data or errors. 88 if data[k] == None or error[k] == None: 89 random[j].append(None) 90 continue 91 92 # Gaussian randomisation. 93 random[j].append(gauss(data[k], error[k])) 94 95 # Dictionary type data. 96 if isinstance(data, dict): 97 # Loop over the Monte Carlo simulations. 98 random = [] 99 for j in xrange(cdp.sim_number): 100 # Randomise the data. 101 random.append({}) 102 for id in data.keys(): 103 # No data or errors. 104 if data[id] == None or error[id] == None: 105 random[j][id] = None 106 continue 107 108 # Gaussian randomisation. 109 random[j][id] = gauss(data[id], error[id]) 110 111 # Pack the simulation data. 112 pack_sim_data(data_index, random)

113 114

115 -def error_analysis(prune=0.0):

116 """Function for calculating errors from the Monte Carlo simulations. 117 118 The standard deviation formula used to calculate the errors is the square root of the 119 bias-corrected variance, given by the formula:: 120 121 __________________________ 122 / 1 123 sd = / ----- * sum({Xi - Xav}^2) 124 \/ n - 1 125 126 where 127 - n is the total number of simulations. 128 - Xi is the parameter value for simulation i. 129 - Xav is the mean parameter value for all simulations. 130 131 132 @keyword prune: The amount to prune the upper and lower tails the distribution. If set to 133 0.0, no simulations will be pruned. If set to 1.0, all simulations will be 134 pruned. This argument should be set 0.0 to avoid meaningless statistics. 135 @type prune: float 136 """ 137 138 # Test if the current data pipe exists. 139 pipes.test() 140 141 # Test if simulations have been set up. 142 if not hasattr(cdp, 'sim_state'): 143 raise RelaxError("Monte Carlo simulations have not been set up.") 144 145 # Model loop, return simulation chi2 array, return selected simulation array, return simulation parameter array, and set error functions. 146 model_loop = get_specific_fn('model_loop', cdp.pipe_type) 147 if prune > 0.0: 148 return_sim_chi2 = get_specific_fn('return_sim_chi2', cdp.pipe_type) 149 return_selected_sim = get_specific_fn('return_selected_sim', cdp.pipe_type) 150 return_sim_param = get_specific_fn('return_sim_param', cdp.pipe_type) 151 set_error = get_specific_fn('set_error', cdp.pipe_type) 152 153 # Loop over the models. 154 for model_info in model_loop(): 155 # Get the selected simulation array. 156 select_sim = return_selected_sim(model_info) 157 158 # Initialise an array of indices to prune (an empty array means no pruning). 159 indices_to_skip = [] 160 161 # Pruning. 162 if prune > 0.0: 163 # Get the array of simulation chi-squared values. 164 chi2_array = return_sim_chi2(model_info) 165 166 # The total number of simulations. 167 n = len(chi2_array) 168 169 # Create a sorted array of chi-squared values. 170 chi2_sorted = sorted(deepcopy(chi2_array)) 171 172 # Number of indices to remove from one side of the chi2 distribution. 173 num = int(float(n) * 0.5 * prune) 174 175 # Remove the lower tail. 176 for i in xrange(num): 177 indices_to_skip.append(chi2_array.index(chi2_sorted[i])) 178 179 # Remove the upper tail. 180 for i in xrange(n-num, n): 181 indices_to_skip.append(chi2_array.index(chi2_sorted[i])) 182 183 # Loop over the parameters. 184 index = 0 185 while True: 186 # Get the array of simulation parameters for the index. 187 param_array = return_sim_param(model_info, index) 188 189 # Break (no more parameters). 190 if param_array == None: 191 break 192 193 # Simulation parameters with values (ie not None). 194 if param_array[0] != None: 195 # The total number of simulations. 196 n = 0 197 for i in xrange(len(param_array)): 198 # Skip deselected simulations. 199 if not select_sim[i]: 200 continue 201 202 # Prune. 203 if i in indices_to_skip: 204 continue 205 206 # Increment n. 207 n = n + 1 208 209 # Calculate the sum of the parameter value for all simulations. 210 Xsum = 0.0 211 for i in xrange(len(param_array)): 212 # Skip deselected simulations. 213 if not select_sim[i]: 214 continue 215 216 # Prune. 217 if i in indices_to_skip: 218 continue 219 220 # Sum. 221 Xsum = Xsum + param_array[i] 222 223 # Calculate the mean parameter value for all simulations. 224 if n == 0: 225 Xav = 0.0 226 else: 227 Xav = Xsum / float(n) 228 229 # Calculate the sum part of the standard deviation. 230 sd = 0.0 231 for i in xrange(len(param_array)): 232 # Skip deselected simulations. 233 if not select_sim[i]: 234 continue 235 236 # Prune. 237 if i in indices_to_skip: 238 continue 239 240 # Sum. 241 sd = sd + (param_array[i] - Xav)**2 242 243 # Calculate the standard deviation. 244 if n <= 1: 245 sd = 0.0 246 else: 247 sd = sqrt(sd / (float(n) - 1.0)) 248 249 # Simulation parameters with the value None. 250 else: 251 sd = None 252 253 # Set the parameter error. 254 set_error(model_info, index, sd) 255 256 # Increment the parameter index. 257 index = index + 1 258 259 # Turn off the Monte Carlo simulation state, as the MC analysis is now finished. 260 cdp.sim_state = False

261 262

263 -def initial_values():

264 """Set the initial simulation parameter values.""" 265 266 # Test if the current data pipe exists. 267 pipes.test() 268 269 # Test if simulations have been set up. 270 if not hasattr(cdp, 'sim_state'): 271 raise RelaxError("Monte Carlo simulations have not been set up.") 272 273 # Specific initial Monte Carlo parameter value function setup. 274 init_sim_values = get_specific_fn('init_sim_values', cdp.pipe_type) 275 276 # Set the initial parameter values. 277 init_sim_values()

278 279

280 -def off():

281 """Turn simulations off.""" 282 283 # Test if the current data pipe exists. 284 pipes.test() 285 286 # Test if simulations have been set up. 287 if not hasattr(cdp, 'sim_state'): 288 raise RelaxError("Monte Carlo simulations have not been set up.") 289 290 # Turn simulations off. 291 cdp.sim_state = False

292 293

294 -def on():

295 """Turn simulations on.""" 296 297 # Test if the current data pipe exists. 298 pipes.test() 299 300 # Test if simulations have been set up. 301 if not hasattr(cdp, 'sim_state'): 302 raise RelaxError("Monte Carlo simulations have not been set up.") 303 304 # Turn simulations on. 305 cdp.sim_state = True

306 307

308 -def select_all_sims(number=None, all_select_sim=None):

309 """Set the select flag of all simulations of all models to one. 310 311 @keyword number: The number of Monte Carlo simulations to set up. 312 @type number: int 313 @keyword all_select_sim: The selection status of the Monte Carlo simulations. The first 314 dimension of this matrix corresponds to the simulation and the 315 second corresponds to the models. 316 @type all_select_sim: list of lists of bool 317 """ 318 319 # Model loop and set the selected simulation array functions. 320 model_loop = get_specific_fn('model_loop', cdp.pipe_type) 321 set_selected_sim = get_specific_fn('set_selected_sim', cdp.pipe_type) 322 323 # Create the selected simulation array with all simulations selected. 324 if all_select_sim == None: 325 select_sim = [True]*number 326 327 # Loop over the models. 328 i = 0 329 for model_info in model_loop(): 330 # Set up the selected simulation array. 331 if all_select_sim != None: 332 select_sim = all_select_sim[i] 333 334 # Set the selected simulation array. 335 set_selected_sim(model_info, select_sim) 336 337 # Model index. 338 i = i + 1

339 340

341 -def setup(number=None, all_select_sim=None):

342 """Function for setting up Monte Carlo simulations. 343 344 @keyword number: The number of Monte Carlo simulations to set up. 345 @type number: int 346 @keyword all_select_sim: The selection status of the Monte Carlo simulations. The first 347 dimension of this matrix corresponds to the simulation and the 348 second corresponds to the instance. 349 @type all_select_sim: list of lists of bool 350 """ 351 352 # Test if the current data pipe exists. 353 pipes.test() 354 355 # Create a number of MC sim data structures. 356 cdp.sim_number = number 357 cdp.sim_state = True 358 359 # Select all simulations. 360 select_all_sims(number=number, all_select_sim=all_select_sim)

361

Source Code for Module generic_fns.monte_carlo