- Timestamp:
- 09/21/12 15:07:30 (4 years ago)
- File:
-
- 1 moved
Legend:
- Unmodified
- Added
- Removed
-
branches/UQ/math/legacy/support_hypercube_scenario.py
r565 r566 1 1 #!/usr/bin/env python 2 2 __doc__ = """ 3 support_hypercube_scenario.py [options] filename datafile4 5 generate scenario support plots from file written with 'write_support_file' 6 and a dataset file 7 8 The options "bounds", " axes", and "iters" all take indicator strings.3 support_hypercube_scenario.py [options] filename (datafile) 4 5 generate scenario support plots from file written with 'write_support_file'; 6 generate legacy data and cones from a dataset file, if provided 7 8 The options "bounds", "dim", and "iters" all take indicator strings. 9 9 The bounds should be given as a quoted list of tuples. For example, using 10 10 bounds = "[(.062,.125),(0,30),(2300,3200)]" will set the lower and upper bounds 11 for x to be (.062,.125), y to be (0,30), and z to be (2300,3200). Similarly, 12 axes also accepts a quoted list of tuples; however, for axes, the first tuple 13 indicates which parameters are along the x direction, the second tuple for 14 the y direction, and the third tuple for the z direction. Thus, axes = 15 "[(2,3),(6,7),(10,11)]" would set the 2nd and 3rd parameters along x. Iters, 16 however, accepts a list of strings instead of a list of tuples. For example, 11 for x to be (.062,.125), y to be (0,30), and z to be (2300,3200). I bounds 12 are to not be strictly enforced, append an asterisk '*' to the bounds. 13 The dim (dimensions of the scenario) should be a quoted tuple. For example, 14 dim = "(1,1,2)" will convert the params to a two-member 3-D dataset. Iters 15 takes a list of strings instead of a tuple or a list of tuples. For example, 17 16 iters = "[':']" will plot all iters in a single plot. Alternatively, 18 17 iters = "[':2','2:']" will split the iters into two plots, while … … 28 27 Required Inputs: 29 28 filename name of the python convergence logfile (e.g. paramlog.py) 29 30 Additional Inputs: 30 31 datafile name of the dataset textfile (e.g. StAlDataset.txt) 31 32 """ … … 33 34 34 35 #### building the cone primitive #### 35 def cone_builder(slope, bounds ):36 def cone_builder(slope, bounds, strict=True): 36 37 """ factory to create a cone primitive 37 38 … … 66 67 lb,ub = zip(*bounds) 67 68 # if False, the cone surface may violate bounds 68 69 #strict = True # always respect the bounds 69 70 70 71 def cone(position, top=True): … … 97 98 98 99 #### plotting the cones #### 99 def plot_cones(ax, data, slope, bounds, color='0.75', axis=None ):100 def plot_cones(ax, data, slope, bounds, color='0.75', axis=None, strict=True): 100 101 """plot double cones for a given dataset 101 102 … … 111 112 bounds = swap(bounds, axis) 112 113 data = [swap(pt,axis) for pt in data] 113 cone = cone_builder(slope, bounds )114 cone = cone_builder(slope, bounds, strict=strict) 114 115 # plot the upper and lower cone 115 116 for datapt in data: … … 117 118 if _cone: 118 119 X,Z,Y = swap(_cone, axis) # 'unswap' the axes 119 ax.plot_surface(X, Y,Z, rstride=1, cstride=1, color=color) 120 ax.plot_surface(X, Y,Z, rstride=1, cstride=1, color=color, \ 121 linewidths=0, alpha=.1) 120 122 _cone = cone(datapt, top=False) 121 123 if _cone: 122 124 X,Z,Y = swap(_cone, axis) # 'unswap' the axes 123 ax.plot_surface(X, Y,Z, rstride=1, cstride=1, color=color) 125 ax.plot_surface(X, Y,Z, rstride=1, cstride=1, color=color, \ 126 linewidths=0, alpha=.1) 124 127 return ax 125 128 126 def plot_data(ax, data, bounds, color='red' ):129 def plot_data(ax, data, bounds, color='red', strict=True): 127 130 """plot datapoints for a given dataset 128 131 … … 132 135 color -- string name (or rbg value) of color to use for datapoints 133 136 """ 134 137 # strict = True # always respect the bounds 135 138 lb,ub = zip(*bounds) 136 139 # plot the datapoints themselves … … 140 143 d_hi = ub - position 141 144 d_lo = lb - position 142 if strict and (any(d_hi < 0) or any(d_lo > 0)): 145 if strict and (any(d_hi < 0) or any(d_lo > 0)): #XXX: any or just in Y ? 143 146 pass # don't plot if datapt is out of bounds 144 147 else: … … 170 173 ax.set_zlabel(labels[2]) # cone "center" axis 171 174 return ax 175 176 def get_slope(data, replace=None): 177 """ replace one slope in a list of slopes with '1.0' 178 (i.e. replace a slope from the dataset with the unit slope) 179 180 data -- dataset object, where coordinates are 3-tuples and values are floats 181 replace -- selected axis (an int) to plot values NOT coords 182 """ 183 slope = data.lipschitz 184 if replace not in range(len(slope)): # don't replace an axis 185 return slope 186 return slope[:replace] + [1.0] + slope[replace+1:] 172 187 173 188 def get_coords(data, replace=None): … … 182 197 values = data.values 183 198 if replace not in range(len(slope)): # don't replace an axis 184 return coords, slope 185 slope = slope[:replace] + [1.0] + slope[replace+1:] 186 coords = [coords[i][:replace] + [values[i]] + coords[i][replace+1:] \ 187 for i in range(len(coords))] 188 return coords, slope 199 return coords 200 return [list(coords[i][:replace]) + [values[i]] + \ 201 list(coords[i][replace+1:]) for i in range(len(coords))] 189 202 190 203 def swap(alist, index=None): … … 198 211 return alist[:index] + alist[index+1:] + alist[index:index+1] 199 212 213 from support_convergence import best_dimensions 214 200 215 201 216 if __name__ == '__main__': … … 208 223 metavar="STR",default="[(0,1),(0,1),(0,1)]", 209 224 help="indicator string to set hypercube bounds") 210 # parser.add_option("-x","--axes",action="store",dest="xyz",\211 # metavar="STR",default="[(0,),(1,),(2,)]",212 # help="indicator string to assign parameter to axis")225 parser.add_option("-i","--iters",action="store",dest="iters",\ 226 metavar="STR",default="['-1']", 227 help="indicator string to select iterations to plot") 213 228 parser.add_option("-l","--label",action="store",dest="label",\ 214 229 metavar="STR",default="['','','']", 215 230 help="string to assign label to axis") 231 parser.add_option("-d","--dim",action="store",dest="dim",\ 232 metavar="STR",default="(1,1,1)", 233 help="indicator string set to scenario dimensions") 216 234 parser.add_option("-r","--replace",action="store",dest="replace",\ 217 235 metavar="INT",default=None, 218 236 help="id # of axis to plot values NOT coords") 219 parser.add_option("-v","--vertical",action="store_true",dest="vertical",\ 220 default=False,help="always plot cones on vertical axis") 237 parser.add_option("-n","--nid",action="store",dest="id",\ 238 metavar="INT",default=None, 239 help="id # of the nth simultaneous points to plot") 240 parser.add_option("-s","--scale",action="store",dest="scale",\ 241 metavar="FLOAT",default=1.0, 242 help="multiplier for scaling color of points in plot") 243 parser.add_option("-o","--datapts",action="store_true",dest="legacy",\ 244 default=False,help="plot legacy data, if provided") 245 parser.add_option("-v","--cones",action="store_true",dest="cones",\ 246 default=False,help="plot cones, if provided") 247 parser.add_option("-z","--vertical",action="store_true",dest="vertical",\ 248 default=False,help="always plot values on vertical axis") 249 parser.add_option("-f","--flat",action="store_true",dest="flatten",\ 250 default=False,help="show selected iterations in a single plot") 221 251 parsed_opts, parsed_args = parser.parse_args() 222 252 223 # try: # get the name of the parameter log file 224 # file = parsed_args[0] 225 # import re 226 # file = re.sub('\.py*.$', '', file) #XXX: strip off .py* extension 227 # except: 228 # raise IOError, "please provide log file name" 229 # try: # read standard logfile 230 # from mystic.munge import logfile_reader, raw_to_support 231 # _step, params, _cost = logfile_reader(file) 232 # params, _cost = raw_to_support(params, _cost) 233 # except: 234 # exec "from %s import params" % file 235 # #exec "from %s import meta" % file 236 # # would be nice to use meta = ['wx','wx2','x','x2','wy',...] 253 try: # get the name of the parameter log file 254 file = parsed_args[0] 255 import re 256 file = re.sub('\.py*.$', '', file) #XXX: strip off .py* extension 257 except: 258 raise IOError, "please provide log file name" 259 try: # read standard logfile 260 from mystic.munge import logfile_reader, raw_to_support 261 _step, params, _cost = logfile_reader(file) 262 params, _cost = raw_to_support(params, _cost) 263 except: 264 exec "from %s import params" % file 265 #exec "from %s import meta" % file 266 # would be nice to use meta = ['wx','wx2','x','x2','wy',...] 267 268 from dirac_measure import scenario 269 from legacydata import dataset 270 try: # select whether to plot the cones 271 cones = parsed_opts.cones 272 except: 273 cones = False 274 275 try: # select whether to plot the legacy data 276 legacy = parsed_opts.legacy 277 except: 278 legacy = False 237 279 238 280 try: # get the name of the dataset file 239 file = parsed_args[ 0]281 file = parsed_args[1] 240 282 from legacydata import load_dataset 241 283 data = load_dataset(file) 242 284 except: 243 raise IOError, "please provide dataset file name" 285 # raise IOError, "please provide dataset file name" 286 data = dataset() 287 cones = False 288 legacy = False 289 290 try: # select the scenario dimensions 291 npts = eval(parsed_opts.dim) # format is "(1,1,1)" 292 except: 293 npts = (1,1,1) #XXX: better in parsed_args ? 244 294 245 295 try: # select the bounds 246 bounds = eval(parsed_opts.bounds) # format is "[(60,105),(0,30),(2.1,2.8)]" 247 except: 296 _bounds = parsed_opts.bounds 297 if _bounds[-1] == "*": #XXX: then bounds are NOT strictly enforced 298 _bounds = _bounds[:-1] 299 strict = False 300 else: 301 strict = True 302 bounds = eval(_bounds) # format is "[(60,105),(0,30),(2.1,2.8)]" 303 except: 304 strict = True 248 305 bounds = [(0,1),(0,1),(0,1)] 249 306 … … 253 310 label = ['','',''] 254 311 312 x = params[-1] 313 try: # select which iterations to plot 314 select = eval(parsed_opts.iters) # format is "[':2','2:4','5','6:']" 315 except: 316 select = ['-1'] 317 #select = [':'] 318 #select = [':1'] 319 #select = [':2','2:'] 320 #select = [':1','1:2','2:3','3:'] 321 #select = ['0','1','2','3'] 322 323 try: # collapse non-consecutive iterations into a single plot... 324 flatten = parsed_opts.flatten 325 except: 326 flatten = False 327 328 try: # select which 'id' to plot results for 329 id = int(parsed_opts.id) 330 except: 331 id = None # i.e. 'all' **or** use id=0, which should be 'best' energy ? 332 333 try: # scale the color in plotting the weights 334 scale = float(parsed_opts.scale) 335 except: 336 scale = 1.0 # color = color**scale 337 255 338 try: # select which axis to plot 'values' on 256 339 xs = int(parsed_opts.replace) … … 263 346 vertical_cones = False 264 347 265 # ensure all terms of bounds a nd xyz are tuples348 # ensure all terms of bounds are tuples 266 349 for bound in bounds: 267 350 if not isinstance(bound, tuple): 268 351 raise TypeError, "bounds should be tuples of (lower_bound,upper_bound)" 269 # for i in range(len(xyz)): 270 # if isinstance(xyz[i], int): 271 # xyz[i] = (xyz[i],) 272 # elif not isinstance(xyz[i], tuple): 273 # raise TypeError, "xyz should be tuples of (param1,param2,param3,...)" 274 275 # get coords (and values) for selected axes 276 coords, slope = get_coords(data, xs) 277 print "bounds: %s" % bounds 278 print "slope: %s" % slope 279 #print "coords: %s" % coords 352 353 # ensure all terms of select are strings that have a ":" 354 for i in range(len(select)): 355 if isinstance(select[i], int): select[i] = str(select[i]) 356 if select[i] == '-1': select[i] = 'len(x)-1:len(x)' 357 elif not select[i].count(':'): 358 select[i] += ':' + str(int(select[i])+1) 359 360 # take only the selected 'id' 361 if id != None: 362 param = [] 363 for j in range(len(params)): 364 param.append([p[id] for p in params[j]]) 365 params = param[:] 366 367 # at this point, we should have: 368 #bounds = [(60,105),(0,30),(2.1,2.8)] or [(None,None),(None,None),(None,None)] 369 #select = ['-1:'] or [':'] or [':1','1:2','2:3','3:'] or similar 370 #id = 0 or None 371 372 # get dataset coords (and values) for selected axes 373 if data: 374 slope = get_slope(data, xs) 375 coords = get_coords(data, xs) 376 #print "bounds: %s" % bounds 377 #print "slope: %s" % slope 378 #print "coords: %s" % coords 379 #else: 380 # slope = [] 381 # coords = [] 382 383 from mpl_toolkits.mplot3d import Axes3D 384 import matplotlib.pyplot as plt 385 from matplotlib.axes import subplot_class_factory 386 Subplot3D = subplot_class_factory(Axes3D) 387 388 plots = len(select) 389 if not flatten: 390 dim1,dim2 = best_dimensions(plots) 391 else: dim1,dim2 = 1,1 280 392 281 393 # use the default bounds where not specified … … 285 397 if bounds[i][1] is None: bounds[i][1] = 1 286 398 287 import matplotlib 288 import matplotlib.pyplot as plt 289 from mpl_toolkits.mplot3d import Axes3D 399 # swap the axes appropriately when plotting cones 400 if vertical_cones and xs in range(len(bounds)): # we are replacing an axis 401 # adjust slope, bounds, and data so cone axis is last 402 if data: 403 slope = swap(slope, xs) 404 coords = [swap(pt,xs) for pt in coords] 405 bounds = swap(bounds, xs) 406 label = swap(label, xs) 407 axis = None 408 else: 409 axis = xs 410 411 # correctly bound the first plot. there must be at least one plot 290 412 fig = plt.figure() 291 ax = Axes3D(fig)#,aspect='equal') 292 if xs in range(len(bounds)): # we are replacing an axis 293 if vertical_cones: 294 # adjust slope, bounds, and data so cone axis is last 295 slope = swap(slope, xs) 296 bounds = swap(bounds, xs) 297 coords = [swap(pt,xs) for pt in coords] 298 label = swap(label, xs) 299 xs = None 300 plot_cones(ax, coords, slope, bounds, axis=xs) 301 plot_data(ax, coords, bounds) 302 #clip_axes(ax, bounds) 303 label_axes(ax, label) 413 ax1 = Subplot3D(fig, dim1,dim2,1) 414 ax1.plot([bounds[0][0]],[bounds[1][0]],[bounds[2][0]]) 415 ax1.plot([bounds[0][1]],[bounds[1][1]],[bounds[2][1]]) 416 if not flatten: 417 exec "plt.title('iterations[%s]')" % select[0] 418 else: 419 exec "plt.title('iterations[*]')" 420 if cones and data and axis in range(len(bounds)): 421 plot_cones(ax1,coords,slope,bounds,axis=axis,strict=strict) 422 if legacy and data: 423 plot_data(ax1,coords,bounds,strict=strict) 424 #clip_axes(ax1,bounds) 425 label_axes(ax1,label) 426 a = [ax1] 427 428 # set up additional plots 429 if not flatten: 430 for i in range(2, plots + 1): 431 exec "ax%d = Subplot3D(fig, dim1,dim2,%d)" % (i,i) 432 exec "ax%d.plot([bounds[0][0]],[bounds[1][0]],[bounds[2][0]])" % i 433 exec "ax%d.plot([bounds[0][1]],[bounds[1][1]],[bounds[2][1]])" % i 434 exec "plt.title('iterations[%s]')" % select[i - 1] 435 if cones and data and axis in range(len(bounds)): 436 exec "plot_cones(ax%d,coords,slope,bounds,axis=axis,strict=strict)" % i 437 if legacy and data: 438 exec "plot_data(ax%d,coords,bounds,strict=strict)" % i 439 #exec "clip_axes(ax%d,bounds)" % i 440 exec "label_axes(ax%d,label)" % i 441 exec "a.append(ax%d)" % i 442 443 # turn each "n:m" in select to a list 444 _select = [] 445 for sel in select: 446 if sel[0] == ':': _select.append("0"+sel) 447 else: _select.append(sel) 448 for i in range(len(_select)): 449 if _select[i][-1] == ':': select[i] = _select[i]+str(len(x)) 450 else: select[i] = _select[i] 451 for i in range(len(select)): 452 p = select[i].split(":") 453 if p[0][0] == '-': p[0] = "len(x)"+p[0] 454 if p[1][0] == '-': p[1] = "len(x)"+p[1] 455 select[i] = p[0]+":"+p[1] 456 steps = [eval("range(%s)" % sel.replace(":",",")) for sel in select] 457 458 # at this point, we should have: 459 #steps = [[0,1],[1,2],[2,3],[3,4,5,6,7,8]] or similar 460 if flatten: 461 from mystic.tools import flatten 462 steps = [list(flatten(steps))] 463 464 # plot all the scenario "data" 465 from numpy import inf, e 466 scale = e**(scale - 1.0) 467 for v in range(len(steps)): 468 if len(steps[v]) > 1: qp = float(max(steps[v])) 469 else: qp = inf 470 for s in steps[v]: 471 par = eval("[params[q][%s][0] for q in range(len(params))]" % s) 472 pm = scenario() 473 pm.load(par, npts) 474 d = dataset() 475 d.load(pm.coords, pm.values) 476 # dot color determined by number of simultaneous iterations 477 t = str((s/qp)**scale) 478 # get and plot dataset coords for selected axes 479 plot_data(a[v], get_coords(d, xs), bounds, color=t, strict=strict) 480 304 481 plt.show() 305 482
Note: See TracChangeset
for help on using the changeset viewer.