The code is out of order, but I might be able to find a spare hour in the next few days to quickly redesign this. So please don't remove the functionality. If you can plot isosurfaces in matplotlib, that would be awesome as not everyone can run OpenDX. What do you think of the previous ideas though? Especially the cdp storage ideas? Could you extend this to help with investigating the anomalies? For example, could you use more options for extending your matplotlib space mapping, that are specific to matplotlib?
Regards, Edward On 13 October 2014 20:14, Troels Emtekær Linnet <tlin...@nmr-relax.com> wrote: > Hi Edward. > > Yes, I know I am on the edge of putting the code here. > > It is a function related to matplotlib, but takes the data from dx map. > > So, I did not know where to put it. > > I am using the code, to back track some issues I have with some single spins > showing some abnormalities in fitting. > > I need to know what happens, and plotting helps with this. > > I do not have time to implement more, as you have suggested. > > If you think the code is out of order, to be put here, I will remove it > again, and make some small scripts on the wiki. > > Best > Troels > > 2014-10-13 18:34 GMT+02:00 Edward d'Auvergne <edw...@nmr-relax.com>: >> >> 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 > > _______________________________________________ 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
