Re: r22952 - in /trunk/test_suite/system_tests: relax_disp.py scripts/relax_disp/cpmg_synthetic.py -- May 05, 2014 - 09:55

Hi Troels,

I have two suggestions.  The first is that it would be useful to
branch the tests into two - one noise-free, the other with the
r2eff_errs_1 = [0.05, -0.05, 0.05, -0.05, 0.05, -0.05, 0.05, -0.05]
errors.  The more tests the better and noise-free is good for testing
how perfect an implementation is and randomised is good for testing
the robustness of the model/implementation.  Though I cannot see how
r2eff_errs_1 are used in the test yet.

A different point is that the system test script is quite large now at
400+ lines.  To avoid tests from becoming too large, it is a good idea
to produce many small tests which each check one feature rather than
one large test covering everything.  This helps in debugging,
especially when large code refactorisations occur in the future, as
then you see a series of failures simultaneously.  This is better than
seeing one test failing, fixing it, seeing it fail at another point,
fixing that, then seeing another failure.  A simple way to do this is
to have a script in the test_suite/shared_data directories which sets
everything up, maybe even preforms some minimisation, and then creates
a relax state file.  That state file is then added to the repository
and the many small tests then load it and run from there.  I've used
this trick for many of the other analysis types.  There is also not
enough relax state testing for relaxation dispersion, meaning that
loading old dispersion state file in newer relax versions might break
in the future.

Cheers,

Edward



On 4 May 2014 22:41,  <tlinnet@xxxxxxxxxxxxx> wrote:
> 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)
>
>
> _______________________________________________
> relax (http://www.nmr-relax.com)
>
> This is the relax-commits mailing list
> relax-commits@xxxxxxx
>
> To unsubscribe from this list, get a password
> reminder, or change your subscription options,
> visit the list information page at
> https://mail.gna.org/listinfo/relax-commits