Author: bugman
Date: Tue May 7 12:29:02 2013
New Revision: 19668
URL: http://svn.gna.org/viewcvs/relax?rev=19668&view=rev
Log:
Added support for handling missing data in the relaxation dispersion analysis.
This support was mentioned in the post
http://thread.gmane.org/gmane.science.nmr.relax.devel/3835.
Modified:
branches/relax_disp/specific_analyses/relax_disp/__init__.py
branches/relax_disp/specific_analyses/relax_disp/disp_data.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=19668&r1=19667&r2=19668&view=diff
==============================================================================
--- branches/relax_disp/specific_analyses/relax_disp/__init__.py (original)
+++ branches/relax_disp/specific_analyses/relax_disp/__init__.py Tue May 7
12:29:02 2013
@@ -131,9 +131,10 @@
# Initialise the data structures for the target function.
values = zeros((1, 1, cdp.dispersion_points), float64)
errors = zeros((1, 1, cdp.dispersion_points), float64)
+ missing = 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_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)
+ 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,
missing=missing, 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)
@@ -976,7 +977,7 @@
# Loop over the spin blocks.
for spins, spin_ids in self.model_loop():
# The R2eff/R1rho data.
- values, errors = return_r2eff_arrays(spins=spins,
spin_ids=spin_ids, fields=fields, field_count=field_count)
+ values, errors, missing = return_r2eff_arrays(spins=spins,
spin_ids=spin_ids, fields=fields, field_count=field_count)
# Create the initial parameter vector.
param_vector = assemble_param_vector(spins=spins)
@@ -1009,7 +1010,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,
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, 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):
Modified: branches/relax_disp/specific_analyses/relax_disp/disp_data.py
URL:
http://svn.gna.org/viewcvs/relax/branches/relax_disp/specific_analyses/relax_disp/disp_data.py?rev=19668&r1=19667&r2=19668&view=diff
==============================================================================
--- branches/relax_disp/specific_analyses/relax_disp/disp_data.py (original)
+++ branches/relax_disp/specific_analyses/relax_disp/disp_data.py Tue May 7
12:29:02 2013
@@ -24,7 +24,7 @@
"""Functions for handling relaxation dispersion data within the relax data
store."""
# Python module imports.
-from numpy import float64, zeros
+from numpy import float64, int32, ones, zeros
# relax module imports.
from lib.errors import RelaxError, RelaxNoSpectraError
@@ -329,16 +329,17 @@
@type field_count: int
@keyword sim_index: The index of the simulation to return the data
of. This should be None if the normal data is required.
@type sim_index: None or int
- @return: The numpy array structure of the R2eff/R1rho
values and the structure for the errors. For each structure, the first
dimension corresponds to the spins of a spin block, the second to the
spectrometer field strength, and the third is the dispersion points.
- @rtype: numpy rank-3 float array, numpy rank-3 float
array
+ @return: The numpy array structures of the R2eff/R1rho
values, errors and missing data. For each structure, the first dimension
corresponds to the spins of a spin block, the second to the spectrometer
field strength, and the third is the dispersion points.
+ @rtype: numpy rank-3 float array, numpy rank-3 float
array, numpy rank-3 int array
"""
# The spin count.
spin_num = len(spins)
- # Initialise the data structures for the target function.
+ # Initialise the data structures for the target function (errors are set
to one to avoid divide by zero for missing data in the chi-squared function).
values = zeros((spin_num, field_count, cdp.dispersion_points), float64)
- errors = zeros((spin_num, field_count, cdp.dispersion_points), float64)
+ errors = ones((spin_num, field_count, cdp.dispersion_points), float64)
+ missing = ones((spin_num, field_count, cdp.dispersion_points), int32)
# Pack the R2eff/R1rho data.
data_flag = False
@@ -384,12 +385,15 @@
# The errors.
errors[spin_index, field_index, disp_pt_index] =
spin.r2eff_err[key]
+ # Flip the missing flag to off.
+ missing[spin_index, field_index, disp_pt_index] = 0
+
# 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)
# Return the structures.
- return values, errors
+ return values, errors, missing
def return_key(frq=None, point=None, time=None):
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=19668&r1=19667&r2=19668&view=diff
==============================================================================
--- branches/relax_disp/target_functions/relax_disp.py (original)
+++ branches/relax_disp/target_functions/relax_disp.py Tue May 7 12:29:02
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, 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,
cpmg_frqs=None, spin_lock_nu1=None, scaling_matrix=None):
"""Relaxation dispersion target functions for optimisation.
Models
@@ -60,6 +60,8 @@
@type values: numpy rank-3 float array
@keyword errors: The R2eff/R1rho errors. The three
dimensions must correspond to those of the values argument.
@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 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.
@@ -71,6 +73,12 @@
# Check the args.
if model not in [MODEL_R2EFF, MODEL_LM63, MODEL_CR72]:
raise RelaxError("The model '%s' is unknown." % model)
+ if values == None:
+ raise RelaxError("No values have been supplied to the target
function.")
+ if errors == None:
+ raise RelaxError("No errors have been supplied to the target
function.")
+ if missing == None:
+ raise RelaxError("No missing data information has been supplied
to the target function.")
# Store the arguments.
self.num_params = num_params
@@ -79,6 +87,7 @@
self.num_disp_points = num_disp_points
self.values = values
self.errors = errors
+ self.missing = missing
self.cpmg_frqs = cpmg_frqs
self.spin_lock_nu1 = spin_lock_nu1
self.scaling_matrix = scaling_matrix
@@ -119,6 +128,11 @@
# 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)
+ # 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):
+ if self.missing[spin_index, frq_index, point_index]:
+ self.back_calc[spin_index, frq_index, point_index] =
self.values[spin_index, frq_index, point_index]
+
# Calculate and return the chi-squared value.
chi2_sum += chi2(self.values[spin_index, frq_index],
self.back_calc[spin_index, frq_index], self.errors[spin_index, frq_index])
r19668 - in /branches/relax_disp: specific_analyses/relax_disp/ target_functions/