Author: tlinnet
Date: Sat Oct 11 17:35:46 2014
New Revision: 26247
URL: http://svn.gna.org/viewcvs/relax?rev=26247&view=rev
Log:
Added scripts to make surface plots of spin independents parameters dw and
r2a.
Added:
trunk/test_suite/shared_data/dispersion/KTeilum_FMPoulsen_MAkke_2006/surface_chi2_clustered_fitting/1_create_surface_data_S65_dw_r2a_FT128.py
trunk/test_suite/shared_data/dispersion/KTeilum_FMPoulsen_MAkke_2006/surface_chi2_clustered_fitting/2_make_surface_plot_S65_dw_r2a_FT128.py
Added:
trunk/test_suite/shared_data/dispersion/KTeilum_FMPoulsen_MAkke_2006/surface_chi2_clustered_fitting/1_create_surface_data_S65_dw_r2a_FT128.py
URL:
http://svn.gna.org/viewcvs/relax/trunk/test_suite/shared_data/dispersion/KTeilum_FMPoulsen_MAkke_2006/surface_chi2_clustered_fitting/1_create_surface_data_S65_dw_r2a_FT128.py?rev=26247&view=auto
==============================================================================
---
trunk/test_suite/shared_data/dispersion/KTeilum_FMPoulsen_MAkke_2006/surface_chi2_clustered_fitting/1_create_surface_data_S65_dw_r2a_FT128.py
(added)
+++
trunk/test_suite/shared_data/dispersion/KTeilum_FMPoulsen_MAkke_2006/surface_chi2_clustered_fitting/1_create_surface_data_S65_dw_r2a_FT128.py
Sat Oct 11 17:35:46 2014
@@ -0,0 +1,121 @@
+# Python imports
+from os import getcwd, sep
+from numpy import array, float64, zeros
+
+# relax module imports.
+from lib.io import open_write_file, write_data
+from pipe_control.mol_res_spin import display_spin, generate_spin_string,
return_spin
+from pipe_control import value
+from specific_analyses.api import return_api
+
+# Create pipe
+pipe.create('relax_disp', 'relax_disp')
+
+# The specific analysis API object.
+api = return_api()
+
+# Variables
+prev_data_path = getcwd()
+result_filename = 'FT_-_TSMFK01_-_min_-_128_-_free_spins.bz2'
+
+# Read data in
+results.read(prev_data_path + sep + result_filename)
+
+# Get residue of interest. S65 is
+cur_spin_id = ":%i@%s"%(65, 'N')
+
+# Get the spin container.
+cur_spin = return_spin(cur_spin_id)
+
+# Get the chi2 value
+pre_chi2 = cur_spin.chi2
+
+# Define surface map settings.
+dx_inc = 50
+
+# Lower bounds
+params = ['dw', 'r2a']
+lower = [0.0, 6.0]
+upper = [20.0, 12.0]
+
+# Get the current point for clustered mininimisation.
+pcm = [cur_spin.dw, cur_spin.r2a['SQ CPMG - 499.86214000 MHz']]
+print("Min cluster point %s=%3.3f, %s=%3.3f, with chi2=%3.3f" % (params[0],
pcm[0], params[1], pcm[1], pre_chi2))
+headings = [params[0], params[1], "chi2"]
+
+# Initialise.
+# Number of parameters
+n = 2
+
+# Get the default map bounds.
+bounds = zeros((n, 2), float64)
+
+# Lower bounds.
+bounds[:, 0] = array(lower, float64)
+
+# Upper bounds.
+bounds[:, 1] = array(upper, float64)
+
+# Setup the step sizes.
+step_size = zeros(n, float64)
+step_size = (bounds[:, 1] - bounds[:, 0]) / dx_inc
+
+# Placeholder to update values.
+values = zeros(n, float64)
+
+# Initial value of the first parameter.
+values[0] = bounds[0, 0]
+percent = 0.0
+percent_inc = 100.0 / (dx_inc + 1.0)**(n)
+print("%-10s%8.3f%-1s" % ("Progress:", percent, "%"))
+
+# Collect all chi2, to help finding a reasobale chi level.
+all_chi = []
+
+# Collect data.
+data = []
+# Append point as first data.
+data.append(["%3.3f"%pcm[0], "%3.3f"%pcm[1], "%3.3f"%pre_chi2 ])
+
+# Loop over the first parameter.
+for i in range((dx_inc + 1)):
+ # Initial value of the second parameter.
+ values[1] = bounds[1, 0]
+
+ # Loop over the second parameter.
+ for j in range((dx_inc + 1)):
+ # Set the value.
+ value.set(val=values, param=params, spin_id=cur_spin_id, force=True)
+
+ # Calculate the function values.
+ api.calculate(spin_id=cur_spin_id, verbosity=0)
+
+ # Get the minimisation statistics for the model.
+ k_stat, n_stat, chi2 = api.model_statistics(spin_id=cur_spin_id)
+ #print(k_stat, n_stat, chi2, "point is %s=%3.3f, %s=%3.3f"%
(params[0], values[0], params[1], values[1]))
+
+ # Progress incrementation and printout.
+ percent = percent + percent_inc
+ print("%-10s%8.3f%-8s%-8g" % ("Progress:", percent, "%, " +
repr(values) + ", f(x): ", chi2))
+
+ # Append to data.
+ data.append(["%3.3f"%values[0], "%3.3f"%values[1], "%3.3f"%chi2 ])
+
+ # Save all values of chi2. To help find reasonale level for the
Innermost, Inner, Middle and Outer Isosurface.
+ all_chi.append(chi2)
+
+ # Increment the value of the second parameter.
+ values[1] = values[1] + step_size[1]
+
+ # Increment the value of the first parameter.
+ values[0] = values[0] + step_size[0]
+
+print("\nMin cluster point %s=%3.3f, %s=%3.3f, with chi2=%3.3f" %
(params[0], pcm[0], params[1], pcm[1], pre_chi2))
+
+# Open file
+file_name = '1_create_surface_data_S65_dw_r2a_FT128.txt'
+surface_file = open_write_file(file_name=file_name, dir=None, force=True)
+write_data(out=surface_file, headings=headings, data=data)
+
+# Close file
+surface_file.close()
Added:
trunk/test_suite/shared_data/dispersion/KTeilum_FMPoulsen_MAkke_2006/surface_chi2_clustered_fitting/2_make_surface_plot_S65_dw_r2a_FT128.py
URL:
http://svn.gna.org/viewcvs/relax/trunk/test_suite/shared_data/dispersion/KTeilum_FMPoulsen_MAkke_2006/surface_chi2_clustered_fitting/2_make_surface_plot_S65_dw_r2a_FT128.py?rev=26247&view=auto
==============================================================================
---
trunk/test_suite/shared_data/dispersion/KTeilum_FMPoulsen_MAkke_2006/surface_chi2_clustered_fitting/2_make_surface_plot_S65_dw_r2a_FT128.py
(added)
+++
trunk/test_suite/shared_data/dispersion/KTeilum_FMPoulsen_MAkke_2006/surface_chi2_clustered_fitting/2_make_surface_plot_S65_dw_r2a_FT128.py
Sat Oct 11 17:35:46 2014
@@ -0,0 +1,115 @@
+import numpy as np
+import scipy.interpolate
+from numpy.ma import masked_where
+
+from mpl_toolkits.mplot3d import axes3d
+import matplotlib.pyplot as plt
+from matplotlib import cm
+
+resfile = open('1_create_surface_data_S65_dw_r2a_FT128.txt', 'r')
+
+lines = resfile.readlines()
+resfile.close()
+
+params = lines[0].split()[1:]
+mp = lines[1].split()
+mp_x = np.array([float(mp[0])])
+mp_y = np.array([float(mp[1])])
+mp_z = np.array([float(mp[2])])
+min_point = np.concatenate((mp_x, mp_y, mp_z))
+
+# Collect data
+x = []
+y = []
+z = []
+
+nr_dp = len(lines[2:])
+
+for line in lines[2:]:
+ x_l, y_l, z_l = line.split()
+ x.append(float(x_l))
+ y.append(float(y_l))
+ z.append(float(z_l))
+
+# Make numpy arrays
+x = np.asarray(x)
+y = np.asarray(y)
+z = np.asarray(z)
+
+x_min = x.min()
+x_max = x.max()
+y_min = y.min()
+y_max = y.max()
+z_min = z.min()
+z_max = z.max()
+
+
+# Expand axis and tile, to make mesh.
+x_tile = np.tile(x[np.newaxis, Ellipsis], (nr_dp, 1) )
+y_tile = np.tile(y[Ellipsis, np.newaxis], (1, nr_dp) )
+
+# Or do it by meshgrid
+x_mesh, y_mesh = np.meshgrid(x, y)
+
+# Test if new axis and tiling is the same
+print np.array_equal(x_tile, x_mesh)
+print np.array_equal(y_tile, y_mesh)
+
+# 2d contour plot from 3 lists : x, y and z?
+#
http://stackoverflow.com/questions/9008370/python-2d-contour-plot-from-3-lists-x-y-and-rho
+
+# Set up a regular grid of interpolation points
+int_points = 300
+xi, yi = np.linspace(x_min, x_max, int_points), np.linspace(y_min, y_max,
int_points)
+xi, yi = np.meshgrid(xi, yi)
+
+# This causes memor problem or a very long time.
+#rbf = scipy.interpolate.Rbf(x, y, z, function='linear')
+#zi = rbf(xi, yi)
+
+zi = scipy.interpolate.griddata((x, y), z, (xi, yi), method='linear')
+#z_mesh = scipy.interpolate.griddata((x, y), z, (x_mesh, y_mesh),
method='linear')
+
+fig = plt.figure()
+ax = fig.gca(projection='3d')
+
+# Set which x, y, z to plot
+x_p = xi
+y_p = yi
+z_p = zi
+
+# Replace out-lier value.
+# First get index os largest values
+out_val = 5*z_min
+z_p_mask = masked_where(z_p >= out_val, z_p)
+z_mask = masked_where(z >= out_val, z)
+
+# Replace with 0.0
+z_p[z_p_mask.mask] = 0.0
+z[z_mask.mask] = 0.0
+# Find new max
+new_max = np.max(z_p)
+z_p[z_p_mask.mask] = new_max
+z[z_mask.mask] = new_max
+
+ax.plot_surface(x_p, y_p, z_p, rstride=8, cstride=8, alpha=0.3)
+
+cset = ax.contour(x_p, y_p, z_p, zdir='z', offset=0, cmap=cm.coolwarm)
+cset = ax.contour(x_p, y_p, z_p, zdir='x', offset=x_min, cmap=cm.coolwarm)
+cset = ax.contour(x_p, y_p, z_p, zdir='y', offset=y_min, cmap=cm.coolwarm)
+
+ax.scatter(mp_x, mp_y, mp_z, c='r', marker='o', s=200)
+ax.scatter(x, y, z, c='b', marker='o', s=5)
+
+ax.set_xlabel('%s'%params[0])
+ax.set_xlim(x_min, x_max)
+
+ax.set_ylabel('%s'%params[1])
+ax.set_ylim(y_min, y_max)
+
+ax.set_zlabel('%s'%params[2])
+ax.set_zlim(0, out_val)
+
+plt.show()
+
+
r26247 - /trunk/test_suite/shared_data/dispersion/KTeilum_FMPoulsen_MAkke_2006/surface_chi2_clustered_fitting/