pipe_control.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 numpy import ndarray 27 from random import gauss 28 29 # relax module imports. 30 from lib import statistics 31 from lib.errors import RelaxError 32 from pipe_control import pipes 33 from specific_analyses.setup import get_specific_fn 34 35

36 -def create_data(method=None):

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

114 115

116 -def error_analysis():

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

207 208

209 -def initial_values():

210 """Set the initial simulation parameter values.""" 211 212 # Test if the current data pipe exists. 213 pipes.test() 214 215 # Test if simulations have been set up. 216 if not hasattr(cdp, 'sim_state'): 217 raise RelaxError("Monte Carlo simulations have not been set up.") 218 219 # Specific initial Monte Carlo parameter value function setup. 220 init_sim_values = get_specific_fn('init_sim_values', cdp.pipe_type) 221 222 # Set the initial parameter values. 223 init_sim_values()

224 225

226 -def off():

227 """Turn simulations off.""" 228 229 # Test if the current data pipe exists. 230 pipes.test() 231 232 # Test if simulations have been set up. 233 if not hasattr(cdp, 'sim_state'): 234 raise RelaxError("Monte Carlo simulations have not been set up.") 235 236 # Turn simulations off. 237 cdp.sim_state = False

238 239

240 -def on():

241 """Turn simulations on.""" 242 243 # Test if the current data pipe exists. 244 pipes.test() 245 246 # Test if simulations have been set up. 247 if not hasattr(cdp, 'sim_state'): 248 raise RelaxError("Monte Carlo simulations have not been set up.") 249 250 # Turn simulations on. 251 cdp.sim_state = True

252 253

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

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

290 291

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

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

312

Source Code for Module pipe_control.monte_carlo