pipe_control.monte

1 ############################################################################### 2 # # 3 # Copyright (C) 2004-2014 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.api import return_api 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 # The specific analysis API object. 56 api = return_api() 57 58 # Loop over the models. 59 for data_index in api.base_data_loop(): 60 # Create the Monte Carlo data. 61 if method == 'back_calc': 62 data = api.create_mc_data(data_index) 63 64 # Get the original data. 65 else: 66 data = api.return_data(data_index) 67 68 # No data, so skip. 69 if data == None: 70 continue 71 72 # Get the errors. 73 error = api.return_error(data_index) 74 75 # List type data. 76 if isinstance(data, list) or isinstance(data, ndarray): 77 # Loop over the Monte Carlo simulations. 78 random = [] 79 for j in range(cdp.sim_number): 80 # Randomise the data. 81 random.append([]) 82 for k in range(len(data)): 83 # No data or errors. 84 if data[k] == None or error[k] == None: 85 random[j].append(None) 86 continue 87 88 # Gaussian randomisation. 89 random[j].append(gauss(data[k], error[k])) 90 91 # Dictionary type data. 92 if isinstance(data, dict): 93 # Loop over the Monte Carlo simulations. 94 random = [] 95 for j in range(cdp.sim_number): 96 # Randomise the data. 97 random.append({}) 98 for id in data.keys(): 99 # No data or errors. 100 if data[id] == None or error[id] == None: 101 random[j][id] = None 102 continue 103 104 # Gaussian randomisation. 105 random[j][id] = gauss(data[id], error[id]) 106 107 # Pack the simulation data. 108 api.sim_pack_data(data_index, random)

109 110

111 -def error_analysis():

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

199 200

201 -def initial_values():

202 """Set the initial simulation parameter values.""" 203 204 # Test if the current data pipe exists. 205 pipes.test() 206 207 # Test if simulations have been set up. 208 if not hasattr(cdp, 'sim_state'): 209 raise RelaxError("Monte Carlo simulations have not been set up.") 210 211 # The specific analysis API object. 212 api = return_api() 213 214 # Set the initial parameter values. 215 api.sim_init_values()

216 217

218 -def off():

219 """Turn simulations off.""" 220 221 # Test if the current data pipe exists. 222 pipes.test() 223 224 # Test if simulations have been set up. 225 if not hasattr(cdp, 'sim_state'): 226 raise RelaxError("Monte Carlo simulations have not been set up.") 227 228 # Turn simulations off. 229 cdp.sim_state = False

230 231

232 -def on():

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

244 245

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

247 """Set the select flag of all simulations of all models to one. 248 249 @keyword number: The number of Monte Carlo simulations to set up. 250 @type number: int 251 @keyword all_select_sim: The selection status of the Monte Carlo simulations. The first 252 dimension of this matrix corresponds to the simulation and the 253 second corresponds to the models. 254 @type all_select_sim: list of lists of bool 255 """ 256 257 # The specific analysis API object. 258 api = return_api() 259 260 # Create the selected simulation array with all simulations selected. 261 if all_select_sim == None: 262 select_sim = [True]*number 263 264 # Loop over the models. 265 i = 0 266 for model_info in api.model_loop(): 267 # Skip function. 268 if api.skip_function(model_info): 269 continue 270 271 # Set up the selected simulation array. 272 if all_select_sim != None: 273 select_sim = all_select_sim[i] 274 275 # Set the selected simulation array. 276 api.set_selected_sim(model_info, select_sim) 277 278 # Model index. 279 i += 1

280 281

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

283 """Function for setting up Monte Carlo simulations. 284 285 @keyword number: The number of Monte Carlo simulations to set up. 286 @type number: int 287 @keyword all_select_sim: The selection status of the Monte Carlo simulations. The first 288 dimension of this matrix corresponds to the simulation and the 289 second corresponds to the instance. 290 @type all_select_sim: list of lists of bool 291 """ 292 293 # Test if the current data pipe exists. 294 pipes.test() 295 296 # Create a number of MC sim data structures. 297 cdp.sim_number = number 298 cdp.sim_state = True 299 300 # Select all simulations. 301 select_all_sims(number=number, all_select_sim=all_select_sim)

302

Source Code for Module pipe_control.monte_carlo