Author: tlinnet
Date: Sun May 4 22:41:30 2014
New Revision: 22952
URL: http://svn.gna.org/viewcvs/relax?rev=22952&view=rev
Log:
Modified system_tests/scripts/relax_disp/cpmg_synthetic.py and the
corresponding systemtests.
Relax_disp.test_cpmg_synthetic_cr72.
Relax_disp.test_cpmg_synthetic_cr72_full_noise_cluster.
sr #3154: (https://gna.org/support/?3154) Implementation of Baldwin (2014)
B14 model - 2-site exact solution model for all time scales.
This follows the tutorial for adding relaxation dispersion models at:
http://wiki.nmr-relax.com/Tutorial_for_adding_relaxation_dispersion_models_to_relax#Debugging
Modified:
trunk/test_suite/system_tests/relax_disp.py
trunk/test_suite/system_tests/scripts/relax_disp/cpmg_synthetic.py
Modified: trunk/test_suite/system_tests/relax_disp.py
URL:
http://svn.gna.org/viewcvs/relax/trunk/test_suite/system_tests/relax_disp.py?rev=22952&r1=22951&r2=22952&view=diff
==============================================================================
--- trunk/test_suite/system_tests/relax_disp.py (original)
+++ trunk/test_suite/system_tests/relax_disp.py Sun May 4 22:41:30 2014
@@ -498,8 +498,8 @@
# The grid search size (the number of increments per dimension).
# If None, use the default values.
- GRID = None
- #GRID = 21
+ #GRID = None
+ GRID = 13
# Perform Grid Search.
if GRID:
# Set the R20 parameters in the default grid search using the
minimum R2eff value.
@@ -514,7 +514,7 @@
for param in MODEL_PARAMS[MODEL]:
self.interpreter.value.set(param=param, index=None)
# Do a grid search, which will store the chi2 value.
- self.interpreter.grid_search(lower=None, upper=None, inc=1,
constraints=True, verbosity=1)
+ self.interpreter.grid_search(lower=None, upper=None, inc=1,
constraints=True, verbosity=1)
# Store result.
for spin, mol_name, resi, resn, spin_id in spin_loop(full_info=True,
return_id=True, skip_desel=True):
@@ -731,7 +731,7 @@
relax_disp.Relax_disp(pipe_name='origin - relax_disp (Sun Feb 23
19:36:51 2014)', pipe_bundle='relax_disp (Sun Feb 23 19:36:51 2014)',
results_dir=self.tmpdir, models=['R2eff', 'No Rex'], grid_inc=11,
mc_sim_num=2, modsel='AIC', pre_run_dir=pre_run_dir, insignificance=1.0,
numeric_only=True, mc_sim_all_models=False, eliminate=True)
- def test_cpmg_synthetic(self):
+ def test_cpmg_synthetic_cr72(self):
"""Test synthetic cpmg data.
This script will produce synthetic CPMG R2eff values according to
the selected model, and the fit the selected model.
@@ -761,7 +761,7 @@
exps = [exp_1, exp_2]
spins = [
- ['Ala', 1, 'N', {'r2': {r20_key_1:2, r20_key_2:2}, 'kex': 1000,
'pA': 0.99, 'dw': 2} ]
+ ['Ala', 1, 'N', {'r2': {r20_key_1:2, r20_key_2:2}, 'r2a':
{r20_key_1:2, r20_key_2:2}, 'r2b': {r20_key_1:2, r20_key_2:2}, 'kex': 1000,
'pA': 0.99, 'dw': 2} ]
]
# Collect the data to be used.
@@ -783,6 +783,9 @@
# The grid search size (the number of increments per dimension).
ds.GRID_INC = 8
+ # The do clustering.
+ ds.do_cluster = False
+
# The function tolerance. This is used to terminate minimisation
once the function value between iterations is less than the tolerance.
# The default value is 1e-25.
ds.set_func_tol = 1e-9
@@ -799,6 +802,12 @@
# The plot_curves.
ds.plot_curves = False
+
+ # The conversion for ShereKhan at http://sherekhan.bionmr.org/.
+ ds.sherekhan_input = False
+
+ # Make a dx map to be opened om OpenDX. To map the hypersurface of
chi2, when altering kex, dw and pA.
+ ds.opendx = False
# The set r2eff err.
ds.r2eff_err = 0.1
@@ -849,12 +858,12 @@
print("WARNING: rel change level is above %.2f, and
is %.4f."%(ds.rel_change, rel_change))
print("###################################")
- ## Make test on parameters.
- self.assertAlmostEqual(set_val, min_val, 10)
-
-
- def test_cpmg_synthetic_cr72_fail(self):
- """Test synthetic cpmg data. Make CR72 model fail, for small noise.
+ ## Make test on parameters.
+ self.assertAlmostEqual(set_val, min_val, 2)
+
+
+ def test_cpmg_synthetic_cr72_full_noise_cluster(self):
+ """Test synthetic cpmg data. For CR72 with small noise and cluster.
This script will produce synthetic CPMG R2eff values according to
the selected model, and the fit the selected model.
"""
@@ -869,22 +878,24 @@
r20_key_1 = generate_r20_key(exp_type=EXP_TYPE_CPMG_SQ, frq=sfrq_1)
time_T2_1 = 0.06
ncycs_1 = [2, 4, 8, 10, 20, 30, 40, 60]
- r2eff_err_1 = [0.05, -0.05, 0.05, -0.05, 0.05, -0.05, 0.05, -0.05]
- exp_1 = [sfrq_1, time_T2_1, ncycs_1, r2eff_err_1]
+ r2eff_errs_1 = [0.05, -0.05, 0.05, -0.05, 0.05, -0.05, 0.05, -0.05]
+ #r2eff_errs_1 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+ exp_1 = [sfrq_1, time_T2_1, ncycs_1, r2eff_errs_1]
sfrq_2 = 499.8908617*1E6
r20_key_2 = generate_r20_key(exp_type=EXP_TYPE_CPMG_SQ, frq=sfrq_2)
time_T2_2 = 0.05
ncycs_2 = [2, 4, 8, 10, 30, 35, 40, 50]
- r2eff_err_2 = [0.05, -0.05, 0.05, -0.05, 0.05, -0.05, 0.05, -0.05]
- exp_2 = [sfrq_2, time_T2_2, ncycs_2, r2eff_err_2]
+ r2eff_errs_2 = [0.05, -0.05, 0.05, -0.05, 0.05, -0.05, 0.05, -0.05]
+ #r2eff_errs_2 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+ exp_2 = [sfrq_2, time_T2_2, ncycs_2, r2eff_errs_2]
# Collect all exps
exps = [exp_1, exp_2]
spins = [
- ['Ala', 1, 'N', {'r2': {r20_key_1:10, r20_key_2:11.5}, 'kex':
1000, 'pA': 0.99, 'dw': 2} ],
- ['Ala', 2, 'N', {'r2': {r20_key_1:13, r20_key_2:14.5}, 'kex':
1000, 'pA': 0.99, 'dw': 1} ]
+ ['Ala', 1, 'N', {'r2': {r20_key_1:10, r20_key_2:11.5}, 'r2a':
{r20_key_1:10, r20_key_2:11.5}, 'r2b': {r20_key_1:10, r20_key_2:11.5}, 'kex':
1000, 'pA': 0.99, 'dw': 2} ],
+ ['Ala', 2, 'N', {'r2': {r20_key_1:13, r20_key_2:14.5}, 'r2a':
{r20_key_1:13, r20_key_2:14.5}, 'r2b': {r20_key_1:13, r20_key_2:14.5}, 'kex':
1000, 'pA': 0.99, 'dw': 1} ]
]
# Collect the data to be used.
@@ -904,15 +915,18 @@
ds.insignificance = 0.0
# The grid search size (the number of increments per dimension).
- ds.GRID_INC = 11
+ ds.GRID_INC = 13
+
+ # The do clustering.
+ ds.do_cluster = True
# The function tolerance. This is used to terminate minimisation
once the function value between iterations is less than the tolerance.
# The default value is 1e-25.
- ds.set_func_tol = 1e-15
+ ds.set_func_tol = 1e-8
# The maximum number of iterations.
# The default value is 1e7.
- ds.set_max_iter = 1000000
+ ds.set_max_iter = 10000
# The verbosity level.
ds.verbosity = 1
@@ -922,6 +936,12 @@
# The plot_curves.
ds.plot_curves = False
+
+ # The conversion for ShereKhan at http://sherekhan.bionmr.org/.
+ ds.sherekhan_input = False
+
+ # Make a dx map to be opened om OpenDX. To map the hypersurface of
chi2, when altering kex, dw and pA.
+ ds.opendx = False
# The set r2eff err.
ds.r2eff_err = 0.1
@@ -940,7 +960,9 @@
cur_spin = return_spin(cur_spin_id)
grid_params = ds.grid_results[i][3]
- min_params = ds.min_results[i][3]
+
+ # Extract the clust results.
+ min_params = ds.clust_results[i][3]
# Now read the parameters.
print("For spin: '%s'"%cur_spin_id)
for mo_param in cur_spin.params:
@@ -961,7 +983,7 @@
print("###################################")
## Make test on R2.
- self.assertAlmostEqual(set_r2_frq/10, min_r2_frq/10,
1)
+ self.assertAlmostEqual(set_r2_frq, min_r2_frq, 1)
else:
grid_val = grid_params[mo_param]
min_val = min_params[mo_param]
@@ -972,7 +994,7 @@
print("WARNING: rel change level is above %.2f, and
is %.4f."%(ds.rel_change, rel_change))
print("###################################")
- ## Make test on parameters.
+ ## Make test on parameters. Only if breaking the
relative change.
self.assertAlmostEqual(set_val, min_val, 1)
Modified: trunk/test_suite/system_tests/scripts/relax_disp/cpmg_synthetic.py
URL:
http://svn.gna.org/viewcvs/relax/trunk/test_suite/system_tests/scripts/relax_disp/cpmg_synthetic.py?rev=22952&r1=22951&r2=22952&view=diff
==============================================================================
--- trunk/test_suite/system_tests/scripts/relax_disp/cpmg_synthetic.py
(original)
+++ trunk/test_suite/system_tests/scripts/relax_disp/cpmg_synthetic.py Sun
May 4 22:41:30 2014
@@ -7,11 +7,11 @@
# relax module imports.
from auto_analyses.relax_disp import Relax_disp
-from lib.io import mkdir_nofail, open_write_file
+from lib.io import open_write_file
from data_store import Relax_data_store; ds = Relax_data_store()
-from pipe_control.mol_res_spin import return_spin, spin_loop
-from specific_analyses.relax_disp.data import generate_r20_key,
loop_exp_frq, loop_offset, loop_frq, loop_offset_point
-from specific_analyses.relax_disp.variables import EXP_TYPE_CPMG_SQ
+from pipe_control.mol_res_spin import return_spin
+from specific_analyses.relax_disp.data import generate_r20_key,
loop_exp_frq, loop_offset_point
+from specific_analyses.relax_disp.variables import EXP_TYPE_CPMG_SQ,
MODEL_PARAMS
from specific_analyses.relax_disp import optimisation
from status import Status; status = Status()
@@ -28,6 +28,7 @@
time_T2_1 = 0.06
ncycs_1 = [2, 4, 8, 10, 20, 30, 40, 60]
r2eff_errs_1 = [0.05, -0.05, 0.05, -0.05, 0.05, -0.05, 0.05, -0.05]
+ #r2eff_errs_1 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
exp_1 = [sfrq_1, time_T2_1, ncycs_1, r2eff_errs_1]
sfrq_2 = 499.8908617*1E6
@@ -35,6 +36,7 @@
time_T2_2 = 0.05
ncycs_2 = [2, 4, 8, 10, 30, 35, 40, 50]
r2eff_errs_2 = [0.05, -0.05, 0.05, -0.05, 0.05, -0.05, 0.05, -0.05]
+ #r2eff_errs_2 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
exp_2 = [sfrq_2, time_T2_2, ncycs_2, r2eff_errs_2]
# Collect all exps
@@ -63,19 +65,24 @@
if not hasattr(ds, 'insignificance'):
ds.insignificance = 0.0
-# The grid search size (the number of increments per dimension).
+# The grid search size (the number of increments per dimension). "None" will
set default values.
if not hasattr(ds, 'GRID_INC'):
- ds.GRID_INC = 12
+ #ds.GRID_INC = None
+ ds.GRID_INC = 13
+
+# The do clustering.
+if not hasattr(ds, 'do_cluster'):
+ ds.do_cluster = False
# The function tolerance. This is used to terminate minimisation once the
function value between iterations is less than the tolerance.
# The default value is 1e-25.
if not hasattr(ds, 'set_func_tol'):
- ds.set_func_tol = 1e-25
+ ds.set_func_tol = 1e-8
# The maximum number of iterations.
# The default value is 1e7.
if not hasattr(ds, 'set_max_iter'):
- ds.set_max_iter = 10000000
+ ds.set_max_iter = 10000
# The verbosity level.
if not hasattr(ds, 'verbosity'):
@@ -88,6 +95,15 @@
# The plot_curves.
if not hasattr(ds, 'plot_curves'):
ds.plot_curves = True
+
+# The conversion for ShereKhan at http://sherekhan.bionmr.org/.
+if not hasattr(ds, 'sherekhan_input'):
+ ds.sherekhan_input = True
+
+# Make a dx map to be opened om OpenDX.
+# To map the hypersurface of chi2, when altering kex, dw and pA.
+if not hasattr(ds, 'opendx'):
+ ds.opendx = True
# The set r2eff err.
if not hasattr(ds, 'r2eff_err'):
@@ -109,7 +125,6 @@
pipe_type = 'relax_disp'
pipe_bundle = 'relax_disp'
pipe_name_r2eff = "%s_%s_R2eff"%(ds.data[0], pipe_name)
-#pipe_name_r2eff_calc = "%s_calc"%(pipe_name_r2eff)
pipe.create(pipe_name=pipe_name, pipe_type=pipe_type, bundle = pipe_bundle)
@@ -256,15 +271,31 @@
# Then select model.
relax_disp.select_model(model=ds.data[0])
-# Set the R20 parameters in the default grid search using the minimum R2eff
value.
-if ds.set_grid_r20_from_min_r2eff:
- relax_disp.set_grid_r20_from_min_r2eff(force=False)
+# Do a dx map.
+# To map the hypersurface of chi2, when altering kex, dw and pA.
+if ds.opendx:
+ dx.map(params=['dw', 'pA', 'kex'], map_type='Iso3D', spin_id=":1@N",
inc=20, lower=None, upper=None, axis_incs=5, file_prefix='map',
dir=ds.resdir, point=None, point_file='point', remap=None)
# Remove insignificant
relax_disp.insignificance(level=ds.insignificance)
# Perform Grid Search.
-grid_search(lower=None, upper=None, inc=ds.GRID_INC, constraints=True,
verbosity=ds.verbosity)
+if ds.GRID_INC:
+ # Set the R20 parameters in the default grid search using the minimum
R2eff value.
+ # This speeds it up considerably.
+ if ds.set_grid_r20_from_min_r2eff:
+ relax_disp.set_grid_r20_from_min_r2eff(force=False)
+
+ # Then do grid search.
+ grid_search(lower=None, upper=None, inc=ds.GRID_INC, constraints=True,
verbosity=ds.verbosity)
+
+# If no Grid search, set the default values.
+else:
+ for param in MODEL_PARAMS[ds.data[0]]:
+ value.set(param=param, index=None)
+ # Do a grid search, which will store the chi2 value.
+ #grid_search(lower=None, upper=None, inc=10, constraints=True,
verbosity=ds.verbosity)
+
# Define function to store grid results.
def save_res(res_spins):
@@ -286,9 +317,38 @@
ds.grid_results = save_res(cur_spins)
## Now do minimisation.
+
minimise(min_algor='simplex', func_tol=ds.set_func_tol,
max_iter=ds.set_max_iter, constraints=True, scaling=True,
verbosity=ds.verbosity)
+# Save results
ds.min_results = save_res(cur_spins)
+
+# Now do clustering
+if ds.do_cluster:
+ # Change pipe.
+ pipe_name_MODEL_CLUSTER = "%s_%s_CLUSTER"%(pipe_name, ds.data[0])
+ pipe.copy(pipe_from=pipe_name, pipe_to=pipe_name_MODEL_CLUSTER)
+ pipe.switch(pipe_name=pipe_name_MODEL_CLUSTER)
+
+ # Copy R2eff, but not the original parameters
+ value.copy(pipe_from=pipe_name_r2eff, pipe_to=pipe_name_MODEL_CLUSTER,
param='r2eff')
+
+ # Then select model.
+ relax_disp.select_model(model=ds.data[0])
+
+ # Then cluster
+ relax_disp.cluster('model_cluster', ":1-100")
+
+ # Copy the parameters from before.
+ relax_disp.parameter_copy(pipe_from=pipe_name_MODEL,
pipe_to=pipe_name_MODEL_CLUSTER)
+
+ # Now minimise.
+ minimise(min_algor='simplex', func_tol=ds.set_func_tol,
max_iter=ds.set_max_iter, constraints=True, scaling=True,
verbosity=ds.verbosity)
+
+ # Save results
+ ds.clust_results = save_res(cur_spins)
+else:
+ ds.clust_results = ds.min_results
# Compare results.
for i in range(len(cur_spins)):
@@ -298,22 +358,26 @@
grid_params = ds.grid_results[i][3]
min_params = ds.min_results[i][3]
+ clust_params = ds.clust_results[i][3]
# Now read the parameters.
- print("For spin: '%s'"%cur_spin_id)
+ if ds.print_res:
+ print("For spin: '%s'"%cur_spin_id)
for mo_param in cur_spin.params:
# The R2 is a dictionary, depending on spectrometer frequency.
if isinstance(getattr(cur_spin, mo_param), dict):
grid_r2 = grid_params[mo_param]
min_r2 = min_params[mo_param]
+ clust_r2 = clust_params[mo_param]
set_r2 = params[mo_param]
for key, val in set_r2.items():
grid_r2_frq = grid_r2[key]
min_r2_frq = min_r2[key]
+ clust_r2_frq = min_r2[key]
set_r2_frq = set_r2[key]
frq = float(key.split(EXP_TYPE_CPMG_SQ+' -
')[-1].split('MHz')[0])
- rel_change = sqrt( (min_r2_frq -
set_r2_frq)**2/(min_r2_frq)**2 )
+ rel_change = sqrt( (clust_r2_frq -
set_r2_frq)**2/(clust_r2_frq)**2 )
if ds.print_res:
- print("%s %s %s %s %.1f GRID=%.3f MIN=%.3f SET=%.3f
RELC=%.3f"%(cur_spin.model, res_name, cur_spin_id, mo_param, frq,
grid_r2_frq, min_r2_frq, set_r2_frq, rel_change) )
+ print("%s %s %s %s %.1f GRID=%.3f MIN=%.3f CLUST=%.3f
SET=%.3f RELC=%.3f"%(cur_spin.model, res_name, cur_spin_id, mo_param, frq,
grid_r2_frq, min_r2_frq, clust_r2_frq, set_r2_frq, rel_change) )
if rel_change > ds.rel_change:
if ds.print_res:
print("WARNING: rel change level is above %.2f, and
is %.4f."%(ds.rel_change, rel_change))
@@ -321,10 +385,11 @@
else:
grid_val = grid_params[mo_param]
min_val = min_params[mo_param]
+ clust_val = clust_params[mo_param]
set_val = params[mo_param]
- rel_change = sqrt( (min_val - set_val)**2/(min_val)**2 )
+ rel_change = sqrt( (clust_val - set_val)**2/(clust_val)**2 )
if ds.print_res:
- print("%s %s %s %s GRID=%.3f MIN=%.3f SET=%.3f
RELC=%.3f"%(cur_spin.model, res_name, cur_spin_id, mo_param, grid_val,
min_val, set_val, rel_change) )
+ print("%s %s %s %s GRID=%.3f MIN=%.3f CLUST=%.3f SET=%.3f
RELC=%.3f"%(cur_spin.model, res_name, cur_spin_id, mo_param, grid_val,
min_val, clust_val, set_val, rel_change) )
if rel_change > ds.rel_change:
if ds.print_res:
print("WARNING: rel change level is above %.2f, and is
%.4f."%(ds.rel_change, rel_change))
@@ -333,3 +398,9 @@
# Plot curves.
if ds.plot_curves:
relax_disp.plot_disp_curves(dir=ds.resdir, force=True)
+
+# The conversion for ShereKhan at http://sherekhan.bionmr.org/.
+if ds.sherekhan_input:
+ relax_disp.cluster('sherekhan', ":1-100")
+ print(cdp.clustering)
+ relax_disp.sherekhan_input(force=True, spin_id=None, dir=ds.resdir)
r22952 - in /trunk/test_suite/system_tests: relax_disp.py scripts/relax_disp/cpmg_synthetic.py