Context Navigation

← Previous Changeset
Next Changeset →

Changeset 804

Timestamp:

07/16/15 16:36:54 (10 months ago)

Author:

mmckerns

Message:

add model_plotter as a callable function

Location:

mystic

Files:

: 1 added
: 1 edited

mystic/model_plotter.py (added)
scripts/mystic_model_plotter.py (modified) (1 diff)

Legend:

: Unmodified
: Added
: Removed

mystic/scripts/mystic_model_plotter.py

-                      r790
+                      r804
+#
 # Author: Mike McKerns (mmckerns @caltech and @uqfoundation)
+# Author: Jean-Christophe Fillion-Robin (jchris.fillionr @kitware.com)
 # Copyright (c) 1997-2015 California Institute of Technology.
 # License: 3-clause BSD.  The full license text is available at:
 #  - http://mmckerns.github.io/project/mystic/browser/mystic/LICENSE
+__doc__ = """
+mystic_model_plotter.py [options] model (filename)
+from mystic.model_plotter import model_plotter
+generate surface contour plots for model, specified by full import path
+generate model trajectory from logfile (or solver restart file), if provided
+The option "bounds" takes an indicator string, where the bounds should
+be given as comma-separated slices. For example, using bounds = "-1:10, 0:20"
+will set the lower and upper bounds for x to be (-1,10) and y to be (0,20).
+The "step" can also be given, to control the number of lines plotted in the
+grid. Thus "-1:10:.1, 0:20" would set the bounds as above, but use increments
+of .1 along x and the default step along y.  For models with > 2D, the bounds
+can be used to specify 2 dimensions plus fixed values for remaining dimensions.
+Thus, "-1:10, 0:20, 1.0" would plot the 2D surface where the z-axis was fixed
+at z=1.0.
+The option "label" takes comma-separated strings. For example, label = "x,y,"
+will place 'x' on the x-axis, 'y' on the y-axis, and nothing on the z-axis.
+LaTeX is also accepted. For example, label = "$ h $, $ {\\alpha}$, $ v$" will
+label the axes with standard LaTeX math formatting. Note that the leading
+space is required, while a trailing space aligns the text with the axis
+instead of the plot frame.
+The option "reduce" can be given to reduce the output of a model to a scalar,
+thus converting 'model(params)' to 'reduce(model(params))'. A reducer is given
+by the import path (e.g. 'numpy.add'). The option "scale" will convert the plot
+to log-scale, and scale the cost by 'z=log(4*z*scale+1)+2'. This is useful for
+visualizing small contour changes around the minimium. If using log-scale
+produces negative numbers, the option "shift" can be used to shift the cost
+by 'z=z+shift'. Both shift and scale are intended to help visualize contours.
+Required Inputs:
+  model               full import path for the model (e.g. mystic.models.rosen)
+Additional Inputs:
+  filename            name of the convergence logfile (e.g. log.txt)
+"""
+from mpl_toolkits.mplot3d import axes3d
+import matplotlib.pyplot as plt
+from matplotlib import cm
+from mystic.munge import read_history
+from mystic.munge import logfile_reader, raw_to_support
+#XXX: better if reads single id only? (e.g. same interface as read_history)
+def get_history(source, ids=None):
+    """get params and cost from the given source
+source is the name of the trajectory logfile (or solver instance)
+if provided, ids are the list of 'run ids' to select
+    """
+    try: # if it's a logfile, it might be multi-id
+        step, param, cost = logfile_reader(source)
+    except: # it's not a logfile, so read and return
+        param, cost = read_history(source)
+        return [param],[cost]
+    # split (i,id) into iteration and id
+    multinode = len(step[0]) - 1  #XXX: what if step = []?
+    if multinode: id = [i[1] for i in step]
+    else: id = [0 for i in step]
+    params = [[] for i in range(max(id) + 1)]
+    costs = [[] for i in range(len(params))]
+    # populate params for each id with the corresponding (param,cost)
+    for i in range(len(id)):
+        if ids is None or id[i] in ids: # take only the selected 'id'
+            params[id[i]].append(param[i])
+            costs[id[i]].append(cost[i])
+    params = [r for r in params if len(r)] # only keep selected 'ids'
+    costs = [r for r in costs if len(r)] # only keep selected 'ids'
+    # convert to support format
+    for i in range(len(params)):
+        params[i], costs[i] = raw_to_support(params[i], costs[i])
+    return params, costs
+def get_instance(location, *args, **kwds):
+    """given the import location of a model or model class, return the model
+args and kwds will be passed to the constructor of the model class
+    """
+    package, target = location.rsplit('.',1)
+    exec "from %s import %s as model" % (package, target)
+    import inspect
+    if inspect.isclass(model):
+        model = model(*args, **kwds)
+    return model
+def parse_input(option):
+    """parse 'option' string into 'select', 'axes', and 'mask'
+select contains the dimension specifications on which to plot
+axes holds the indicies of the parameters selected to plot
+mask is a dictionary of the parameter indicies and fixed values
+For example,
+    >>> select, axes, mask = parse_input("-1:10:.1, 0.0, 5.0, -50:50:.5")
+    >>> select
+    [0, 3]
+    >>> axes
+    "-1:10:.1, -50:50:.5"
+    >>> mask
+    {1: 0.0, 2: 5.0}
+    """
+    option = option.split(',')
+    select = []
+    axes = []
+    mask = {}
+    for index,value in enumerate(option):
+        if ":" in value:
+            select.append(index)
+            axes.append(value)
+        else:
+            mask.update({index:float(value)})
+    axes = ','.join(axes)
+    return select, axes, mask
+def parse_axes(option, grid=True):
+    """parse option string into grid axes; using modified numpy.ogrid notation
+For example:
+  option='-1:10:.1, 0:10:.1' yields x,y=ogrid[-1:10:.1,0:10:.1],
+If grid is False, accept options suitable for line plotting.
+For example:
+  option='-1:10' yields x=ogrid[-1:10] and y=0,
+  option='-1:10, 2' yields x=ogrid[-1:10] and y=2,
+Returns tuple (x,y) with 'x,y' defined above.
+    """
+    import numpy
+    option = option.split(',')
+    opt = dict(zip(['x','y','z'],option))
+    if len(option) > 2 or len(option) < 1:
+        raise ValueError("invalid format string: '%s'" % ','.join(option))
+    z = bool(grid)
+    if len(option) == 1: opt['y'] = '0'
+    xd = True if ':' in opt['x'] else False
+    yd = True if ':' in opt['y'] else False
+    #XXX: accepts option='3:1', '1:1', and '1:2:10' (try to catch?)
+    if xd and yd:
+        try: # x,y form a 2D grid
+            exec('x,y = numpy.ogrid[%s,%s]' % (opt['x'],opt['y']))
+        except: # AttributeError:
+            raise ValueError("invalid format string: '%s'" % ','.join(option))
+    elif xd and not z:
+        try:
+            exec('x = numpy.ogrid[%s]' % opt['x'])
+            y = float(opt['y'])
+        except: # (AttributeError, SyntaxError, ValueError):
+            raise ValueError("invalid format string: '%s'" % ','.join(option))
+    elif yd and not z:
+        try:
+            x = float(opt['x'])
+            exec('y = numpy.ogrid[%s]' % opt['y'])
+        except: # (AttributeError, SyntaxError, ValueError):
+            raise ValueError("invalid format string: '%s'" % ','.join(option))
+    else:
+        raise ValueError("invalid format string: '%s'" % ','.join(option))
+    if not x.size or not y.size:
+        raise ValueError("invalid format string: '%s'" % ','.join(option))
+    return x,y
+def draw_projection(x, cost, scale=True, shift=False, style=None, figure=None):
+    """draw a solution trajectory (for overlay on a 1D plot)
+x is the sequence of values for one parameter (i.e. a parameter trajectory)
+cost is the sequence of costs (i.e. the solution trajectory)
+if scale is provided, scale the intensity as 'z = log(4*z*scale+1)+2'
+if shift is provided, shift the intensity as 'z = z+shift' (useful for -z's)
+if style is provided, set the line style (e.g. 'w-o', 'k-', 'ro')
+if figure is provided, plot to an existing figure
+    """
+    if not figure: figure = plt.figure()
+    ax = figure.gca()
+    ax.autoscale(tight=True)
+    if style in [None, False]:
+        style = 'k-o'
+    import numpy
+    if shift:
+        if shift is True: #NOTE: MAY NOT be the exact minimum
+            shift = max(-numpy.min(cost), 0.0) + 0.5 # a good guess
+        cost = numpy.asarray(cost)+shift
+    cost = numpy.asarray(cost)
+    if scale:
+        cost = numpy.log(4*cost*scale+1)+2
+    ax.plot(x,cost, style, linewidth=2, markersize=4)
+    #XXX: need to 'correct' the z-axis (or provide easy conversion)
+    return figure
+def draw_trajectory(x, y, cost=None, scale=True, shift=False, style=None, figure=None):
+    """draw a solution trajectory (for overlay on a contour plot)
+x is a sequence of values for one parameter (i.e. a parameter trajectory)
+y is a sequence of values for one parameter (i.e. a parameter trajectory)
+cost is the solution trajectory (i.e. costs); if provided, plot a 3D contour
+if scale is provided, scale the intensity as 'z = log(4*z*scale+1)+2'
+if shift is provided, shift the intensity as 'z = z+shift' (useful for -z's)
+if style is provided, set the line style (e.g. 'w-o', 'k-', 'ro')
+if figure is provided, plot to an existing figure
+    """
+    if not figure: figure = plt.figure()
+    if cost: kwds = {'projection':'3d'} # 3D
+    else: kwds = {}                     # 2D
+    ax = figure.gca(**kwds)
+    if style in [None, False]:
+        style = 'w-o' #if not scale else 'k-o'
+    if cost: # is 3D, cost is needed
+        import numpy
+        if shift:
+            if shift is True: #NOTE: MAY NOT be the exact minimum
+                shift = max(-numpy.min(cost), 0.0) + 0.5 # a good guess
+            cost = numpy.asarray(cost)+shift
+        if scale:
+            cost = numpy.asarray(cost)
+            cost = numpy.log(4*cost*scale+1)+2
+        ax.plot(x,y,cost, style, linewidth=2, markersize=4)
+        #XXX: need to 'correct' the z-axis (or provide easy conversion)
+    else:    # is 2D, cost not needed
+        ax.plot(x,y, style, linewidth=2, markersize=4)
+    return figure
+def draw_slice(f, x, y=None, scale=True, shift=False):
+    """plot a slice of a 2D function 'f' in 1D
+x is an array used to set up the axis
+y is a fixed value for the 2nd axis
+if scale is provided, scale the intensity as 'z = log(4*z*scale+1)+2'
+if shift is provided, shift the intensity as 'z = z+shift' (useful for -z's)
+NOTE: when plotting the 'y-axis' at fixed 'x',
+pass the array to 'y' and the fixed value to 'x'
+    """
+    import numpy
+    if y is None:
+        y = 0.0
+    x, y = numpy.meshgrid(x, y)
+    plotx = True if numpy.all(y == y[0,0]) else False
+    z = 0*x
+    s,t = x.shape
+    for i in range(s):
+        for j in range(t):
+            xx,yy = x[i,j], y[i,j]
+            z[i,j] = f([xx,yy])
+    if shift:
+        if shift is True: shift = max(-numpy.min(z), 0.0) + 0.5 # exact minimum
+        z = z+shift
+    if scale: z = numpy.log(4*z*scale+1)+2
+    #XXX: need to 'correct' the z-axis (or provide easy conversion)
+    fig = plt.figure()
+    ax = fig.gca()
+    ax.autoscale(tight=True)
+    if plotx:
+        ax.plot(x.reshape(-1), z.reshape(-1))
+    else:
+        ax.plot(y.reshape(-1), z.reshape(-1))
+    return fig
+def draw_contour(f, x, y=None, surface=False, fill=True, scale=True, shift=False, density=5):
+    """draw a contour plot for a given 2D function 'f'
+x and y are arrays used to set up a 2D mesh grid
+if fill is True, color fill the contours
+if surface is True, plot the contours as a 3D projection
+if scale is provided, scale the intensity as 'z = log(4*z*scale+1)+2'
+if shift is provided, shift the intensity as 'z = z+shift' (useful for -z's)
+use density to adjust the number of contour lines
+    """
+    import numpy
+    if y is None:
+        y = x
+    x, y = numpy.meshgrid(x, y)
+    z = 0*x
+    s,t = x.shape
+    for i in range(s):
+        for j in range(t):
+            xx,yy = x[i,j], y[i,j]
+            z[i,j] = f([xx,yy])
+    if shift:
+        if shift is True: shift = max(-numpy.min(z), 0.0) + 0.5 # exact minimum
+        z = z+shift
+    if scale: z = numpy.log(4*z*scale+1)+2
+    #XXX: need to 'correct' the z-axis (or provide easy conversion)
+    fig = plt.figure()
+    if surface and fill is None: # 'hidden' option; full 3D surface plot
+        ax = fig.gca(projection='3d')
+        d = max(11 - density, 1) # or 1/density ?
+        kwds = {'rstride':d,'cstride':d,'cmap':cm.jet,'linewidth':0}
+        ax.plot_surface(x, y, z, **kwds)
+    else:
+        if surface: kwds = {'projection':'3d'} # 3D
+        elif surface is None: # 1D
+            raise NotImplementedError('need to add an option string parser')
+        else: kwds = {}                        # 2D
+        ax = fig.gca(**kwds)
+        density = 10*density
+        if fill: plotf = ax.contourf  # filled contours
+        else: plotf = ax.contour      # wire contours
+        plotf(x, y, z, density, cmap=cm.jet)
+    return fig
+__doc__ = model_plotter.__doc__
 if __name__ == '__main__':
-   #FIXME: should be able to:
-   # - apply a constraint as a region of NaN -- apply when 'xx,yy=x[ij],y[ij]'
-   # - apply a penalty by shifting the surface (plot w/alpha?) -- as above
-   # - build an appropriately-sized default grid (from logfile info)
-   # - move all mulit-id param/cost reading into read_history
-   #FIXME: current issues:
-   # - 1D slice and projection work for 2D function, but aren't "pretty"
-   # - 1D slice and projection for 1D function, is it meaningful and correct?
-   # - should be able to plot from solver.genealogy (multi-monitor?) [1D,2D,3D?]
-   # - should be able to scale 'z-axis' instead of scaling 'z' itself
-   #   (see https://github.com/matplotlib/matplotlib/issues/209)
-   # - if trajectory outside contour grid, will increase bounds
-   #   (see support_hypercube.py for how to fix bounds)
+    #XXX: note that 'argparse' is new as of python2.7
+    from optparse import OptionParser
+    parser = OptionParser(usage=__doc__)
+    parser.add_option("-b","--bounds",action="store",dest="bounds",\
+                      metavar="STR",default="-5:5:.1, -5:5:.1",
+                      help="indicator string to set plot bounds and density")
+    parser.add_option("-l","--label",action="store",dest="label",\
+                      metavar="STR",default=",,",
+                      help="string to assign label to axis")
+    parser.add_option("-n","--nid",action="store",dest="id",\
+                      metavar="INT",default=None,
+                      help="id # of the nth simultaneous points to plot")
+    parser.add_option("-i","--iter",action="store",dest="stop",\
+                      metavar="STR",default=":",
+                      help="string for smallest:largest iterations to plot")
+    parser.add_option("-r","--reduce",action="store",dest="reducer",\
+                      metavar="STR",default="None",
+                      help="import path of output reducer function")
+    parser.add_option("-x","--scale",action="store",dest="zscale",\
+                      metavar="INT",default=0.0,
+                      help="scale plotted cost by z=log(4*z*scale+1)+2")
+    parser.add_option("-z","--shift",action="store",dest="zshift",\
+                      metavar="INT",default=0.0,
+                      help="shift plotted cost by z=z+shift")
+    parser.add_option("-f","--fill",action="store_true",dest="fill",\
+                      default=False,help="plot using filled contours")
+    parser.add_option("-d","--depth",action="store_true",dest="surface",\
+                      default=False,help="plot contours showing depth in 3D")
+    parser.add_option("-o","--dots",action="store_true",dest="dots",\
+                      default=False,help="show trajectory points in plot")
+    parser.add_option("-j","--join",action="store_true",dest="line",\
+                      default=False,help="connect trajectory points in plot")
+    parsed_opts, parsed_args = parser.parse_args()
+    # get the import path for the model
+    model = parsed_args[0]  # e.g. 'mystic.models.rosen'
+    if "None" == model: model = None #XXX: 'required'... allow this?
+    try: # get the name of the parameter log file
+      source = parsed_args[1]  # e.g. 'log.txt'
+    except:
+      source = None
+    try: # select the bounds
+      options = parsed_opts.bounds  # format is "-1:10:.1, -1:10:.1, 1.0"
+    except:
+      options = "-5:5:.1, -5:5:.1"
+    try: # plot using filled contours
+      fill = parsed_opts.fill
+    except:
+      fill = False
+    try: # plot contours showing depth in 3D
+      surface = parsed_opts.surface
+    except:
+      surface = False
+    #XXX: can't do '-x' with no argument given  (use T/F instead?)
+    try: # scale plotted cost by z=log(4*z*scale+1)+2
+      scale = float(parsed_opts.zscale)
+      if not scale: scale = False
+    except:
+      scale = False
+    #XXX: can't do '-z' with no argument given
+    try: # shift plotted cost by z=z+shift
+      shift = float(parsed_opts.zshift)
+      if not shift: shift = False
+    except:
+      shift = False
+    try: # import path of output reducer function
+      reducer = parsed_opts.reducer  # e.g. 'numpy.add'
+      if "None" == reducer: reducer = None
+    except:
+      reducer = None
+    style = '-' # default linestyle
+    if parsed_opts.dots:
+      mark = 'o' # marker=mark
+      # when using 'dots', also can turn off 'line'
+      if not parsed_opts.line:
+        style = '' # linestyle='None'
+    else:
+      mark = ''
+    color = 'w' if fill else 'k'
+    style = color + style + mark
+    try: # select labels for the axes
+      label = parsed_opts.label.split(',')  # format is "x, y, z"
+    except:
+      label = ['','','']
+    try: # select which 'id' to plot results for
+      ids = (int(parsed_opts.id),) #XXX: allow selecting more than one id ?
+    except:
+      ids = None # i.e. 'all'
+    try: # select which iteration to stop plotting at
+      stop = parsed_opts.stop  # format is "1:10:1"
+      stop = stop if ":" in stop else ":"+stop
+    except:
+      stop = ":"
+    #################################################
+    solver = None  # set to 'mystic.solvers.fmin' (or similar) for 'live' fits
+    #NOTE: 'live' runs constrain params explicitly in the solver, then reduce
+    #      dimensions appropriately so results can be 2D contour plotted.
+    #      When working with legacy results that have more than 2 params,
+    #      the trajectory WILL NOT follow the masked surface generated
+    #      because the masked params were NOT fixed when the solver was run.
+    #################################################
+    from mystic.tools import reduced, masked, partial
+    # process inputs
+    select, spec, mask = parse_input(options)
+    x,y = parse_axes(spec, grid=True) # grid=False for 1D plots
+    #FIXME: does grid=False still make sense here...?
+    if reducer: reducer = get_instance(reducer)
+    if solver and (not source or not model):
+        raise RuntimeError('a model and results filename are required')
+    elif not source and not model:
+        raise RuntimeError('a model or a results file is required')
+    if model:
+        model = get_instance(model)
+        # need a reducer if model returns an array
+        if reducer: model = reduced(reducer, arraylike=False)(model)
+    if solver:
+        # if 'live'... pick a solver
+        solver = 'mystic.solvers.fmin'
+        solver = get_instance(solver)
+        xlen = len(select)+len(mask)
+        if solver.__name__.startswith('diffev'):
+            initial = [(-1,1)]*xlen
+        else:
+            initial = [0]*xlen
+        from mystic.monitors import VerboseLoggingMonitor
+        itermon = VerboseLoggingMonitor(filename=source, new=True)
+        # explicitly constrain parameters
+        model = partial(mask)(model)
+        # solve
+        sol = solver(model, x0=initial, itermon=itermon)
+        #-OVERRIDE-INPUTS-#
+        import numpy
+        # read trajectories from monitor (comment out to use logfile)
+        source = itermon
+        # if negative minimum, shift by the 'solved minimum' plus an epsilon
+        shift = max(-numpy.min(itermon.y), 0.0) + 0.5 # a good guess
+        #-----------------#
+    if model: # for plotting, implicitly constrain by reduction
+        model = masked(mask)(model)
+       ## plot the surface in 1D
+       #if solver: v=sol[-1]
+       #elif source: v=cost[-1]
+       #else: v=None
+       #fig0 = draw_slice(model, x=x, y=v, scale=scale, shift=shift)
+        # plot the surface in 2D or 3D
+        fig = draw_contour(model, x, y, surface=surface, fill=fill, scale=scale, shift=shift)
+    else:
+       #fig0 = None
+        fig = None
+    if source:
+        # params are the parameter trajectories
+        # cost is the solution trajectory
+        params, cost = get_history(source, ids)
+        if len(cost) > 1: style = style[1:] # 'auto-color' #XXX: or grayscale?
+        for p,c in zip(params, cost):
+           ## project trajectory on a 1D slice of model surface #XXX: useful?
+           #s = select[0] if len(select) else 0
+           #px = p[int(s)] # draw_projection requires one parameter
+           ## ignore everything after 'stop'
+           #_c = eval('c[%s]' % stop)
+           #_x = eval('px[%s]' % stop)
+           #fig0 = draw_projection(_x,_c, style=style, scale=scale, shift=shift, figure=fig0)
+            # plot the trajectory on the model surface (2D or 3D)
+            # get two selected params #XXX: what if len(select)<2? or len(p)<2?
+            p = [p[int(i)] for i in select[:2]]
+            px,py = p # draw_trajectory requires two parameters
+            # ignore everything after 'stop'
+            _x = eval('px[%s]' % stop)
+            _y = eval('py[%s]' % stop)
+            _c = eval('c[%s]' % stop) if surface else None
+            fig = draw_trajectory(_x,_y,_c, style=style, scale=scale, shift=shift, figure=fig)
+    # add labels to the axes
+    if surface: kwds = {'projection':'3d'} # 3D
+    else: kwds = {}                        # 2D
+    ax = fig.gca(**kwds)
+    ax.set_xlabel(label[0])
+    ax.set_ylabel(label[1])
+    if surface: ax.set_zlabel(label[2])
+    plt.show()
+    import sys
+    model_plotter(cmdargs=sys.argv[1:])

Note: See TracChangeset for help on using the changeset viewer.

Context Navigation

Changeset 804

Legend:

mystic/scripts/mystic_model_plotter.py

Download in other formats: