Author: bugman
Date: Thu May 23 19:18:00 2013
New Revision: 19712
URL: http://svn.gna.org/viewcvs/relax?rev=19712&view=rev
Log:
Fixes for the LM63 dispersion CPMG model.
The 'r2' model parameter is now an array as there is one R2 value per
magnetic field strength. And
the 'rex' parameter has been renamed to 'phi_ex' and is scaled quadratically
with the field strength
within the optimisation target function.
Modified:
branches/relax_disp/specific_analyses/relax_disp/__init__.py
branches/relax_disp/specific_analyses/relax_disp/parameters.py
branches/relax_disp/target_functions/relax_disp.py
Modified: branches/relax_disp/specific_analyses/relax_disp/__init__.py
URL:
http://svn.gna.org/viewcvs/relax/branches/relax_disp/specific_analyses/relax_disp/__init__.py?rev=19712&r1=19711&r2=19712&view=diff
==============================================================================
--- branches/relax_disp/specific_analyses/relax_disp/__init__.py (original)
+++ branches/relax_disp/specific_analyses/relax_disp/__init__.py Thu May 23
19:18:00 2013
@@ -90,8 +90,8 @@
self.PARAMS.add('spin_lock_nu1', scope='spin', py_type=dict,
grace_string='\\qSpin-lock field strength (Hz)\\Q')
self.PARAMS.add('r2eff', scope='spin', default=15.0, desc='The
effective transversal relaxation rate', set='params', py_type=dict,
grace_string='\\qR\\s2,eff\\N\\Q (rad.s\\S-1\\N)', err=True, sim=True)
self.PARAMS.add('i0', scope='spin', default=10000.0, desc='The
initial intensity', py_type=dict, set='params', grace_string='\\qI\\s0\\Q',
err=True, sim=True)
- self.PARAMS.add('r2', scope='spin', default=15.0, desc='The
transversal relaxation rate', set='params', py_type=float,
grace_string='\\qR\\s2\\N\\Q (rad.s\\S-1\\N)', err=True, sim=True)
- self.PARAMS.add('rex', scope='spin', default=5.0, desc='The chemical
exchange contribution to R2 (sigma_ex = Rex / omega**2)', set='params',
py_type=float, grace_string='\\xF\\B\\sex\\N\\q (s.rad\\S-1\\N)', err=True,
sim=True)
+ self.PARAMS.add('r2', scope='spin', default=15.0, desc='The
transversal relaxation rate', set='params', py_type=list,
grace_string='\\qR\\s2\\N\\Q (rad.s\\S-1\\N)', err=True, sim=True)
+ self.PARAMS.add('phi_ex', scope='spin', default=5.0, desc='The
pA.pB.dw**2 value scaled by wH (phi_ex = pA * pB * Delta_omega**2 /
omega_H**2)', set='params', py_type=float, grace_string='\\xF\\B\\sex\\N\\q
(p\\sA\\N.p\\sB\\N.\\xDw\\B\\S2\\N / \\xw\\B\\sH\\N\\S2\\N)', err=True,
sim=True)
self.PARAMS.add('kex', scope='spin', default=10000.0, desc='The
exchange rate', set='params', py_type=float, grace_string='\\qk\\sex\\N\\Q
(rad.s\\S-1\\N)', err=True, sim=True)
self.PARAMS.add('r2a', scope='spin', default=15.0, desc='The
transversal relaxation rate for state A', set='params', py_type=float,
grace_string='\\qR\\s2,A\\N\\Q (rad.s\\S-1\\N)', err=True, sim=True)
self.PARAMS.add('ka', scope='spin', default=10000.0, desc='The
exchange rate from state A to state B', set='params', py_type=float,
grace_string='\\qk\\sA\\N\\Q (rad.s\\S-1\\N)', err=True, sim=True)
@@ -131,7 +131,7 @@
values, errors, missing = return_r2eff_arrays(spins=[spin],
spin_ids=[spin_id], fields=fields, field_count=field_count)
# Initialise the relaxation dispersion fit functions.
- model = Dispersion(model=cdp.model,
num_params=param_num(spins=[spin]), num_spins=1, num_frq=field_count,
num_disp_points=cdp.dispersion_points, values=values, errors=errors,
missing=missing, cpmg_frqs=return_cpmg_frqs(ref_flag=False),
spin_lock_nu1=return_spin_lock_nu1(ref_flag=False),
scaling_matrix=scaling_matrix)
+ model = Dispersion(model=cdp.model,
num_params=param_num(spins=[spin]), num_spins=1, num_frq=field_count,
num_disp_points=cdp.dispersion_points, values=values, errors=errors,
missing=missing, frqs=fields, cpmg_frqs=return_cpmg_frqs(ref_flag=False),
spin_lock_nu1=return_spin_lock_nu1(ref_flag=False),
scaling_matrix=scaling_matrix)
# Make a single function call. This will cause back calculation and
the data will be stored in the class instance.
model.func(param_vector)
@@ -382,10 +382,10 @@
lower.append(0.0)
upper.append(40.0)
- # Chemical exchange contribution to 'R2'.
- elif spin.params[i] == 'rex':
+ # The pA.pB.dw**2/wH**2 parameter.
+ elif spin.params[i] == 'phi_ex':
lower.append(0.0)
- upper.append(20.0)
+ upper.append(1e-17)
# Exchange rate.
elif spin.params[i] == 'kex':
@@ -727,12 +727,18 @@
# LM63 model.
elif model == MODEL_LM63:
print("The Luz and Meiboom (1963) 2-site fast exchange model.")
- params = ['r2', 'rex', 'kex']
+ params = []
+ for i in range(cdp.spectro_frq_count):
+ params.append('r2')
+ params += ['phi_ex', 'kex']
# CR72 model.
elif model == MODEL_CR72:
print("The Carver and Richards (1972) 2-site model for all time
scales.")
- params = ['r2', 'r2a', 'ka', 'dw']
+ params = []
+ for i in range(cdp.spectro_frq_count):
+ params.append('r2')
+ params += ['r2a', 'ka', 'dw']
# Invalid model.
else:
@@ -1056,7 +1062,7 @@
print("Unconstrained grid search size: %s (constraints
may decrease this size).\n" % grid_size)
# Initialise the function to minimise.
- model = Dispersion(model=cdp.model,
num_params=param_num(spins=spins), num_spins=len(spins), num_frq=field_count,
num_disp_points=cdp.dispersion_points, values=values, errors=errors,
missing=missing, cpmg_frqs=return_cpmg_frqs(ref_flag=False),
spin_lock_nu1=return_spin_lock_nu1(ref_flag=False),
scaling_matrix=scaling_matrix)
+ model = Dispersion(model=cdp.model,
num_params=param_num(spins=spins), num_spins=len(spins), num_frq=field_count,
num_disp_points=cdp.dispersion_points, values=values, errors=errors,
missing=missing, frqs=fields, cpmg_frqs=return_cpmg_frqs(ref_flag=False),
spin_lock_nu1=return_spin_lock_nu1(ref_flag=False),
scaling_matrix=scaling_matrix)
# Grid search.
if search('^[Gg]rid', min_algor):
@@ -1217,7 +1223,7 @@
_table = uf_tables.add_table(label="table: dispersion curve-fit data
type patterns", caption="Relaxation dispersion curve fitting data type string
matching patterns.")
_table.add_headings(["Data type", "Object name"])
_table.add_row(["Transversal relaxation rate", "'r2'"])
- _table.add_row(["Chemical exchange contribution to 'R2'", "'rex'"])
+ _table.add_row(["The pA.pB.dw**2/wH**2 parameter", "'phi_ex'"])
_table.add_row(["Exchange rate", "'kex'"])
_table.add_row(["Transversal relaxation rate for state A", "'r2a'"])
_table.add_row(["Exchange rate from state A to state B", "'ka'"])
Modified: branches/relax_disp/specific_analyses/relax_disp/parameters.py
URL:
http://svn.gna.org/viewcvs/relax/branches/relax_disp/specific_analyses/relax_disp/parameters.py?rev=19712&r1=19711&r2=19712&view=diff
==============================================================================
--- branches/relax_disp/specific_analyses/relax_disp/parameters.py (original)
+++ branches/relax_disp/specific_analyses/relax_disp/parameters.py Thu May 23
19:18:00 2013
@@ -110,20 +110,20 @@
# Transversal relaxation rate.
if spin.params[i] == 'r2':
if sim_index != None:
- param_vector.append(spin.r2_sim[sim_index])
- elif spin.r2 == None:
+ param_vector.append(spin.r2_sim[sim_index][i])
+ elif spin.r2 == [] or spin.r2[i] == None:
param_vector.append(0.0)
else:
- param_vector.append(spin.r2)
-
- # Chemical exchange contribution to 'R2'.
- if spin.params[i] == 'rex':
- if sim_index != None:
- param_vector.append(spin.rex_sim[sim_index])
- elif spin.rex == None:
+ param_vector.append(spin.r2[i])
+
+ # The pA.pB.dw**2/wH**2 parameter.
+ if spin.params[i] == 'phi_ex':
+ if sim_index != None:
+ param_vector.append(spin.phi_ex_sim[sim_index])
+ elif spin.phi_ex == None:
param_vector.append(0.0)
else:
- param_vector.append(spin.rex)
+ param_vector.append(spin.phi_ex)
# Exchange rate.
elif spin.params[i] == 'kex':
@@ -223,9 +223,9 @@
if spin.params[i] == 'r2':
scaling_matrix[param_index, param_index] = 10
- # Chemical exchange contribution to 'R2' scaling.
- elif spin.params[i] == 'rex':
- scaling_matrix[param_index, param_index] = 10
+ # The pA.pB.dw**2/wH**2 parameter.
+ elif spin.params[i] == 'phi_ex':
+ scaling_matrix[param_index, param_index] = 1e18
# Exchange rate scaling.
elif spin.params[i] == 'kex':
@@ -326,21 +326,32 @@
# Reset the parameter index.
param_index = 0
+ # Initialise the parameter if needed.
+ if 'r2' in spin.params:
+ if sim_index != None:
+ spin.r2_sim[sim_index] = []
+ for i in range(cdp.spectro_frq_count):
+ spin.r2_sim[sim_index].append(None)
+ else:
+ spin.r2 = []
+ for i in range(cdp.spectro_frq_count):
+ spin.r2.append(None)
+
# Loop over each parameter.
for i in range(len(spin.params)):
# Transversal relaxation rate.
if spin.params[i] == 'r2':
if sim_index != None:
- spin.r2_sim[sim_index] = param_vector[param_index]
+ spin.r2_sim[sim_index][i] = param_vector[param_index]
else:
- spin.r2 = param_vector[param_index]
-
- # Chemical exchange contribution to 'R2'.
- if spin.params[i] == 'rex':
+ spin.r2[i] = param_vector[param_index]
+
+ # The pA.pB.dw**2/wH**2 parameter.
+ if spin.params[i] == 'phi_ex':
if sim_index != None:
- spin.rex_sim[sim_index] = param_vector[param_index]
+ spin.phi_ex_sim[sim_index] =
param_vector[param_index]
else:
- spin.rex = param_vector[param_index]
+ spin.phi_ex = param_vector[param_index]
# Exchange rate.
elif spin.params[i] == 'kex':
@@ -439,9 +450,9 @@
b.append(0.0)
j += 1
- # Relaxation rates and Rex.
- elif search('^r', spin.params[k]):
- # Rex, R2A >= 0.
+ # Relaxation rates and phi_ex.
+ elif spin.params[k] in ['r2a', 'phi_ex']:
+ # phi_ex, R2A >= 0.
A.append(zero_array * 0.0)
A[j][i] = 1.0
b.append(0.0)
Modified: branches/relax_disp/target_functions/relax_disp.py
URL:
http://svn.gna.org/viewcvs/relax/branches/relax_disp/target_functions/relax_disp.py?rev=19712&r1=19711&r2=19712&view=diff
==============================================================================
--- branches/relax_disp/target_functions/relax_disp.py (original)
+++ branches/relax_disp/target_functions/relax_disp.py Thu May 23 19:18:00
2013
@@ -34,7 +34,7 @@
class Dispersion:
- def __init__(self, model=None, num_params=None, num_spins=None,
num_frq=None, num_disp_points=None, values=None, errors=None, missing=None,
cpmg_frqs=None, spin_lock_nu1=None, scaling_matrix=None):
+ def __init__(self, model=None, num_params=None, num_spins=None,
num_frq=None, num_disp_points=None, values=None, errors=None, missing=None,
frqs=None, cpmg_frqs=None, spin_lock_nu1=None, scaling_matrix=None):
"""Relaxation dispersion target functions for optimisation.
Models
@@ -62,6 +62,8 @@
@type errors: numpy rank-3 float array
@keyword missing: The data structure indicating missing
R2eff/R1rho data. The three dimensions must correspond to those of the
values argument.
@type missing: numpy rank-3 int array
+ @keyword frqs: The spectrometer frequencies in Hz.
+ @type frqs: numpy rank-1 float array
@keyword cpmg_frqs: The CPMG frequencies in Hertz for each
separate dispersion point. This will be ignored for R1rho experiments.
@type cpmg_frqs: numpy rank-1 float array
@keyword spin_lock_nu1: The spin-lock field strengths in Hertz
for each separate dispersion point. This will be ignored for CPMG
experiments.
@@ -88,6 +90,7 @@
self.values = values
self.errors = errors
self.missing = missing
+ self.frqs = frqs
self.cpmg_frqs = cpmg_frqs
self.spin_lock_nu1 = spin_lock_nu1
self.scaling_matrix = scaling_matrix
@@ -118,6 +121,12 @@
if self.scaling_flag:
params = dot(params, self.scaling_matrix)
+ # Unpack the parameter values.
+ index = self.num_frq - 1
+ R20 = params[:index+1]
+ phi_ex = params[index+1]
+ kex = params[index+2]
+
# Initialise.
chi2_sum = 0.0
@@ -125,8 +134,11 @@
for spin_index in range(self.num_spins):
# Loop over the spectrometer frequencies.
for frq_index in range(self.num_frq):
+ # Spectrometer field strength scaling.
+ phi_ex_scaled = phi_ex * self.frqs[frq_index]**2
+
# Back calculate the R2eff values.
- r2eff_LM63(r20=params[0], phi_ex=params[1], kex=params[2],
cpmg_frqs=self.cpmg_frqs, back_calc=self.back_calc[spin_index, frq_index],
num_points=self.num_disp_points)
+ r2eff_LM63(r20=R20[frq_index], phi_ex=phi_ex, kex=kex,
cpmg_frqs=self.cpmg_frqs, back_calc=self.back_calc[spin_index, frq_index],
num_points=self.num_disp_points)
# For all missing data points, set the back-calculated value
to the measured values so that it has no effect on the chi-squared value.
for point_index in range(self.num_disp_points):
r19712 - in /branches/relax_disp: specific_analyses/relax_disp/ target_functions/