mailRe: r22933 - /trunk/test_suite/system_tests/relax_disp.py


Others Months | Index by Date | Thread Index
>>   [Date Prev] [Date Next] [Thread Prev] [Thread Next]

Header


Content

Posted by Edward d'Auvergne on May 02, 2014 - 15:16:
Hi,

From the graphs, the fits nevertheless look reasonable.  The residuals
are all within 2 standard deviations from zero.  This is therefore a
statistically significant, hence reasonable fit.  This tells me that
the problem is one which is inherent to relaxation dispersion itself.
For your own sanity, I would recommend you try a few things:

- Use this synthetic data with other dispersion software and see if
the same result is obtained.  As this is the CR72 model, you should
try ShereKhan.  Maybe also try the 'MMQ CR72' model and then compare
to cpmg_fit.  Doing this, I think you will find that this is not a
problem with relax.

- Map the space to see if two or more minima are present.  Try the
dx.map user function to obtain 3D isosurfaces for 3 of the parameters,
fixing the other parameters to the correct values.  It could be that
there is one minimum with pA = 0.5.  Or more likely, the errors you
have added have crushed a shallow minimum at the real parameter values
so much that it no longer exists.  This 3D visualisation of the space
is incredibly powerful, you will learn a lot about the problem.  See
my model-free model elimination paper
(http://dx.doi.org/10.1007/s10858-006-9007-z) and model-free
minimisation paper (http://dx.doi.org/10.1007/s10858-007-9214-2) for
examples.  Some Linux distributions bundle opendx
(https://packages.debian.org/search?keywords=dx), so setting up a
VirtualBox image might be the easiest way to run this on a Mac.

Does this fail if you do not add noise?

Regards,

Edward

On 2 May 2014 14:57, Troels Emtekær Linnet <tlinnet@xxxxxxxxxxxxx> wrote:
Hi Edward.

I just added this system test:

relax -s Relax_disp.test_cpmg_synthetic_cr72_fail -d

I am quite concerned about this.

It proves that small R2eff errors of 0.05, will make the minimisation
go terrible wrong.

For a spin with dw=2, it goes okay, but with dw=1, it finds pA=0.5
instead of 0.99.
For spin: ':1@N'
CR72 Ala :1@N r2 599.9 GRID=10.067 MIN=9.985 SET=10.000 RELC=0.002
CR72 Ala :1@N r2 499.9 GRID=11.531 MIN=11.496 SET=11.500 RELC=0.000
CR72 Ala :1@N pA GRID=0.955 MIN=0.989 SET=0.990 RELC=0.001
CR72 Ala :1@N dw GRID=0.909 MIN=1.931 SET=2.000 RELC=0.036
CR72 Ala :1@N kex GRID=1819.000 MIN=1027.289 SET=1000.000 RELC=0.027

For spin: ':2@N'
CR72 Ala :2@N r2 599.9 GRID=12.979 MIN=13.073 SET=13.000 RELC=0.006
CR72 Ala :2@N r2 499.9 GRID=14.470 MIN=14.477 SET=14.500 RELC=0.002
CR72 Ala :2@N pA GRID=0.500 MIN=0.500 SET=0.990 RELC=0.980
WARNING: rel change level is above 0.05, and is 0.9800.
###################################

CR72 Ala :2@N dw GRID=6.364 MIN=3.505 SET=1.000 RELC=0.715
WARNING: rel change level is above 0.05, and is 0.7147.
###################################

CR72 Ala :2@N kex GRID=1.000 MIN=1.885 SET=1000.000 RELC=529.556
WARNING: rel change level is above 0.05, and is 529.5564.
###################################

Try also this, to get the graphs in a test folder.

   relax test_suite/system_tests/scripts/relax_disp/cpmg_synthetic.py

This really looks like my problem for my SOD1 WT dataset.

We think that all of the dw values are 2x to big.
So, I am quite suspicious for the minimisation of dw values.


Ughhh....


Best
Troels




---------- Forwarded message ----------
From:  <tlinnet@xxxxxxxxxxxxx>
Date: 2014-05-02 14:47 GMT+02:00
Subject: r22933 - /trunk/test_suite/system_tests/relax_disp.py
To: relax-commits@xxxxxxx


Author: tlinnet
Date: Fri May  2 14:47:19 2014
New Revision: 22933

URL: http://svn.gna.org/viewcvs/relax?rev=22933&view=rev
Log:
Added a systemtest, which proves that small dw values of 1, makes the
mimisation goes wrong.

This is for synthetic data with R2eff values of +/- 0.05, which is to
be expected for real data.

Modified:
    trunk/test_suite/system_tests/relax_disp.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=22933&r1=22932&r2=22933&view=diff
==============================================================================
--- trunk/test_suite/system_tests/relax_disp.py (original)
+++ trunk/test_suite/system_tests/relax_disp.py Fri May  2 14:47:19 2014
@@ -760,6 +760,129 @@

                         ## 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.
+
+        This script will produce synthetic CPMG R2eff values
according to the selected model, and the fit the selected model.
+        """
+
+        # Reset.
+        #self.interpreter.reset()
+
+        ## Set Experiments.
+        model = "CR72"
+        # Exp 1
+        sfrq_1 = 599.8908617*1E6
+        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]
+
+        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]
+
+        # 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} ]
+            ]
+
+        # Collect the data to be used.
+        ds.data = [model, spins, exps]
+
+        # The tmp directory. None is the local directory.
+        ds.tmpdir = ds.tmpdir
+
+        # The results directory. None is the local directory.
+        #ds.resdir = None
+        ds.resdir = ds.tmpdir
+
+        # Do set_grid_r20_from_min_r2eff ?.
+        ds.set_grid_r20_from_min_r2eff = True
+
+        # Remove insignificant level.
+        ds.insignificance = 0.0
+
+        # The grid search size (the number of increments per dimension).
+        ds.GRID_INC = 11
+
+        # 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
+
+        # The maximum number of iterations.
+        # The default value is 1e7.
+        ds.set_max_iter = 1000000
+
+        # The verbosity level.
+        ds.verbosity = 1
+
+        # The rel_change WARNING level.
+        ds.rel_change = 0.05
+
+        # The plot_curves.
+        ds.plot_curves = False
+
+        # The set r2eff err.
+        ds.r2eff_err = 0.1
+
+        # The print result info.
+        ds.print_res = False
+
+        # Execute the script.
+        self.interpreter.run(script_file=status.install_path +
sep+'test_suite'+sep+'system_tests'+sep+'scripts'+sep+'relax_disp'+sep+'cpmg_synthetic.py')
+
+        cur_spins = ds.data[1]
+        # Compare results.
+        for i in range(len(cur_spins)):
+            res_name, res_num, spin_name, params = cur_spins[i]
+            cur_spin_id = ":%i@%s"%(res_num, spin_name)
+            cur_spin = return_spin(cur_spin_id)
+
+            grid_params = ds.grid_results[i][3]
+            min_params = ds.min_results[i][3]
+            # Now read the parameters.
+            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]
+                    set_r2 = params[mo_param]
+                    for key, val in set_r2.items():
+                        grid_r2_frq = grid_r2[key]
+                        min_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 = math.sqrt( (min_r2_frq -
set_r2_frq)**2/(min_r2_frq)**2 )
+                        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) )
+                        if rel_change > ds.rel_change:
+                            print("WARNING: rel change level is above
%.2f, and is %.4f."%(ds.rel_change, rel_change))
+                            print("###################################")
+
+                        ## Make test on R2.
+                        self.assertAlmostEqual(set_r2_frq/10, 
min_r2_frq/10, 1)
+                else:
+                    grid_val = grid_params[mo_param]
+                    min_val = min_params[mo_param]
+                    set_val = params[mo_param]
+                    rel_change = math.sqrt( (min_val -
set_val)**2/(min_val)**2 )
+                    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) )
+                    if rel_change > ds.rel_change:
+                        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, 1)


     def test_curve_type_cpmg_fixed_time(self):


_______________________________________________
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

_______________________________________________
relax (http://www.nmr-relax.com)

This is the relax-devel mailing list
relax-devel@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-devel



Related Messages


Powered by MHonArc, Updated Fri May 02 16:40:08 2014