Changeset 804
- Timestamp:
- 07/16/15 16:36:54 (10 months ago)
- Location:
- mystic
- Files:
-
- 1 added
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
mystic/scripts/mystic_model_plotter.py
r790 r804 2 2 # 3 3 # Author: Mike McKerns (mmckerns @caltech and @uqfoundation) 4 # Author: Jean-Christophe Fillion-Robin (jchris.fillionr @kitware.com) 4 5 # Copyright (c) 1997-2015 California Institute of Technology. 5 6 # License: 3-clause BSD. The full license text is available at: 6 7 # - http://mmckerns.github.io/project/mystic/browser/mystic/LICENSE 7 8 8 __doc__ = """ 9 mystic_model_plotter.py [options] model (filename) 9 from mystic.model_plotter import model_plotter 10 10 11 generate surface contour plots for model, specified by full import path 12 generate model trajectory from logfile (or solver restart file), if provided 13 14 The option "bounds" takes an indicator string, where the bounds should 15 be given as comma-separated slices. For example, using bounds = "-1:10, 0:20" 16 will set the lower and upper bounds for x to be (-1,10) and y to be (0,20). 17 The "step" can also be given, to control the number of lines plotted in the 18 grid. Thus "-1:10:.1, 0:20" would set the bounds as above, but use increments 19 of .1 along x and the default step along y. For models with > 2D, the bounds 20 can be used to specify 2 dimensions plus fixed values for remaining dimensions. 21 Thus, "-1:10, 0:20, 1.0" would plot the 2D surface where the z-axis was fixed 22 at z=1.0. 23 24 The option "label" takes comma-separated strings. For example, label = "x,y," 25 will place 'x' on the x-axis, 'y' on the y-axis, and nothing on the z-axis. 26 LaTeX is also accepted. For example, label = "$ h $, $ {\\alpha}$, $ v$" will 27 label the axes with standard LaTeX math formatting. Note that the leading 28 space is required, while a trailing space aligns the text with the axis 29 instead of the plot frame. 30 31 The option "reduce" can be given to reduce the output of a model to a scalar, 32 thus converting 'model(params)' to 'reduce(model(params))'. A reducer is given 33 by the import path (e.g. 'numpy.add'). The option "scale" will convert the plot 34 to log-scale, and scale the cost by 'z=log(4*z*scale+1)+2'. This is useful for 35 visualizing small contour changes around the minimium. If using log-scale 36 produces negative numbers, the option "shift" can be used to shift the cost 37 by 'z=z+shift'. Both shift and scale are intended to help visualize contours. 38 39 Required Inputs: 40 model full import path for the model (e.g. mystic.models.rosen) 41 42 Additional Inputs: 43 filename name of the convergence logfile (e.g. log.txt) 44 """ 45 46 from mpl_toolkits.mplot3d import axes3d 47 import matplotlib.pyplot as plt 48 from matplotlib import cm 49 50 from mystic.munge import read_history 51 from mystic.munge import logfile_reader, raw_to_support 52 53 #XXX: better if reads single id only? (e.g. same interface as read_history) 54 def get_history(source, ids=None): 55 """get params and cost from the given source 56 57 source is the name of the trajectory logfile (or solver instance) 58 if provided, ids are the list of 'run ids' to select 59 """ 60 try: # if it's a logfile, it might be multi-id 61 step, param, cost = logfile_reader(source) 62 except: # it's not a logfile, so read and return 63 param, cost = read_history(source) 64 return [param],[cost] 65 66 # split (i,id) into iteration and id 67 multinode = len(step[0]) - 1 #XXX: what if step = []? 68 if multinode: id = [i[1] for i in step] 69 else: id = [0 for i in step] 70 71 params = [[] for i in range(max(id) + 1)] 72 costs = [[] for i in range(len(params))] 73 # populate params for each id with the corresponding (param,cost) 74 for i in range(len(id)): 75 if ids is None or id[i] in ids: # take only the selected 'id' 76 params[id[i]].append(param[i]) 77 costs[id[i]].append(cost[i]) 78 params = [r for r in params if len(r)] # only keep selected 'ids' 79 costs = [r for r in costs if len(r)] # only keep selected 'ids' 80 81 # convert to support format 82 for i in range(len(params)): 83 params[i], costs[i] = raw_to_support(params[i], costs[i]) 84 return params, costs 85 86 87 def get_instance(location, *args, **kwds): 88 """given the import location of a model or model class, return the model 89 90 args and kwds will be passed to the constructor of the model class 91 """ 92 package, target = location.rsplit('.',1) 93 exec "from %s import %s as model" % (package, target) 94 import inspect 95 if inspect.isclass(model): 96 model = model(*args, **kwds) 97 return model 98 99 100 def parse_input(option): 101 """parse 'option' string into 'select', 'axes', and 'mask' 102 103 select contains the dimension specifications on which to plot 104 axes holds the indicies of the parameters selected to plot 105 mask is a dictionary of the parameter indicies and fixed values 106 107 For example, 108 >>> select, axes, mask = parse_input("-1:10:.1, 0.0, 5.0, -50:50:.5") 109 >>> select 110 [0, 3] 111 >>> axes 112 "-1:10:.1, -50:50:.5" 113 >>> mask 114 {1: 0.0, 2: 5.0} 115 """ 116 option = option.split(',') 117 select = [] 118 axes = [] 119 mask = {} 120 for index,value in enumerate(option): 121 if ":" in value: 122 select.append(index) 123 axes.append(value) 124 else: 125 mask.update({index:float(value)}) 126 axes = ','.join(axes) 127 return select, axes, mask 128 129 130 def parse_axes(option, grid=True): 131 """parse option string into grid axes; using modified numpy.ogrid notation 132 133 For example: 134 option='-1:10:.1, 0:10:.1' yields x,y=ogrid[-1:10:.1,0:10:.1], 135 136 If grid is False, accept options suitable for line plotting. 137 For example: 138 option='-1:10' yields x=ogrid[-1:10] and y=0, 139 option='-1:10, 2' yields x=ogrid[-1:10] and y=2, 140 141 Returns tuple (x,y) with 'x,y' defined above. 142 """ 143 import numpy 144 option = option.split(',') 145 opt = dict(zip(['x','y','z'],option)) 146 if len(option) > 2 or len(option) < 1: 147 raise ValueError("invalid format string: '%s'" % ','.join(option)) 148 z = bool(grid) 149 if len(option) == 1: opt['y'] = '0' 150 xd = True if ':' in opt['x'] else False 151 yd = True if ':' in opt['y'] else False 152 #XXX: accepts option='3:1', '1:1', and '1:2:10' (try to catch?) 153 if xd and yd: 154 try: # x,y form a 2D grid 155 exec('x,y = numpy.ogrid[%s,%s]' % (opt['x'],opt['y'])) 156 except: # AttributeError: 157 raise ValueError("invalid format string: '%s'" % ','.join(option)) 158 elif xd and not z: 159 try: 160 exec('x = numpy.ogrid[%s]' % opt['x']) 161 y = float(opt['y']) 162 except: # (AttributeError, SyntaxError, ValueError): 163 raise ValueError("invalid format string: '%s'" % ','.join(option)) 164 elif yd and not z: 165 try: 166 x = float(opt['x']) 167 exec('y = numpy.ogrid[%s]' % opt['y']) 168 except: # (AttributeError, SyntaxError, ValueError): 169 raise ValueError("invalid format string: '%s'" % ','.join(option)) 170 else: 171 raise ValueError("invalid format string: '%s'" % ','.join(option)) 172 if not x.size or not y.size: 173 raise ValueError("invalid format string: '%s'" % ','.join(option)) 174 return x,y 175 176 177 def draw_projection(x, cost, scale=True, shift=False, style=None, figure=None): 178 """draw a solution trajectory (for overlay on a 1D plot) 179 180 x is the sequence of values for one parameter (i.e. a parameter trajectory) 181 cost is the sequence of costs (i.e. the solution trajectory) 182 if scale is provided, scale the intensity as 'z = log(4*z*scale+1)+2' 183 if shift is provided, shift the intensity as 'z = z+shift' (useful for -z's) 184 if style is provided, set the line style (e.g. 'w-o', 'k-', 'ro') 185 if figure is provided, plot to an existing figure 186 """ 187 if not figure: figure = plt.figure() 188 ax = figure.gca() 189 ax.autoscale(tight=True) 190 191 if style in [None, False]: 192 style = 'k-o' 193 import numpy 194 if shift: 195 if shift is True: #NOTE: MAY NOT be the exact minimum 196 shift = max(-numpy.min(cost), 0.0) + 0.5 # a good guess 197 cost = numpy.asarray(cost)+shift 198 cost = numpy.asarray(cost) 199 if scale: 200 cost = numpy.log(4*cost*scale+1)+2 201 202 ax.plot(x,cost, style, linewidth=2, markersize=4) 203 #XXX: need to 'correct' the z-axis (or provide easy conversion) 204 return figure 205 206 207 def draw_trajectory(x, y, cost=None, scale=True, shift=False, style=None, figure=None): 208 """draw a solution trajectory (for overlay on a contour plot) 209 210 x is a sequence of values for one parameter (i.e. a parameter trajectory) 211 y is a sequence of values for one parameter (i.e. a parameter trajectory) 212 cost is the solution trajectory (i.e. costs); if provided, plot a 3D contour 213 if scale is provided, scale the intensity as 'z = log(4*z*scale+1)+2' 214 if shift is provided, shift the intensity as 'z = z+shift' (useful for -z's) 215 if style is provided, set the line style (e.g. 'w-o', 'k-', 'ro') 216 if figure is provided, plot to an existing figure 217 """ 218 if not figure: figure = plt.figure() 219 220 if cost: kwds = {'projection':'3d'} # 3D 221 else: kwds = {} # 2D 222 ax = figure.gca(**kwds) 223 224 if style in [None, False]: 225 style = 'w-o' #if not scale else 'k-o' 226 if cost: # is 3D, cost is needed 227 import numpy 228 if shift: 229 if shift is True: #NOTE: MAY NOT be the exact minimum 230 shift = max(-numpy.min(cost), 0.0) + 0.5 # a good guess 231 cost = numpy.asarray(cost)+shift 232 if scale: 233 cost = numpy.asarray(cost) 234 cost = numpy.log(4*cost*scale+1)+2 235 ax.plot(x,y,cost, style, linewidth=2, markersize=4) 236 #XXX: need to 'correct' the z-axis (or provide easy conversion) 237 else: # is 2D, cost not needed 238 ax.plot(x,y, style, linewidth=2, markersize=4) 239 return figure 240 241 242 def draw_slice(f, x, y=None, scale=True, shift=False): 243 """plot a slice of a 2D function 'f' in 1D 244 245 x is an array used to set up the axis 246 y is a fixed value for the 2nd axis 247 if scale is provided, scale the intensity as 'z = log(4*z*scale+1)+2' 248 if shift is provided, shift the intensity as 'z = z+shift' (useful for -z's) 249 250 NOTE: when plotting the 'y-axis' at fixed 'x', 251 pass the array to 'y' and the fixed value to 'x' 252 """ 253 import numpy 254 255 if y is None: 256 y = 0.0 257 x, y = numpy.meshgrid(x, y) 258 plotx = True if numpy.all(y == y[0,0]) else False 259 260 z = 0*x 261 s,t = x.shape 262 for i in range(s): 263 for j in range(t): 264 xx,yy = x[i,j], y[i,j] 265 z[i,j] = f([xx,yy]) 266 if shift: 267 if shift is True: shift = max(-numpy.min(z), 0.0) + 0.5 # exact minimum 268 z = z+shift 269 if scale: z = numpy.log(4*z*scale+1)+2 270 #XXX: need to 'correct' the z-axis (or provide easy conversion) 271 272 fig = plt.figure() 273 ax = fig.gca() 274 ax.autoscale(tight=True) 275 if plotx: 276 ax.plot(x.reshape(-1), z.reshape(-1)) 277 else: 278 ax.plot(y.reshape(-1), z.reshape(-1)) 279 return fig 280 281 282 def draw_contour(f, x, y=None, surface=False, fill=True, scale=True, shift=False, density=5): 283 """draw a contour plot for a given 2D function 'f' 284 285 x and y are arrays used to set up a 2D mesh grid 286 if fill is True, color fill the contours 287 if surface is True, plot the contours as a 3D projection 288 if scale is provided, scale the intensity as 'z = log(4*z*scale+1)+2' 289 if shift is provided, shift the intensity as 'z = z+shift' (useful for -z's) 290 use density to adjust the number of contour lines 291 """ 292 import numpy 293 294 if y is None: 295 y = x 296 x, y = numpy.meshgrid(x, y) 297 298 z = 0*x 299 s,t = x.shape 300 for i in range(s): 301 for j in range(t): 302 xx,yy = x[i,j], y[i,j] 303 z[i,j] = f([xx,yy]) 304 if shift: 305 if shift is True: shift = max(-numpy.min(z), 0.0) + 0.5 # exact minimum 306 z = z+shift 307 if scale: z = numpy.log(4*z*scale+1)+2 308 #XXX: need to 'correct' the z-axis (or provide easy conversion) 309 310 fig = plt.figure() 311 if surface and fill is None: # 'hidden' option; full 3D surface plot 312 ax = fig.gca(projection='3d') 313 d = max(11 - density, 1) # or 1/density ? 314 kwds = {'rstride':d,'cstride':d,'cmap':cm.jet,'linewidth':0} 315 ax.plot_surface(x, y, z, **kwds) 316 else: 317 if surface: kwds = {'projection':'3d'} # 3D 318 elif surface is None: # 1D 319 raise NotImplementedError('need to add an option string parser') 320 else: kwds = {} # 2D 321 ax = fig.gca(**kwds) 322 density = 10*density 323 if fill: plotf = ax.contourf # filled contours 324 else: plotf = ax.contour # wire contours 325 plotf(x, y, z, density, cmap=cm.jet) 326 return fig 327 11 __doc__ = model_plotter.__doc__ 328 12 329 13 if __name__ == '__main__': 330 #FIXME: should be able to:331 # - apply a constraint as a region of NaN -- apply when 'xx,yy=x[ij],y[ij]'332 # - apply a penalty by shifting the surface (plot w/alpha?) -- as above333 # - build an appropriately-sized default grid (from logfile info)334 # - move all mulit-id param/cost reading into read_history335 #FIXME: current issues:336 # - 1D slice and projection work for 2D function, but aren't "pretty"337 # - 1D slice and projection for 1D function, is it meaningful and correct?338 # - should be able to plot from solver.genealogy (multi-monitor?) [1D,2D,3D?]339 # - should be able to scale 'z-axis' instead of scaling 'z' itself340 # (see https://github.com/matplotlib/matplotlib/issues/209)341 # - if trajectory outside contour grid, will increase bounds342 # (see support_hypercube.py for how to fix bounds)343 14 344 #XXX: note that 'argparse' is new as of python2.7 345 from optparse import OptionParser 346 parser = OptionParser(usage=__doc__) 347 parser.add_option("-b","--bounds",action="store",dest="bounds",\ 348 metavar="STR",default="-5:5:.1, -5:5:.1", 349 help="indicator string to set plot bounds and density") 350 parser.add_option("-l","--label",action="store",dest="label",\ 351 metavar="STR",default=",,", 352 help="string to assign label to axis") 353 parser.add_option("-n","--nid",action="store",dest="id",\ 354 metavar="INT",default=None, 355 help="id # of the nth simultaneous points to plot") 356 parser.add_option("-i","--iter",action="store",dest="stop",\ 357 metavar="STR",default=":", 358 help="string for smallest:largest iterations to plot") 359 parser.add_option("-r","--reduce",action="store",dest="reducer",\ 360 metavar="STR",default="None", 361 help="import path of output reducer function") 362 parser.add_option("-x","--scale",action="store",dest="zscale",\ 363 metavar="INT",default=0.0, 364 help="scale plotted cost by z=log(4*z*scale+1)+2") 365 parser.add_option("-z","--shift",action="store",dest="zshift",\ 366 metavar="INT",default=0.0, 367 help="shift plotted cost by z=z+shift") 368 parser.add_option("-f","--fill",action="store_true",dest="fill",\ 369 default=False,help="plot using filled contours") 370 parser.add_option("-d","--depth",action="store_true",dest="surface",\ 371 default=False,help="plot contours showing depth in 3D") 372 parser.add_option("-o","--dots",action="store_true",dest="dots",\ 373 default=False,help="show trajectory points in plot") 374 parser.add_option("-j","--join",action="store_true",dest="line",\ 375 default=False,help="connect trajectory points in plot") 376 parsed_opts, parsed_args = parser.parse_args() 377 378 # get the import path for the model 379 model = parsed_args[0] # e.g. 'mystic.models.rosen' 380 if "None" == model: model = None #XXX: 'required'... allow this? 381 382 try: # get the name of the parameter log file 383 source = parsed_args[1] # e.g. 'log.txt' 384 except: 385 source = None 386 387 try: # select the bounds 388 options = parsed_opts.bounds # format is "-1:10:.1, -1:10:.1, 1.0" 389 except: 390 options = "-5:5:.1, -5:5:.1" 391 392 try: # plot using filled contours 393 fill = parsed_opts.fill 394 except: 395 fill = False 396 397 try: # plot contours showing depth in 3D 398 surface = parsed_opts.surface 399 except: 400 surface = False 401 402 #XXX: can't do '-x' with no argument given (use T/F instead?) 403 try: # scale plotted cost by z=log(4*z*scale+1)+2 404 scale = float(parsed_opts.zscale) 405 if not scale: scale = False 406 except: 407 scale = False 408 409 #XXX: can't do '-z' with no argument given 410 try: # shift plotted cost by z=z+shift 411 shift = float(parsed_opts.zshift) 412 if not shift: shift = False 413 except: 414 shift = False 415 416 try: # import path of output reducer function 417 reducer = parsed_opts.reducer # e.g. 'numpy.add' 418 if "None" == reducer: reducer = None 419 except: 420 reducer = None 421 422 style = '-' # default linestyle 423 if parsed_opts.dots: 424 mark = 'o' # marker=mark 425 # when using 'dots', also can turn off 'line' 426 if not parsed_opts.line: 427 style = '' # linestyle='None' 428 else: 429 mark = '' 430 color = 'w' if fill else 'k' 431 style = color + style + mark 432 433 try: # select labels for the axes 434 label = parsed_opts.label.split(',') # format is "x, y, z" 435 except: 436 label = ['','',''] 437 438 try: # select which 'id' to plot results for 439 ids = (int(parsed_opts.id),) #XXX: allow selecting more than one id ? 440 except: 441 ids = None # i.e. 'all' 442 443 try: # select which iteration to stop plotting at 444 stop = parsed_opts.stop # format is "1:10:1" 445 stop = stop if ":" in stop else ":"+stop 446 except: 447 stop = ":" 448 449 ################################################# 450 solver = None # set to 'mystic.solvers.fmin' (or similar) for 'live' fits 451 #NOTE: 'live' runs constrain params explicitly in the solver, then reduce 452 # dimensions appropriately so results can be 2D contour plotted. 453 # When working with legacy results that have more than 2 params, 454 # the trajectory WILL NOT follow the masked surface generated 455 # because the masked params were NOT fixed when the solver was run. 456 ################################################# 457 458 from mystic.tools import reduced, masked, partial 459 460 # process inputs 461 select, spec, mask = parse_input(options) 462 x,y = parse_axes(spec, grid=True) # grid=False for 1D plots 463 #FIXME: does grid=False still make sense here...? 464 if reducer: reducer = get_instance(reducer) 465 if solver and (not source or not model): 466 raise RuntimeError('a model and results filename are required') 467 elif not source and not model: 468 raise RuntimeError('a model or a results file is required') 469 if model: 470 model = get_instance(model) 471 # need a reducer if model returns an array 472 if reducer: model = reduced(reducer, arraylike=False)(model) 473 474 if solver: 475 # if 'live'... pick a solver 476 solver = 'mystic.solvers.fmin' 477 solver = get_instance(solver) 478 xlen = len(select)+len(mask) 479 if solver.__name__.startswith('diffev'): 480 initial = [(-1,1)]*xlen 481 else: 482 initial = [0]*xlen 483 from mystic.monitors import VerboseLoggingMonitor 484 itermon = VerboseLoggingMonitor(filename=source, new=True) 485 # explicitly constrain parameters 486 model = partial(mask)(model) 487 # solve 488 sol = solver(model, x0=initial, itermon=itermon) 489 490 #-OVERRIDE-INPUTS-# 491 import numpy 492 # read trajectories from monitor (comment out to use logfile) 493 source = itermon 494 # if negative minimum, shift by the 'solved minimum' plus an epsilon 495 shift = max(-numpy.min(itermon.y), 0.0) + 0.5 # a good guess 496 #-----------------# 497 498 if model: # for plotting, implicitly constrain by reduction 499 model = masked(mask)(model) 500 501 ## plot the surface in 1D 502 #if solver: v=sol[-1] 503 #elif source: v=cost[-1] 504 #else: v=None 505 #fig0 = draw_slice(model, x=x, y=v, scale=scale, shift=shift) 506 # plot the surface in 2D or 3D 507 fig = draw_contour(model, x, y, surface=surface, fill=fill, scale=scale, shift=shift) 508 else: 509 #fig0 = None 510 fig = None 511 512 if source: 513 # params are the parameter trajectories 514 # cost is the solution trajectory 515 params, cost = get_history(source, ids) 516 if len(cost) > 1: style = style[1:] # 'auto-color' #XXX: or grayscale? 517 518 for p,c in zip(params, cost): 519 ## project trajectory on a 1D slice of model surface #XXX: useful? 520 #s = select[0] if len(select) else 0 521 #px = p[int(s)] # draw_projection requires one parameter 522 ## ignore everything after 'stop' 523 #_c = eval('c[%s]' % stop) 524 #_x = eval('px[%s]' % stop) 525 #fig0 = draw_projection(_x,_c, style=style, scale=scale, shift=shift, figure=fig0) 526 527 # plot the trajectory on the model surface (2D or 3D) 528 # get two selected params #XXX: what if len(select)<2? or len(p)<2? 529 p = [p[int(i)] for i in select[:2]] 530 px,py = p # draw_trajectory requires two parameters 531 # ignore everything after 'stop' 532 _x = eval('px[%s]' % stop) 533 _y = eval('py[%s]' % stop) 534 _c = eval('c[%s]' % stop) if surface else None 535 fig = draw_trajectory(_x,_y,_c, style=style, scale=scale, shift=shift, figure=fig) 536 537 # add labels to the axes 538 if surface: kwds = {'projection':'3d'} # 3D 539 else: kwds = {} # 2D 540 ax = fig.gca(**kwds) 541 ax.set_xlabel(label[0]) 542 ax.set_ylabel(label[1]) 543 if surface: ax.set_zlabel(label[2]) 544 545 plt.show() 15 import sys 16 model_plotter(cmdargs=sys.argv[1:]) 546 17 547 18
Note: See TracChangeset
for help on using the changeset viewer.