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.)
 
(36 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'''.
* Click on 'Edit Visual Programs...',
+
* Click on <code>Edit Visual Programs...</code>.
* Select the map.net program created by relax,
+
* Select the <code>map.net</code> program created by relax,
* Select the menu entry 'Execute->Execute on change'.
+
 
 +
Now in the <code>Visual Program Editor</code> or <code>VPE</code>.
 +
* 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.
+
* 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 -> 100,000.
+
* In the second, 100 -> 100000.
* In the third, 20 -> 20,000.
+
* In the third, 20 -> 20000.
* In the last, 7 -> 7,000.
+
* In the last, 7 -> 7000.
  
 
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 ==
 +
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.]]
 +
|}
 +
 
 +
{|style="margin: 0 auto;"
 +
| [[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 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|{{gna bug url|22024}} Bug 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 52: Line 99:
 
relax cpmg_synthetic.py
 
relax cpmg_synthetic.py
 
</source>
 
</source>
Here is the script '''cpmg_synthetic.py'''
+
 
<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 =
 +
###############################################################################
 +
#                                                                            #
 +
# Copyright (C) 2013-2014 Troels E. Linnet                                    #
 +
#                                                                            #
 +
# This file is part of the program relax (http://www.nmr-relax.com).          #
 +
#                                                                            #
 +
# This program is free software: you can redistribute it and/or modify        #
 +
# it under the terms of the GNU General Public License as published by        #
 +
# the Free Software Foundation, either version 3 of the License, or          #
 +
# (at your option) any later version.                                        #
 +
#                                                                            #
 +
# This program is distributed in the hope that it will be useful,            #
 +
# but WITHOUT ANY WARRANTY; without even the implied warranty of              #
 +
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              #
 +
# GNU General Public License for more details.                                #
 +
#                                                                            #
 +
# You should have received a copy of the GNU General Public License          #
 +
# along with this program.  If not, see <http://www.gnu.org/licenses/>.      #
 +
#                                                                            #
 +
###############################################################################
 
# Script for calculating synthetics CPMG data.
 
# Script for calculating synthetics CPMG data.
  
Line 60: 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 81: 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 121: Line 195:
 
     ncycs_1 = [28, 4, 32, 60, 2, 10, 16, 8, 20, 50, 18, 40, 6, 12, 24]
 
     ncycs_1 = [28, 4, 32, 60, 2, 10, 16, 8, 20, 50, 18, 40, 6, 12, 24]
 
     # Here you define the direct R2eff errors (rad/s), as being added or subtracted for the created R2eff point in the corresponding ncyc cpmg frequence.
 
     # Here you define the direct R2eff errors (rad/s), as being added or subtracted for the created R2eff point in the corresponding ncyc cpmg frequence.
     r2eff_errs_1 = [0.05, -0.05, 0.05, -0.05, 0.05, -0.05, 0.05, -0.05, 0.05, -0.05, 0.05, -0.05, 0.05, -0.05, 0.05]
+
     #r2eff_errs_1 = [0.05, -0.05, 0.05, -0.05, 0.05, -0.05, 0.05, -0.05, 0.05, -0.05, 0.05, -0.05, 0.05, -0.05, 0.05]
     #r2eff_errs_1 = [0.0] * len(ncycs_1)
+
     r2eff_errs_1 = [0.0] * len(ncycs_1)
 
     exp_1 = [sfrq_1, time_T2_1, ncycs_1, r2eff_errs_1]
 
     exp_1 = [sfrq_1, time_T2_1, ncycs_1, r2eff_errs_1]
  
Line 135: Line 209:
  
 
     # Collect all exps
 
     # Collect all exps
     #exps = [exp_1, exp_2]
+
     exps = [exp_1, exp_2]
     exps = [exp_1]
+
     #exps = [exp_1]
  
 
     # Add more spins here.
 
     # Add more spins here.
 
     spins = [
 
     spins = [
             ['Ala', 1, 'N', {'r2': {r20_key_1: 20.0, r20_key_2: 20.0}, 'r2a': {r20_key_1: 20.0, r20_key_2: 20.0}, 'r2b': {r20_key_1: 20.0, r20_key_2: 20.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: 13.0, r20_key_2: 14.5}, 'r2a': {r20_key_1: 13.0, r20_key_2: 14.5}, 'r2b': {r20_key_1: 13.0, r20_key_2: 14.5}, '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 154: 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 191: 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
  
 +
# Which parameters to map on chi2 surface.
 +
if not hasattr(ds, 'dx_params'):
 +
    ds.dx_params = ['dw', 'pA', 'kex']
 +
 +
# How many increements to map for each of the parameters.
 +
# Above 20 is good. 70 - 100 is a quite good resolution.
 +
# 4 is for speed tests.
 
if not hasattr(ds, 'dx_inc'):
 
if not hasattr(ds, 'dx_inc'):
     ds.dx_inc = 20
+
     #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]
 +
 
 +
##################################################################################
  
# The set r2eff err.
+
# Define repeating functions.
if not hasattr(ds, 'r2eff_err'):
+
##################################################################################
    ds.r2eff_err = 0.1
+
# 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 225: Line 413:
  
 
# Create the data pipe.
 
# Create the data pipe.
pipe_name = 'base pipe'
+
ds.pipe_name = 'base pipe'
pipe_type = 'relax_disp'
+
ds.pipe_type = 'relax_disp'
pipe_bundle = 'relax_disp'
+
ds.pipe_bundle = 'relax_disp'
pipe_name_r2eff = "%s_%s_R2eff"%(model_create, pipe_name)
+
ds.pipe_name_r2eff = "%s_%s_R2eff"%(model_create, ds.pipe_name )
pipe.create(pipe_name=pipe_name, pipe_type=pipe_type, bundle = pipe_bundle)
+
pipe.create(pipe_name=ds.pipe_name , pipe_type=ds.pipe_type, bundle = ds.pipe_bundle)
  
 
# Generate the sequence.
 
# Generate the sequence.
Line 242: 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 253: Line 441:
 
         nu_cpmg = ncyc / time_T2
 
         nu_cpmg = ncyc / time_T2
 
         cur_id = '%s_%.1f' % (exp_id, nu_cpmg)
 
         cur_id = '%s_%.1f' % (exp_id, nu_cpmg)
        print cur_id
 
 
         ids.append(cur_id)
 
         ids.append(cur_id)
  
Line 271: Line 458:
  
 
## Now prepare to calculate the synthetic R2eff values.
 
## Now prepare to calculate the synthetic R2eff values.
pipe.copy(pipe_from=pipe_name, pipe_to=pipe_name_r2eff, bundle_to = pipe_bundle)
+
pipe.copy(pipe_from=ds.pipe_name , pipe_to=ds.pipe_name_r2eff, bundle_to = ds.pipe_bundle)
pipe.switch(pipe_name=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)
  
# First loop over the spins and set the model parameters.
+
# First loop over the defined spins and set the model parameters.
for res_name, res_num, spin_name, params in cur_spins:
+
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_id = ":%i@%s"%(res_num, spin_name)
 
     cur_spin = return_spin(cur_spin_id)
 
     cur_spin = return_spin(cur_spin_id)
    #print cur_spin.model, cur_spin.name, cur_spin.isotope
 
  
     #print as
+
     if ds.print_res:
    # Now set the parameters.
+
        print("For spin: '%s'"%cur_spin_id)
 
     for mo_param in cur_spin.params:
 
     for mo_param in cur_spin.params:
 
         # The R2 is a dictionary, depending on spectrometer frequency.
 
         # The R2 is a dictionary, depending on spectrometer frequency.
 
         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 set_r2.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) })
                 print cur_spin.model, res_name, cur_spin_id, mo_param, key, float(val)
+
                 print(cur_spin.model, res_name, cur_spin_id, mo_param, key, float(val))
 +
            # Set it back
 
             setattr(cur_spin, mo_param, set_r2)
 
             setattr(cur_spin, mo_param, set_r2)
 
         else:
 
         else:
Line 298: Line 487:
 
             setattr(cur_spin, mo_param, float(params[mo_param]))
 
             setattr(cur_spin, mo_param, float(params[mo_param]))
 
             after = getattr(cur_spin, mo_param)
 
             after = getattr(cur_spin, mo_param)
             print cur_spin.model, res_name, cur_spin_id, mo_param, before
+
             print(cur_spin.model, res_name, cur_spin_id, mo_param, before)
 
 
  
 
## 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 328: 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 357: Line 544:
 
         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.
+
### Now do fitting.
        i += 1
 
 
 
print("Did following number of iterations: %i"%i)
 
 
 
# Now do fitting.
 
 
 
 
# Change pipe.
 
# Change pipe.
pipe_name_MODEL = "%s_%s"%(pipe_name, model_analyse)
+
ds.pipe_name_MODEL = "%s_%s"%(ds.pipe_name , model_analyse)
pipe.copy(pipe_from=pipe_name, pipe_to=pipe_name_MODEL, bundle_to = pipe_bundle)
+
pipe.copy(pipe_from=ds.pipe_name , pipe_to=ds.pipe_name_MODEL, bundle_to = ds.pipe_bundle)
pipe.switch(pipe_name=pipe_name_MODEL)
+
pipe.switch(pipe_name=ds.pipe_name_MODEL)
  
 
# Copy R2eff, but not the original parameters
 
# Copy R2eff, but not the original parameters
value.copy(pipe_from=pipe_name_r2eff, pipe_to=pipe_name_MODEL, param='r2eff')
+
value.copy(pipe_from=ds.pipe_name_r2eff, pipe_to=ds.pipe_name_MODEL, param='r2eff')
  
 
# Then select model.
 
# Then select model.
Line 376: Line 557:
  
 
print("Analysing with MODEL:%s."%(model_analyse))
 
print("Analysing with MODEL:%s."%(model_analyse))
 
# Do a dx map.
 
# To map the hypersurface of chi2, when altering kex, dw and pA.
 
if ds.opendx:
 
    for cur_spin, mol_name, resi, resn, spin_id in spin_loop(full_info=True, return_id=True, skip_desel=True):
 
        cur_spin_id = spin_id
 
        file_name = "map%s" % (cur_spin_id .replace('#', '_').replace(':', '_').replace('@', '_'))
 
        dx.map(params=['dw', 'pA', 'kex'], map_type='Iso3D', spin_id=":1@N", inc=ds.dx_inc, 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)
 
  
 
# Remove insignificant
 
# Remove insignificant
Line 407: 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 437: Line 590:
 
if ds.do_cluster:
 
if ds.do_cluster:
 
     # Change pipe.
 
     # Change pipe.
     pipe_name_MODEL_CLUSTER = "%s_%s_CLUSTER"%(pipe_name, model_create)
+
     ds.pipe_name_MODEL_CLUSTER = "%s_%s_CLUSTER"%(ds.pipe_name , model_create)
     pipe.copy(pipe_from=pipe_name, pipe_to=pipe_name_MODEL_CLUSTER)
+
     pipe.copy(pipe_from=ds.pipe_name , pipe_to=ds.pipe_name_MODEL_CLUSTER)
     pipe.switch(pipe_name=pipe_name_MODEL_CLUSTER)
+
     pipe.switch(pipe_name=ds.pipe_name_MODEL_CLUSTER)
  
 
     # Copy R2eff, but not the original parameters
 
     # Copy R2eff, but not the original parameters
     value.copy(pipe_from=pipe_name_r2eff, pipe_to=pipe_name_MODEL_CLUSTER, param='r2eff')
+
     value.copy(pipe_from=ds.pipe_name_r2eff, pipe_to=ds.pipe_name_MODEL_CLUSTER, param='r2eff')
  
 
     # Then select model.
 
     # Then select model.
Line 451: Line 604:
  
 
     # Copy the parameters from before.
 
     # Copy the parameters from before.
     relax_disp.parameter_copy(pipe_from=pipe_name_MODEL, pipe_to=pipe_name_MODEL_CLUSTER)
+
     relax_disp.parameter_copy(pipe_from=ds.pipe_name_MODEL, pipe_to=ds.pipe_name_MODEL_CLUSTER)
  
 
     # Now minimise.
 
     # Now minimise.
Line 471: Line 624:
 
     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.
+
# Looping over data, to collect the results.
if ds.print_res:
+
# Define which dx_params to collect for.
    print("\n########################")
+
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)
    print("Generated data with MODEL:%s"%(model_create))
 
    print("Analysing with MODEL:%s."%(model_analyse))
 
    print("########################\n")
 
  
for i in range(len(cur_spins)):
+
## Do a dx map.
    res_name, res_num, spin_name, params = cur_spins[i]
+
# To map the hypersurface of chi2, when altering kex, dw and pA.
    cur_spin_id = ":%i@%s"%(res_num, spin_name)
+
if ds.opendx:
    cur_spin = return_spin(cur_spin_id)
+
    # First switch pipe, since dx.map will go through parameters and end up a "bad" place. :-)
 +
    ds.pipe_name_MODEL_MAP = "%s_%s_map"%(ds.pipe_name, model_analyse)
 +
    pipe.copy(pipe_from=ds.pipe_name , pipe_to=ds.pipe_name_MODEL_MAP, bundle_to = ds.pipe_bundle)
 +
    pipe.switch(pipe_name=ds.pipe_name_MODEL_MAP)
 +
 
 +
    # Copy R2eff, but not the original parameters
 +
    value.copy(pipe_from=ds.pipe_name_r2eff, pipe_to=ds.pipe_name_MODEL_MAP, param='r2eff')
 +
 
 +
    # Then select model.
 +
    relax_disp.select_model(model=model_analyse)
 +
 
 +
    # 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)):
 +
        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)
 +
 
 +
        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)
  
    grid_params = ds.grid_results[i][3]
+
        if ds.print_res:
    min_params = ds.min_results[i][3]
+
            print("For spin: '%s'"%cur_spin_id)
    clust_params = ds.clust_results[i][3]
+
            print("Params for dx map is")
    # Now read the parameters.
+
            print(ds.dx_params)
 +
            print("Point creating dx map is")
 +
            print(ds.dx_set_val)
 +
            print("Minimises point in dx map is")
 +
            print(ds.dx_clust_val)
  
    if ds.print_res:
+
# Print result again after dx map.
        print("For spin: '%s'"%cur_spin_id)
+
if ds.opendx:
     for mo_param in cur_spin.params:
+
     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])
        # 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")
 
</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