Source Code for Module generic_fns.monte_carlo

  1  ############################################################################### 
  2  #                                                                             # 
  3  # Copyright (C) 2004 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  from copy import deepcopy 
 24  from math import sqrt 
 25  from random import gauss 
 26   
 27   
 28 -class Monte_carlo: 
 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          # Test if the run exists. 
 42          if not run in self.relax.data.run_names: 
 43              raise RelaxNoRunError, run 
 44   
 45          # Test if simulations have been set up. 
 46          if not hasattr(self.relax.data, 'sim_state'): 
 47              raise RelaxError, "Monte Carlo simulations for the run " + `run` + " have not been set up." 
 48   
 49          # Test if sequence data is loaded. 
 50          if not self.relax.data.res.has_key(run): 
 51              raise RelaxNoSequenceError, run 
 52   
 53          # Test the method argument. 
 54          valid_methods = ['back_calc', 'direct'] 
 55          if method not in valid_methods: 
 56              raise RelaxError, "The method " + `method` + " is not valid." 
 57   
 58          # Function type. 
 59          function_type = self.relax.data.run_types[self.relax.data.run_names.index(run)] 
 60   
 61          # Specific Monte Carlo data creation, data return, and error return function setup. 
 62          create_mc_data = self.relax.specific_setup.setup('create_mc_data', function_type) 
 63          return_data = self.relax.specific_setup.setup('return_data', function_type) 
 64          return_error = self.relax.specific_setup.setup('return_error', function_type) 
 65          pack_sim_data = self.relax.specific_setup.setup('pack_sim_data', function_type) 
 66   
 67          # Loop over the sequence. 
 68          for i in xrange(len(self.relax.data.res[run])): 
 69              # Skip unselected residues. 
 70              if not self.relax.data.res[run][i].select: 
 71                  continue 
 72   
 73              # Create the Monte Carlo data. 
 74              if method == 'back_calc': 
 75                  data = create_mc_data(run, i) 
 76   
 77              # Get the original data. 
 78              else: 
 79                  data = return_data(run, i) 
 80   
 81              # Get the errors. 
 82              error = return_error(run, i) 
 83   
 84              # Loop over the Monte Carlo simulations. 
 85              random = [] 
 86              for j in xrange(self.relax.data.sim_number[run]): 
 87                  # Randomise the data. 
 88                  random.append([]) 
 89                  for k in xrange(len(data)): 
 90                      # No data or errors. 
 91                      if data[k] == None or error[k] == None: 
 92                          random[j].append(None) 
 93                          continue 
 94   
 95                      # Gaussian randomisation. 
 96                      random[j].append(gauss(data[k], error[k])) 
 97   
 98              # Pack the simulation data. 
 99              pack_sim_data(run, i, random) 
100   
101   
102 -    def error_analysis(self, run=None, prune=0): 
103          """Function for calculating errors from the Monte Carlo simulations. 
104   
105          The standard deviation formula used to calculate the errors is the square root of the 
106          bias-corrected variance, given by the formula: 
107   
108                         ____________________________ 
109                        /   1 
110              sd  =    /  ----- * sum({Xi - Xav}^2)] 
111                     \/   n - 1 
112   
113          where: 
114              n is the total number of simulations. 
115              Xi is the parameter value for simulation i. 
116              Xav is the mean parameter value for all simulations. 
117          """ 
118   
119          # Test if the run exists. 
120          if not run in self.relax.data.run_names: 
121              raise RelaxNoRunError, run 
122   
123          # Test if simulations have been set up. 
124          if not hasattr(self.relax.data, 'sim_state'): 
125              raise RelaxError, "Monte Carlo simulations for the run " + `run` + " have not been set up." 
126   
127          # Function type. 
128          function_type = self.relax.data.run_types[self.relax.data.run_names.index(run)] 
129   
130          # Specific number of instances, return simulation chi2 array, return simulation parameter array, and set error functions. 
131          count_num_instances = self.relax.specific_setup.setup('num_instances', function_type) 
132          if prune > 0.0: 
133              return_sim_chi2 = self.relax.specific_setup.setup('return_sim_chi2', function_type) 
134          return_sim_param = self.relax.specific_setup.setup('return_sim_param', function_type) 
135          set_error = self.relax.specific_setup.setup('set_error', function_type) 
136   
137          # Count the number of instances. 
138          num_instances = count_num_instances(run) 
139   
140          # Loop over the instances. 
141          for instance in xrange(num_instances): 
142              # Initialise an array of indecies to prune. 
143              indecies_to_skip = [] 
144   
145              # Pruning. 
146              if prune > 0.0: 
147                  # Get the array of simulation chi-squared values. 
148                  chi2_array = return_sim_chi2(run, instance) 
149   
150                  # The total number of simulations. 
151                  n = len(chi2_array) 
152   
153                  # Create a sorted array of chi-squared values. 
154                  chi2_sorted = deepcopy(chi2_array) 
155                  chi2_sorted.sort() 
156   
157                  # Number of indecies to remove from one side of the chi2 distribution. 
158                  num = int(float(n) * 0.5 * prune) 
159   
160                  # Remove the lower tail. 
161                  for i in xrange(num): 
162                      indecies_to_skip.append(chi2_array.index(chi2_sorted[i])) 
163   
164                  # Remove the upper tail. 
165                  for i in xrange(n-num, n): 
166                      indecies_to_skip.append(chi2_array.index(chi2_sorted[i])) 
167   
168              # Loop over the parameters. 
169              index = 0 
170              while 1: 
171                  # Get the array of simulation parameters for the index. 
172                  param_array = return_sim_param(run, instance, index) 
173   
174                  # Break (no more parameters). 
175                  if param_array == None: 
176                      break 
177   
178                  # Simulation parameters with values (ie not None). 
179                  if param_array[0] != None: 
180                      # The total number of simulations. 
181                      n = len(param_array) 
182   
183                      # Calculate the mean parameter value for all simulations. 
184                      Xav = 0.0 
185                      for i in xrange(n): 
186                          # Prune. 
187                          if i in indecies_to_skip: 
188                              continue 
189   
190                          # Sum. 
191                          Xav = Xav + param_array[i] 
192                      Xav = Xav / (float(n) - float(len(indecies_to_skip))) 
193   
194                      # Calculate the sum part of the standard deviation. 
195                      sd = 0.0 
196                      for i in xrange(n): 
197                          # Prune. 
198                          if i in indecies_to_skip: 
199                              continue 
200   
201                          # Sum. 
202                          sd = sd + (param_array[i] - Xav)**2 
203   
204                      # Calculate the standard deviation. 
205                      sd = sqrt(sd / (float(n) - float(len(indecies_to_skip)) - 1.0)) 
206   
207                  # Simulation parameters with the value None. 
208                  else: 
209                      sd = None 
210   
211                  # Set the parameter error. 
212                  set_error(run, instance, index, sd) 
213   
214                  # Increment the parameter index. 
215                  index = index + 1 
216   
217   
218 -    def initial_values(self, run=None): 
219          """Function for setting the initial simulation parameter values.""" 
220   
221          # Test if the run exists. 
222          if not run in self.relax.data.run_names: 
223              raise RelaxNoRunError, run 
224   
225          # Test if simulations have been set up. 
226          if not hasattr(self.relax.data, 'sim_state'): 
227              raise RelaxError, "Monte Carlo simulations for the run " + `run` + " have not been set up." 
228   
229          # Function type. 
230          function_type = self.relax.data.run_types[self.relax.data.run_names.index(run)] 
231   
232          # Specific initial Monte Carlo parameter value function setup. 
233          init_sim_values = self.relax.specific_setup.setup('init_sim_values', function_type) 
234   
235          # Set the initial parameter values. 
236          init_sim_values(run) 
237   
238   
239 -    def off(self, run=None): 
240          """Function for turning simulations off.""" 
241   
242          # Test if the run exists. 
243          if not run in self.relax.data.run_names: 
244              raise RelaxNoRunError, run 
245   
246          # Test if simulations have been set up. 
247          if not hasattr(self.relax.data, 'sim_state'): 
248              raise RelaxError, "Monte Carlo simulations for the run " + `run` + " have not been set up." 
249   
250          # Turn simulations off. 
251          self.relax.data.sim_state[run] = 0 
252   
253   
254 -    def on(self, run=None): 
255          """Function for turning simulations on.""" 
256   
257          # Test if the run exists. 
258          if not run in self.relax.data.run_names: 
259              raise RelaxNoRunError, run 
260   
261          # Test if simulations have been set up. 
262          if not hasattr(self.relax.data, 'sim_state'): 
263              raise RelaxError, "Monte Carlo simulations for the run " + `run` + " have not been set up." 
264   
265          # Turn simulations on. 
266          self.relax.data.sim_state[run] = 1 
267   
268   
269 -    def setup(self, run=None, number=0): 
270          """Function for setting up Monte Carlo simulations.""" 
271   
272          # Test if the run exists. 
273          if not run in self.relax.data.run_names: 
274              raise RelaxNoRunError, run 
275   
276          # Test if Monte Carlo simulations have already been set up for the given run. 
277          if hasattr(self.relax.data, 'sim_number') and self.relax.data.sim_number.has_key(run): 
278              raise RelaxError, "Monte Carlo simulations for the run " + `run` + " have already been set up." 
279   
280          # Create the data structure 'sim_number' if it doesn't exist. 
281          if not hasattr(self.relax.data, 'sim_number'): 
282              self.relax.data.sim_number = {} 
283   
284          # Add the simulation number. 
285          self.relax.data.sim_number[run] = number 
286   
287          # Create the data structure 'sim_state'. 
288          if not hasattr(self.relax.data, 'sim_state'): 
289              self.relax.data.sim_state = {} 
290   
291          # Turn simulations on. 
292          self.relax.data.sim_state[run] = 1 
293