Re: r26258 - /trunk/pipe_control/opendx.py -- October 13, 2014 - 18:34

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,  <tlinnet@xxxxxxxxxxxxx> 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-commits@xxxxxxx
>
> 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