specific_analyses.frame_order.parameters

1 ############################################################################### 2 # # 3 # Copyright (C) 2009-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 specific_analyses.frame_order.data import pivot_fixed, translation_fixed 31 32

33 -def assemble_param_vector(sim_index=None):

34 """Assemble and return the parameter vector. 35 36 @return: The parameter vector. 37 @rtype: numpy rank-1 array 38 @keyword sim_index: The Monte Carlo simulation index. 39 @type sim_index: int 40 """ 41 42 # Initialise. 43 param_vect = [] 44 45 # Parameter name extension. 46 ext = '' 47 if sim_index != None: 48 ext = '_sim' 49 50 # Loop over all model parameters. 51 for param_name in cdp.params: 52 # Get the object. 53 obj = getattr(cdp, param_name+ext) 54 55 # Add it to the parameter vector. 56 if sim_index == None: 57 param_vect.append(obj) 58 else: 59 param_vect.append(obj[sim_index]) 60 61 # Return as a numpy array. 62 return array(param_vect, float64)

63 64

65 -def linear_constraints(scaling_matrix=None):

66 """Create the linear constraint matrices A and b. 67 68 Standard notation 69 ================= 70 71 The parameter constraints for the motional amplitude parameters are:: 72 73 0 <= theta <= pi, 74 0 <= theta_x <= theta_y <= pi, 75 -0.125 <= S <= 1, 76 0 <= sigma_max <= pi, 77 78 The pivot point parameter, domain position parameters, and eigenframe parameters are unconstrained. 79 80 81 Matrix notation 82 =============== 83 84 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:: 85 86 | 1 0 0 0 0 | | 0 | 87 | | | | 88 |-1 0 0 0 0 | | -pi | 89 | | | | 90 | 0 1 0 0 0 | | 0 | 91 | | | | 92 | 0 -1 0 0 0 | | theta | | -pi | 93 | | | | | | 94 | 0 -1 1 0 0 | | theta_x | | 0 | 95 | | | | | | 96 | 0 0 1 0 0 | . | theta_y | >= | 0 | 97 | | | | | | 98 | 0 0 -1 0 0 | | S | | -pi | 99 | | | | | | 100 | 0 0 0 1 0 | | sigma_max | | -0.125 | 101 | | | | 102 | 0 0 0 -1 0 | | -1 | 103 | | | | 104 | 0 0 0 0 1 | | 0 | 105 | | | | 106 | 0 0 0 0 -1 | | -pi | 107 108 109 @keyword scaling_matrix: The diagonal, square scaling matrix. 110 @type scaling_matrix: numpy rank-2 square matrix 111 @return: The matrices A and b. 112 @rtype: numpy rank-2 NxM array, numpy rank-1 N array 113 """ 114 115 # Initialisation (0..j..m). 116 A = [] 117 b = [] 118 n = param_num() 119 zero_array = zeros(n, float64) 120 i = 0 121 j = 0 122 123 # Loop over the parameters of the model. 124 for i in range(n): 125 # The cone opening angles and sigma_max. 126 if cdp.params[i] in ['cone_theta', 'cone_theta_x', 'cone_theta_y', 'cone_sigma_max']: 127 # 0 <= theta <= pi. 128 A.append(zero_array * 0.0) 129 A.append(zero_array * 0.0) 130 A[j][i] = 1.0 131 A[j+1][i] = -1.0 132 b.append(0.0) 133 b.append(-pi / scaling_matrix[i, i]) 134 j = j + 2 135 136 # The pseudo-ellipse restriction (theta_x <= theta_y). 137 if cdp.params[i] == 'cone_theta_y': 138 for m in range(n): 139 if cdp.params[m] == 'cone_theta_x': 140 A.append(zero_array * 0.0) 141 A[j][i] = 1.0 142 A[j][m] = -1.0 143 b.append(0.0) 144 j = j + 1 145 146 147 # The order parameter. 148 if cdp.params[i] == 'cone_s1': 149 # -0.125 <= S <= 1. 150 A.append(zero_array * 0.0) 151 A.append(zero_array * 0.0) 152 A[j][i] = 1.0 153 A[j+1][i] = -1.0 154 b.append(-0.125 / scaling_matrix[i, i]) 155 b.append(-1 / scaling_matrix[i, i]) 156 j = j + 2 157 158 # Convert to numpy data structures. 159 A = array(A, float64) 160 b = array(b, float64) 161 162 # Return the constraint objects. 163 return A, b

164 165

166 -def param_num():

167 """Determine the number of parameters in the model. 168 169 @return: The number of model parameters. 170 @rtype: int 171 """ 172 173 # Update the model, just in case. 174 update_model() 175 176 # Simple parameter list count. 177 return len(cdp.params)

178 179

180 -def update_model():

181 """Update the model parameters as necessary.""" 182 183 # Re-initialise the list of model parameters. 184 cdp.params = [] 185 186 # The pivot parameters. 187 if not pivot_fixed(): 188 cdp.params.append('pivot_x') 189 cdp.params.append('pivot_y') 190 cdp.params.append('pivot_z') 191 192 # Double rotor. 193 if cdp.model == 'double rotor': 194 cdp.params.append('pivot_x_2') 195 cdp.params.append('pivot_y_2') 196 cdp.params.append('pivot_z_2') 197 198 # The average domain position translation parameters. 199 if not translation_fixed(): 200 cdp.params.append('ave_pos_x') 201 cdp.params.append('ave_pos_y') 202 cdp.params.append('ave_pos_z') 203 204 # The tensor rotation, or average domain position. 205 if cdp.model not in ['free rotor', 'iso cone, free rotor']: 206 cdp.params.append('ave_pos_alpha') 207 cdp.params.append('ave_pos_beta') 208 cdp.params.append('ave_pos_gamma') 209 210 # Frame order eigenframe - the full frame. 211 if cdp.model in ['pseudo-ellipse', 'pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor']: 212 cdp.params.append('eigen_alpha') 213 cdp.params.append('eigen_beta') 214 cdp.params.append('eigen_gamma') 215 216 # Frame order eigenframe - the isotropic cone axis. 217 if cdp.model in ['iso cone', 'free rotor', 'iso cone, torsionless', 'iso cone, free rotor', 'double rotor']: 218 cdp.params.append('axis_theta') 219 cdp.params.append('axis_phi') 220 221 # Frame order eigenframe - the second rotation axis. 222 if cdp.model in ['double rotor']: 223 cdp.params.append('axis_theta_2') 224 cdp.params.append('axis_phi_2') 225 226 # Frame order eigenframe - the rotor axis alpha angle. 227 if cdp.model in ['rotor']: 228 cdp.params.append('axis_alpha') 229 230 # Cone parameters - pseudo-elliptic cone parameters. 231 if cdp.model in ['pseudo-ellipse', 'pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor']: 232 cdp.params.append('cone_theta_x') 233 cdp.params.append('cone_theta_y') 234 235 # Cone parameters - single isotropic angle or order parameter. 236 if cdp.model in ['iso cone', 'iso cone, torsionless']: 237 cdp.params.append('cone_theta') 238 if cdp.model in ['iso cone, free rotor']: 239 cdp.params.append('cone_s1') 240 241 # Cone parameters - torsion angle. 242 if cdp.model in ['double rotor', 'rotor', 'line', 'iso cone', 'pseudo-ellipse']: 243 cdp.params.append('cone_sigma_max') 244 245 # Cone parameters - 2nd torsion angle. 246 if cdp.model in ['double rotor']: 247 cdp.params.append('cone_sigma_max_2') 248 249 # Initialise the parameters in the current data pipe. 250 for param in cdp.params: 251 if not param in ['pivot_x', 'pivot_y', 'pivot_z', 'pivot_x_2', 'pivot_y_2', 'pivot_z_2'] and not hasattr(cdp, param): 252 setattr(cdp, param, 0.0)

253

Source Code for Module specific_analyses.frame_order.parameters