Author: bugman
Date: Fri May 3 21:57:50 2013
New Revision: 19661
URL: http://svn.gna.org/viewcvs/relax?rev=19661&view=rev
Log:
Redesigned the optimisation code of the dispersion analysis specific class
for the new target functions.
This includes the assembling of R2eff/R1rho values instead of peak heights,
and a number of small
fixes.
Modified:
branches/relax_disp/specific_analyses/relax_disp/__init__.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=19661&r1=19660&r2=19661&view=diff
==============================================================================
--- branches/relax_disp/specific_analyses/relax_disp/__init__.py (original)
+++ branches/relax_disp/specific_analyses/relax_disp/__init__.py Fri May 3
21:57:50 2013
@@ -55,7 +55,7 @@
from pipe_control.result_files import add_result_file
from specific_analyses.api_base import API_base
from specific_analyses.api_common import API_common
-from specific_analyses.relax_disp.disp_data import exp_curve_index_from_key,
exp_curve_key_from_index, intensity_key, loop_all_data,
loop_dispersion_point, loop_exp_curve, loop_spectrometer, relax_time,
return_key
+from specific_analyses.relax_disp.disp_data import exp_curve_index_from_key,
exp_curve_key_from_index, intensity_key, loop_all_data,
loop_dispersion_point, loop_exp_curve, loop_spectrometer, relax_time,
return_cpmg_frqs, return_key, return_spin_lock_nu1
from specific_analyses.relax_disp.parameters import assemble_param_vector,
assemble_scaling_matrix, disassemble_param_vector, linear_constraints,
param_index_to_param_info, param_num
from specific_analyses.relax_disp.variables import CPMG_EXP, FIXED_TIME_EXP,
MODEL_R2EFF, MODEL_LM63, MODEL_CR72, R1RHO_EXP, VAR_TIME_EXP
from target_functions.relax_disp import Dispersion
@@ -129,11 +129,11 @@
scaling_matrix = assemble_scaling_matrix(spins=[spin], scaling=False)
# Initialise the data structures for the target function.
- values = zeros((1, cdp.dispersion_points, cdp.num_time_pts), float64)
- errors = zeros((1, cdp.dispersion_points, cdp.num_time_pts), float64)
+ values = zeros((1, 1, cdp.dispersion_points), float64)
+ errors = zeros((1, 1, cdp.dispersion_points), float64)
# Initialise the relaxation dispersion fit functions.
- model = Dispersion(model=cdp.model,
num_params=param_num(spins=[spin]), num_spins=1,
num_exp_curves=cdp.dispersion_points, num_times=cdp.num_time_pts,
values=values, errors=errors, cpmg_frqs=cdp.cpmg_frqs_list,
spin_lock_nu1=cdp.spin_lock_nu1_list, relax_times=cdp.relax_time_list,
scaling_matrix=scaling_matrix)
+ model = Dispersion(model=cdp.model,
num_params=param_num(spins=[spin]), num_spins=1, num_frq=1,
num_disp_points=cdp.dispersion_points, values=values, errors=errors,
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)
@@ -337,57 +337,60 @@
lower = []
upper = []
- # First the spin specific parameters.
- for spin_index in range(len(spins)):
- # Alias the spin.
- spin = spins[spin_index]
-
- # Loop over each exponential curve.
- for exp_i, key in loop_exp_curve():
- # Loop over the parameters.
- for i in range(len(spin.params)):
- # R2eff relaxation rate (from 0 to 40 s^-1).
- if spin.params[i] == 'r2eff':
- lower.append(0.0)
- upper.append(40.0)
-
- # Intensity.
- elif spin.params[i] == 'i0':
- lower.append(0.0)
- upper.append(max(spin.intensities.values()))
-
- # Then the spin block specific parameters.
- spin = spins[0]
- for i in range(len(spin.params)):
- # R2 relaxation rate (from 0 to 40 s^-1).
- if spin.params[i] == 'r2':
- lower.append(0.0)
- upper.append(40.0)
-
- # Chemical exchange contribution to 'R2'.
- elif spin.params[i] == 'rex':
- lower.append(0.0)
- upper.append(20.0)
-
- # Exchange rate.
- elif spin.params[i] == 'kex':
- lower.append(0.0)
- upper.append(100000.0)
-
- # Transversal relaxation rate for state A.
- elif spin.params[i] == 'r2a':
- lower.append(0.0)
- upper.append(20.0)
-
- # Exchange rate from state A to state B.
- elif spin.params[i] == 'ka':
- lower.append(0.0)
- upper.append(100000.0)
-
- # Chemical shift difference between states A and B.
- elif spin.params[i] == 'dw':
- lower.append(0.0)
- upper.append(10000.0)
+ # The R2eff model.
+ if cdp.model == MODEL_R2EFF:
+ for spin_index in range(len(spins)):
+ # Alias the spin.
+ spin = spins[spin_index]
+
+ # Loop over each exponential curve.
+ for exp_i, key in loop_exp_curve():
+ # Loop over the parameters.
+ for i in range(len(spin.params)):
+ # R2eff relaxation rate (from 0 to 40 s^-1).
+ if spin.params[i] == 'r2eff':
+ lower.append(0.0)
+ upper.append(40.0)
+
+ # Intensity.
+ elif spin.params[i] == 'i0':
+ lower.append(0.0)
+ upper.append(max(spin.intensities.values()))
+
+ # All other models.
+ else:
+ # Only use the parameters of the first spin of the cluster.
+ spin = spins[0]
+ for i in range(len(spin.params)):
+ # R2 relaxation rate (from 0 to 40 s^-1).
+ if spin.params[i] == 'r2':
+ lower.append(0.0)
+ upper.append(40.0)
+
+ # Chemical exchange contribution to 'R2'.
+ elif spin.params[i] == 'rex':
+ lower.append(0.0)
+ upper.append(20.0)
+
+ # Exchange rate.
+ elif spin.params[i] == 'kex':
+ lower.append(0.0)
+ upper.append(100000.0)
+
+ # Transversal relaxation rate for state A.
+ elif spin.params[i] == 'r2a':
+ lower.append(0.0)
+ upper.append(20.0)
+
+ # Exchange rate from state A to state B.
+ elif spin.params[i] == 'ka':
+ lower.append(0.0)
+ upper.append(100000.0)
+
+ # Chemical shift difference between states A and B.
+ elif spin.params[i] == 'dw':
+ lower.append(0.0)
+ upper.append(10000.0)
# The full grid size.
grid_size = 1
@@ -460,7 +463,7 @@
# Get the grid search minimisation options.
lower_new, upper_new = None, None
if match('^[Gg]rid', min_algor):
- grid_size, inc_new, lower_new, upper_new, =
self._grid_search_setup(spins=[spin], param_vector=param_vector, lower=lower,
upper=upper, inc=inc, scaling_matrix=scaling_matrix)
+ grid_size, inc_new, lower_new, upper_new =
self._grid_search_setup(spins=[spin], param_vector=param_vector, lower=lower,
upper=upper, inc=inc, scaling_matrix=scaling_matrix)
# Linear constraints.
A, b = None, None
@@ -976,18 +979,24 @@
spin_num = len(spins)
# Initialise the data structures for the target function.
- values = zeros((spin_num, field_count, cdp.dispersion_points,
num_time_pts), float64)
- errors = zeros((spin_num, field_count, cdp.dispersion_points,
num_time_pts), float64)
-
- # Pack the peak intensity data.
+ values = zeros((spin_num, field_count, cdp.dispersion_points),
float64)
+ errors = zeros((spin_num, field_count, cdp.dispersion_points),
float64)
+
+ # Pack the R2eff/R1rho data.
+ data_flag = False
for spin_index in range(spin_num):
# Alias the spin.
spin = spins[spin_index]
+ # No data.
+ if not hasattr(spin, 'r2eff'):
+ continue
+ data_flag = True
+
# The keys.
- keys = list(spin.intensities.keys())
-
- # Loop over the spectral data.
+ keys = list(spin.r2eff.keys())
+
+ # Loop over the R2eff data.
for key in keys:
# The indices.
disp_pt_index = 0
@@ -995,22 +1004,22 @@
disp_pt_index =
cdp.cpmg_frqs_list.index(cdp.cpmg_frqs[key])
elif cdp.exp_type == 'r1rho':
disp_pt_index =
cdp.spin_lock_nu1_list.index(cdp.spin_lock_nu1[key])
- time_index = 0
- if hasattr(cdp, 'relax_time_list'):
- time_index =
cdp.relax_time_list.index(cdp.relax_times[key])
field_index = 0
if hasattr(cdp, 'frqs'):
field_index = fields.index(cdp.frqs[keys])
# The values.
if sim_index == None:
- values[spin_index, field_index, disp_pt_index,
time_index] = spin.intensities[key]
+ values[spin_index, field_index, disp_pt_index] =
spin.r2eff[key]
else:
- values[spin_index, field_index, disp_pt_index,
time_index] = spin.intensity_sim[key][sim_index]
+ values[spin_index, field_index, disp_pt_index] =
spin.r2eff_sim[key][sim_index]
# The errors.
- errors[spin_index, field_index, disp_pt_index,
time_index] = spin.intensity_err[key]
-
+ errors[spin_index, field_index, disp_pt_index] =
spin.r2eff_err[key]
+
+ # No R2eff/R1rho data for the spin cluster.
+ if not data_flag:
+ raise RelaxError("No R2eff/R1rho data could be found for the
spin cluster %s." % spin_ids)
# Create the initial parameter vector.
param_vector = assemble_param_vector(spins=spins)
@@ -1023,7 +1032,7 @@
# Get the grid search minimisation options.
lower_new, upper_new = None, None
if match('^[Gg]rid', min_algor):
- grid_size, inc_new, lower_new, upper_new, sparseness =
self._grid_search_setup(spins=spins, param_vector=param_vector, lower=lower,
upper=upper, inc=inc, scaling_matrix=scaling_matrix)
+ grid_size, inc_new, lower_new, upper_new =
self._grid_search_setup(spins=spins, param_vector=param_vector, lower=lower,
upper=upper, inc=inc, scaling_matrix=scaling_matrix)
# Linear constraints.
A, b = None, None
@@ -1043,11 +1052,11 @@
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=spin_num,
num_exp_curves=cdp.dispersion_points, num_times=cdp.num_time_pts,
values=values, errors=errors, cpmg_frqs=cdp.cpmg_frqs_list,
spin_lock_nu1=cdp.spin_lock_nu1_list, relax_times=cdp.relax_time_list,
scaling_matrix=scaling_matrix)
+ model = Dispersion(model=cdp.model,
num_params=param_num(spins=spins), num_spins=spin_num, num_frq=field_count,
num_disp_points=cdp.dispersion_points, values=values, errors=errors,
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):
- results = grid(func=model.func, args=(), num_incs=inc_new,
lower=lower_new, upper=upper_new, A=A, b=b, sparseness=sparseness,
verbosity=verbosity)
+ results = grid(func=model.func, args=(), num_incs=inc_new,
lower=lower_new, upper=upper_new, A=A, b=b, verbosity=verbosity)
# Unpack the results.
param_vector, chi2, iter_count, warning = results
r19661 - /branches/relax_disp/specific_analyses/relax_disp/__init__.py