Difference between revisions of "Dx map"

From relax wiki
Jump to navigation Jump to search
m (Switch to the {{gna bug url}} and {{gna bug url}} templates to remove dead Gna! links.)
 
(19 intermediate revisions by 3 users not shown)
Line 1: Line 1:
 +
__TOC__
 +
 +
== Sample data to try dx ==
 +
<source lang="bash">
 +
cd $HOME/Desktop
 +
mkdir test_dx
 +
cd test_dx
 +
 +
# See sample scripts
 +
ls -la $HOME/software/relax/sample_scripts/model_free
 +
cp $HOME/software/relax/sample_scripts/model_free/map.py .
 +
 +
# Get data
 +
ls -la $HOME/software/relax/test_suite/shared_data/model_free/sphere
 +
cp $HOME/software/relax/test_suite/shared_data/model_free/sphere/*.out .
 +
cp $HOME/software/relax/test_suite/shared_data/model_free/sphere/*.pdb .
 +
 +
# Modify map script
 +
cp map.py map_mod.py
 +
sed -i.bak "s/res_num_col=1)/mol_name_col=1, res_num_col=2, res_name_col=3, spin_num_col=4, spin_name_col=5)/g" map_mod.py
 +
sed -i.bak "s/600/900/g" map_mod.py
 +
sed -i.bak "s/spin\.name/#spin\.name/g" map_mod.py
 +
sed -i.bak "s/res_num_col=1, data_col=3, error_col=4)/mol_name_col=1, res_num_col=2, res_name_col=3, spin_num_col=4, spin_name_col=5, data_col=6, error_col=7)/g" map_mod.py
 +
sed -i.bak "s/sequence\.attach_protons/#sequence\.attach_protons/g" map_mod.py
 +
sed -i.bak "spin/a spin" map_mod.py
 +
 +
# Add these two lines to map_mod.py
 +
spin.element(element='H', spin_id='@H')
 +
spin.isotope('1H', spin_id='@H')
 +
 +
# Run
 +
relax map_mod.py
 +
</source>
 +
 
== Code to generate map ==
 
== Code to generate map ==
 
<source lang="Python">
 
<source lang="Python">
Line 8: Line 42:
 
for cur_spin, mol_name, resi, resn, spin_id in spin_loop(full_info=True, return_id=True, skip_desel=True):
 
for cur_spin, mol_name, resi, resn, spin_id in spin_loop(full_info=True, return_id=True, skip_desel=True):
 
     file_name = "map%s" % (cur_spin_id .replace('#', '_').replace(':', '_').replace('@', '_'))
 
     file_name = "map%s" % (cur_spin_id .replace('#', '_').replace(':', '_').replace('@', '_'))
     dx.map(params=['dw', 'pA', 'kex'], map_type='Iso3D', spin_id=":1@N", inc=70, lower=None, upper=None, axis_incs=5, file_prefix=file_name, dir=ds.resdir, point=None, point_file='point', remap=None)
+
     dx.map(params=['dw', 'pA', 'kex'], map_type='Iso3D', spin_id=":1@N", inc=70, lower=None, upper=None, axis_incs=5,
 +
        file_prefix=file_name, dir=ds.resdir, point=None, point_file='point', remap=None)
 
     #vp_exec:  A flag specifying whether to execute the visual program automatically at start-up.
 
     #vp_exec:  A flag specifying whether to execute the visual program automatically at start-up.
 
     dx.execute(file_prefix=file_name, dir=ds.resdir, dx_exe='dx', vp_exec=True)
 
     dx.execute(file_prefix=file_name, dir=ds.resdir, dx_exe='dx', vp_exec=True)
Line 17: Line 52:
  
 
== How to use dx ==
 
== How to use dx ==
* Run 'dx',
+
Run <code>dx</code>, and then in the '''Data explorer''' or '''DE'''.
In teh '''Data explorer''' or '''DE'''.
+
* Click on <code>Edit Visual Programs...</code>.
* Click on '''Edit Visual Programs...'''.
+
* Select the <code>map.net</code> program created by relax,
* Select the map.net program created by relax,
 
  
Now in the '''Visual Program Editor''' or '''VPE'''.
+
Now in the <code>Visual Program Editor</code> or <code>VPE</code>.
* Select the menu entry 'Execute->Execute on change'.
+
* Select the menu entry <code>Execute->Execute on change</code>.
  
 
That's it.
 
That's it.
  
You now have a 3D frame, but nothing in it.  <br>
+
You now have a 3D frame, but nothing in it.  Therefore the contour levels must be too low or high. From the map file, the values are in the hundreds of thousands.
Therefore the contour levels must be too low or high. <br>
 
From the map file, the values are in the hundreds of thousands.
 
  
 
Then:
 
Then:
* In the main program window, double click on the 'Isosurface elements'.
+
* In the main program window, double click on the <code>Isosurface elements</code>.
 
* Change the values until you see surfaces.  In the first the value is 500.  I changed this to 500,000. Multiply all by 1000.
 
* Change the values until you see surfaces.  In the first the value is 500.  I changed this to 500,000. Multiply all by 1000.
 
* In the second, 100 -> 100000.
 
* In the second, 100 -> 100000.
Line 40: Line 72:
 
This should maybe be performed by the dx.map user function, determining reasonable contour levels.   
 
This should maybe be performed by the dx.map user function, determining reasonable contour levels.   
  
With a bit of zooming, clicking on 'File -> Save image' in the "Surface" window, "allowing rendering", and outputting to a large TIFF file, "save current", then "apply".
+
With a bit of zooming, clicking on <code>File -> Save image</code> in the <code>Surface</code> window, <code>allowing rendering</code>, and outputting to a large TIFF file, <code>save current</code>, then <code>apply</code>.
  
An example image cropped and converted to PNG in the GIMP at <br>
+
An example image cropped and converted to PNG in the GIMP at https://gna.org/bugs/download.php?file_id=20641.  
https://gna.org/bugs/download.php?file_id=20641.  
 
  
Note that for a good resolution plot, you will need many more increments. <br>
+
Note that for a good resolution plot, you will need many more increments. Using the lower and upper dx.map arguments will be useful to zoom into the space.
Using the lower and upper dx.map arguments will be useful to zoom into the space.
 
  
 
== Example images ==
 
== Example images ==
These images is from [https://gna.org/bugs/index.php?22024 bug 22024Minimisation space for CR72 is catastrophic. The chi2 surface over dw and pA is bounded.]
+
These images is from {{gna bug link|22024|text=bug 22024.  Minimisation space for CR72 is catastrophic. The chi2 surface over dw and pA is bounded.}}
 +
 
 +
{|style="margin: 0 auto;"
 +
| [[File:Bug 22024 R2eff.png|thumb|center|upright=2|{{gna bug url|22024}} Bug 22024: Minimisation space for CR72 is catastrophic. The chi2 surface over dw and pA is bounded.]]
 +
| [[File:Bug 22024 Dw kex.png|thumb|upright=2|{{gna bug url|22024}} Bug 22024: Minimisation space for CR72 is catastrophic. The chi2 surface over dw and pA is bounded.]]
 +
|}
  
[[File:Bug 22024 R2eff.png|thumb|upright=2|https://gna.org/bugs/index.php?22024Bug 22024: Minimisation space for CR72 is catastrophic. The chi2 surface over dw and pA is bounded.]]
+
{|style="margin: 0 auto;"
[[File:Bug 22024 Dw kex.png|thumb|upright=2|https://gna.org/bugs/index.php?22024Bug 22024: Minimisation space for CR72 is catastrophic. The chi2 surface over dw and pA is bounded.]]
+
| [[File:Bug 22024 Dw kex pA.png|thumb|upright=2|{{gna bug url|22024}} Bug 22024: Minimisation space for CR72 is catastrophic. The chi2 surface over dw and pA is bounded.]]
[[File:Bug 22024 Dw kex pA.png|thumb|upright=2|https://gna.org/bugs/index.php?22024Bug 22024: Minimisation space for CR72 is catastrophic. The chi2 surface over dw and pA is bounded.]]
+
| [[File:Bug 22024 Dw pA.png|thumb|upright=2|{{gna bug url|22024}} Bug 22024: Minimisation space for CR72 is catastrophic. The chi2 surface over dw and pA is bounded.]]
[[File:Bug 22024 Dw pA.png|thumb|upright=2|https://gna.org/bugs/index.php?22024Bug 22024: Minimisation space for CR72 is catastrophic. The chi2 surface over dw and pA is bounded.]]
+
| [[File:Bug 22024 Kex pA.png|thumb|upright=2|{{gna bug url|22024}} Bug 22024: Minimisation space for CR72 is catastrophic. The chi2 surface over dw and pA is bounded.]]
[[File:Bug 22024 Kex pA.png|thumb|upright=2|https://gna.org/bugs/index.php?22024Bug 22024: Minimisation space for CR72 is catastrophic. The chi2 surface over dw and pA is bounded.]]
+
|}
  
 
== Code to generate map ==
 
== Code to generate map ==
Line 64: Line 99:
 
relax cpmg_synthetic.py
 
relax cpmg_synthetic.py
 
</source>
 
</source>
Here it the code
+
 
<source lang="Python">
+
Here is the synthetic data generator code:
 +
 
 +
{{collapsible script
 +
| type  = relax script
 +
| title  = Script for calculating synthetics CPMG data.
 +
| lang   = python
 +
| script =
 
###############################################################################
 
###############################################################################
 
#                                                                            #
 
#                                                                            #
Line 92: Line 133:
 
from tempfile import mkdtemp
 
from tempfile import mkdtemp
 
from math import sqrt
 
from math import sqrt
 +
import sys
  
 
# relax module imports.
 
# relax module imports.
Line 113: Line 155:
  
 
# Analysis variables.
 
# Analysis variables.
#####################
+
##################################################################################
 
# The dispersion model to test.
 
# The dispersion model to test.
 
if not hasattr(ds, 'data'):
 
if not hasattr(ds, 'data'):
Line 172: Line 214:
 
     # Add more spins here.
 
     # Add more spins here.
 
     spins = [
 
     spins = [
             ['Ala', 1, 'N', {'r2': {r20_key_1: 25.0, r20_key_2: 24.0}, 'r2a': {r20_key_1: 25.0, r20_key_2: 24.0}, 'r2b': {r20_key_1: 25.0, r20_key_2: 24.0}, 'kex': 2200.0, 'pA': 0.993, 'dw': 1.5} ]
+
             ['Ala', 1, 'N', {'r2': {r20_key_1: 25.0, r20_key_2: 24.0}, 'r2a': {r20_key_1: 25.0, r20_key_2: 24.0}, 'r2b': {r20_key_1: 25.0, r20_key_2: 24.0}, 'kex': 2200.0, 'pA': 0.993, 'dw': 0.5} ]
 
             #['Ala', 2, 'N', {'r2': {r20_key_1: 15.0, r20_key_2: 14.0}, 'r2a': {r20_key_1: 15.0, r20_key_2: 14.0}, 'r2b': {r20_key_1: 15.0, r20_key_2: 14.0}, 'kex': 2200.0, 'pA': 0.993, 'dw': 1.5} ]
 
             #['Ala', 2, 'N', {'r2': {r20_key_1: 15.0, r20_key_2: 14.0}, 'r2a': {r20_key_1: 15.0, r20_key_2: 14.0}, 'r2b': {r20_key_1: 15.0, r20_key_2: 14.0}, 'kex': 2200.0, 'pA': 0.993, 'dw': 1.5} ]
 
             ]
 
             ]
  
 +
    # Collect the eksperiment to simulate.
 
     ds.data = [model_create, model_analyse, spins, exps]
 
     ds.data = [model_create, model_analyse, spins, exps]
  
 +
# Setup variables to be used
 +
##################################################################################
 
# The tmp directory.
 
# The tmp directory.
 
if not hasattr(ds, 'tmpdir'):
 
if not hasattr(ds, 'tmpdir'):
Line 186: Line 231:
 
     ds.resdir = None
 
     ds.resdir = None
  
# Do set_grid_r20_from_min_r2eff ?.
+
# Do set_grid_r20_from_min_r2eff.
 
if not hasattr(ds, 'set_grid_r20_from_min_r2eff'):
 
if not hasattr(ds, 'set_grid_r20_from_min_r2eff'):
 
     ds.set_grid_r20_from_min_r2eff = True
 
     ds.set_grid_r20_from_min_r2eff = True
Line 223: Line 268:
 
# The plot_curves.
 
# The plot_curves.
 
if not hasattr(ds, 'plot_curves'):
 
if not hasattr(ds, 'plot_curves'):
     ds.plot_curves = True
+
     ds.plot_curves = False
  
 
# The conversion for ShereKhan at http://sherekhan.bionmr.org/.
 
# The conversion for ShereKhan at http://sherekhan.bionmr.org/.
 
if not hasattr(ds, 'sherekhan_input'):
 
if not hasattr(ds, 'sherekhan_input'):
 
     ds.sherekhan_input = False
 
     ds.sherekhan_input = False
 +
 +
# The set r2eff err, for defining the error on the graphs and in the fitting weight.
 +
# We set the error to be the same for all ncyc eksperiments.
 +
if not hasattr(ds, 'r2eff_err'):
 +
    ds.r2eff_err = 0.1
 +
 +
# The print result info.
 +
if not hasattr(ds, 'print_res'):
 +
    ds.print_res = True
  
 
# Make a dx map to be opened om OpenDX.
 
# Make a dx map to be opened om OpenDX.
 
# To map the hypersurface of chi2, when altering kex, dw and pA.
 
# To map the hypersurface of chi2, when altering kex, dw and pA.
 
if not hasattr(ds, 'opendx'):
 
if not hasattr(ds, 'opendx'):
 +
    #ds.opendx = False
 
     ds.opendx = True
 
     ds.opendx = True
if not hasattr(ds, 'dx_inc'):
+
 
    ds.dx_inc = 20
+
# Which parameters to map on chi2 surface.
 
if not hasattr(ds, 'dx_params'):
 
if not hasattr(ds, 'dx_params'):
 
     ds.dx_params = ['dw', 'pA', 'kex']
 
     ds.dx_params = ['dw', 'pA', 'kex']
  
# The set r2eff err.
+
# How many increements to map for each of the parameters.
if not hasattr(ds, 'r2eff_err'):
+
# Above 20 is good. 70 - 100 is a quite good resolution.
     ds.r2eff_err = 0.1
+
# 4 is for speed tests.
 +
if not hasattr(ds, 'dx_inc'):
 +
    #ds.dx_inc = 70
 +
    ds.dx_inc = 4
 +
 
 +
# If bounds set to None, uses the normal Grid search bounds.
 +
if not hasattr(ds, 'dx_lower_bounds'):
 +
     ds.dx_lower_bounds = None
 +
    #ds.dx_lower_bounds = [0.0, 0.5, 1.0 ]
 +
 
 +
# If bounds set to None, uses the normal Grid search bounds.
 +
if not hasattr(ds, 'dx_upper_bounds'):
 +
    #ds.dx_upper_bounds = None
 +
    ds.dx_upper_bounds = [10.0, 1.0, 3000.0]
 +
 
 +
# If dx_chi_surface is None, it will find Innermost, Inner, Middle and Outer Isosurface
 +
# at 10, 20, 50 and 90 percentile of all chi2 values.
 +
if not hasattr(ds, 'dx_chi_surface'):
 +
    ds.dx_chi_surface = None
 +
    #ds.dx_chi_surface = [1000.0, 20000.0, 30000., 30000.00]
 +
 
 +
##################################################################################
 +
 
 +
# Define repeating functions.
 +
##################################################################################
 +
# Define function to store grid results.
 +
def save_res(res_spins):
 +
    res_list = []
 +
    for res_name, res_num, spin_name, params in res_spins:
 +
        cur_spin_id = ":%i@%s"%(res_num, spin_name)
 +
        cur_spin = return_spin(cur_spin_id)
 +
 
 +
        # Save all the paramers to loop throgh later.
 +
        ds.model_analyse_params = cur_spin.params
 +
 
 +
        par_dic = {}
 +
        # Now read the parameters.
 +
        for mo_param in cur_spin.params:
 +
            par_dic.update({mo_param : getattr(cur_spin, mo_param) })
 +
 
 +
        # Append result.
 +
        res_list.append([res_name, res_num, spin_name, par_dic])
 +
 
 +
    return res_list
 +
 
 +
# Define function print relative change
 +
def print_res(cur_spins=None, grid_params=None, min_params=None, clust_params=None, dx_params=[]):
 +
    # Define list to hold data.
 +
    dx_set_val = list(range(len(dx_params)))
 +
    dx_clust_val = list(range(len(dx_params)))
 +
 
 +
    # Compare results.
 +
    if ds.print_res:
 +
        print("\n########################")
 +
        print("Generated data with MODEL:%s"%(model_create))
 +
        print("Analysing with MODEL:%s."%(model_analyse))
 +
        print("########################\n")
 +
 
 +
    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)
 +
 
 +
        # Now read the parameters.
 +
        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 min_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( (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 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("###################################")
 +
                            print("WARNING: %s Have relative change above %.2f, and is %.4f."%(key, ds.rel_change, rel_change))
 +
                            print("###################################\n")
 +
            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( (clust_val - set_val)**2/(clust_val)**2 )
 +
                if ds.print_res:
 +
                    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("###################################")
 +
                        print("WARNING: %s Have relative change above %.2f, and is %.4f."%(mo_param, ds.rel_change, rel_change))
 +
                        print("###################################\n")
  
# The number of Monte Carlo simulations to be used for the error analyses.
+
                # Store to dx map.
if not hasattr(ds, 'MC_NUM'):
+
                if mo_param in dx_params:
    ds.MC_NUM = 3
+
                    dx_set_val[ds.dx_params.index(mo_param)] = set_val
 +
                    dx_clust_val[ds.dx_params.index(mo_param)] = clust_val
  
# The print result info.
+
    return dx_set_val, dx_clust_val
if not hasattr(ds, 'print_res'):
 
    ds.print_res = True
 
  
 +
# Now starting
 +
##################################################################################
 
# Set up the data pipe.
 
# Set up the data pipe.
#######################
+
##################################################################################
  
 
# Extract the models
 
# Extract the models
Line 275: Line 430:
 
exps = ds.data[3]
 
exps = ds.data[3]
  
# Now loop over the experiments
+
# Now loop over the experiments, to set the variables in relax.
 
exp_ids = []
 
exp_ids = []
 
for exp in exps:
 
for exp in exps:
Line 306: Line 461:
 
pipe.switch(pipe_name=ds.pipe_name_r2eff)
 
pipe.switch(pipe_name=ds.pipe_name_r2eff)
  
# Then select model.
+
# Then select the model to create data.
 
relax_disp.select_model(model=model_create)
 
relax_disp.select_model(model=model_create)
  
Line 321: Line 476:
 
         if isinstance(getattr(cur_spin, mo_param), dict):
 
         if isinstance(getattr(cur_spin, mo_param), dict):
 
             set_r2 = params[mo_param]
 
             set_r2 = params[mo_param]
             for key, val in getattr(cur_spin, mo_param).items():
+
 
 +
             for key, val in set_r2.items():
 
                 # Update value to float
 
                 # Update value to float
 
                 set_r2.update({ key : float(val) })
 
                 set_r2.update({ key : float(val) })
Line 334: Line 490:
  
 
## Now doing the back calculation of R2eff values.
 
## Now doing the back calculation of R2eff values.
 
 
# First loop over the frequencies.
 
# First loop over the frequencies.
i = 0
 
 
for exp_type, frq, ei, mi in loop_exp_frq(return_indices=True):
 
for exp_type, frq, ei, mi in loop_exp_frq(return_indices=True):
 
     exp_id = exp_ids[mi]
 
     exp_id = exp_ids[mi]
Line 360: Line 514:
 
         file.close()
 
         file.close()
  
         # Read in the R2eff file to create the structure
+
         # Read in the R2eff file to create the structure.
 +
        # This is a trick, or else relax complains.
 
         relax_disp.r2eff_read_spin(id=exp_id, spin_id=cur_spin_id, file=file_name, dir=ds.tmpdir, disp_point_col=1, data_col=2, error_col=3)
 
         relax_disp.r2eff_read_spin(id=exp_id, spin_id=cur_spin_id, file=file_name, dir=ds.tmpdir, disp_point_col=1, data_col=2, error_col=3)
  
         ###Now back calculate, and stuff it back.
+
         ### Now back calculate values from parameters, and stuff R2eff it back.
 
         print("Generating data with MODEL:%s, for spin id:%s"%(model_create, cur_spin_id))
 
         print("Generating data with MODEL:%s, for spin id:%s"%(model_create, cur_spin_id))
 
         r2effs = optimisation.back_calc_r2eff(spin=cur_spin, spin_id=cur_spin_id)
 
         r2effs = optimisation.back_calc_r2eff(spin=cur_spin, spin_id=cur_spin_id)
Line 388: Line 543:
 
         # Read in the R2eff file to put into spin structure.
 
         # Read in the R2eff file to put into spin structure.
 
         relax_disp.r2eff_read_spin(id=exp_id, spin_id=cur_spin_id, file=file_name, dir=ds.resdir, disp_point_col=1, data_col=2, error_col=3)
 
         relax_disp.r2eff_read_spin(id=exp_id, spin_id=cur_spin_id, file=file_name, dir=ds.resdir, disp_point_col=1, data_col=2, error_col=3)
 
        # Add to counter.
 
        i += 1
 
 
print("Did following number of iterations: %i"%i)
 
  
 
### Now do fitting.
 
### Now do fitting.
Line 428: Line 578:
 
     #grid_search(lower=None, upper=None, inc=10, constraints=True, verbosity=ds.verbosity)
 
     #grid_search(lower=None, upper=None, inc=10, constraints=True, verbosity=ds.verbosity)
  
 
+
# Save result.
# Define function to store grid results.
 
def save_res(res_spins):
 
    res_list = []
 
    for res_name, res_num, spin_name, params in res_spins:
 
        cur_spin_id = ":%i@%s"%(res_num, spin_name)
 
        cur_spin = return_spin(cur_spin_id)
 
 
 
        par_dic = {}
 
        # Now read the parameters.
 
        for mo_param in cur_spin.params:
 
            par_dic.update({mo_param : getattr(cur_spin, mo_param) })
 
 
 
        # Append result.
 
        res_list.append([res_name, res_num, spin_name, par_dic])
 
 
 
    return res_list
 
 
 
 
ds.grid_results = save_res(cur_spins)
 
ds.grid_results = save_res(cur_spins)
  
 
## Now do minimisation.
 
## 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)
 
minimise(min_algor='simplex', func_tol=ds.set_func_tol, max_iter=ds.set_max_iter, constraints=True, scaling=True, verbosity=ds.verbosity)
  
Line 491: Line 623:
 
     print(cdp.clustering)
 
     print(cdp.clustering)
 
     relax_disp.sherekhan_input(force=True, spin_id=None, dir=ds.resdir)
 
     relax_disp.sherekhan_input(force=True, spin_id=None, dir=ds.resdir)
 
# Compare results.
 
if ds.print_res:
 
    print("\n########################")
 
    print("Generated data with MODEL:%s"%(model_create))
 
    print("Analysing with MODEL:%s."%(model_analyse))
 
    print("########################\n")
 
  
 
# Looping over data, to collect the results.
 
# Looping over data, to collect the results.
 
# Define which dx_params to collect for.
 
# Define which dx_params to collect for.
ds.dx_set_val = list(range(len(ds.dx_params)))
+
ds.dx_set_val, ds.dx_clust_val = print_res(cur_spins=cur_spins, grid_params=ds.grid_results[i][3], min_params=ds.min_results[i][3], clust_params=ds.clust_results[i][3], dx_params=ds.dx_params)
ds.dx_clust_val = list(range(len(ds.dx_params)))
 
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)
 
 
 
    # Fetch data.
 
    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.
 
    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 getattr(cur_spin, mo_param).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( (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 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("###################################")
 
                        print("WARNING: %s Have relative change above %.2f, and is %.4f."%(key, ds.rel_change, rel_change))
 
                        print("###################################\n")
 
        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( (clust_val - set_val)**2/(clust_val)**2 )
 
            if ds.print_res:
 
                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("###################################")
 
                    print("WARNING: %s Have relative change above %.2f, and is %.4f."%(mo_param, ds.rel_change, rel_change))
 
                    print("###################################\n")
 
 
 
            # Store to dx map.
 
            if mo_param in ds.dx_params:
 
                ds.dx_set_val[ds.dx_params.index(mo_param)] = set_val
 
                ds.dx_clust_val[ds.dx_params.index(mo_param)] = clust_val
 
  
 
## Do a dx map.
 
## Do a dx map.
Line 570: Line 642:
 
     relax_disp.select_model(model=model_analyse)
 
     relax_disp.select_model(model=model_analyse)
  
     # First loop over the defined spins and set the model parameters.
+
     # First loop over the defined spins and set the model parameters which is not in the dx.dx_params.
 
     for i in range(len(cur_spins)):
 
     for i in range(len(cur_spins)):
 
         res_name, res_num, spin_name, params = cur_spins[i]
 
         res_name, res_num, spin_name, params = cur_spins[i]
 
         cur_spin_id = ":%i@%s"%(res_num, spin_name)
 
         cur_spin_id = ":%i@%s"%(res_num, spin_name)
 
         cur_spin = return_spin(cur_spin_id)
 
         cur_spin = return_spin(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):
 +
                set_r2 = params[mo_param]
 +
                if mo_param not in ds.dx_params:
 +
                    for key, val in set_r2.items():
 +
                        # Update value to float
 +
                        set_r2.update({ key : float(val) })
 +
                        print("Setting param:%s to :%f"%(key, float(val)))
 +
                    # Set it back
 +
                    setattr(cur_spin, mo_param, set_r2)
 +
 +
            # For non dict values.
 +
            else:
 +
                if mo_param not in ds.dx_params:
 +
                    before = getattr(cur_spin, mo_param)
 +
                    setattr(cur_spin, mo_param, float(params[mo_param]))
 +
                    after = getattr(cur_spin, mo_param)
 +
                    print("Setting param:%s to :%f"%(mo_param, after))
 +
 +
        cur_model = model_analyse.replace(' ', '_')
 +
        file_name_map = "%s_map%s" % (cur_model, cur_spin_id.replace('#', '_').replace(':', '_').replace('@', '_'))
 +
        file_name_point = "%s_point%s" % (cur_model, cur_spin_id .replace('#', '_').replace(':', '_').replace('@', '_'))
 +
        dx.map(params=ds.dx_params, map_type='Iso3D', spin_id=cur_spin_id, inc=ds.dx_inc, lower=ds.dx_lower_bounds, upper=ds.dx_upper_bounds, axis_incs=10, file_prefix=file_name_map, dir=ds.resdir, point=[ds.dx_set_val, ds.dx_clust_val], point_file=file_name_point, chi_surface=ds.dx_chi_surface)
 +
        # vp_exec:  A flag specifying whether to execute the visual program automatically at start-up.
 +
        #dx.execute(file_prefix=file_name, dir=ds.resdir, dx_exe='dx', vp_exec=True)
  
 
         if ds.print_res:
 
         if ds.print_res:
 +
            print("For spin: '%s'"%cur_spin_id)
 
             print("Params for dx map is")
 
             print("Params for dx map is")
 
             print(ds.dx_params)
 
             print(ds.dx_params)
             print("Point param for dx map is")
+
             print("Point creating dx map is")
 
             print(ds.dx_set_val)
 
             print(ds.dx_set_val)
        cur_model = model_analyse.replace(' ', '_')
+
            print("Minimises point in dx map is")
        file_name_map = "%s_map%s" % (cur_model, cur_spin_id.replace('#', '_').replace(':', '_').replace('@', '_'))
+
            print(ds.dx_clust_val)
        file_name_point = "%s_point%s" % (cur_model, cur_spin_id .replace('#', '_').replace(':', '_').replace('@', '_'))
+
 
        dx.map(params=ds.dx_params, map_type='Iso3D', spin_id=cur_spin_id, inc=ds.dx_inc, lower=None, upper=None, axis_incs=10, file_prefix=file_name_map, dir=ds.resdir, point=[ds.dx_set_val, ds.dx_clust_val], point_file=file_name_point, remap=None)
+
# Print result again after dx map.
        #vp_exec:  A flag specifying whether to execute the visual program automatically at start-up.
+
if ds.opendx:
        #dx.execute(file_prefix=file_name, dir=ds.resdir, dx_exe='dx', vp_exec=True)
+
    print_res(cur_spins=cur_spins, grid_params=ds.grid_results[i][3], min_params=ds.min_results[i][3], clust_params=ds.clust_results[i][3])
</source>
+
}}
  
 
= See also =
 
= See also =
 
[[Install dx]]
 
[[Install dx]]
 
[[Category:User_functions]]
 
[[Category:User_functions]]

Latest revision as of 19:44, 16 October 2020

Sample data to try dx

cd $HOME/Desktop
mkdir test_dx
cd test_dx

# See sample scripts
ls -la $HOME/software/relax/sample_scripts/model_free
cp $HOME/software/relax/sample_scripts/model_free/map.py .

# Get data
ls -la $HOME/software/relax/test_suite/shared_data/model_free/sphere
cp $HOME/software/relax/test_suite/shared_data/model_free/sphere/*.out .
cp $HOME/software/relax/test_suite/shared_data/model_free/sphere/*.pdb .

# Modify map script
cp map.py map_mod.py
sed -i.bak "s/res_num_col=1)/mol_name_col=1, res_num_col=2, res_name_col=3, spin_num_col=4, spin_name_col=5)/g" map_mod.py
sed -i.bak "s/600/900/g" map_mod.py
sed -i.bak "s/spin\.name/#spin\.name/g" map_mod.py
sed -i.bak "s/res_num_col=1, data_col=3, error_col=4)/mol_name_col=1, res_num_col=2, res_name_col=3, spin_num_col=4, spin_name_col=5, data_col=6, error_col=7)/g" map_mod.py
sed -i.bak "s/sequence\.attach_protons/#sequence\.attach_protons/g" map_mod.py
sed -i.bak "spin/a spin" map_mod.py

# Add these two lines to map_mod.py
spin.element(element='H', spin_id='@H')
spin.isotope('1H', spin_id='@H')

# Run
relax map_mod.py

Code to generate map

# See relax help in prompt
help(dx.map)

from pipe_control.mol_res_spin import return_spin, spin_loop

for cur_spin, mol_name, resi, resn, spin_id in spin_loop(full_info=True, return_id=True, skip_desel=True):
    file_name = "map%s" % (cur_spin_id .replace('#', '_').replace(':', '_').replace('@', '_'))
    dx.map(params=['dw', 'pA', 'kex'], map_type='Iso3D', spin_id=":1@N", inc=70, lower=None, upper=None, axis_incs=5,
        file_prefix=file_name, dir=ds.resdir, point=None, point_file='point', remap=None)
    #vp_exec:  A flag specifying whether to execute the visual program automatically at start-up.
    dx.execute(file_prefix=file_name, dir=ds.resdir, dx_exe='dx', vp_exec=True)

How to install dx

See install dx

How to use dx

Run dx, and then in the Data explorer or DE.

  • Click on Edit Visual Programs....
  • Select the map.net program created by relax,

Now in the Visual Program Editor or VPE.

  • Select the menu entry Execute->Execute on change.

That's it.

You now have a 3D frame, but nothing in it. Therefore the contour levels must be too low or high. From the map file, the values are in the hundreds of thousands.

Then:

  • In the main program window, double click on the Isosurface elements.
  • Change the values until you see surfaces. In the first the value is 500. I changed this to 500,000. Multiply all by 1000.
  • In the second, 100 -> 100000.
  • In the third, 20 -> 20000.
  • In the last, 7 -> 7000.

This should maybe be performed by the dx.map user function, determining reasonable contour levels.

With a bit of zooming, clicking on File -> Save image in the Surface window, allowing rendering, and outputting to a large TIFF file, save current, then apply.

An example image cropped and converted to PNG in the GIMP at https://gna.org/bugs/download.php?file_id=20641.

Note that for a good resolution plot, you will need many more increments. Using the lower and upper dx.map arguments will be useful to zoom into the space.

Example images

These images is from bug 22024. Minimisation space for CR72 is catastrophic. The chi2 surface over dw and pA is bounded.

https://web.archive.org/web/gna.org/bugs/?22024 Bug 22024: Minimisation space for CR72 is catastrophic. The chi2 surface over dw and pA is bounded.
https://web.archive.org/web/gna.org/bugs/?22024 Bug 22024: Minimisation space for CR72 is catastrophic. The chi2 surface over dw and pA is bounded.
https://web.archive.org/web/gna.org/bugs/?22024 Bug 22024: Minimisation space for CR72 is catastrophic. The chi2 surface over dw and pA is bounded.
https://web.archive.org/web/gna.org/bugs/?22024 Bug 22024: Minimisation space for CR72 is catastrophic. The chi2 surface over dw and pA is bounded.
https://web.archive.org/web/gna.org/bugs/?22024 Bug 22024: Minimisation space for CR72 is catastrophic. The chi2 surface over dw and pA is bounded.

Code to generate map

This script will generate data that can help visualize the different models and the minimisation algorithms in relax.

Call the script cpmg_synthetic.py or similar. Remember .py ending

relax cpmg_synthetic.py

Here is the synthetic data generator code:

See also

Install dx