Hi Troels, This is clearly pushing the boundaries of this design, far beyond its scope of OpenDX specific functionality. Do you have other ideas of how you would like to extend this? Would you like relax to launch the matplotlib viewing code as well? A small redesign is clearly now required to prevent the dx user functions from becoming over complicated. I can help with this. From the code changes, I am guessing that you would like to try to generalise the mapping functionality to be software independent. Therefore I would suggest we do the following:
1) Shift the Map class maybe into pipe_control.mapping.base. This was planned for the future anyway, but I'll wait for your response before making changes, in case you are still working on this code. 2) Covert the Map class into a base class (this doesn't require any changes). 3) Create one OpenDX specific class and one matplotlib specific class (in pipe_control.mapping.opendx and pipe_control.mapping.matplotlib). These would then have a run() method which replaces what the current __init__() method is doing. 4) Rename the user_functions.dx module to user_functions.map. 5) Rename the dx.map user function to map.dx_create. 6) Rename the dx.execute user function to map.dx_execute. 7) Create the new map.matplotlib_create user function. A second idea would be: 1) Shift the Map class maybe into pipe_control.mapping.isosurface. 2) Rename the user_functions.dx module to user_functions.map. 3) Rename the dx.map user function to map.create. This would then have options for creating the different map types, having values of 'OpenDX', 'matplotlib', and 'text' (the last is for the 4 column text file you have already implemented). This could be allowed to be a list as well, if you would really like to have all maps created simultaneously. 4) Rename the dx.execute user function to map.dx_execute. 5) Create the new map.matplotlib_execute user function. The first design would be simpler, in that many more OpenDX or matplotlib specific arguments can be added without killing the GUI user function window. What do you think? Does either of these match your aims? An alternative idea that would be much more powerful would be to to shift a lot of the current dx.map functionality into a new map.setup user function. We then store things in the current data pipe such as cdp.map_parameters, cdp.map_type, cdp.map_inc, cdp.map_lower, cdp.map_upper, cdp.map_axis_inc, cdp.map_points, cdp.map_isosurface_levels. I could also quickly create a special Python object called cdp.map with XML methods which has the variables stored there (e.g. cdp.map.parameters and cdp.map.type). The new map.setup user function would also reset the cdp.map_data or cdp.map.data variable to an empty list. You may like this next idea - the Map class run() method (from the first idea above) returns all of the chi2 values as a list. We then store this in cdp.map.data. Then the next time the Map class is set up, we pass in the cdp.map.data values. If this list contains values, the the target function calls via the specific analysis API are skipped. Then if you call map.setup, map.dx_create, map.text_create, and map.matplotlib_create, the last two user functions will be super fast as map.dx_create would have already generated the chi2 values. We also then have the added benefit of having all the base mapping data stored in relax results and save files. Regards, Edward On 13 October 2014 17:18, <tlin...@nmr-relax.com> wrote: > 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() > > > _______________________________________________ > relax (http://www.nmr-relax.com) > > This is the relax-commits mailing list > relax-comm...@gna.org > > To unsubscribe from this list, get a password > reminder, or change your subscription options, > visit the list information page at > https://mail.gna.org/listinfo/relax-commits _______________________________________________ relax (http://www.nmr-relax.com) This is the relax-devel mailing list relax-devel@gna.org To unsubscribe from this list, get a password reminder, or change your subscription options, visit the list information page at https://mail.gna.org/listinfo/relax-devel
