Author: tlinnet
Date: Mon Oct 13 17:18:56 2014
New Revision: 26258
URL: http://svn.gna.org/viewcvs/relax?rev=26258&view=rev
Log:
Added the write out of a matplotlib command file, to plot surfaces of a dx
map.
It uses the minimum chi2 value in the map space, to define surface defitions.
It creates a X,Y; X,Z; Y,Z map, where the values in the missing dimension has
been cut at the minimum chi2 value.
For each map, it creates a projected 3d map of the parameters and the chi2
value, and a heat map for the contours.
It also scatters the minimum chi2 value, the 4 smallest chi2 values, and maps
any points in the point file, to a scatter point.
Mapping the points from file to map points, is done by finding the shortest
Euclidean distance in the space from the points to any map points.
Modified:
trunk/pipe_control/opendx.py
Modified: trunk/pipe_control/opendx.py
URL:
http://svn.gna.org/viewcvs/relax/trunk/pipe_control/opendx.py?rev=26258&r1=26257&r2=26258&view=diff
==============================================================================
--- trunk/pipe_control/opendx.py (original)
+++ trunk/pipe_control/opendx.py Mon Oct 13 17:18:56 2014
@@ -174,6 +174,9 @@
if create_par_file:
self.create_par_chi2(file_prefix=self.file_prefix,
par_chi2_vals=self.par_chi2_vals)
+ # Generate the matplotlib script file, to plot surfaces
+ self.matplotlib_surface_plot()
+
## Generate the file with parameters and associated chi2 value for
the points send to dx.
if self.num_points >= 1 and create_par_file:
# Calculate the parameter and associated chi2 values for the
points.
@@ -449,3 +452,235 @@
string = string + " " + repr(val)
val = val + loc_inc
self.tick_locations.append("{" + string + " }")
+
+
+ def matplotlib_surface_plot(self):
+ """Function to write matplotlib script to plot surfaces of
parameters."""
+
+ # Add ".par" to file prefix
+ mapfile_name = '"%s.par"' % self.file_prefix
+
+ # If point file_file is different from None
+ if self.point_file != None:
+ pointfile_name = '"%s.par"' % self.point_file
+ else:
+ pointfile_name = "None"
+
+ # Open the file.
+ plot_file = open_write_file(file_name=self.file_prefix+'.py',
dir=self.dir, force=True)
+
+ matplotlib_file = [
+ 'import numpy as np'+"\n",
+ 'import scipy.interpolate'+"\n",
+ 'from numpy.ma import masked_where'+"\n",
+ ''+"\n",
+ 'from mpl_toolkits.mplot3d import axes3d'+"\n",
+ 'import matplotlib.pyplot as plt'+"\n",
+ 'from matplotlib import cm'+"\n",
+ ''+"\n",
+ '# Open file and get header.'+"\n",
+ 'mapfile_name = %s'%mapfile_name+"\n",
+ 'pointfile_name = %s'%pointfile_name+"\n",
+ ''+"\n",
+ 'mapfile = open(mapfile_name, "r")'+"\n",
+ 'lines = mapfile.readlines()'+"\n",
+ 'mapfile.close()'+"\n",
+ 'header = lines[0].split()[1:]'+"\n",
+ ''+"\n",
+ '# Prepare the dtype for reading file.'+"\n",
+ 'dtype_str = "i8,f8,f8,f8,f8,i8,f8,f8,f8,f8"'+"\n",
+ ''+"\n",
+ 'print("Fileheader is: %s"%header)'+"\n",
+ 'print("Value types are: %s"%dtype_str)'+"\n",
+ ''+"\n",
+ '# Load the data.'+"\n",
+ 'data = np.genfromtxt(fname=mapfile_name, dtype=dtype_str,
names=header)'+"\n",
+ ''+"\n",
+ '# Load the point data'+"\n",
+ 'if pointfile_name:'+"\n",
+ ' # Load the point data.'+"\n",
+ ' data_p = np.genfromtxt(fname=pointfile_name,
dtype=dtype_str, names=header)'+"\n",
+ ' '+"\n",
+ ''+"\n",
+ '# Define where to cut the data, as the minimum.'+"\n",
+ 'header_min = header[6:10]'+"\n",
+ ''+"\n",
+ '# Define to cut at min map point.'+"\n",
+ 'map_min_par0 = data[header_min[0]][0]'+"\n",
+ 'map_min_par1 = data[header_min[1]][0]'+"\n",
+ 'map_min_par2 = data[header_min[2]][0]'+"\n",
+ 'map_min_chi2 = data[header_min[3]][0]'+"\n",
+ ''+"\n",
+ '# Now get the headers for the data.'+"\n",
+ 'header_val = header[1:5]'+"\n",
+ ''+"\n",
+ '# Define which 2D maps to create, as a list of 2 parameters,
and at which third parameter to cut the values.'+"\n",
+ 'maps_xy = [header_val[0], header_val[1], header_val[2],
map_min_par2]'+"\n",
+ 'maps_xz = [header_val[0], header_val[2], header_val[1],
map_min_par1]'+"\n",
+ 'maps_yz = [header_val[1], header_val[2], header_val[0],
map_min_par0]'+"\n",
+ ''+"\n",
+ 'maps = [maps_xy, maps_xz, maps_yz]'+"\n",
+ ''+"\n",
+ '# Nr of columns is number of maps.'+"\n",
+ 'nr_cols = 1'+"\n",
+ '# Nr of rows, is 2, for 3d projection and imshow'+"\n",
+ 'nr_rows = 2'+"\n",
+ ''+"\n",
+ '# Loop over the maps:'+"\n",
+ 'for x_par, y_par, z_par, z_cut in maps:'+"\n",
+ ' # Define figure'+"\n",
+ ' fig = plt.figure()'+"\n",
+ ''+"\n",
+ ' # Define c_par'+"\n",
+ ' c_par = header_val[3]'+"\n",
+ ''+"\n",
+ ' # Now get the values for the map.'+"\n",
+ ' map_x = data[x_par]'+"\n",
+ ' map_y = data[y_par]'+"\n",
+ ' map_z = data[z_par]'+"\n",
+ ' map_c = data[c_par]'+"\n",
+ ''+"\n",
+ ' # Now define which map to create.'+"\n",
+ ' mask_xy = masked_where(map_z == z_cut, map_z)'+"\n",
+ ' map_mask_x = map_x[mask_xy.mask]'+"\n",
+ ' map_mask_y = map_y[mask_xy.mask]'+"\n",
+ ' map_mask_c = map_c[mask_xy.mask]'+"\n",
+ ''+"\n",
+ ' # Define min and max values.'+"\n",
+ ' map_mask_x_min = map_mask_x.min()'+"\n",
+ ' map_mask_x_max = map_mask_x.max()'+"\n",
+ ' map_mask_y_min = map_mask_y.min()'+"\n",
+ ' map_mask_y_max = map_mask_y.max()'+"\n",
+ ' map_mask_c_min = map_mask_c.min()'+"\n",
+ ' map_mask_c_max = map_mask_c.max()'+"\n",
+ ''+"\n",
+ ' # Set up a regular grid of interpolation points'+"\n",
+ ' int_points = 300'+"\n",
+ ' xi, yi = np.linspace(map_mask_x_min, map_mask_x_max,
int_points), np.linspace(map_mask_y_min, map_mask_y_max, int_points)'+"\n",
+ ' xi, yi = np.meshgrid(xi, yi)'+"\n",
+ ''+"\n",
+ ' # Interpolate to create grid'+"\n",
+ ' ci = scipy.interpolate.griddata((map_mask_x, map_mask_y),
map_mask_c, (xi, yi), method="linear")'+"\n",
+ ''+"\n",
+ ' # Set which x, y, z to plot'+"\n",
+ ' x_p = xi'+"\n",
+ ' y_p = yi'+"\n",
+ ' c_p = ci'+"\n",
+ ''+"\n",
+ ' # Cut map at a certain height.'+"\n",
+ ' # First get index os largest values'+"\n",
+ ' #out_val = 5*map_mask_c_min'+"\n",
+ ' out_val = map_mask_c_max'+"\n",
+ ' ci_mask = masked_where(ci >= out_val, ci)'+"\n",
+ ''+"\n",
+ ' # Replace with 0.0'+"\n",
+ ' ci[ci_mask.mask] = 0.0'+"\n",
+ ' # Find new max'+"\n",
+ ' new_max = np.max(ci)'+"\n",
+ ''+"\n",
+ ' # Insert values in array.'+"\n",
+ ' ci[ci_mask.mask] = new_max'+"\n",
+ ''+"\n",
+ ' # Create figure and plot'+"\n",
+ ' ax = fig.add_subplot(nr_rows, nr_cols, 1,
projection="3d")'+"\n",
+ ' ax.plot_surface(x_p, y_p, c_p, rstride=8, cstride=8,
alpha=0.3)'+"\n",
+ ''+"\n",
+ ' # Possible add scatter points for map.'+"\n",
+ ' #ax.scatter(map_x, map_y, map_c, c="b", marker="o",
s=5)'+"\n",
+ ''+"\n",
+ ' # One could also make the mesh just from the values, but
this require much memory.'+"\n",
+ ' ##ax.scatter(x_p, y_p, c_p, c="y", marker="o", s=5)'+"\n",
+ ''+"\n",
+ ' # Add contour levels on sides.'+"\n",
+ ' ax.contour(x_p, y_p, c_p, zdir="z", offset=0,
cmap=cm.coolwarm)'+"\n",
+ ' ax.contour(x_p, y_p, c_p, zdir="x", offset=map_mask_x_min,
cmap=cm.coolwarm)'+"\n",
+ ' ax.contour(x_p, y_p, c_p, zdir="y", offset=map_mask_y_min,
cmap=cm.coolwarm)'+"\n",
+ ''+"\n",
+ ' # Add scatter values, for 5 lowest values.'+"\n",
+ ' x_par_min = x_par + "_sort"'+"\n",
+ ' y_par_min = y_par + "_sort"'+"\n",
+ ' c_par_min = c_par + "_sort"'+"\n",
+ ' mp_x = data[x_par_min][0:5]'+"\n",
+ ' mp_y = data[y_par_min][0:5]'+"\n",
+ ' mp_c = data[c_par_min][0:5]'+"\n",
+ ' ax.scatter(mp_x[0], mp_y[0], mp_c[0], c="r", marker="o",
s=200)'+"\n",
+ ' ax.scatter(mp_x[1:], mp_y[1:], mp_c[1:], c="g", marker="o",
s=100)'+"\n",
+ ''+"\n",
+ ' # Add points from file, as the closest point in map.'+"\n",
+ ' if pointfile_name:'+"\n",
+ ' if data_p[x_par].ndim == 0:'+"\n",
+ ' points_x = np.asarray([data_p[x_par]])'+"\n",
+ ' points_y = np.asarray([data_p[y_par]])'+"\n",
+ ' else:'+"\n",
+ ' points_x = data_p[x_par]'+"\n",
+ ' points_y = data_p[y_par]'+"\n",
+ ''+"\n",
+ ' # Normalize, by division of largest number of map'+"\n",
+ ' points_x_norm = points_x / map_mask_x_max'+"\n",
+ ' points_y_norm = points_y / map_mask_y_max'+"\n",
+ ' map_mask_x_norm = map_mask_x / map_mask_x_max'+"\n",
+ ' map_mask_y_norm = map_mask_y / map_mask_y_max'+"\n",
+ ''+"\n",
+ ' p_x = []'+"\n",
+ ' p_y = []'+"\n",
+ ' p_c = []'+"\n",
+ ' # Now calculate the Euclidean distance in the space, to
find map point best represents the point.'+"\n",
+ ' for i, point_x_norm in enumerate(points_x_norm):'+"\n",
+ ' point_y_norm = points_y_norm[i]'+"\n",
+ ''+"\n",
+ ' # Get the distance.'+"\n",
+ ' dist = np.sqrt( (map_mask_x_norm - point_x_norm)**2
+ (map_mask_y_norm - point_y_norm)**2)'+"\n",
+ ''+"\n",
+ ' # Return the indices of the minimum values along an
axis.'+"\n",
+ ' min_index = np.argmin(dist)'+"\n",
+ ' p_x.append(map_mask_x[min_index])'+"\n",
+ ' p_y.append(map_mask_y[min_index])'+"\n",
+ ' p_c.append(map_mask_c[min_index])'+"\n",
+ ''+"\n",
+ ' # Convert to numpy array'+"\n",
+ ' p_x = np.asarray(p_x)'+"\n",
+ ' p_y = np.asarray(p_y)'+"\n",
+ ' p_c = np.asarray(p_c)'+"\n",
+ ''+"\n",
+ ' # Plot points'+"\n",
+ ' ax.scatter(p_x, p_y, p_c, c="m", marker="o",
s=50)'+"\n",
+ ''+"\n",
+ ''+"\n",
+ ' # Set label'+"\n",
+ ' ax.set_xlabel("%s"%x_par)'+"\n",
+ ' ax.set_ylabel("%s"%y_par)'+"\n",
+ ' ax.set_zlabel("%s"%c_par)'+"\n",
+ ''+"\n",
+ ' # Create figure and plot'+"\n",
+ ' ax = fig.add_subplot(nr_rows, nr_cols, 2)'+"\n",
+ ' fig_imshow = ax.imshow(ci, vmin=map_mask_c_min,
vmax=map_mask_c_max, origin="lower", extent=[map_mask_x_min, map_mask_x_max,
map_mask_y_min, map_mask_y_max])'+"\n",
+ ''+"\n",
+ ' # Add scatter values, for 5 lowest values.'+"\n",
+ ' ax.scatter(mp_x[0], mp_y[0], c=mp_c[0], marker="o",
s=200)'+"\n",
+ ' ax.scatter(mp_x[1:], mp_y[1:], c="g", marker="o",
s=100)'+"\n",
+ ''+"\n",
+ ' # Also add point to this map.'+"\n",
+ ' if pointfile_name:'+"\n",
+ ' # Plot points'+"\n",
+ ' ax.scatter(p_x, p_y, c="m", marker="o", s=50)'+"\n",
+ ''+"\n",
+ ' # Set label'+"\n",
+ ' ax.set_xlabel("%s"%x_par)'+"\n",
+ ' ax.set_ylabel("%s"%y_par)'+"\n",
+ ''+"\n",
+ ' # Add colorbar.'+"\n",
+ ' fig.subplots_adjust(right=0.8)'+"\n",
+ ' cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.3])'+"\n",
+ ' fig.colorbar(fig_imshow, cax=cbar_ax)'+"\n",
+ ''+"\n",
+ '# Show plot.'+"\n",
+ 'plt.show()'+"\n",
+ ''+"\n",
+ ]
+
+ # Loop over the lines and write.
+ for line in matplotlib_file:
+ plot_file.write(line)
+
+ # Close the file.
+ plot_file.close()
r26258 - /trunk/pipe_control/opendx.py