"Fossies" - the Fresh Open Source Software Archive

Member "relax-5.0.0/specific_analyses/relax_disp/api.py" (18 Nov 2019, 63612 Bytes) of package /linux/privat/relax-5.0.0.src.tar.bz2:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) Python source code syntax highlighting (style: standard) with prefixed line numbers. Alternatively you can here view or download the uninterpreted source code file. See also the latest Fossies "Diffs" side-by-side code changes report for "api.py": 4.1.3_vs_5.0.0.

    1 ###############################################################################
    2 #                                                                             #
    3 # Copyright (C) 2003-2008,2013-2014,2016,2019 Edward d'Auvergne               #
    4 # Copyright (C) 2006 Chris MacRaild                                           #
    5 # Copyright (C) 2008-2009 Sebastien Morin                                     #
    6 # Copyright (C) 2013-2015 Troels E. Linnet                                    #
    7 #                                                                             #
    8 # This file is part of the program relax (http://www.nmr-relax.com).          #
    9 #                                                                             #
   10 # This program is free software: you can redistribute it and/or modify        #
   11 # it under the terms of the GNU General Public License as published by        #
   12 # the Free Software Foundation, either version 3 of the License, or           #
   13 # (at your option) any later version.                                         #
   14 #                                                                             #
   15 # This program is distributed in the hope that it will be useful,             #
   16 # but WITHOUT ANY WARRANTY; without even the implied warranty of              #
   17 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               #
   18 # GNU General Public License for more details.                                #
   19 #                                                                             #
   20 # You should have received a copy of the GNU General Public License           #
   21 # along with this program.  If not, see <http://www.gnu.org/licenses/>.       #
   22 #                                                                             #
   23 ###############################################################################
   24 
   25 # Module docstring.
   26 """The relaxation dispersion API object."""
   27 
   28 # Python module imports.
   29 import dep_check
   30 if dep_check.bmrblib_module:
   31     import bmrblib
   32 from copy import deepcopy
   33 from numpy import int32, sqrt, zeros
   34 from re import match, search
   35 import string
   36 import sys
   37 from types import MethodType
   38 
   39 # relax module imports.
   40 from lib.arg_check import is_list, is_str_list
   41 from lib.dispersion.variables import EXP_TYPE_CPMG_PROTON_MQ, EXP_TYPE_CPMG_PROTON_SQ, MODEL_LIST_MMQ, MODEL_R2EFF, PARAMS_R20
   42 from lib.errors import RelaxError, RelaxImplementError
   43 from lib.text.sectioning import subsection
   44 from multi import Processor_box
   45 from pipe_control import pipes, sequence
   46 from pipe_control.exp_info import bmrb_write_citations, bmrb_write_methods, bmrb_write_software
   47 from pipe_control.mol_res_spin import bmrb_write_entity, check_mol_res_spin_data, get_molecule_names, return_spin, spin_loop
   48 from pipe_control.pipes import check_pipe
   49 from pipe_control.spectrometer import check_spectrometer_setup
   50 from pipe_control.sequence import return_attached_protons
   51 from specific_analyses.api_base import API_base
   52 from specific_analyses.api_common import API_common
   53 from specific_analyses.relax_disp.checks import check_model_type
   54 from specific_analyses.relax_disp.data import average_intensity, calc_rotating_frame_params, find_intensity_keys, generate_r20_key, has_exponential_exp_type, has_proton_mmq_cpmg, loop_cluster, loop_exp_frq, loop_exp_frq_offset_point, loop_time, pack_back_calc_r2eff, return_param_key_from_data, spin_ids_to_containers
   55 from specific_analyses.relax_disp.optimisation import Disp_memo, Disp_minimise_command, back_calc_peak_intensities, back_calc_r2eff, calculate_r2eff, minimise_r2eff
   56 from specific_analyses.relax_disp.parameter_object import Relax_disp_params
   57 from specific_analyses.relax_disp.parameters import get_param_names, get_value, loop_parameters, param_conversion, param_index_to_param_info, param_num, r1_setup
   58 
   59 
   60 class Relax_disp(API_base, API_common):
   61     """Class containing functions for relaxation dispersion curve fitting."""
   62 
   63     # Class variable for storing the class instance (for the singleton design pattern).
   64     instance = None
   65 
   66     def __init__(self):
   67         """Initialise the class by placing API_common methods into the API."""
   68 
   69         # Place methods into the API.
   70         self.data_init = self._data_init_spin
   71         self.model_type = self._model_type_local
   72         self.return_conversion_factor = self._return_no_conversion_factor
   73         self.return_value = self.return_value
   74 
   75         # Place a copy of the parameter list object in the instance namespace.
   76         self._PARAMS = Relax_disp_params()
   77 
   78 
   79     def base_data_loop(self):
   80         """Custom generator method for looping over the base data.
   81 
   82         For the R2eff model, the base data is the peak intensity data defining a single exponential curve.  For all other models, the base data is the R2eff/R1rho values for individual spins.
   83 
   84 
   85         @return:    For the R2eff model, a tuple of the spin container and the exponential curve identifying key (the CPMG frequency or R1rho spin-lock field strength).  For all other models, the spin container and ID from the spin loop.
   86         @rtype:     (tuple of SpinContainer instance and float) or (SpinContainer instance and str)
   87         """
   88 
   89         # The R2eff model data (the base data is peak intensities).
   90         if cdp.model_type == MODEL_R2EFF:
   91             # Loop over the sequence.
   92             for spin, spin_id in spin_loop(return_id=True):
   93                 # Skip deselected spins.
   94                 if not spin.select:
   95                     continue
   96 
   97                 # Skip spins with no peak intensity data.
   98                 if not hasattr(spin, 'peak_intensity'):
   99                     continue
  100 
  101                 # Loop over each spectrometer frequency and dispersion point.
  102                 for exp_type, frq, offset, point in loop_exp_frq_offset_point():
  103                     yield spin, spin_id, exp_type, frq, offset, point
  104 
  105         # All other models (the base data is the R2eff/R1rho values).
  106         else:
  107             # 1H MMQ flag.
  108             proton_mmq_flag = has_proton_mmq_cpmg()
  109 
  110             # Loop over the sequence.
  111             for spin, spin_id in spin_loop(return_id=True):
  112                 # Skip deselected spins.
  113                 if not spin.select:
  114                     continue
  115 
  116                 # Skip protons for MMQ data.
  117                 if spin.model in MODEL_LIST_MMQ and spin.isotope == '1H':
  118                     continue
  119 
  120                 # Get the attached proton.
  121                 proton = None
  122                 if proton_mmq_flag:
  123                     proton = return_attached_protons(spin_hash=spin._hash)[0]
  124 
  125                 # Skip spins with no R2eff/R1rho values.
  126                 if not hasattr(spin, 'r2eff') and not hasattr(proton, 'r2eff'):
  127                     continue
  128 
  129                 # Yield the spin container and ID.
  130                 yield spin, spin_id
  131 
  132 
  133     def bmrb_write(self, file_path, version=None):
  134         """Write the model-free results to a BMRB NMR-STAR v3.1 formatted file.
  135 
  136         @param file_path:   The full file path.
  137         @type file_path:    str
  138         @keyword version:   The BMRB NMR-STAR dictionary format to output to.
  139         @type version:      str
  140         """
  141 
  142         # This function is not yet implemented, as it would require a re-write of the relax_data bmrb_write(star) function, and proper handling of cdp.ri_ids.
  143         # It was also not readily possible to find examples of submitted CPMG data in the BMRB database.
  144 
  145         # Not implemented.
  146         raise RelaxImplementError('bmrb_write')
  147 
  148         # Checks.
  149         check_spectrometer_setup(escalate=2)
  150 
  151         # Alias the current data pipe.
  152         cdp = pipes.get_pipe()
  153 
  154         # Initialise the NMR-STAR data object.
  155         star = bmrblib.create_nmr_star('relax_relaxation_dispersion_results', file_path, version)
  156 
  157         # Initialise the spin specific data lists.
  158         mol_name_list = []
  159         res_num_list = []
  160         res_name_list = []
  161         atom_name_list = []
  162 
  163         isotope_list = []
  164         element_list = []
  165 
  166         chi2_list = []
  167         model_list = []
  168 
  169         # Store the spin specific data in lists for later use.
  170         for spin, mol_name, res_num, res_name, spin_id in spin_loop(full_info=True, return_id=True):
  171             # Skip the protons.
  172             if spin.name == 'H' or (hasattr(spin, 'element') and spin.element == 'H'):
  173                 warn(RelaxWarning("Skipping the proton spin '%s'." % spin_id))
  174                 continue
  175 
  176             # Check the data for None (not allowed in BMRB!).
  177             if res_num == None:
  178                 raise RelaxError("For the BMRB, the residue of spin '%s' must be numbered." % spin_id)
  179             if res_name == None:
  180                 raise RelaxError("For the BMRB, the residue of spin '%s' must be named." % spin_id)
  181             if spin.name == None:
  182                 raise RelaxError("For the BMRB, the spin '%s' must be named." % spin_id)
  183             if not hasattr(spin, 'isotope') or spin.isotope == None:
  184                 raise RelaxError("For the BMRB, the spin isotope type of '%s' must be specified." % spin_id)
  185             if not hasattr(spin, 'element') or spin.element == None:
  186                 raise RelaxError("For the BMRB, the spin element type of '%s' must be specified.  Please use the spin user function for setting the element type." % spin_id)
  187 
  188             # The molecule/residue/spin info.
  189             mol_name_list.append(mol_name)
  190             res_num_list.append(res_num)
  191             res_name_list.append(res_name)
  192             atom_name_list.append(spin.name)
  193 
  194             # The nuclear isotope.
  195             if hasattr(spin, 'isotope'):
  196                 isotope_list.append(int(spin.isotope.strip(string.ascii_letters)))
  197             else:
  198                 isotope_list.append(None)
  199 
  200             # The element.
  201             if hasattr(spin, 'element'):
  202                 element_list.append(spin.element)
  203             else:
  204                 element_list.append(None)
  205 
  206             # Opt stats.
  207             if hasattr(spin, 'chi2'):
  208                 chi2_list.append(spin.chi2)
  209             else:
  210                 chi2_list.append(None)
  211 
  212             # Model-free model.
  213             model_list.append(spin.model)
  214 
  215         # Convert the molecule names into the entity IDs.
  216         entity_ids = zeros(len(mol_name_list), int32)
  217         mol_names = get_molecule_names()
  218         for i in range(len(mol_name_list)):
  219             for j in range(len(mol_names)):
  220                 if mol_name_list[i] == mol_names[j]:
  221                     entity_ids[i] = j+1
  222 
  223 
  224         # Create Supergroup 2 : The citations.
  225         ######################################
  226 
  227         # Generate the citations saveframe.
  228         bmrb_write_citations(star)
  229 
  230 
  231         # Create Supergroup 3 : The molecular assembly saveframes.
  232         ##########################################################
  233 
  234         # Generate the entity saveframe.
  235         bmrb_write_entity(star)
  236 
  237 
  238         # Create Supergroup 4:  The experimental descriptions saveframes.
  239         #################################################################
  240 
  241         # Generate the method saveframes.
  242         bmrb_write_methods(star)
  243 
  244         # Generate the software saveframe.
  245         software_ids, software_labels = bmrb_write_software(star)
  246 
  247 
  248     def calculate(self, spin_id=None, scaling_matrix=None, verbosity=1, sim_index=None):
  249         """Calculate the model chi-squared value or the R2eff values for fixed time period data.
  250 
  251         @keyword spin_id:           The spin identification string.
  252         @type spin_id:              None or str
  253         @keyword scaling_matrix:    The per-model list of diagonal and square scaling matrices.
  254         @type scaling_matrix:       list of numpy rank-2, float64 array or list of None
  255         @keyword verbosity:         The amount of information to print.  The higher the value, the greater the verbosity.
  256         @type verbosity:            int
  257         @keyword sim_index:         The MC simulation index (unused).
  258         @type sim_index:            None
  259         """
  260 
  261         # Data checks.
  262         check_pipe()
  263         check_mol_res_spin_data()
  264         check_model_type()
  265 
  266         # Special exponential curve-fitting for the R2eff model.
  267         if cdp.model_type == MODEL_R2EFF:
  268             calculate_r2eff()
  269 
  270         # Calculate the chi-squared value.
  271         else:
  272             # 1H MMQ flag.
  273             proton_mmq_flag = has_proton_mmq_cpmg()
  274 
  275             # Loop over the spin blocks.
  276             model_index = -1
  277             for spin_ids in self.model_loop():
  278                 # Increment the model index.
  279                 model_index += 1
  280 
  281                 # The spin containers.
  282                 spins = spin_ids_to_containers(spin_ids)
  283 
  284                 # Skip deselected clusters.
  285                 skip = True
  286                 for spin in spins:
  287                     if spin.select:
  288                         skip = False
  289                 if skip:
  290                     continue
  291 
  292                 # The back calculated values.
  293                 back_calc = back_calc_r2eff(spins=spins, spin_ids=spin_ids, store_chi2=True)
  294 
  295                 # Pack the data.
  296                 for i, spin_id in enumerate(spin_ids):
  297                     spin = spins[i]
  298                     pack_back_calc_r2eff(spin=spin, spin_id=spin_id, si=i, back_calc=back_calc, proton_mmq_flag=proton_mmq_flag)
  299 
  300 
  301     def constraint_algorithm(self):
  302         """Return the 'Log barrier' optimisation constraint algorithm.
  303 
  304         @return:    The 'Log barrier' constraint algorithm.
  305         @rtype:     str
  306         """
  307 
  308         # The log barrier algorithm, as required by minfx.
  309         return 'Log barrier'
  310 
  311 
  312     def create_mc_data(self, data_id):
  313         """Create the Monte Carlo peak intensity data.
  314 
  315         @param data_id: The tuple of the spin container and the exponential curve identifying key, as yielded by the base_data_loop() generator method.
  316         @type data_id:  SpinContainer instance and float
  317         @return:        The Monte Carlo simulation data.
  318         @rtype:         list of floats
  319         """
  320 
  321         # The R2eff model (with peak intensity base data).
  322         if cdp.model_type == MODEL_R2EFF:
  323             # Unpack the data.
  324             spin, spin_id, exp_type, frq, offset, point = data_id
  325 
  326             # Back calculate the peak intensities.
  327             values = back_calc_peak_intensities(spin=spin, spin_id=spin_id, exp_type=exp_type, frq=frq, offset=offset, point=point)
  328 
  329         # All other models (with R2eff/R1rho base data).
  330         else:
  331             # 1H MMQ flag.
  332             proton_mmq_flag = has_proton_mmq_cpmg()
  333 
  334             # Unpack the data.
  335             spin, spin_id = data_id
  336 
  337             # Back calculate the R2eff/R1rho data.
  338             back_calc = back_calc_r2eff(spins=[spin], spin_ids=[spin_id])
  339 
  340             # Get the attached proton data.
  341             if proton_mmq_flag:
  342                 proton = return_attached_protons(spin_hash=spin._hash)[0]
  343                 proton_back_calc = back_calc_r2eff(spins=[proton], spin_ids=[spin_id])
  344 
  345             # Convert to a dictionary matching the R2eff data structure.
  346             values = {}
  347             for exp_type, frq, offset, point, ei, mi, oi, di in loop_exp_frq_offset_point(return_indices=True):
  348                 # Alias the correct data.
  349                 current_bc = back_calc
  350                 current_spin = spin
  351                 if exp_type in [EXP_TYPE_CPMG_PROTON_SQ, EXP_TYPE_CPMG_PROTON_MQ]:
  352                     current_spin = proton
  353                     current_bc = proton_back_calc
  354 
  355                 # The parameter key.
  356                 param_key = return_param_key_from_data(exp_type=exp_type, frq=frq, offset=offset, point=point)
  357 
  358                 # Skip missing data.
  359                 if not hasattr(current_spin, 'r2eff') or param_key not in current_spin.r2eff:
  360                     continue
  361 
  362                 # Store the result.
  363                 values[param_key] = back_calc[ei][0][mi][oi][di]
  364 
  365         # Return the MC data.
  366         return values
  367 
  368 
  369     def deselect(self, sim_index=None, model_info=None):
  370         """Deselect models or simulations.
  371 
  372         @keyword sim_index:     The optional Monte Carlo simulation index.  If None, then models will be deselected, otherwise the given simulation will.
  373         @type sim_index:        None or int
  374         @keyword model_info:    The list of spins and spin IDs per cluster originating from model_loop().
  375         @type model_info:       list of SpinContainer instances, list of str
  376         """
  377 
  378         # Loop over all the spins and deselect them.
  379         for spin_id in model_info:
  380             # Get the spin.
  381             spin = return_spin(spin_id=spin_id)
  382 
  383             # Spin deselection.
  384             if sim_index == None:
  385                 spin.select = False
  386 
  387             # Simulation deselection.
  388             else:
  389                 spin.select_sim[sim_index] = False
  390 
  391 
  392     def duplicate_data(self, pipe_from=None, pipe_to=None, model_info=None, global_stats=False, verbose=True):
  393         """Duplicate the data specific to a single model.
  394 
  395         @keyword pipe_from:     The data pipe to copy the data from.
  396         @type pipe_from:        str
  397         @keyword pipe_to:       The data pipe to copy the data to.
  398         @type pipe_to:          str
  399         @keyword model_info:    The list of spins and spin IDs per cluster originating from model_loop().
  400         @type model_info:       list of SpinContainer instances, list of str
  401         @keyword global_stats:  The global statistics flag.
  402         @type global_stats:     bool
  403         @keyword verbose:       A flag which if True will cause info to be printed out.
  404         @type verbose:          bool
  405         """
  406 
  407         # First create the pipe_to data pipe, if it doesn't exist, but don't switch to it.
  408         if not pipes.has_pipe(pipe_to):
  409             pipes.create(pipe_to, pipe_type='relax_disp', switch=False)
  410 
  411         # Get the data pipes.
  412         dp_from = pipes.get_pipe(pipe_from)
  413         dp_to = pipes.get_pipe(pipe_to)
  414 
  415         # Duplicate all non-sequence specific data.
  416         for data_name in dir(dp_from):
  417             # Skip the container objects.
  418             if data_name in ['mol', 'interatomic', 'structure', 'exp_info', 'result_files']:
  419                 continue
  420 
  421             # Skip dispersion specific parameters.
  422             if data_name in ['model']:
  423                 continue
  424 
  425             # Skip special objects.
  426             if search('^_', data_name) or data_name in dp_from.__class__.__dict__:
  427                 continue
  428 
  429             # Get the original object.
  430             data_from = getattr(dp_from, data_name)
  431 
  432             # The data already exists.
  433             if hasattr(dp_to, data_name):
  434                 # Get the object in the target pipe.
  435                 data_to = getattr(dp_to, data_name)
  436 
  437                 # The data must match!
  438                 if data_from != data_to:
  439                     raise RelaxError("The object " + repr(data_name) + " is not consistent between the pipes " + repr(pipe_from) + " and " + repr(pipe_to) + ".")
  440 
  441                 # Skip the data.
  442                 continue
  443 
  444             # Duplicate the data.
  445             setattr(dp_to, data_name, deepcopy(data_from))
  446 
  447         # No sequence data, so skip the rest.
  448         if dp_from.mol.is_empty():
  449             return
  450 
  451         # Duplicate the sequence data if it doesn't exist.
  452         if dp_to.mol.is_empty():
  453             sequence.copy(pipe_from=pipe_from, pipe_to=pipe_to, preserve_select=True, verbose=verbose)
  454 
  455         # Loop over the cluster.
  456         for id in model_info:
  457             # The original spin container.
  458             spin = return_spin(spin_id=id, pipe=pipe_from)
  459 
  460             # Duplicate the spin specific data.
  461             for name in dir(spin):
  462                 # Skip special objects.
  463                 if search('^__', name):
  464                     continue
  465 
  466                 # Get the object.
  467                 obj = getattr(spin, name)
  468 
  469                 # Skip methods.
  470                 if isinstance(obj, MethodType):
  471                     continue
  472 
  473                 # Duplicate the object.
  474                 new_obj = deepcopy(getattr(spin, name))
  475                 setattr(dp_to.mol[spin._mol_index].res[spin._res_index].spin[spin._spin_index], name, new_obj)
  476 
  477 
  478     def eliminate(self, name, value, args, sim=None, model_info=None):
  479         """Relaxation dispersion model elimination, parameter by parameter.
  480 
  481         @param name:            The parameter name.
  482         @type name:             str
  483         @param value:           The parameter value.
  484         @type value:            float
  485         @param args:            The c1 and c2 elimination constant overrides.
  486         @type args:             None or tuple of float
  487         @keyword sim:           The Monte Carlo simulation index.
  488         @type sim:              int
  489         @keyword model_info:    The list of spins and spin IDs per cluster originating from model_loop().
  490         @type model_info:       list of SpinContainer instances, list of str
  491         @return:                True if the model is to be eliminated, False otherwise.
  492         @rtype:                 bool
  493         """
  494 
  495         # Skip the R2eff model parameters.
  496         if name in ['r2eff', 'i0']:
  497             return False
  498 
  499         # Default limits.
  500         c1 = 0.501
  501         c2 = 0.999
  502         c3 = 1.0
  503 
  504         # Depack the arguments.
  505         if args != None:
  506             c1, c2, c3 = args
  507 
  508         # Elimination text.
  509         elim_text = "Data pipe '%s':  The %s parameter of %.5f is %s, eliminating " % (pipes.cdp_name(), name, value, "%s")
  510         if sim == None:
  511             elim_text += "the spin cluster %s." % model_info
  512         else:
  513             elim_text += "simulation %i of the spin cluster %s." % (sim, model_info)
  514 
  515         # The pA parameter.
  516         if name == 'pA':
  517             if value < c1:
  518                 print(elim_text % ("less than  %.5f" % c1))
  519                 return True
  520             if value > c2:
  521                 print(elim_text % ("greater than  %.5f" % c2))
  522                 return True
  523 
  524         # The tex parameter.
  525         if name == 'tex':
  526             if value > c3:
  527                 print(elim_text % ("greater than  %.5f" % c3))
  528                 return True
  529 
  530         # Accept model.
  531         return False
  532 
  533 
  534     def get_param_names(self, model_info=None):
  535         """Return a vector of parameter names.
  536 
  537         @keyword model_info:    The list of spins and spin IDs per cluster originating from model_loop().
  538         @type model_info:       list of SpinContainer instances, list of str
  539         @return:                The vector of parameter names.
  540         @rtype:                 list of str
  541         """
  542 
  543         # Get the spin containers.
  544         spins = []
  545         for spin_id in model_info:
  546             # Get the spin.
  547             spin = return_spin(spin_id=spin_id)
  548 
  549             # Skip deselected spins.
  550             if not spin.select:
  551                 continue
  552 
  553             # Add the spin.
  554             spins.append(spin)
  555 
  556         # No spins.
  557         if not len(spins):
  558             return None
  559 
  560         # Call the get_param_names() function.
  561         return get_param_names(spins)
  562 
  563 
  564     def get_param_values(self, model_info=None, sim_index=None):
  565         """Return a vector of parameter values.
  566 
  567         @keyword model_info:    The list of spins and spin IDs per cluster originating from model_loop().
  568         @type model_info:       list of SpinContainer instances, list of str
  569         @keyword sim_index:     The Monte Carlo simulation index.
  570         @type sim_index:        int
  571         @return:                The vector of parameter values.
  572         @rtype:                 list of str
  573         """
  574 
  575         # Get the spin containers.
  576         spins = []
  577         for spin_id in model_info:
  578             # Get the spin.
  579             spin = return_spin(spin_id=spin_id)
  580 
  581             # Skip deselected spins.
  582             if not spin.select:
  583                 continue
  584 
  585             # Add the spin.
  586             spins.append(spin)
  587 
  588         # No spins.
  589         if not len(spins):
  590             return None
  591 
  592         # Set up the R1 parameter, if needed.
  593         r1_setup()
  594 
  595         # Loop over the parameters of the cluster, fetching their values.
  596         values = []
  597         for param_name, param_index, si, r20_key in loop_parameters(spins=spins):
  598             values.append(get_value(spins=spins, sim_index=sim_index, param_name=param_name, spin_index=si, r20_key=r20_key))
  599 
  600         # Return all values.
  601         return values
  602 
  603 
  604     def grid_search(self, lower=None, upper=None, inc=None, scaling_matrix=None, constraints=True, verbosity=1, sim_index=None):
  605         """The relaxation dispersion curve fitting grid search function.
  606 
  607         @keyword lower:             The per-model lower bounds of the grid search which must be equal to the number of parameters in the model.
  608         @type lower:                list of list of numbers
  609         @keyword upper:             The per-model upper bounds of the grid search which must be equal to the number of parameters in the model.
  610         @type upper:                list of list of numbers
  611         @keyword inc:               The per-model increments for each dimension of the space for the grid search. The number of elements in the array must equal to the number of parameters in the model.
  612         @type inc:                  list of list of int
  613         @keyword scaling_matrix:    The per-model list of diagonal and square scaling matrices.
  614         @type scaling_matrix:       list of numpy rank-2, float64 array or list of None
  615         @keyword constraints:       If True, constraints are applied during the grid search (eliminating parts of the grid).  If False, no constraints are used.
  616         @type constraints:          bool
  617         @keyword verbosity:         A flag specifying the amount of information to print.  The higher the value, the greater the verbosity.
  618         @type verbosity:            int
  619         @keyword sim_index:         The index of the simulation to apply the grid search to.  If None, the normal model is optimised.
  620         @type sim_index:            int
  621         """
  622 
  623         # Minimisation.
  624         self.minimise(min_algor='grid', lower=lower, upper=upper, inc=inc, scaling_matrix=scaling_matrix, constraints=constraints, verbosity=verbosity, sim_index=sim_index)
  625 
  626 
  627     def map_bounds(self, param, spin_id=None):
  628         """Create bounds for the OpenDX mapping function.
  629 
  630         @param param:       The name of the parameter to return the lower and upper bounds of.
  631         @type param:        str
  632         @param spin_id:     The spin identification string (unused).
  633         @type spin_id:      None
  634         @return:            The upper and lower bounds of the parameter.
  635         @rtype:             list of float
  636         """
  637 
  638         # Is the parameter is valid?
  639         if not self._PARAMS.contains(param):
  640             raise RelaxError("The parameter '%s' is not valid for this data pipe type." % param)
  641 
  642         # Return the spin.
  643         spin = return_spin(spin_id=spin_id)
  644 
  645         # Loop over each spectrometer frequency and dispersion point to collect param_keys.
  646         param_keys = []
  647         for exp_type, frq, offset, point in loop_exp_frq_offset_point():
  648             # The parameter key.
  649             param_key = return_param_key_from_data(exp_type=exp_type, frq=frq, offset=offset, point=point)
  650 
  651             # Collect the key.
  652             param_keys.append(param_key)
  653 
  654         # The initial parameter vector.
  655         param_vector = []
  656 
  657         # Collect param_names.
  658         param_names = []
  659         for param_name, param_index, si, r20_key in loop_parameters(spins=[spin]):
  660             # Add to the param vector.
  661             param_vector.append([0.0])
  662 
  663             # Collect parameter names.
  664             param_names.append(param_name)
  665 
  666         # Loop over the parameter names.
  667         for i in range(len(param_names)):
  668             # Test if the parameter is in the list:
  669 
  670             if param_names[i] == param:
  671                 return [self._PARAMS.grid_lower(param, incs=0, model_info=[spin_id]), self._PARAMS.grid_upper(param, incs=0, model_info=[spin_id])]
  672 
  673 
  674     def minimise(self, min_algor=None, min_options=None, func_tol=None, grad_tol=None, max_iterations=None, constraints=False, scaling_matrix=None, verbosity=0, sim_index=None, lower=None, upper=None, inc=None):
  675         """Relaxation dispersion curve fitting function.
  676 
  677         @keyword min_algor:         The minimisation algorithm to use.
  678         @type min_algor:            str
  679         @keyword min_options:       An array of options to be used by the minimisation algorithm.
  680         @type min_options:          array of str
  681         @keyword func_tol:          The function tolerance which, when reached, terminates optimisation.  Setting this to None turns of the check.
  682         @type func_tol:             None or float
  683         @keyword grad_tol:          The gradient tolerance which, when reached, terminates optimisation.  Setting this to None turns of the check.
  684         @type grad_tol:             None or float
  685         @keyword max_iterations:    The maximum number of iterations for the algorithm.
  686         @type max_iterations:       int
  687         @keyword constraints:       If True, constraints are used during optimisation.
  688         @type constraints:          bool
  689         @keyword scaling_matrix:    The per-model list of diagonal and square scaling matrices.
  690         @type scaling_matrix:       list of numpy rank-2, float64 array or list of None
  691         @keyword verbosity:         The amount of information to print.  The higher the value, the greater the verbosity.
  692         @type verbosity:            int
  693         @keyword sim_index:         The index of the simulation to optimise.  This should be None if normal optimisation is desired.
  694         @type sim_index:            None or int
  695         @keyword lower:             The per-model lower bounds of the grid search which must be equal to the number of parameters in the model.  This optional argument is only used when doing a grid search.
  696         @type lower:                list of list of numbers
  697         @keyword upper:             The per-model upper bounds of the grid search which must be equal to the number of parameters in the model.  This optional argument is only used when doing a grid search.
  698         @type upper:                list of list of numbers
  699         @keyword inc:               The per-model increments for each dimension of the space for the grid search. The number of elements in the array must equal to the number of parameters in the model.  This argument is only used when doing a grid search.
  700         @type inc:                  list of list of int
  701         """
  702 
  703         # Data checks.
  704         check_mol_res_spin_data()
  705         check_model_type()
  706 
  707         # Check the optimisation algorithm.
  708         algor = min_algor
  709         if min_algor == 'Log barrier':
  710             algor = min_options[0]
  711 
  712         allow = False
  713         # Check the model type.
  714         if hasattr(cdp, 'model_type'):
  715             # Set the model type:
  716             model_type = cdp.model_type
  717 
  718             if model_type == MODEL_R2EFF:
  719                 if match('^[Gg]rid$', algor):
  720                     allow = True
  721 
  722                 elif match('^[Ss]implex$', algor):
  723                     allow = True
  724 
  725                 # Quasi-Newton BFGS minimisation.
  726                 elif match('^[Bb][Ff][Gg][Ss]$', algor):
  727                     allow = True
  728 
  729                 # Newton minimisation.
  730                 elif match('^[Nn]ewton$', algor):
  731                     allow = True
  732 
  733                 # Newton minimisation.
  734                 elif match('^[Nn]ewton$', algor):
  735                     allow = True
  736 
  737                 # Constrained method, Method of Multipliers.
  738                 elif match('^[Mm][Oo][Mm]$', algor) or match('[Mm]ethod of [Mm]ultipliers$', algor):
  739                     allow = True
  740                     
  741                 # Constrained method, Logarithmic barrier function.
  742                 elif match('^[Ll]og [Bb]arrier$', algor):
  743                     allow = True
  744 
  745             # If the Jacobian and Hessian matrix have not been specified for fitting, 'simplex' should be used.
  746             else:
  747                 if match('^[Gg]rid$', algor):
  748                     allow = True
  749 
  750                 elif match('^[Ss]implex$', algor):
  751                     allow = True
  752 
  753         # Do not allow, if no model has been specified.
  754         else:
  755             model_type = 'None'
  756             # Do not allow.
  757             allow = False
  758 
  759         if not allow:
  760             raise RelaxError("Minimisation algorithm '%s' is not allowed, since function gradients for model '%s' is not implemented.  Only the 'simplex' minimisation algorithm is supported for the relaxation dispersion analysis of this model."%(algor, model_type))
  761 
  762         # Initialise some empty data pipe structures so that the target function set up does not fail.
  763         if not hasattr(cdp, 'cpmg_frqs_list'):
  764             cdp.cpmg_frqs_list = []
  765         if not hasattr(cdp, 'spin_lock_nu1_list'):
  766             cdp.spin_lock_nu1_list = []
  767 
  768         # Get the Processor box singleton (it contains the Processor instance) and alias the Processor.
  769         processor_box = Processor_box() 
  770         processor = processor_box.processor
  771 
  772         # The number of time points for the exponential curves (if present).
  773         num_time_pts = 1
  774         if hasattr(cdp, 'num_time_pts'):
  775             num_time_pts = cdp.num_time_pts
  776 
  777         # Number of spectrometer fields.
  778         fields = [None]
  779         field_count = 1
  780         if hasattr(cdp, 'spectrometer_frq'):
  781             fields = cdp.spectrometer_frq_list
  782             field_count = cdp.spectrometer_frq_count
  783 
  784         # Loop over the spin blocks.
  785         model_index = -1
  786         for spin_ids in self.model_loop():
  787             # Increment the model index.
  788             model_index += 1
  789 
  790             # The spin containers.
  791             spins = spin_ids_to_containers(spin_ids)
  792 
  793             # Skip deselected clusters.
  794             skip = True
  795             for spin in spins:
  796                 if spin.select:
  797                     skip = False
  798             if skip:
  799                 continue
  800 
  801             # Alias the grid options.
  802             lower_i, upper_i, inc_i = None, None, None
  803             if min_algor == 'grid':
  804                 lower_i = lower[model_index]
  805                 upper_i = upper[model_index]
  806                 inc_i = inc[model_index]
  807 
  808             # Special exponential curve-fitting for the R2eff model.
  809             if cdp.model_type == MODEL_R2EFF:
  810                 # Sanity checks.
  811                 if not has_exponential_exp_type():
  812                     raise RelaxError("The R2eff model with the fixed time period dispersion experiments cannot be optimised.")
  813 
  814                 # Optimisation.
  815                 minimise_r2eff(spins=spins, spin_ids=spin_ids, min_algor=min_algor, min_options=min_options, func_tol=func_tol, grad_tol=grad_tol, max_iterations=max_iterations, constraints=constraints, scaling_matrix=scaling_matrix[model_index], verbosity=verbosity, sim_index=sim_index, lower=lower_i, upper=upper_i, inc=inc_i)
  816 
  817                 # Skip the rest.
  818                 continue
  819 
  820             # Set up the slave command object.
  821             command = Disp_minimise_command(spins=spins, spin_ids=spin_ids, sim_index=sim_index, scaling_matrix=scaling_matrix[model_index], min_algor=min_algor, min_options=min_options, func_tol=func_tol, grad_tol=grad_tol, max_iterations=max_iterations, constraints=constraints, verbosity=verbosity, lower=lower_i, upper=upper_i, inc=inc_i, fields=fields, param_names=get_param_names(spins=spins, full=True))
  822 
  823             # Set up the memo.
  824             memo = Disp_memo(spins=spins, spin_ids=spin_ids, sim_index=sim_index, scaling_matrix=scaling_matrix[model_index], verbosity=verbosity)
  825 
  826             # Add the slave command and memo to the processor queue.
  827             processor.add_to_queue(command, memo)
  828 
  829 
  830     def model_desc(self, model_info=None):
  831         """Return a description of the model.
  832 
  833         @keyword model_info:    The list of spins and spin IDs per cluster originating from model_loop().
  834         @type model_info:       list of SpinContainer instances, list of str
  835         @return:                The model description.
  836         @rtype:                 str
  837         """
  838 
  839         # The model loop is over the spin clusters, so return a description of the cluster.
  840         return "The spin cluster %s." % model_info
  841 
  842 
  843     def model_loop(self):
  844         """Loop over the spin groupings for one model applied to multiple spins.
  845 
  846         @return:    The list of spins per block will be yielded, as well as the list of spin IDs.
  847         @rtype:     tuple of list of SpinContainer instances and list of str
  848         """
  849 
  850         # Loop over individual spins for the R2eff model.
  851         if not hasattr(cdp, 'model_type') or cdp.model_type == MODEL_R2EFF:
  852             # The spin loop.
  853             for spin, spin_id in spin_loop(return_id=True):
  854                 # Skip deselected spins
  855                 if not spin.select:
  856                     continue
  857 
  858                 # Yield the spin ID as a list.
  859                 yield [spin_id]
  860 
  861          # The cluster loop.
  862         else:
  863             for spin_ids in loop_cluster(skip_desel=False):
  864                 yield spin_ids
  865 
  866 
  867     def model_statistics(self, model_info=None, spin_id=None, global_stats=None):
  868         """Return the k, n, and chi2 model statistics.
  869 
  870         k - number of parameters.
  871         n - number of data points.
  872         chi2 - the chi-squared value.
  873 
  874 
  875         @keyword model_info:    The list of spins and spin IDs per cluster originating from model_loop().
  876         @type model_info:       list of SpinContainer instances, list of str
  877         @keyword spin_id:       The spin ID string to override the model_info argument.  This is ignored in the N-state model.
  878         @type spin_id:          None or str
  879         @keyword global_stats:  A parameter which determines if global or local statistics are returned.  For the N-state model, this argument is ignored.
  880         @type global_stats:     None or bool
  881         @return:                The optimisation statistics, in tuple format, of the number of parameters (k), the number of data points (n), and the chi-squared value (chi2).
  882         @rtype:                 tuple of (int, int, float)
  883         """
  884 
  885         # Get the spin containers (the spin ID overrides the model info).
  886         if spin_id != None:
  887             spins = [return_spin(spin_id=spin_id)]
  888         else:
  889             spins = spin_ids_to_containers(model_info)
  890 
  891         # The number of parameters for the cluster.
  892         k = param_num(spins=spins)
  893 
  894         # The number of points from all spins.
  895         n = 0
  896         for spin in spins:
  897             # Skip deselected spins.
  898             if not spin.select:
  899                 continue
  900 
  901             n += len(spin.r2eff)
  902 
  903         # Take the chi-squared from the first spin of the cluster (which has a value).
  904         chi2 = None
  905         for spin in spins:
  906             # Skip deselected spins.
  907             if not spin.select:
  908                 continue
  909 
  910             if hasattr(spin, 'chi2'):
  911                 chi2 = spin.chi2
  912                 break
  913 
  914         # Return the values.
  915         return k, n, chi2
  916 
  917 
  918     def overfit_deselect(self, data_check=True, verbose=True):
  919         """Deselect spins which have insufficient data to support minimisation.
  920 
  921         @keyword data_check:    A flag to signal if the presence of base data is to be checked for.
  922         @type data_check:       bool
  923         @keyword verbose:       A flag which if True will allow printouts.
  924         @type verbose:          bool
  925         """
  926 
  927         # Test the sequence data exists and the model is setup.
  928         check_mol_res_spin_data()
  929         check_model_type()
  930 
  931         # 1H MMQ flag.
  932         proton_mmq_flag = has_proton_mmq_cpmg()
  933 
  934         # Loop over spin data.
  935         for spin, spin_id in spin_loop(return_id=True, skip_desel=True):
  936             # Skip protons for MMQ data.
  937             if spin.model in MODEL_LIST_MMQ and spin.isotope == '1H':
  938                 continue
  939 
  940             # Get the attached proton.
  941             proton = None
  942             if proton_mmq_flag:
  943                 # Get all protons.
  944                 proton_spins = return_attached_protons(spin_hash=spin._hash)
  945 
  946                 # Only one allowed.
  947                 if len(proton_spins) > 1:
  948                     print("Multiple protons attached to the spin '%s', but one one attached proton is supported for the MMQ-type models." % spin_id)
  949                     spin.select = False
  950                     continue
  951 
  952                 # Alias the single proton.
  953                 if len(proton_spins):
  954                     proton = proton_spins[0]
  955 
  956             # Check if data exists.
  957             if not hasattr(spin, 'r2eff') and not hasattr(proton, 'r2eff'):
  958                 print("No R2eff data could be found, deselecting the '%s' spin." % spin_id)
  959                 spin.select = False
  960                 continue
  961 
  962 
  963     def print_model_title(self, prefix=None, model_info=None):
  964         """Print out the model title.
  965 
  966         @keyword prefix:        The starting text of the title.  This should be printed out first, followed by the model information text.
  967         @type prefix:           str
  968         @keyword model_info:    The list of spins and spin IDs per cluster originating from model_loop().
  969         @type model_info:       list of SpinContainer instances, list of str
  970         """
  971 
  972         # The printout.
  973         subsection(file=sys.stdout, text=prefix+"the spin block %s"%model_info, prespace=2)
  974 
  975 
  976     def return_data(self, data_id=None):
  977         """Return the peak intensity data structure.
  978 
  979         @param data_id: The spin ID string, as yielded by the base_data_loop() generator method.
  980         @type data_id:  str
  981         @return:        The peak intensity data structure.
  982         @rtype:         list of float
  983         """
  984 
  985         # The R2eff model.
  986         if cdp.model_type == MODEL_R2EFF:
  987             # Unpack the data.
  988             spin, spin_id, exp_type, frq, offset, point = data_id
  989 
  990             # Return the data.
  991             return spin.peak_intensity
  992 
  993         # All other models.
  994         else:
  995             raise RelaxImplementError
  996 
  997 
  998     def return_error(self, data_id=None):
  999         """Return the standard deviation data structure.
 1000 
 1001         @param data_id: The tuple of the spin container and the exponential curve identifying key, as yielded by the base_data_loop() generator method.
 1002         @type data_id:  SpinContainer instance and float
 1003         @return:        The standard deviation data structure.
 1004         @rtype:         list of float
 1005         """
 1006 
 1007         # The R2eff model.
 1008         if cdp.model_type == MODEL_R2EFF:
 1009             # Unpack the data.
 1010             spin, spin_id, exp_type, frq, offset, point = data_id
 1011 
 1012             # Generate the data structure to return.
 1013             errors = []
 1014             for time in loop_time(exp_type=exp_type, frq=frq, offset=offset, point=point):
 1015                 errors.append(average_intensity(spin=spin, exp_type=exp_type, frq=frq, offset=offset, point=point, time=time, error=True))
 1016 
 1017         # All other models.
 1018         else:
 1019             # Unpack the data.
 1020             spin, spin_id = data_id
 1021 
 1022             # 1H MMQ flag.
 1023             proton_mmq_flag = has_proton_mmq_cpmg()
 1024 
 1025             # Get the attached proton.
 1026             proton = None
 1027             if proton_mmq_flag:
 1028                 proton = return_attached_protons(spin_hash=spin._hash)[0]
 1029 
 1030             # The errors.
 1031             r2eff_err = {}
 1032             if hasattr(spin, 'r2eff_err'):
 1033                 r2eff_err.update(spin.r2eff_err)
 1034             if hasattr(proton, 'r2eff_err'):
 1035                 r2eff_err.update(proton.r2eff_err)
 1036             return r2eff_err
 1037 
 1038         # Return the error list.
 1039         return errors
 1040 
 1041 
 1042     def return_error_red_chi2(self, data_id=None):
 1043         """Return the standard deviation data structure, where standard deviation is from the overall gauss distribution described by the STD_fit of the goodness of fit, where STD_fit = sqrt(chi2/(N-p))
 1044 
 1045         @param data_id: The tuple of the spin container and the exponential curve identifying key, as yielded by the base_data_loop() generator method.
 1046         @type data_id:  SpinContainer instance and float
 1047         @return:        The standard deviation data structure.
 1048         @rtype:         list of float
 1049         """
 1050 
 1051         # The R2eff model.
 1052         if cdp.model_type == MODEL_R2EFF:
 1053             raise RelaxError("Drawing errors from the gauss distribution described by the STD_fit of the goodness of fit, is not possible for the '%s' model."%MODEL_R2EFF)
 1054 
 1055         # Get the errors structure as above.
 1056         errors = self.return_error(data_id=data_id)
 1057 
 1058         # Unpack the data.
 1059         spin, spin_id = data_id
 1060 
 1061         # Loop over the spin groupings for the model.
 1062         for spin_ids in self.model_loop():
 1063             # If the spin of interest is in the returned spin cluster.
 1064             if spin_id in spin_ids:
 1065                 # Get the statistics
 1066                 k, n, chi2 = self.model_statistics(model_info=spin_ids)
 1067 
 1068                 # Calculate degrees of freedom.
 1069                 dof = n - k
 1070 
 1071                 # Calculate reduced chi2, or named as the variance of the squared residuals.
 1072                 red_chi2 = chi2 / float(dof)
 1073 
 1074                 # Calculated the standard deviation.
 1075                 std_red_chi2 = sqrt(red_chi2)
 1076 
 1077                 # Replace values with the stored value.
 1078                 for id in errors:
 1079                     errors[id] = std_red_chi2
 1080 
 1081                 return errors
 1082 
 1083 
 1084     def return_value(self, spin, param, sim=None, bc=False):
 1085         """Return the value and error corresponding to the parameter.
 1086 
 1087         If sim is set to an integer, return the value of the simulation and None.
 1088 
 1089 
 1090         @param spin:    The SpinContainer object.
 1091         @type spin:     SpinContainer
 1092         @param param:   The name of the parameter to return values for.
 1093         @type param:    str
 1094         @keyword sim:   The Monte Carlo simulation index.
 1095         @type sim:      None or int
 1096         @keyword bc:    The back-calculated data flag.  If True, then the back-calculated data will be returned rather than the actual data.
 1097         @type bc:       bool
 1098         @return:        The value and error corresponding to
 1099         @rtype:         tuple of length 2 of floats or None
 1100         """
 1101 
 1102         # Define list of special parameters.
 1103         special_parameters = ['theta', 'w_eff']
 1104 
 1105         # Use api_common function for paramets not defined in special_parameters.
 1106         if param not in special_parameters:
 1107             returnval = self._return_value_general(spin=spin, param=param, sim=sim, bc=bc)
 1108             return returnval
 1109 
 1110         # If parameter in special_parameters, the do the following.
 1111 
 1112         # Initial values.
 1113         value = None
 1114         error = None
 1115 
 1116         # Return the data
 1117         theta_spin_dic, Domega_spin_dic, w_eff_spin_dic, dic_key_list = calc_rotating_frame_params(spin=spin)
 1118 
 1119         # Return for parameter theta
 1120         if param == "theta":
 1121             value = theta_spin_dic
 1122 
 1123         # Return for parameter theta
 1124         elif param == "w_eff":
 1125             value = w_eff_spin_dic
 1126 
 1127         # Return the data.
 1128         return value, error
 1129 
 1130 
 1131     def set_error(self, index, error, model_info=None):
 1132         """Set the parameter errors.
 1133 
 1134         @param index:           The index of the parameter to set the errors for.
 1135         @type index:            int
 1136         @param error:           The error value.
 1137         @type error:            float
 1138         @keyword model_info:    The list of spins and spin IDs per cluster originating from model_loop().
 1139         @type model_info:       list of SpinContainer instances, list of str
 1140         """
 1141 
 1142         # Unpack the data.
 1143         spin_ids = model_info
 1144         spins = spin_ids_to_containers(spin_ids)
 1145 
 1146         # The number of parameters.
 1147         total_param_num = param_num(spins=spins)
 1148 
 1149         # No more model parameters.
 1150         model_param = True
 1151         if index >= total_param_num:
 1152             model_param = False
 1153 
 1154         # The auxiliary cluster parameters.
 1155         aux_params = []
 1156         if 'pA' in spins[0].params:
 1157             aux_params.append('pB')
 1158         if 'pB' in spins[0].params:
 1159             aux_params.append('pA')
 1160         if 'kex' in spins[0].params:
 1161             aux_params.append('tex')
 1162         if 'tex' in spins[0].params:
 1163             aux_params.append('kex')
 1164 
 1165         # Convert the parameter index.
 1166         if model_param:
 1167             param_name, si, mi = param_index_to_param_info(index=index, spins=spins)
 1168         else:
 1169             param_name = aux_params[index - total_param_num]
 1170             si = None
 1171             mi = None
 1172 
 1173         # The parameter error name.
 1174         err_name = param_name + "_err"
 1175 
 1176         # The exponential curve parameters.
 1177         if param_name in ['r2eff', 'i0']:
 1178             # Initialise if needed.
 1179             if not hasattr(spins[si], err_name):
 1180                 setattr(spins[si], err_name, {})
 1181 
 1182             # Set the value.
 1183             setattr(spins[si], err_name, error)
 1184 
 1185         # Model and auxiliary parameters.
 1186         else:
 1187             # If clustered paramater:
 1188             if si == None:
 1189                 for spin in spins:
 1190                     setattr(spin, err_name, error)
 1191 
 1192             # If independent value.
 1193             else:
 1194                 # Set the value.
 1195                 setattr(spins[si], err_name, error)
 1196 
 1197 
 1198     def set_param_values(self, param=None, value=None, index=None, spin_id=None, error=False, force=True):
 1199         """Set the spin specific parameter values.
 1200 
 1201         @keyword param:     The parameter name list.
 1202         @type param:        list of str
 1203         @keyword value:     The parameter value list.
 1204         @type value:        list
 1205         @keyword index:     The index for parameters which are of the list-type.
 1206         @type index:        None or int
 1207         @keyword spin_id:   The spin identification string, only used for spin specific parameters.
 1208         @type spin_id:      None or str
 1209         @keyword error:     A flag which if True will allow the parameter errors to be set instead of the values.
 1210         @type error:        bool
 1211         @keyword force:     A flag which if True will cause current values to be overwritten.  If False, a RelaxError will raised if the parameter value is already set.
 1212         @type force:        bool
 1213         """
 1214 
 1215         # Checks.
 1216         is_str_list(param, 'parameter name')
 1217         is_list(value, 'parameter value')
 1218 
 1219         # Loop over the spin blocks.
 1220         model_index = -1
 1221         for spin_ids in self.model_loop():
 1222             # Increment the model index.
 1223             model_index += 1
 1224 
 1225             # The spin containers.
 1226             spins = spin_ids_to_containers(spin_ids)
 1227 
 1228             # Skip deselected clusters.
 1229             skip = True
 1230             for spin in spins:
 1231                 if spin.select:
 1232                     skip = False
 1233             if skip:
 1234                 continue
 1235 
 1236             # Loop over the parameters.
 1237             for i in range(len(param)):
 1238                 param_i = param[i]
 1239                 value_i = value[i]
 1240 
 1241                 # Is the parameter is valid?
 1242                 if not self._PARAMS.contains(param_i):
 1243                     raise RelaxError("The parameter '%s' is not valid for this data pipe type." % param_i)
 1244 
 1245                 # If the parameter is a global parameter, then change for all spins part of the cluster.
 1246                 if param_i in ['pA', 'pB', 'pC', 'kex', 'k_AB', 'k_BC', 'k_AC', 'tex', 'kB', 'kC', 'kex_AB', 'kex_BC', 'kex_AC']:
 1247                     selection_list = spin_ids
 1248                 else:
 1249                     selection_list = [spin_id]
 1250 
 1251                 # Now loop over selections in the list.
 1252                 for selection in selection_list:
 1253                     # Spin loop.
 1254                     for spin in spin_loop(selection=selection):
 1255                         # The object name.
 1256                         obj_name = param_i
 1257                         if error:
 1258                             obj_name += '_err'
 1259         
 1260                         # Handle the R10 parameters.
 1261                         if param_i in ['r1']:
 1262                             # Loop over the current keys.
 1263                             for exp_type, frq, ei, mi in loop_exp_frq(return_indices=True):
 1264                                 # The parameter key.
 1265                                 key = generate_r20_key(exp_type=exp_type, frq=frq)
 1266         
 1267                                 # Initialise the structure if needed.
 1268                                 if not hasattr(spin, obj_name):
 1269                                     setattr(spin, obj_name, {})
 1270         
 1271                                 # Set the value.
 1272                                 if index == None:
 1273                                     obj = getattr(spin, obj_name)
 1274                                     obj[key] = value_i
 1275         
 1276                                 # If the index is specified, let it match the frequency index
 1277                                 elif mi == index:
 1278                                     obj = getattr(spin, obj_name)
 1279                                     obj[key] = value_i
 1280         
 1281                         # Handle the R20 parameters.
 1282                         elif param_i in PARAMS_R20:
 1283                             # Loop over the current keys.
 1284                             for exp_type, frq, ei, mi in loop_exp_frq(return_indices=True):
 1285                                 # The parameter key.
 1286                                 key = generate_r20_key(exp_type=exp_type, frq=frq)
 1287         
 1288                                 # Initialise the structure if needed.
 1289                                 if not hasattr(spin, obj_name):
 1290                                     setattr(spin, obj_name, {})
 1291         
 1292                                 # Set the value.
 1293                                 if index == None:
 1294                                     obj = getattr(spin, obj_name)
 1295                                     obj[key] = value_i
 1296         
 1297                                 # If the index is specified, let it match the frequency index
 1298                                 elif mi == index:
 1299                                     obj = getattr(spin, obj_name)
 1300                                     obj[key] = value_i
 1301         
 1302                         # Handle the R2eff and I0 parameters.
 1303                         elif param_i in ['r2eff', 'i0'] and not isinstance(value_i, dict):
 1304                             # Loop over all the data.
 1305                             for exp_type, frq, offset, point in loop_exp_frq_offset_point():
 1306                                 # The parameter key.
 1307                                 key = return_param_key_from_data(exp_type=exp_type, frq=frq, offset=offset, point=point)
 1308         
 1309                                 # Initialise the structure if needed.
 1310                                 if not hasattr(spin, obj_name):
 1311                                     setattr(spin, obj_name, {})
 1312         
 1313                                 # Set the value.
 1314                                 obj = getattr(spin, obj_name)
 1315                                 obj[key] = value_i
 1316         
 1317                         # Set the other parameters.
 1318                         else:
 1319                             setattr(spin, obj_name, value_i)
 1320 
 1321 
 1322     def set_selected_sim(self, select_sim, model_info=None):
 1323         """Set the simulation selection flag.
 1324 
 1325         @param select_sim:      The selection flag for the simulations.
 1326         @type select_sim:       bool
 1327         @keyword model_info:    The list of spins and spin IDs per cluster originating from model_loop().
 1328         @type model_info:       list of SpinContainer instances, list of str
 1329         """
 1330 
 1331         # Unpack the data.
 1332         spin_ids = model_info
 1333         spins = spin_ids_to_containers(spin_ids)
 1334 
 1335         # Loop over the spins, storing the structure for each spin.
 1336         for spin in spins:
 1337             spin.select_sim = deepcopy(select_sim)
 1338 
 1339 
 1340     def sim_init_values(self):
 1341         """Initialise the Monte Carlo parameter values."""
 1342 
 1343         # Get the minimisation statistic object names.
 1344         min_names = self.data_names(set='min')
 1345 
 1346         # Set the Monte Carlo parameter values.
 1347         for spin_ids in self.model_loop():
 1348             spins = spin_ids_to_containers(spin_ids)
 1349 
 1350             # Firstly perform any required parameter conversions.
 1351             param_conversion(key=None, spins=spins, sim_index=None)
 1352 
 1353             # Loop over the spins.
 1354             for spin in spins:
 1355                 # Skip deselected spins.
 1356                 if not spin.select:
 1357                     continue
 1358 
 1359                 # Skip protons for MMQ data.
 1360                 if spin.model in MODEL_LIST_MMQ and spin.isotope == '1H':
 1361                     continue
 1362 
 1363                 # Include all parameter conversions.
 1364                 conversion_names = []
 1365                 if 'pB' in spin.params:
 1366                     conversion_names.append('pC')
 1367                 elif 'pA' in spin.params:
 1368                     conversion_names.append('pB')
 1369                 if 'kex' in spin.params:
 1370                     conversion_names.append('tex')
 1371                 elif 'tex' in spin.params:
 1372                     conversion_names.append('kex')
 1373                 if 'kex' in spin.params and 'pA' in spin.params:
 1374                     conversion_names.append('k_AB')
 1375                     conversion_names.append('k_BA')
 1376 
 1377                 # Loop over all the data names.
 1378                 for object_name in spin.params + min_names + conversion_names:
 1379                     # Name for the simulation object.
 1380                     sim_object_name = object_name + '_sim'
 1381 
 1382                     # Create the simulation object.
 1383                     setattr(spin, sim_object_name, [])
 1384 
 1385                     # Get the simulation object.
 1386                     sim_object = getattr(spin, sim_object_name)
 1387 
 1388                     # Copy and append the data.
 1389                     for i in range(cdp.sim_number):
 1390                         sim_object.append(deepcopy(getattr(spin, object_name)))
 1391 
 1392             # Perform the required simulation parameter conversions.
 1393             for i in range(cdp.sim_number):
 1394                 param_conversion(key=None, spins=spins, sim_index=i)
 1395 
 1396 
 1397     def sim_pack_data(self, data_id, sim_data):
 1398         """Pack the Monte Carlo simulation data.
 1399 
 1400         @param data_id:     The tuple of the spin container and the exponential curve identifying key, as yielded by the base_data_loop() generator method.
 1401         @type data_id:      SpinContainer instance and float
 1402         @param sim_data:    The Monte Carlo simulation data.
 1403         @type sim_data:     list of float
 1404         """
 1405 
 1406         # The R2eff model (with peak intensity base data).
 1407         if cdp.model_type == MODEL_R2EFF:
 1408             # Unpack the data.
 1409             spin, spin_id, exp_type, frq, offset, point = data_id
 1410 
 1411             # Initialise the data structure if needed.
 1412             if not hasattr(spin, 'peak_intensity_sim'):
 1413                 spin.peak_intensity_sim = {}
 1414 
 1415             # Loop over each time point.
 1416             ti = 0
 1417             for time in loop_time(exp_type=exp_type, frq=frq, offset=offset, point=point):
 1418                 # Get the intensity keys.
 1419                 int_keys = find_intensity_keys(exp_type=exp_type, frq=frq, offset=offset, point=point, time=time)
 1420 
 1421                 # Loop over the intensity keys.
 1422                 for int_key in int_keys:
 1423                     # Test if the simulation data point already exists.
 1424                     if int_key in spin.peak_intensity_sim:
 1425                         raise RelaxError("Monte Carlo simulation data for the key '%s' already exists." % int_key)
 1426 
 1427                     # Initialise the list.
 1428                     spin.peak_intensity_sim[int_key] = []
 1429 
 1430                     # Loop over the simulations, appending the corresponding data.
 1431                     for i in range(cdp.sim_number):
 1432                         spin.peak_intensity_sim[int_key].append(sim_data[i][ti])
 1433 
 1434                 # Increment the time index.
 1435                 ti += 1
 1436 
 1437         # All other models (with R2eff/R1rho base data).
 1438         else:
 1439             # Unpack the data.
 1440             spin, spin_id = data_id
 1441 
 1442             # 1H MMQ flag.
 1443             proton_mmq_flag = has_proton_mmq_cpmg()
 1444 
 1445             # Get the attached proton.
 1446             proton = None
 1447             if proton_mmq_flag:
 1448                 proton = return_attached_protons(spin_hash=spin._hash)[0]
 1449 
 1450             # Pack the data.
 1451             spin.r2eff_sim = sim_data
 1452             if proton != None:
 1453                 proton.r2eff_sim = sim_data
 1454 
 1455 
 1456     def sim_return_param(self, index, model_info=None):
 1457         """Return the array of simulation parameter values.
 1458 
 1459         @param index:           The index of the parameter to return the array of values for.
 1460         @type index:            int
 1461         @keyword model_info:    The list of spins and spin IDs per cluster originating from model_loop().
 1462         @type model_info:       list of SpinContainer instances, list of str
 1463         @return:                The array of simulation parameter values.
 1464         @rtype:                 list of float
 1465         """
 1466 
 1467         # Unpack the data.
 1468         spin_ids = model_info
 1469         spins = spin_ids_to_containers(spin_ids)
 1470 
 1471         # The number of parameters.
 1472         total_param_num = param_num(spins=spins)
 1473 
 1474         # No more model parameters.
 1475         model_param = True
 1476         if index >= total_param_num:
 1477             model_param = False
 1478 
 1479         # The auxiliary cluster parameters.
 1480         aux_params = []
 1481         for spin in spins:
 1482             if not spin.select:
 1483                 continue
 1484             if 'pA' in spin.params:
 1485                 aux_params.append('pB')
 1486             if 'pB' in spin.params:
 1487                 aux_params.append('pA')
 1488             if 'kex' in spin.params:
 1489                 aux_params.append('tex')
 1490             if 'tex' in spin.params:
 1491                 aux_params.append('kex')
 1492             break
 1493 
 1494         # No more auxiliary parameters.
 1495         total_aux_num = total_param_num + len(aux_params)
 1496         if index >= total_aux_num:
 1497             return
 1498 
 1499         # Convert the parameter index.
 1500         if model_param:
 1501             param_name, si, mi = param_index_to_param_info(index=index, spins=spins)
 1502             if not param_name in ['r2eff', 'i0']:
 1503                 if si == None:
 1504                     si = 0
 1505 
 1506         else:
 1507             param_name = aux_params[index - total_param_num]
 1508             si = 0
 1509             mi = 0
 1510 
 1511         # The exponential curve parameters.
 1512         sim_data = []
 1513         if param_name == 'r2eff':
 1514             for i in range(cdp.sim_number):
 1515                 sim_data.append(spins[si].r2eff_sim[i])
 1516         elif param_name == 'i0':
 1517             for i in range(cdp.sim_number):
 1518                 sim_data.append(spins[si].i0_sim[i])
 1519 
 1520         # Model and auxiliary parameters.
 1521         else:
 1522             sim_data = getattr(spins[si], param_name + "_sim")
 1523 
 1524         # Set the sim data to None if empty.
 1525         if sim_data == []:
 1526             sim_data = None
 1527 
 1528         # Return the object.
 1529         return sim_data
 1530 
 1531 
 1532     def sim_return_selected(self, model_info=None):
 1533         """Return the array of selected simulation flags.
 1534 
 1535         @keyword model_info:    The list of spins and spin IDs per cluster originating from model_loop().
 1536         @type model_info:       list of SpinContainer instances, list of str
 1537         @return:                The array of selected simulation flags.
 1538         @rtype:                 list of int
 1539         """
 1540 
 1541         # Unpack the data.
 1542         spin_ids = model_info
 1543         spins = spin_ids_to_containers(spin_ids)
 1544 
 1545         # Return the array from the first spin, as this array will be identical for all spins in the cluster.
 1546         return spins[0].select_sim