generic_fns.monte

1 ############################################################################### 2 # # 3 # Copyright (C) 2004-2013 Edward d'Auvergne # 4 # # 5 # This file is part of the program relax (http://www.nmr-relax.com). # 6 # # 7 # This program 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 3 of the License, or # 10 # (at your option) any later version. # 11 # # 12 # This program 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 this program. If not, see <http://www.gnu.org/licenses/>. # 19 # # 20 ############################################################################### 21 22 # Module docstring. 23 """Module for performing Monte Carlo simulations for error analysis.""" 24 25 # Python module imports. 26 from copy import deepcopy 27 from math import sqrt 28 from numpy import ndarray 29 from random import gauss 30 31 # relax module imports. 32 from generic_fns.mol_res_spin import exists_mol_res_spin_data 33 from generic_fns import pipes 34 from relax_errors import RelaxError, RelaxNoSequenceError 35 from specific_fns.setup import get_specific_fn 36 37

38 -def create_data(method=None):

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

112 113

114 -def error_analysis():

115 """Function for calculating errors from the Monte Carlo simulations. 116 117 The standard deviation formula used to calculate the errors is the square root of the 118 bias-corrected variance, given by the formula:: 119 120 __________________________ 121 / 1 122 sd = / ----- * sum({Xi - Xav}^2) 123 \/ n - 1 124 125 where 126 - n is the total number of simulations. 127 - Xi is the parameter value for simulation i. 128 - Xav is the mean parameter value for all simulations. 129 """ 130 131 # Test if the current data pipe exists. 132 pipes.test() 133 134 # Test if simulations have been set up. 135 if not hasattr(cdp, 'sim_state'): 136 raise RelaxError("Monte Carlo simulations have not been set up.") 137 138 # Model loop, return simulation chi2 array, return selected simulation array, return simulation parameter array, and set error functions. 139 model_loop = get_specific_fn('model_loop', cdp.pipe_type) 140 return_selected_sim = get_specific_fn('return_selected_sim', cdp.pipe_type) 141 return_sim_param = get_specific_fn('return_sim_param', cdp.pipe_type) 142 set_error = get_specific_fn('set_error', cdp.pipe_type) 143 144 # Loop over the models. 145 for model_info in model_loop(): 146 # Get the selected simulation array. 147 select_sim = return_selected_sim(model_info) 148 149 # Loop over the parameters. 150 index = 0 151 while True: 152 # Get the array of simulation parameters for the index. 153 param_array = return_sim_param(model_info, index) 154 155 # Break (no more parameters). 156 if param_array == None: 157 break 158 159 # Simulation parameters with values (ie not None). 160 if param_array[0] != None: 161 # The total number of simulations. 162 n = 0 163 for i in range(len(param_array)): 164 # Skip deselected simulations. 165 if not select_sim[i]: 166 continue 167 168 # Increment n. 169 n = n + 1 170 171 # Calculate the sum of the parameter value for all simulations. 172 Xsum = 0.0 173 for i in range(len(param_array)): 174 # Skip deselected simulations. 175 if not select_sim[i]: 176 continue 177 178 # Sum. 179 Xsum = Xsum + param_array[i] 180 181 # Calculate the mean parameter value for all simulations. 182 if n == 0: 183 Xav = 0.0 184 else: 185 Xav = Xsum / float(n) 186 187 # Calculate the sum part of the standard deviation. 188 sd = 0.0 189 for i in range(len(param_array)): 190 # Skip deselected simulations. 191 if not select_sim[i]: 192 continue 193 194 # Sum. 195 sd = sd + (param_array[i] - Xav)**2 196 197 # Calculate the standard deviation. 198 if n <= 1: 199 sd = 0.0 200 else: 201 sd = sqrt(sd / (float(n) - 1.0)) 202 203 # Simulation parameters with the value None. 204 else: 205 sd = None 206 207 # Set the parameter error. 208 set_error(model_info, index, sd) 209 210 # Increment the parameter index. 211 index = index + 1 212 213 # Turn off the Monte Carlo simulation state, as the MC analysis is now finished. 214 cdp.sim_state = False

215 216

217 -def initial_values():

218 """Set the initial simulation parameter values.""" 219 220 # Test if the current data pipe exists. 221 pipes.test() 222 223 # Test if simulations have been set up. 224 if not hasattr(cdp, 'sim_state'): 225 raise RelaxError("Monte Carlo simulations have not been set up.") 226 227 # Specific initial Monte Carlo parameter value function setup. 228 init_sim_values = get_specific_fn('init_sim_values', cdp.pipe_type) 229 230 # Set the initial parameter values. 231 init_sim_values()

232 233

234 -def off():

235 """Turn simulations off.""" 236 237 # Test if the current data pipe exists. 238 pipes.test() 239 240 # Test if simulations have been set up. 241 if not hasattr(cdp, 'sim_state'): 242 raise RelaxError("Monte Carlo simulations have not been set up.") 243 244 # Turn simulations off. 245 cdp.sim_state = False

246 247

248 -def on():

249 """Turn simulations on.""" 250 251 # Test if the current data pipe exists. 252 pipes.test() 253 254 # Test if simulations have been set up. 255 if not hasattr(cdp, 'sim_state'): 256 raise RelaxError("Monte Carlo simulations have not been set up.") 257 258 # Turn simulations on. 259 cdp.sim_state = True

260 261

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

263 """Set the select flag of all simulations of all models to one. 264 265 @keyword number: The number of Monte Carlo simulations to set up. 266 @type number: int 267 @keyword all_select_sim: The selection status of the Monte Carlo simulations. The first 268 dimension of this matrix corresponds to the simulation and the 269 second corresponds to the models. 270 @type all_select_sim: list of lists of bool 271 """ 272 273 # Model loop and set the selected simulation array functions. 274 model_loop = get_specific_fn('model_loop', cdp.pipe_type) 275 set_selected_sim = get_specific_fn('set_selected_sim', cdp.pipe_type) 276 skip_function = get_specific_fn('skip_function', cdp.pipe_type) 277 278 # Create the selected simulation array with all simulations selected. 279 if all_select_sim == None: 280 select_sim = [True]*number 281 282 # Loop over the models. 283 i = 0 284 for model_info in model_loop(): 285 # Skip function. 286 if skip_function(model_info): 287 continue 288 289 # Set up the selected simulation array. 290 if all_select_sim != None: 291 select_sim = all_select_sim[i] 292 293 # Set the selected simulation array. 294 set_selected_sim(model_info, select_sim) 295 296 # Model index. 297 i += 1

298 299

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

301 """Function for setting up Monte Carlo simulations. 302 303 @keyword number: The number of Monte Carlo simulations to set up. 304 @type number: int 305 @keyword all_select_sim: The selection status of the Monte Carlo simulations. The first 306 dimension of this matrix corresponds to the simulation and the 307 second corresponds to the instance. 308 @type all_select_sim: list of lists of bool 309 """ 310 311 # Test if the current data pipe exists. 312 pipes.test() 313 314 # Create a number of MC sim data structures. 315 cdp.sim_number = number 316 cdp.sim_state = True 317 318 # Select all simulations. 319 select_all_sims(number=number, all_select_sim=all_select_sim)

320

Source Code for Module generic_fns.monte_carlo