specific_analyses.frame_order.parameters

1 ############################################################################### 2 # # 3 # Copyright (C) 2009-2011,2013-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 handling the frame order model parameters.""" 24 25 # Python module imports. 26 from math import pi 27 from numpy import array, float64, zeros 28 29 # relax module imports. 30 from lib.errors import RelaxError 31 from lib.frame_order.variables import MODEL_DOUBLE_ROTOR, MODEL_FREE_ROTOR, MODEL_ISO_CONE, MODEL_ISO_CONE_FREE_ROTOR, MODEL_ISO_CONE_TORSIONLESS, MODEL_LIST_FREE_ROTORS, MODEL_LIST_ISO_CONE, MODEL_LIST_PSEUDO_ELLIPSE, MODEL_PSEUDO_ELLIPSE, MODEL_ROTOR 32 from lib.list import unique_elements 33 from specific_analyses.frame_order.data import pivot_fixed 34 35

36 -def assemble_param_vector(sim_index=None, unset_fail=False):

37 """Assemble and return the parameter vector. 38 39 @keyword sim_index: The Monte Carlo simulation index. 40 @type sim_index: int 41 @keyword unset_fail: A flag which if True will cause a RelaxError to be raised if the parameter is not set yet. 42 @type unset_fail: bool 43 @return: The parameter vector. 44 @rtype: numpy rank-1 array 45 """ 46 47 # Initialise. 48 param_vect = [] 49 50 # Parameter name extension. 51 ext = '' 52 if sim_index != None: 53 ext = '_sim' 54 55 # Loop over all model parameters. 56 for param_name in cdp.params: 57 # Add the extension. 58 param_name += ext 59 60 # The parameter does not exist yet. 61 if not hasattr(cdp, param_name): 62 if unset_fail: 63 raise RelaxError("The parameter '%s' has not been set." % param_name) 64 else: 65 param_vect.append(None) 66 continue 67 68 # Get the object. 69 obj = getattr(cdp, param_name) 70 71 # Add it to the parameter vector. 72 if sim_index == None: 73 param_vect.append(obj) 74 else: 75 param_vect.append(obj[sim_index]) 76 77 # Return as a numpy array. 78 return array(param_vect, float64)

79 80

81 -def linear_constraints(scaling_matrix=None):

82 """Create the linear constraint matrices A and b. 83 84 Standard notation 85 ================= 86 87 The parameter constraints for the motional amplitude parameters are:: 88 89 0 <= theta <= pi, 90 0 <= theta_x <= theta_y <= pi, 91 -0.125 <= S <= 1, 92 0 <= sigma_max <= pi, 93 94 The pivot point and average domain position parameter constraints are:: 95 96 -999 <= pivot_x <= 999 97 -999 <= pivot_y <= 999 98 -999 <= pivot_z <= 999 99 -500 <= ave_pos_x <= 500 100 -500 <= ave_pos_y <= 500 101 -500 <= ave_pos_y <= 500 102 103 These are necessary to allow for valid PDB representations to be created. The eigenframe parameters are unconstrained. 104 105 106 Matrix notation 107 =============== 108 109 In the notation A.x >= b, where A is an matrix of coefficients, x is an array of parameter values, and b is a vector of scalars, these inequality constraints are:: 110 111 | 1 0 0 0 | | 0 | 112 | | | | 113 |-1 0 0 0 | | -pi | 114 | | | | 115 | 0 1 0 0 | | 0 | 116 | | | theta | | | 117 | 0 -1 0 0 | | | | -pi | 118 | | | theta_x | | | 119 | 0 -1 1 0 | . | | >= | 0 | 120 | | | theta_y | | | 121 | 0 0 1 0 | | | | 0 | 122 | | | sigma_max | | | 123 | 0 0 -1 0 | | -pi | 124 | | | | 125 | 0 0 0 1 | | 0 | 126 | | | | 127 | 0 0 0 -1 | | -pi | 128 129 The pivot and average position constraints in the A.x >= b notation are:: 130 131 | 1 0 0 0 0 0 | | -999.0 | 132 | | | | 133 |-1 0 0 0 0 0 | | -999.0 | 134 | | | | 135 | 0 1 0 0 0 0 | | -999.0 | 136 | | | | 137 | 0 -1 0 0 0 0 | | pivot_x | | -999.0 | 138 | | | | | | 139 | 0 0 1 0 0 0 | | pivot_y | | -999.0 | 140 | | | | | | 141 | 0 0 -1 0 0 0 | | pivot_z | | -999.0 | 142 | | . | | >= | | 143 | 0 0 0 1 0 0 | | ave_pos_x | | -500.0 | 144 | | | | | | 145 | 0 0 0 -1 0 0 | | ave_pos_y | | -500.0 | 146 | | | | | | 147 | 0 0 0 0 1 0 | | ave_pos_z | | -500.0 | 148 | | | | 149 | 0 0 0 0 -1 0 | | -500.0 | 150 | | | | 151 | 0 0 0 0 0 1 | | -500.0 | 152 | | | | 153 | 0 0 0 0 0 -1 | | -500.0 | 154 155 156 @keyword scaling_matrix: The diagonal, square scaling matrix. 157 @type scaling_matrix: numpy rank-2 square matrix 158 @return: The matrices A and b. 159 @rtype: numpy rank-2 NxM array, numpy rank-1 N array 160 """ 161 162 # Initialisation (0..j..m). 163 A = [] 164 b = [] 165 n = param_num() 166 zero_array = zeros(n, float64) 167 i = 0 168 j = 0 169 170 # Loop over the parameters of the model. 171 for i in range(n): 172 # The pivot parameters. 173 if cdp.params[i] in ['pivot_x', 'pivot_y', 'pivot_z']: 174 # -999 <= pivot_i <= 999. 175 A.append(zero_array * 0.0) 176 A.append(zero_array * 0.0) 177 A[j][i] = 1.0 178 A[j+1][i] = -1.0 179 b.append(-999.0 / scaling_matrix[i, i]) 180 b.append(-999.0 / scaling_matrix[i, i]) 181 j = j + 2 182 183 # The average domain translation parameters. 184 if cdp.params[i] in ['ave_pos_x', 'ave_pos_y', 'ave_pos_z']: 185 # -500 <= ave_pos_i <= 500. 186 A.append(zero_array * 0.0) 187 A.append(zero_array * 0.0) 188 A[j][i] = 1.0 189 A[j+1][i] = -1.0 190 b.append(-500.0 / scaling_matrix[i, i]) 191 b.append(-500.0 / scaling_matrix[i, i]) 192 j = j + 2 193 194 # The cone opening angles and sigma_max. 195 if cdp.params[i] in ['cone_theta', 'cone_theta_x', 'cone_theta_y', 'cone_sigma_max']: 196 # 0 <= theta <= pi. 197 A.append(zero_array * 0.0) 198 A.append(zero_array * 0.0) 199 A[j][i] = 1.0 200 A[j+1][i] = -1.0 201 b.append(0.0) 202 b.append(-pi / scaling_matrix[i, i]) 203 j = j + 2 204 205 # The pseudo-ellipse restriction (theta_x <= theta_y). 206 if cdp.params[i] == 'cone_theta_y': 207 for m in range(n): 208 if cdp.params[m] == 'cone_theta_x': 209 A.append(zero_array * 0.0) 210 A[j][i] = 1.0 211 A[j][m] = -1.0 212 b.append(0.0) 213 j = j + 1 214 215 216 # Convert to numpy data structures. 217 A = array(A, float64) 218 b = array(b, float64) 219 220 # No constraints are present. 221 if len(A) == 0: 222 A = None 223 b = None 224 225 # Return the constraint objects. 226 return A, b

227 228

229 -def param_num():

230 """Determine the number of parameters in the model. 231 232 @return: The number of model parameters. 233 @rtype: int 234 """ 235 236 # Update the model, just in case. 237 update_model(verbosity=0) 238 239 # Simple parameter list count. 240 return len(cdp.params)

241 242

243 -def update_model(verbosity=1):

244 """Update the model parameters as necessary. 245 246 @keyword verbosity: If this value is greater than zero, then a number of printouts will occur. 247 @type verbosity: int 248 """ 249 250 # Re-initialise the list of model parameters. 251 cdp.params = [] 252 updated = [] 253 254 # The pivot parameters. 255 if not pivot_fixed(): 256 updated.append("pivot parameters") 257 cdp.params.append('pivot_x') 258 cdp.params.append('pivot_y') 259 cdp.params.append('pivot_z') 260 261 # The 2nd pivot point parameters - the minimum inter rotor axis distance. 262 if cdp.model in [MODEL_DOUBLE_ROTOR]: 263 updated.append("pivot parameters") 264 cdp.params.append('pivot_disp') 265 266 # The average domain position translation parameters. 267 updated.append("average domain position") 268 cdp.params.append('ave_pos_x') 269 cdp.params.append('ave_pos_y') 270 cdp.params.append('ave_pos_z') 271 272 # The tensor rotation, or average domain position. 273 if cdp.model not in MODEL_LIST_FREE_ROTORS: 274 cdp.params.append('ave_pos_alpha') 275 cdp.params.append('ave_pos_beta') 276 cdp.params.append('ave_pos_gamma') 277 278 # Frame order eigenframe - the full frame. 279 if cdp.model in MODEL_LIST_PSEUDO_ELLIPSE + [MODEL_DOUBLE_ROTOR]: 280 updated.append("frame order eigenframe") 281 cdp.params.append('eigen_alpha') 282 cdp.params.append('eigen_beta') 283 cdp.params.append('eigen_gamma') 284 285 # Frame order eigenframe - the isotropic cone axis. 286 if cdp.model in MODEL_LIST_ISO_CONE: 287 updated.append("frame order eigenframe") 288 cdp.params.append('axis_theta') 289 cdp.params.append('axis_phi') 290 291 # Frame order eigenframe - the rotor axis alpha angle. 292 if cdp.model in [MODEL_ROTOR, MODEL_FREE_ROTOR]: 293 updated.append("frame order eigenframe") 294 cdp.params.append('axis_alpha') 295 296 # Cone parameters - pseudo-elliptic cone parameters. 297 if cdp.model in MODEL_LIST_PSEUDO_ELLIPSE: 298 updated.append("cone opening half-angles") 299 cdp.params.append('cone_theta_x') 300 cdp.params.append('cone_theta_y') 301 302 # Cone parameters - single isotropic angle or order parameter. 303 if cdp.model in [MODEL_ISO_CONE, MODEL_ISO_CONE_FREE_ROTOR, MODEL_ISO_CONE_TORSIONLESS]: 304 updated.append("cone opening half-angles") 305 cdp.params.append('cone_theta') 306 307 # Cone parameters - torsion angle. 308 if cdp.model in [MODEL_DOUBLE_ROTOR, MODEL_ROTOR, MODEL_ISO_CONE, MODEL_PSEUDO_ELLIPSE]: 309 updated.append("cone torsion half-angles") 310 cdp.params.append('cone_sigma_max') 311 312 # Cone parameters - 2nd torsion angle. 313 if cdp.model in [MODEL_DOUBLE_ROTOR]: 314 updated.append("cone torsion half-angles") 315 cdp.params.append('cone_sigma_max_2') 316 317 # Printout. 318 if verbosity: 319 print("Reinitialising the list of model parameters: %s" % unique_elements(updated))

320

Source Code for Module specific_analyses.frame_order.parameters