Changeset 809 for mystic


Ignore:
Timestamp:
07/25/15 15:16:13 (10 months ago)
Author:
mmckerns
Message:

added functional interface for scenario plotter

Location:
mystic
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • mystic/mystic/support.py

    r808 r809  
    99""" 
    1010 
    11 __all__ = ['convergence', 'hypercube', 'hypercube_measures'] 
     11__all__ = ['convergence', 'hypercube', 'hypercube_measures', \ 
     12           'hypercube_scenario', 'best_dimensions', 'swap'] 
    1213 
    1314from mpl_toolkits.mplot3d import Axes3D as _Axes3D 
     
    2526# globals 
    2627__quit = False 
     28ZERO = 1.0e-6  # zero-ish 
    2729 
    2830 
     
    3941#   return cand[len(cand)/2], cand[len(cand)/2] 
    4042# return cand[len(cand)/2], cand[len(cand)/2 - 1] 
     43 
     44 
     45#### building the cone primitive #### 
     46def _cone_builder(slope, bounds, strict=True): 
     47  """ factory to create a cone primitive 
     48 
     49  slope -- slope multiplier for cone on the X,Y,Z axes (for mesh construction) 
     50  bounds -- list of tuples of bounds for the plot; (lower,upper) for each axis 
     51""" 
     52  from mystic.math import almostEqual 
     53  import numpy as np 
     54  def cone_mesh(length): 
     55    """ construct a conical mesh for a given length of cone """ 
     56    L1,L2,L3 = slope 
     57    radius = length / L3 #XXX: * 0.5 
     58    r0 = ZERO 
     59 
     60    if almostEqual(radius, r0, tol=r0): radius = r0 
     61    r = np.linspace(radius,radius,6)  
     62    r[0]= np.zeros(r[0].shape)  
     63    r[1] *= r0/radius 
     64    r[5] *= r0/radius 
     65    r[3]= np.zeros(r[3].shape)  
     66 
     67    p = np.linspace(0,2*np.pi,50)  
     68    R,P = np.meshgrid(r,p)  
     69    X,Y = L1 * R*np.cos(P), L2 * R*np.sin(P)  
     70 
     71    tmp=list()  
     72    for i in range(np.size(p)):  
     73      tmp.append([0,0,length,length,length,0]) # len = size(r) 
     74    Z = np.array(tmp)  
     75    return X,Z,Y 
     76 
     77  lb,ub = zip(*bounds) 
     78  # if False, the cone surface may violate bounds 
     79 #strict = True # always respect the bounds 
     80 
     81  def cone(position, top=True): 
     82    """ construct a cone primitive, with vertex at given position 
     83 
     84    position -- tuple of vertex position 
     85    top -- boolean for if 'top' or 'bottom' cone 
     86""" 
     87    from numpy import asarray 
     88    position = asarray(position) 
     89    d_hi = ub - position 
     90    d_lo = lb - position 
     91    if strict and (any(d_hi < 0) or any(d_lo > 0)): 
     92      return None  # don't plot cone if vertex is out of bounds 
     93 
     94    if top: # distance of vertex from upper edge 
     95      l = d_hi[2] 
     96      if not l: l = ZERO 
     97    else:   # distance of vertex from lower edge 
     98      l = d_lo[2] 
     99      if not l: l = -ZERO 
     100    X,Z,Y = cone_mesh(l) 
     101    X = X + position[0] 
     102    Y = Y + position[1] 
     103    Z = Z + position[2] 
     104    return X,Z,Y 
     105 
     106  return cone     
     107 
     108 
     109#### plotting the cones #### 
     110def _plot_bowtie(ax, data, slope, bounds, color='0.75', axis=None, tol=0.0): 
     111  """plot (2D) double cones for a given dataset 
     112 
     113  ax -- matplotlib 'Axes3D' plot object 
     114  data -- list of datapoints, where datapoints are 3-tuples (i.e. x,y,z) 
     115  slope -- slope multiplier for cone on the X,Y,Z axes (for mesh construction) 
     116  bounds -- list of tuples of bounds for the plot; (lower,upper) for each axis 
     117  color -- string name (or rbg value) of color to use for datapoints 
     118  axis -- the axis of the cone 
     119  tol -- distance between center of mass of the double cones and a cone vertex 
     120""" 
     121  if axis not in range(len(bounds)-1): return ax 
     122  from numpy import asarray, inf 
     123  data = asarray(data) 
     124  sl = slope[axis] 
     125  gap = max(0., tol) 
     126  # find the endpoints of the cone:  y = slope * (x0 - x) + y0 
     127  ylo = sl * (bounds[0][0] - data.T[0]) + data.T[1] 
     128  yhi = sl * (bounds[0][1] - data.T[0]) + data.T[1] 
     129  zhi = -sl * (bounds[0][0] - data.T[0]) + data.T[1] 
     130  zlo = -sl * (bounds[0][1] - data.T[0]) + data.T[1] 
     131  xlb = bounds[0][0] 
     132  xub = bounds[0][1] 
     133  ylb = bounds[1][0] 
     134  yub = bounds[1][1] 
     135 
     136  for pt,yl,yu,zl,zu in zip(data,ylo,yhi,zlo,zhi): 
     137   #ax.plot([xlb, pt[0], xub], [yl, pt[1], yu], color='black') 
     138   #ax.plot([xlb, pt[0], xub], [zu, pt[1], zl], color='black') 
     139    ax.fill_between([xlb, pt[0], xub], [ylb]*3, [yl-gap, pt[1]-gap, zl-gap], \ 
     140                                     facecolor=color, alpha=0.2) 
     141    ax.fill_between([xlb, pt[0], xub], [yub]*3, [zu+gap, pt[1]+gap, yu+gap], \ 
     142                                     facecolor=color, alpha=0.2) 
     143  return ax 
     144 
     145 
     146def _plot_cones(ax, data, slope, bounds, color='0.75', axis=None, strict=True, tol=0.0): 
     147  """plot (3D) double cones for a given dataset 
     148 
     149  ax -- matplotlib 'Axes3D' plot object 
     150  data -- list of datapoints, where datapoints are 3-tuples (i.e. x,y,z) 
     151  slope -- slope multiplier for cone on the X,Y,Z axes (for mesh construction) 
     152  bounds -- list of tuples of bounds for the plot; (lower,upper) for each axis 
     153  color -- string name (or rbg value) of color to use for datapoints 
     154  axis -- the axis of the cone 
     155  tol -- distance between center of mass of the double cones and a cone vertex 
     156""" 
     157  from numpy import asarray, zeros 
     158  # adjust slope, bounds, and data so cone axis is last  
     159  slope = swap(slope, axis)  
     160  bounds = swap(bounds, axis)  
     161  data = [swap(pt,axis) for pt in data] 
     162  cone = _cone_builder(slope, bounds, strict=strict) 
     163  # prepare to apply the 'gap' tolerance 
     164  gap = zeros(len(slope)) 
     165  gap[-1:] = [max(0., tol)] #XXX: bad behavior if len(slope) is 0 
     166  # plot the upper and lower cone 
     167  for datapt in data: 
     168    datapt = asarray(datapt) 
     169    _cone = cone(datapt + gap, top=True) 
     170    if _cone: 
     171      X,Z,Y = swap(_cone, axis) # 'unswap' the axes 
     172      ax.plot_surface(X, Y,Z, rstride=1, cstride=1, color=color, \ 
     173                                         linewidths=0, alpha=.1) 
     174    _cone = cone(datapt - gap, top=False) 
     175    if _cone: 
     176      X,Z,Y = swap(_cone, axis) # 'unswap' the axes 
     177      ax.plot_surface(X, Y,Z, rstride=1, cstride=1, color=color, \ 
     178                                         linewidths=0, alpha=.1) 
     179  return ax 
     180 
     181 
     182def _plot_data(ax, data, bounds, color='red', strict=True, **kwds): 
     183  """plot datapoints for a given dataset 
     184 
     185  ax -- matplotlib 'Axes3D' plot object 
     186  data -- list of datapoints, where datapoints are 3-tuples (i.e. x,y,z) 
     187  bounds -- list of tuples of bounds for the plot; (lower,upper) for each axis 
     188  color -- string name (or rbg value) of color to use for datapoints 
     189""" 
     190  _2D = kwds.get('_2D', False) 
     191# strict = True # always respect the bounds 
     192  lb,ub = zip(*bounds) 
     193  # plot the datapoints themselves 
     194  from numpy import asarray 
     195  for datapt in data: 
     196    position = asarray(datapt) 
     197    d_hi = ub - position 
     198    d_lo = lb - position 
     199    if strict and (any(d_hi < 0) or any(d_lo > 0)): #XXX: any or just in Y ? 
     200      pass  # don't plot if datapt is out of bounds 
     201    else: 
     202      if _2D: 
     203        ax.plot([datapt[0]], [datapt[1]], \ 
     204                color=color, marker='o', markersize=10) 
     205      else: 
     206        ax.plot([datapt[0]], [datapt[1]], [datapt[2]], \ 
     207                color=color, marker='o', markersize=10) 
     208  return ax 
     209 
     210 
     211def _clip_axes(ax, bounds): 
     212  """ clip plots to be within given lower and upper bounds 
     213 
     214  ax -- matplotlib 'Axes3D' plot object 
     215  bounds -- list of tuples of bounds for the plot; (lower,upper) for each axis 
     216""" 
     217  lb,ub = zip(*bounds) 
     218  # plot only within [lb,ub] 
     219  ax.set_xlim3d(lb[0], ub[0]) 
     220  ax.set_ylim3d(lb[1], ub[1]) 
     221  ax.set_zlim3d(lb[2], ub[2]) # cone "center" axis 
     222  return ax 
     223 
     224 
     225def _label_axes(ax, labels, **kwds): 
     226  """ label plots with given string labels 
     227 
     228  ax -- matplotlib 'Axes3D' plot object 
     229  labels -- list of string labels for the plot 
     230""" 
     231  _2D = kwds.get('_2D', False) 
     232  ax.set_xlabel(labels[0]) 
     233  ax.set_ylabel(labels[1]) 
     234  if not _2D: 
     235    ax.set_zlabel(labels[2]) # cone "center" axis 
     236  return ax 
     237 
     238 
     239def _get_slope(data, replace=None, mask=None): 
     240  """ replace one slope in a list of slopes with '1.0' 
     241  (i.e. replace a slope from the dataset with the unit slope) 
     242 
     243  data -- dataset object, where coordinates are 3-tuples and values are floats 
     244  replace -- selected axis (an int) to plot values NOT coords 
     245  """ 
     246  slope = data.lipschitz 
     247  if mask in range(len(slope)): 
     248    slope = swap(slope, mask) 
     249  if replace not in range(len(slope)):  # don't replace an axis 
     250    return slope 
     251  return slope[:replace] + [1.0] + slope[replace+1:] 
     252 
     253 
     254def _get_coords(data, replace=None, mask=None): 
     255  """ replace one coordiate axis in a 3-D data set with 'data values' 
     256  (i.e. replace an 'x' axis with the 'y' values of the data) 
     257 
     258  data -- dataset object, where coordinates are 3-tuples and values are floats 
     259  replace -- selected axis (an int) to plot values NOT coords 
     260  """ 
     261  slope = data.lipschitz 
     262  coords = data.coords 
     263  values = data.values 
     264  if mask in range(len(slope)): 
     265    coords = [swap(pt,mask) for pt in coords] 
     266  if replace not in range(len(slope)):  # don't replace an axis 
     267    return coords 
     268  return [list(coords[i][:replace]) + [values[i]] + \ 
     269          list(coords[i][replace+1:]) for i in range(len(coords))] 
     270 
     271 
     272def swap(alist, index=None): 
     273  """ swap the selected list element with the last element in alist 
     274 
     275  alist -- a list of objects 
     276  index -- the selected element 
     277  """ 
     278  if index not in range(len(alist)):  # don't swap an element 
     279    return alist  
     280  return alist[:index] + alist[index+1:] + alist[index:index+1] 
    41281 
    42282 
     
    8391079 
    8401080 
     1081def hypercube_scenario(filename, datafile=None, **kwds): 
     1082    """ 
     1083generate scenario support plots from file written with 'write_support_file'; 
     1084generate legacy data and cones from a dataset file, if provided 
     1085 
     1086Available from the command shell as: 
     1087  support_hypercube_scenario.py filename (datafile) [options] 
     1088 
     1089or as a function call as: 
     1090  mystic.support.hypercube_scenario(filename, datafile=None, **options) 
     1091 
     1092The options "bounds", "dim", and "iters" all take indicator strings. 
     1093The bounds should be given as comma-separated slices. For example, using 
     1094bounds = ".062:.125, 0:30, 2300:3200" will set the lower and upper bounds 
     1095for x to be (.062,.125), y to be (0,30), and z to be (2300,3200). If all 
     1096bounds are to not be strictly enforced, append an asterisk '*' to the string.  
     1097The dim (dimensions of the scenario) should comma-separated ints. For example, 
     1098dim = "1, 1, 2" will convert the params to a two-member 3-D dataset. Iters 
     1099accepts a string built from comma-separated array slices. For example, 
     1100iters = ":" will plot all iters in a single plot. Alternatively, 
     1101iters = ":2, 2:" will split the iters into two plots, while iters = "0" 
     1102will only plot the first iteration. 
     1103 
     1104The option "label" takes comma-separated strings. For example, label = "x,y," 
     1105will place 'x' on the x-axis, 'y' on the y-axis, and nothing on the z-axis. 
     1106LaTeX is also accepted. For example, label = "$ h $, $ {\\alpha}$, $ v$" will 
     1107label the axes with standard LaTeX math formatting. Note that the leading 
     1108space is required, while a trailing space aligns the text with the axis 
     1109instead of the plot frame. The option "filter" is used to select datapoints 
     1110from a given dataset, and takes comma-separated ints. A "mask" is given as 
     1111comma-separated ints; when the mask has more than one int, the plot will be 
     11122D. The option "vertical" will plot the dataset values on the vertical axis; 
     1113for 2D plots, cones are always plotted on the vertical axis. 
     1114 
     1115Required Inputs: 
     1116  filename            name of the python convergence logfile (e.g. paramlog.py) 
     1117 
     1118Additional Inputs: 
     1119  datafile            name of the dataset textfile (e.g. StAlDataset.txt) 
     1120""" 
     1121    import shlex 
     1122    global __quit 
     1123    __quit = False 
     1124 
     1125    # handle the special case where list is provided by sys.argv 
     1126    if isinstance(filename, (list,tuple)) and not kwds: 
     1127        cmdargs = filename # (above is used by script to parse command line) 
     1128    elif isinstance(filename, basestring) and not kwds: 
     1129        cmdargs = shlex.split(filename) 
     1130    # 'everything else' is essentially the functional interface 
     1131    else: 
     1132        out = kwds.get('out', None) 
     1133        bounds = kwds.get('bounds', None) 
     1134        iters = kwds.get('iters', None) 
     1135        label = kwds.get('label', None) 
     1136        dim = kwds.get('dim', None) 
     1137        filter = kwds.get('filter', None) 
     1138        mask = kwds.get('mask', None) 
     1139        nid = kwds.get('nid', None) 
     1140        scale = kwds.get('scale', None) 
     1141        gap = kwds.get('gap', None) 
     1142        data = kwds.get('data', False) 
     1143        cones = kwds.get('cones', False) 
     1144        vertical = kwds.get('vertical', False) 
     1145        flat = kwds.get('flat', False) 
     1146 
     1147        # process "commandline" arguments 
     1148        cmdargs = '' 
     1149        cmdargs += '' if out is None else '--out={} '.format(out) 
     1150        cmdargs += '' if bounds is None else '--bounds="{}" '.format(bounds) 
     1151        cmdargs += '' if iters is None else '--iters="{}" '.format(iters) 
     1152        cmdargs += '' if label is None else '--label="{}" '.format(label) 
     1153        cmdargs += '' if dim is None else '--dim="{}" '.format(dim) 
     1154        cmdargs += '' if filter is None else '--filter="{}" '.format(filter) 
     1155        cmdargs += '' if mask is None else '--mask={} '.format(mask) 
     1156        cmdargs += '' if nid is None else '--nid={} '.format(nid) 
     1157        cmdargs += '' if scale is None else '--scale={} '.format(scale) 
     1158        cmdargs += '' if gap is None else '--gap={} '.format(gap) 
     1159        cmdargs += '' if data == False else '--data ' 
     1160        cmdargs += '' if cones == False else '--cones ' 
     1161        cmdargs += '' if vertical == False else '--vertical ' 
     1162        cmdargs += '' if flat == False else '--flat ' 
     1163        cmdargs = filename.split() + shlex.split(cmdargs) 
     1164 
     1165    #XXX: note that 'argparse' is new as of python2.7 
     1166    from optparse import OptionParser 
     1167    def _exit(self, **kwds): 
     1168        global __quit 
     1169        __quit = True 
     1170    OptionParser.exit = _exit 
     1171 
     1172    parser = OptionParser(usage=hypercube_scenario.__doc__.split('\n\nOptions:')[0]) 
     1173    parser.add_option("-u","--out",action="store",dest="out",\ 
     1174                      metavar="STR",default=None, 
     1175                      help="filepath to save generated plot") 
     1176    parser.add_option("-b","--bounds",action="store",dest="bounds",\ 
     1177                      metavar="STR",default="0:1, 0:1, 0:1", 
     1178                      help="indicator string to set hypercube bounds") 
     1179    parser.add_option("-i","--iters",action="store",dest="iters",\ 
     1180                      metavar="STR",default="-1", 
     1181                      help="indicator string to select iterations to plot") 
     1182    parser.add_option("-l","--label",action="store",dest="label",\ 
     1183                      metavar="STR",default=",,", 
     1184                      help="string to assign label to axis") 
     1185    parser.add_option("-d","--dim",action="store",dest="dim",\ 
     1186                      metavar="STR",default="None", 
     1187                      help="indicator string set to scenario dimensions") 
     1188    parser.add_option("-p","--filter",action="store",dest="filter",\ 
     1189                      metavar="STR",default="None", 
     1190                      help="filter to select datapoints to plot") 
     1191    parser.add_option("-m","--mask",action="store",dest="replace",\ 
     1192                      metavar="INT",default=None, 
     1193                      help="id # of the coordinate axis not to be plotted") 
     1194    parser.add_option("-n","--nid",action="store",dest="id",\ 
     1195                      metavar="INT",default=None, 
     1196                      help="id # of the nth simultaneous points to plot") 
     1197    parser.add_option("-s","--scale",action="store",dest="scale",\ 
     1198                      metavar="INT",default=1.0, 
     1199                      help="grayscale contrast multiplier for points in plot") 
     1200    parser.add_option("-g","--gap",action="store",dest="gap",\ 
     1201                      metavar="INT",default=0.0, 
     1202                      help="distance from double cone center to cone vertex") 
     1203    parser.add_option("-o","--data",action="store_true",dest="legacy",\ 
     1204                      default=False,help="plot legacy data, if provided") 
     1205    parser.add_option("-v","--cones",action="store_true",dest="cones",\ 
     1206                      default=False,help="plot cones, if provided") 
     1207    parser.add_option("-z","--vertical",action="store_true",dest="vertical",\ 
     1208                      default=False,help="always plot values on vertical axis") 
     1209    parser.add_option("-f","--flat",action="store_true",dest="flatten",\ 
     1210                      default=False,help="show selected iterations in a single plot") 
     1211    parsed_opts, parsed_args = parser.parse_args(cmdargs) 
     1212 
     1213    from StringIO import StringIO 
     1214    f = StringIO() 
     1215    parser.print_help(file=f) 
     1216    f.seek(0) 
     1217    if 'Options:' not in hypercube_scenario.__doc__: 
     1218        hypercube_scenario.__doc__ += '\nOptions:%s' % f.read().split('Options:')[-1] 
     1219    f.close() 
     1220 
     1221    if __quit: return 
     1222 
     1223    # get the name of the parameter log file 
     1224    params, _cost = read_history(parsed_args[0]) 
     1225    # would be nice to use meta = ['wx','wx2','x','x2','wy',...] 
     1226    # exec "from %s import meta" % file 
     1227 
     1228    from mystic.math.discrete import scenario 
     1229    from mystic.math.legacydata import dataset 
     1230    try: # select whether to plot the cones 
     1231        cones = parsed_opts.cones 
     1232    except: 
     1233        cones = False 
     1234 
     1235    try: # select whether to plot the legacy data 
     1236        legacy = parsed_opts.legacy 
     1237    except: 
     1238        legacy = False 
     1239 
     1240    try: # get dataset filter 
     1241        filter = parsed_opts.filter 
     1242        if "None" == filter: filter = None 
     1243        else: filter = [int(i) for i in filter.split(",")] # format is "1,5,9" 
     1244    except: 
     1245        filter = None 
     1246 
     1247    try: # select the scenario dimensions 
     1248        npts = parsed_opts.dim 
     1249        if "None" == npts: # npts may have been logged 
     1250            from mystic.munge import read_import 
     1251            npts = read_import(parsed_args[0], 'npts') 
     1252        else: npts = tuple(int(i) for i in dim.split(",")) # format is "1,1,1" 
     1253    except: 
     1254        npts = (1,1,1) #XXX: better in parsed_args ? 
     1255 
     1256    try: # get the name of the dataset file 
     1257        file = parsed_args[1] 
     1258        from mystic.math.legacydata import load_dataset 
     1259        data = load_dataset(file, filter) 
     1260    except: 
     1261#       raise IOError, "please provide dataset file name" 
     1262        data = dataset() 
     1263        cones = False 
     1264        legacy = False 
     1265 
     1266    try: # select the bounds 
     1267        _bounds = parsed_opts.bounds 
     1268        if _bounds[-1] == "*": #XXX: then bounds are NOT strictly enforced 
     1269            _bounds = _bounds[:-1] 
     1270            strict = False 
     1271        else: 
     1272            strict = True 
     1273        bounds = _bounds.split(",")  # format is "60:105, 0:30, 2.1:2.8" 
     1274        bounds = [tuple(float(j) for j in i.split(':')) for i in bounds] 
     1275    except: 
     1276        strict = True 
     1277        bounds = [(0,1),(0,1),(0,1)] 
     1278 
     1279    try: # select labels for the axes 
     1280        label = parsed_opts.label.split(',')  # format is "x, y, z" 
     1281    except: 
     1282        label = ['','',''] 
     1283 
     1284    x = params[-1] 
     1285    try: # select which iterations to plot 
     1286        select = parsed_opts.iters.split(',')  # format is ":2, 2:4, 5, 6:" 
     1287    except: 
     1288        select = ['-1'] 
     1289       #select = [':'] 
     1290       #select = [':1'] 
     1291       #select = [':2','2:'] 
     1292       #select = [':1','1:2','2:3','3:'] 
     1293       #select = ['0','1','2','3'] 
     1294 
     1295    try: # collapse non-consecutive iterations into a single plot... 
     1296        flatten = parsed_opts.flatten 
     1297    except: 
     1298        flatten = False 
     1299 
     1300    try: # select which 'id' to plot results for 
     1301        id = int(parsed_opts.id) 
     1302    except: 
     1303        id = None # i.e. 'all' **or** use id=0, which should be 'best' energy ? 
     1304 
     1305    try: # scale the color in plotting the weights 
     1306        scale = float(parsed_opts.scale) 
     1307    except: 
     1308        scale = 1.0 # color = color**scale 
     1309 
     1310    try: # gap between double cone center and cone vertex 
     1311        gap = float(parsed_opts.gap) 
     1312    except: 
     1313        gap = 0.0 
     1314 
     1315    _2D = False # if False, use 3D plots; if True, use 3D plots 
     1316    cs = None 
     1317    try: # select which axis to plot 'values' on  (3D plot) 
     1318        xs = int(parsed_opts.replace) 
     1319    except: 
     1320        try: # select which axes to mask (2D plot) 
     1321            xs = (int(i) for i in parsed_opts.replace.split(",")) # format is "1,2" 
     1322            xs = list(reversed(sorted(set(xs)))) 
     1323            cs = int(xs[-1]) if xs[-1] != xs[0] else None 
     1324            xs = int(xs[0]) 
     1325            xs,cs = cs,xs # cs will swap coord axes; xs will swap with value axis 
     1326            _2D = True #NOTE: always apply cs swap before xs swap 
     1327        except: 
     1328            xs = None # don't plot values; plot values on 'x' axis with xs = 0 
     1329 
     1330    try: # always plot cones on vertical axis 
     1331        vertical_cones = parsed_opts.vertical 
     1332    except: 
     1333        vertical_cones = False 
     1334    if _2D: # always plot cones on vertical axis 
     1335        vertical_cones = True 
     1336 
     1337    # ensure all terms of bounds are tuples 
     1338    for bound in bounds: 
     1339        if not isinstance(bound, tuple): 
     1340            raise TypeError, "bounds should be tuples of (lower_bound,upper_bound)" 
     1341 
     1342    # ensure all terms of select are strings that have a ":" 
     1343    for i in range(len(select)): 
     1344        if isinstance(select[i], int): select[i] = str(select[i]) 
     1345        if select[i] == '-1': select[i] = 'len(x)-1:len(x)' 
     1346        elif not select[i].count(':'): 
     1347            select[i] += ':' + str(int(select[i])+1) 
     1348 
     1349    # take only the selected 'id' 
     1350    if id != None: 
     1351        param = [] 
     1352        for j in range(len(params)): 
     1353            param.append([p[id] for p in params[j]]) 
     1354        params = param[:] 
     1355 
     1356    # at this point, we should have: 
     1357    #bounds = [(60,105),(0,30),(2.1,2.8)] or [(None,None)]*3 
     1358    #select = ['-1:'] or [':'] or [':1','1:2','2:3','3:'] or similar 
     1359    #id = 0 or None 
     1360 
     1361    # get dataset coords (and values) for selected axes 
     1362    if data: 
     1363        slope = _get_slope(data, xs, cs) 
     1364        coords = _get_coords(data, xs, cs) 
     1365        #print "bounds: %s" % bounds 
     1366        #print "slope: %s" % slope 
     1367        #print "coords: %s" % coords 
     1368   #else: 
     1369   #    slope = [] 
     1370   #    coords = [] 
     1371 
     1372    plots = len(select) 
     1373    if not flatten: 
     1374        dim1,dim2 = best_dimensions(plots) 
     1375    else: dim1,dim2 = 1,1 
     1376 
     1377    # use the default bounds where not specified 
     1378    bounds = [list(i) for i in bounds] 
     1379    for i in range(len(bounds)): 
     1380        if bounds[i][0] is None: bounds[i][0] = 0 
     1381        if bounds[i][1] is None: bounds[i][1] = 1 
     1382 
     1383    # swap the axes appropriately when plotting cones (when replacing an axis) 
     1384    if _2D and xs == 0: 
     1385        if data: 
     1386            slope[0],slope[1] = slope[1],slope[0] 
     1387            coords = [list(reversed(pt[:2]))+pt[2:] for pt in coords] 
     1388       #bounds[0],bounds[1] = bounds[1],bounds[0] 
     1389       #label[0],label[1] = label[1],label[0] 
     1390        axis = xs #None 
     1391    elif not _2D and vertical_cones and xs in range(len(bounds)): 
     1392        # adjust slope, bounds, and data so cone axis is last  
     1393        if data: 
     1394            slope = swap(slope, xs)  
     1395            coords = [swap(pt,xs) for pt in coords] 
     1396        bounds = swap(bounds, xs)  
     1397        label = swap(label, xs) 
     1398        axis = None 
     1399    else: 
     1400        axis = xs 
     1401 
     1402    # correctly bound the first plot.  there must be at least one plot 
     1403    fig = plt.figure()  
     1404    if _2D: 
     1405        ax1 = fig.add_subplot(dim1,dim2,1) 
     1406        ax1.plot([bounds[0][0]],[bounds[1][0]]) 
     1407        ax1.plot([bounds[0][1]],[bounds[1][1]]) 
     1408    else: 
     1409        ax1 = Subplot3D(fig, dim1,dim2,1) 
     1410        ax1.plot([bounds[0][0]],[bounds[1][0]],[bounds[2][0]]) 
     1411        ax1.plot([bounds[0][1]],[bounds[1][1]],[bounds[2][1]]) 
     1412    globals = {'plt':plt,'fig':fig,'dim1':dim1,'dim2':dim2,\ 
     1413               'Subplot3D':Subplot3D,'bounds':bounds,'label':label} 
     1414    if not flatten: 
     1415        code = "plt.title('iterations[%s]');" % select[0] 
     1416    else:  
     1417        code = "plt.title('iterations[*]');" 
     1418    code = compile(code, '<string>', 'exec') 
     1419    exec code in globals 
     1420    if cones and data and xs in range(len(bounds)): 
     1421        if _2D: 
     1422            _plot_bowtie(ax1,coords,slope,bounds,axis=axis,tol=gap) 
     1423        else: 
     1424            _plot_cones(ax1,coords,slope,bounds,axis=axis,strict=strict,tol=gap) 
     1425    if legacy and data: 
     1426        _plot_data(ax1,coords,bounds,strict=strict,_2D=_2D) 
     1427   #_clip_axes(ax1,bounds) 
     1428    _label_axes(ax1,label,_2D=_2D) 
     1429    a = [ax1] 
     1430 
     1431    # set up additional plots 
     1432    if not flatten: 
     1433        for i in range(2, plots + 1): 
     1434            if _2D: 
     1435                code = "ax%d = ax = fig.add_subplot(dim1,dim2,%d);" % (i,i) 
     1436                code += "ax.plot([bounds[0][0]],[bounds[1][0]]);" 
     1437                code += "ax.plot([bounds[0][1]],[bounds[1][1]]);" 
     1438            else: 
     1439                code = "ax%d = ax = Subplot3D(fig, dim1,dim2,%d);" % (i,i) 
     1440                code += "ax.plot([bounds[0][0]],[bounds[1][0]],[bounds[2][0]]);" 
     1441                code += "ax.plot([bounds[0][1]],[bounds[1][1]],[bounds[2][1]]);" 
     1442            code += "plt.title('iterations[%s]');" % select[i - 1] 
     1443            code = compile(code, '<string>', 'exec') 
     1444            exec code in globals 
     1445            ax = globals['ax'] 
     1446            if cones and data and xs in range(len(bounds)): 
     1447                if _2D: 
     1448                    _plot_bowtie(ax,coords,slope,bounds,axis=axis,tol=gap) 
     1449                else: 
     1450                    _plot_cones(ax,coords,slope,bounds,axis=axis,strict=strict,tol=gap) 
     1451            if legacy and data: 
     1452                _plot_data(ax,coords,bounds,strict=strict,_2D=_2D) 
     1453           #_clip_axes(ax,bounds) 
     1454            _label_axes(ax,label,_2D=_2D) 
     1455            a.append(ax) 
     1456 
     1457    # turn each "n:m" in select to a list 
     1458    _select = [] 
     1459    for sel in select: 
     1460        if sel[0] == ':': _select.append("0"+sel) 
     1461        else: _select.append(sel) 
     1462    for i in range(len(_select)): 
     1463        if _select[i][-1] == ':': select[i] = _select[i]+str(len(x)) 
     1464        else: select[i] = _select[i] 
     1465    for i in range(len(select)): 
     1466        p = select[i].split(":") 
     1467        if p[0][0] == '-': p[0] = "len(x)"+p[0] 
     1468        if p[1][0] == '-': p[1] = "len(x)"+p[1] 
     1469        select[i] = p[0]+":"+p[1] 
     1470    steps = [eval("range(%s)" % sel.replace(":",",")) for sel in select] 
     1471 
     1472    # at this point, we should have: 
     1473    #steps = [[0,1],[1,2],[2,3],[3,4,5,6,7,8]] or similar 
     1474    if flatten: 
     1475        from mystic.tools import flatten 
     1476        steps = [list(flatten(steps))] 
     1477 
     1478    # plot all the scenario "data" 
     1479    from numpy import inf, e 
     1480    scale = e**(scale - 1.0) 
     1481    for v in range(len(steps)): 
     1482        if len(steps[v]) > 1: qp = float(max(steps[v])) 
     1483        else: qp = inf 
     1484        for s in steps[v]: 
     1485            par = eval("[params[q][%s][0] for q in range(len(params))]" % s) 
     1486            pm = scenario() 
     1487            pm.load(par, npts) 
     1488            d = dataset() 
     1489            d.load(pm.coords, pm.values) 
     1490            # dot color determined by number of simultaneous iterations 
     1491            t = str((s/qp)**scale) 
     1492            # get and plot dataset coords for selected axes       
     1493            _coords = _get_coords(d, xs, cs) 
     1494            # check if we are replacing an axis 
     1495            if _2D and xs == 0: 
     1496                if data: # adjust data so cone axis is last  
     1497                    _coords = [list(reversed(pt[:2]))+pt[2:] for pt in _coords] 
     1498            elif not _2D and vertical_cones and xs in range(len(bounds)): 
     1499                if data: # adjust data so cone axis is last  
     1500                    _coords = [swap(pt,xs) for pt in _coords] 
     1501            _plot_data(a[v], _coords, bounds, color=t, strict=strict, _2D=_2D) 
     1502 
     1503    if _2D: # strict = True 
     1504        plt.xlim((bounds[0][0],bounds[0][1])) 
     1505        plt.ylim((bounds[1][0],bounds[1][1])) 
     1506 
     1507    if not parsed_opts.out: 
     1508        plt.show() 
     1509    else: 
     1510        fig.savefig(parsed_opts.out) 
     1511 
     1512 
    8411513if __name__ == '__main__': 
    8421514    pass 
  • mystic/scripts/support_hypercube_scenario.py

    r801 r809  
    66#  - http://mmckerns.github.io/project/mystic/browser/mystic/LICENSE 
    77 
    8 __doc__ = """ 
    9 support_hypercube_scenario.py [options] filename (datafile) 
     8from mystic.support import hypercube_scenario 
    109 
    11 generate scenario support plots from file written with 'write_support_file'; 
    12 generate legacy data and cones from a dataset file, if provided 
     10__doc__ = hypercube_scenario.__doc__ 
    1311 
    14 The options "bounds", "dim", and "iters" all take indicator strings. 
    15 The bounds should be given as comma-separated slices. For example, using 
    16 bounds = ".062:.125, 0:30, 2300:3200" will set the lower and upper bounds 
    17 for x to be (.062,.125), y to be (0,30), and z to be (2300,3200). If all 
    18 bounds are to not be strictly enforced, append an asterisk '*' to the string.  
    19 The dim (dimensions of the scenario) should comma-separated ints. For example, 
    20 dim = "1, 1, 2" will convert the params to a two-member 3-D dataset. Iters 
    21 accepts a string built from comma-separated array slices. For example, 
    22 iters = ":" will plot all iters in a single plot. Alternatively, 
    23 iters = ":2, 2:" will split the iters into two plots, while iters = "0" 
    24 will only plot the first iteration. 
     12if __name__ == '__main__': 
    2513 
    26 The option "label" takes comma-separated strings. For example, label = "x,y," 
    27 will place 'x' on the x-axis, 'y' on the y-axis, and nothing on the z-axis. 
    28 LaTeX is also accepted. For example, label = "$ h $, $ {\\alpha}$, $ v$" will 
    29 label the axes with standard LaTeX math formatting. Note that the leading 
    30 space is required, while a trailing space aligns the text with the axis 
    31 instead of the plot frame. The option "filter" is used to select datapoints 
    32 from a given dataset, and takes comma-separated ints. A "mask" is given as 
    33 comma-separated ints; when the mask has more than one int, the plot will be 
    34 2D. The option "vertical" will plot the dataset values on the vertical axis; 
    35 for 2D plots, cones are always plotted on the vertical axis. 
     14    import sys 
    3615 
    37 Required Inputs: 
    38   filename            name of the python convergence logfile (e.g. paramlog.py) 
    39  
    40 Additional Inputs: 
    41   datafile            name of the dataset textfile (e.g. StAlDataset.txt) 
    42 """ 
    43 ZERO = 1.0e-6  # zero-ish 
    44  
    45 #### building the cone primitive #### 
    46 def cone_builder(slope, bounds, strict=True): 
    47   """ factory to create a cone primitive 
    48  
    49   slope -- slope multiplier for cone on the X,Y,Z axes (for mesh construction) 
    50   bounds -- list of tuples of bounds for the plot; (lower,upper) for each axis 
    51 """ 
    52   from mystic.math import almostEqual 
    53   import numpy as np 
    54   def cone_mesh(length): 
    55     """ construct a conical mesh for a given length of cone """ 
    56     L1,L2,L3 = slope 
    57     radius = length / L3 #XXX: * 0.5 
    58     r0 = ZERO 
    59  
    60     if almostEqual(radius, r0, tol=r0): radius = r0 
    61     r = np.linspace(radius,radius,6)  
    62     r[0]= np.zeros(r[0].shape)  
    63     r[1] *= r0/radius 
    64     r[5] *= r0/radius 
    65     r[3]= np.zeros(r[3].shape)  
    66  
    67     p = np.linspace(0,2*np.pi,50)  
    68     R,P = np.meshgrid(r,p)  
    69     X,Y = L1 * R*np.cos(P), L2 * R*np.sin(P)  
    70  
    71     tmp=list()  
    72     for i in range(np.size(p)):  
    73       tmp.append([0,0,length,length,length,0]) # len = size(r) 
    74     Z = np.array(tmp)  
    75     return X,Z,Y 
    76  
    77   lb,ub = zip(*bounds) 
    78   # if False, the cone surface may violate bounds 
    79  #strict = True # always respect the bounds 
    80  
    81   def cone(position, top=True): 
    82     """ construct a cone primitive, with vertex at given position 
    83  
    84     position -- tuple of vertex position 
    85     top -- boolean for if 'top' or 'bottom' cone 
    86 """ 
    87     from numpy import asarray 
    88     position = asarray(position) 
    89     d_hi = ub - position 
    90     d_lo = lb - position 
    91     if strict and (any(d_hi < 0) or any(d_lo > 0)): 
    92       return None  # don't plot cone if vertex is out of bounds 
    93  
    94     if top: # distance of vertex from upper edge 
    95       l = d_hi[2] 
    96       if not l: l = ZERO 
    97     else:   # distance of vertex from lower edge 
    98       l = d_lo[2] 
    99       if not l: l = -ZERO 
    100     X,Z,Y = cone_mesh(l) 
    101     X = X + position[0] 
    102     Y = Y + position[1] 
    103     Z = Z + position[2] 
    104     return X,Z,Y 
    105  
    106   return cone     
     16    hypercube_scenario(sys.argv[1:]) 
    10717 
    10818 
    109 #### plotting the cones #### 
    110 def plot_bowtie(ax, data, slope, bounds, color='0.75', axis=None, tol=0.0): 
    111   """plot (2D) double cones for a given dataset 
    112  
    113   ax -- matplotlib 'Axes3D' plot object 
    114   data -- list of datapoints, where datapoints are 3-tuples (i.e. x,y,z) 
    115   slope -- slope multiplier for cone on the X,Y,Z axes (for mesh construction) 
    116   bounds -- list of tuples of bounds for the plot; (lower,upper) for each axis 
    117   color -- string name (or rbg value) of color to use for datapoints 
    118   axis -- the axis of the cone 
    119   tol -- distance between center of mass of the double cones and a cone vertex 
    120 """ 
    121   if axis not in range(len(bounds)-1): return ax 
    122   from numpy import asarray, inf 
    123   data = asarray(data) 
    124   sl = slope[axis] 
    125   gap = max(0., tol) 
    126   # find the endpoints of the cone:  y = slope * (x0 - x) + y0 
    127   ylo = sl * (bounds[0][0] - data.T[0]) + data.T[1] 
    128   yhi = sl * (bounds[0][1] - data.T[0]) + data.T[1] 
    129   zhi = -sl * (bounds[0][0] - data.T[0]) + data.T[1] 
    130   zlo = -sl * (bounds[0][1] - data.T[0]) + data.T[1] 
    131   xlb = bounds[0][0] 
    132   xub = bounds[0][1] 
    133   ylb = bounds[1][0] 
    134   yub = bounds[1][1] 
    135  
    136   for pt,yl,yu,zl,zu in zip(data,ylo,yhi,zlo,zhi): 
    137    #ax.plot([xlb, pt[0], xub], [yl, pt[1], yu], color='black') 
    138    #ax.plot([xlb, pt[0], xub], [zu, pt[1], zl], color='black') 
    139     ax.fill_between([xlb, pt[0], xub], [ylb]*3, [yl-gap, pt[1]-gap, zl-gap], \ 
    140                                      facecolor=color, alpha=0.2) 
    141     ax.fill_between([xlb, pt[0], xub], [yub]*3, [zu+gap, pt[1]+gap, yu+gap], \ 
    142                                      facecolor=color, alpha=0.2) 
    143   return ax 
    144  
    145 def plot_cones(ax, data, slope, bounds, color='0.75', axis=None, strict=True, tol=0.0): 
    146   """plot (3D) double cones for a given dataset 
    147  
    148   ax -- matplotlib 'Axes3D' plot object 
    149   data -- list of datapoints, where datapoints are 3-tuples (i.e. x,y,z) 
    150   slope -- slope multiplier for cone on the X,Y,Z axes (for mesh construction) 
    151   bounds -- list of tuples of bounds for the plot; (lower,upper) for each axis 
    152   color -- string name (or rbg value) of color to use for datapoints 
    153   axis -- the axis of the cone 
    154   tol -- distance between center of mass of the double cones and a cone vertex 
    155 """ 
    156   from numpy import asarray, zeros 
    157   # adjust slope, bounds, and data so cone axis is last  
    158   slope = swap(slope, axis)  
    159   bounds = swap(bounds, axis)  
    160   data = [swap(pt,axis) for pt in data] 
    161   cone = cone_builder(slope, bounds, strict=strict) 
    162   # prepare to apply the 'gap' tolerance 
    163   gap = zeros(len(slope)) 
    164   gap[-1:] = [max(0., tol)] #XXX: bad behavior if len(slope) is 0 
    165   # plot the upper and lower cone 
    166   for datapt in data: 
    167     datapt = asarray(datapt) 
    168     _cone = cone(datapt + gap, top=True) 
    169     if _cone: 
    170       X,Z,Y = swap(_cone, axis) # 'unswap' the axes 
    171       ax.plot_surface(X, Y,Z, rstride=1, cstride=1, color=color, \ 
    172                                          linewidths=0, alpha=.1) 
    173     _cone = cone(datapt - gap, top=False) 
    174     if _cone: 
    175       X,Z,Y = swap(_cone, axis) # 'unswap' the axes 
    176       ax.plot_surface(X, Y,Z, rstride=1, cstride=1, color=color, \ 
    177                                          linewidths=0, alpha=.1) 
    178   return ax 
    179  
    180 def plot_data(ax, data, bounds, color='red', strict=True): 
    181   """plot datapoints for a given dataset 
    182  
    183   ax -- matplotlib 'Axes3D' plot object 
    184   data -- list of datapoints, where datapoints are 3-tuples (i.e. x,y,z) 
    185   bounds -- list of tuples of bounds for the plot; (lower,upper) for each axis 
    186   color -- string name (or rbg value) of color to use for datapoints 
    187 """ 
    188 # strict = True # always respect the bounds 
    189   lb,ub = zip(*bounds) 
    190   # plot the datapoints themselves 
    191   from numpy import asarray 
    192   for datapt in data: 
    193     position = asarray(datapt) 
    194     d_hi = ub - position 
    195     d_lo = lb - position 
    196     if strict and (any(d_hi < 0) or any(d_lo > 0)): #XXX: any or just in Y ? 
    197       pass  # don't plot if datapt is out of bounds 
    198     else: 
    199       if _2D: 
    200         ax.plot([datapt[0]], [datapt[1]], \ 
    201                 color=color, marker='o', markersize=10) 
    202       else: 
    203         ax.plot([datapt[0]], [datapt[1]], [datapt[2]], \ 
    204                 color=color, marker='o', markersize=10) 
    205   return ax 
    206  
    207 def clip_axes(ax, bounds): 
    208   """ clip plots to be within given lower and upper bounds 
    209  
    210   ax -- matplotlib 'Axes3D' plot object 
    211   bounds -- list of tuples of bounds for the plot; (lower,upper) for each axis 
    212 """ 
    213   lb,ub = zip(*bounds) 
    214   # plot only within [lb,ub] 
    215   ax.set_xlim3d(lb[0], ub[0]) 
    216   ax.set_ylim3d(lb[1], ub[1]) 
    217   ax.set_zlim3d(lb[2], ub[2]) # cone "center" axis 
    218   return ax 
    219  
    220 def label_axes(ax, labels): 
    221   """ label plots with given string labels 
    222  
    223   ax -- matplotlib 'Axes3D' plot object 
    224   labels -- list of string labels for the plot 
    225 """ 
    226   ax.set_xlabel(labels[0]) 
    227   ax.set_ylabel(labels[1]) 
    228   if not _2D: 
    229     ax.set_zlabel(labels[2]) # cone "center" axis 
    230   return ax 
    231  
    232 def get_slope(data, replace=None, mask=None): 
    233   """ replace one slope in a list of slopes with '1.0' 
    234   (i.e. replace a slope from the dataset with the unit slope) 
    235  
    236   data -- dataset object, where coordinates are 3-tuples and values are floats 
    237   replace -- selected axis (an int) to plot values NOT coords 
    238   """ 
    239   slope = data.lipschitz 
    240   if mask in range(len(slope)): 
    241     slope = swap(slope, mask) 
    242   if replace not in range(len(slope)):  # don't replace an axis 
    243     return slope 
    244   return slope[:replace] + [1.0] + slope[replace+1:] 
    245  
    246 def get_coords(data, replace=None, mask=None): 
    247   """ replace one coordiate axis in a 3-D data set with 'data values' 
    248   (i.e. replace an 'x' axis with the 'y' values of the data) 
    249  
    250   data -- dataset object, where coordinates are 3-tuples and values are floats 
    251   replace -- selected axis (an int) to plot values NOT coords 
    252   """ 
    253   slope = data.lipschitz 
    254   coords = data.coords 
    255   values = data.values 
    256   if mask in range(len(slope)): 
    257     coords = [swap(pt,mask) for pt in coords] 
    258   if replace not in range(len(slope)):  # don't replace an axis 
    259     return coords 
    260   return [list(coords[i][:replace]) + [values[i]] + \ 
    261           list(coords[i][replace+1:]) for i in range(len(coords))] 
    262  
    263 def swap(alist, index=None): 
    264   """ swap the selected list element with the last element in alist 
    265  
    266   alist -- a list of objects 
    267   index -- the selected element 
    268   """ 
    269   if index not in range(len(alist)):  # don't swap an element 
    270     return alist  
    271   return alist[:index] + alist[index+1:] + alist[index:index+1] 
    272  
    273 from support_convergence import best_dimensions 
    274  
    275  
    276 if __name__ == '__main__': 
    277   #XXX: note that 'argparse' is new as of python2.7 
    278   from optparse import OptionParser 
    279   parser = OptionParser(usage=__doc__) 
    280   parser.add_option("-b","--bounds",action="store",dest="bounds",\ 
    281                     metavar="STR",default="0:1, 0:1, 0:1", 
    282                     help="indicator string to set hypercube bounds") 
    283   parser.add_option("-i","--iters",action="store",dest="iters",\ 
    284                     metavar="STR",default="-1", 
    285                     help="indicator string to select iterations to plot") 
    286   parser.add_option("-l","--label",action="store",dest="label",\ 
    287                     metavar="STR",default=",,", 
    288                     help="string to assign label to axis") 
    289   parser.add_option("-d","--dim",action="store",dest="dim",\ 
    290                     metavar="STR",default="None", 
    291                     help="indicator string set to scenario dimensions") 
    292   parser.add_option("-p","--filter",action="store",dest="filter",\ 
    293                     metavar="STR",default="None", 
    294                     help="filter to select datapoints to plot") 
    295   parser.add_option("-m","--mask",action="store",dest="replace",\ 
    296                     metavar="INT",default=None, 
    297                     help="id # of the coordinate axis not to be plotted") 
    298   parser.add_option("-n","--nid",action="store",dest="id",\ 
    299                     metavar="INT",default=None, 
    300                     help="id # of the nth simultaneous points to plot") 
    301   parser.add_option("-s","--scale",action="store",dest="scale",\ 
    302                     metavar="INT",default=1.0, 
    303                     help="grayscale contrast multiplier for points in plot") 
    304   parser.add_option("-g","--gap",action="store",dest="gap",\ 
    305                     metavar="INT",default=0.0, 
    306                     help="distance from double cone center to cone vertex") 
    307   parser.add_option("-o","--data",action="store_true",dest="legacy",\ 
    308                     default=False,help="plot legacy data, if provided") 
    309   parser.add_option("-v","--cones",action="store_true",dest="cones",\ 
    310                     default=False,help="plot cones, if provided") 
    311   parser.add_option("-z","--vertical",action="store_true",dest="vertical",\ 
    312                     default=False,help="always plot values on vertical axis") 
    313   parser.add_option("-f","--flat",action="store_true",dest="flatten",\ 
    314                     default=False,help="show selected iterations in a single plot") 
    315   parsed_opts, parsed_args = parser.parse_args() 
    316  
    317   # get the name of the parameter log file 
    318   from mystic.munge import read_history 
    319   params, _cost = read_history(parsed_args[0]) 
    320   # would be nice to use meta = ['wx','wx2','x','x2','wy',...] 
    321   # exec "from %s import meta" % file 
    322  
    323   from mystic.math.discrete import scenario 
    324   from mystic.math.legacydata import dataset 
    325   try: # select whether to plot the cones 
    326     cones = parsed_opts.cones 
    327   except: 
    328     cones = False 
    329  
    330   try: # select whether to plot the legacy data 
    331     legacy = parsed_opts.legacy 
    332   except: 
    333     legacy = False 
    334  
    335   try: # get dataset filter 
    336     filter = parsed_opts.filter 
    337     if "None" == filter: filter = None 
    338     else: filter = [int(i) for i in filter.split(",")] # format is "1,5,9" 
    339   except: 
    340     filter = None 
    341  
    342   try: # select the scenario dimensions 
    343     npts = parsed_opts.dim 
    344     if "None" == npts: # npts may have been logged 
    345       from mystic.munge import read_import 
    346       npts = read_import(parsed_args[0], 'npts') 
    347     else: npts = tuple(int(i) for i in dim.split(",")) # format is "1,1,1" 
    348   except: 
    349     npts = (1,1,1) #XXX: better in parsed_args ? 
    350  
    351   try: # get the name of the dataset file 
    352     file = parsed_args[1] 
    353     from mystic.math.legacydata import load_dataset 
    354     data = load_dataset(file, filter) 
    355   except: 
    356 #   raise IOError, "please provide dataset file name" 
    357     data = dataset() 
    358     cones = False 
    359     legacy = False 
    360  
    361   try: # select the bounds 
    362     _bounds = parsed_opts.bounds 
    363     if _bounds[-1] == "*": #XXX: then bounds are NOT strictly enforced 
    364       _bounds = _bounds[:-1] 
    365       strict = False 
    366     else: 
    367       strict = True 
    368     bounds = _bounds.split(",")  # format is "60:105, 0:30, 2.1:2.8" 
    369     bounds = [tuple(float(j) for j in i.split(':')) for i in bounds] 
    370   except: 
    371     strict = True 
    372     bounds = [(0,1),(0,1),(0,1)] 
    373  
    374   try: # select labels for the axes 
    375     label = parsed_opts.label.split(',')  # format is "x, y, z" 
    376   except: 
    377     label = ['','',''] 
    378  
    379   x = params[-1] 
    380   try: # select which iterations to plot 
    381     select = parsed_opts.iters.split(',')  # format is ":2, 2:4, 5, 6:" 
    382   except: 
    383     select = ['-1'] 
    384    #select = [':'] 
    385    #select = [':1'] 
    386    #select = [':2','2:'] 
    387    #select = [':1','1:2','2:3','3:'] 
    388    #select = ['0','1','2','3'] 
    389  
    390   try: # collapse non-consecutive iterations into a single plot... 
    391     flatten = parsed_opts.flatten 
    392   except: 
    393     flatten = False 
    394  
    395   try: # select which 'id' to plot results for 
    396     id = int(parsed_opts.id) 
    397   except: 
    398     id = None # i.e. 'all' **or** use id=0, which should be 'best' energy ? 
    399  
    400   try: # scale the color in plotting the weights 
    401     scale = float(parsed_opts.scale) 
    402   except: 
    403     scale = 1.0 # color = color**scale 
    404  
    405   try: # gap between double cone center and cone vertex 
    406     gap = float(parsed_opts.gap) 
    407   except: 
    408     gap = 0.0 
    409  
    410   _2D = False # if False, use 3D plots; if True, use 3D plots 
    411   _2D = False # if False, use 3D plots; if True, use 3D plots 
    412   cs = None 
    413   try: # select which axis to plot 'values' on  (3D plot) 
    414     xs = int(parsed_opts.replace) 
    415   except: 
    416     try: # select which axes to mask (2D plot) 
    417       xs = (int(i) for i in parsed_opts.replace.split(",")) # format is "1,2" 
    418       xs = list(reversed(sorted(set(xs)))) 
    419       cs = int(xs[-1]) if xs[-1] != xs[0] else None 
    420       xs = int(xs[0]) 
    421       xs,cs = cs,xs # cs will swap coord axes; xs will swap with value axis 
    422       _2D = True #NOTE: always apply cs swap before xs swap 
    423     except: 
    424       xs = None # don't plot values; plot values on 'x' axis with xs = 0 
    425  
    426   try: # always plot cones on vertical axis 
    427     vertical_cones = parsed_opts.vertical 
    428   except: 
    429     vertical_cones = False 
    430   if _2D: # always plot cones on vertical axis 
    431     vertical_cones = True 
    432  
    433   # ensure all terms of bounds are tuples 
    434   for bound in bounds: 
    435     if not isinstance(bound, tuple): 
    436       raise TypeError, "bounds should be tuples of (lower_bound,upper_bound)" 
    437  
    438   # ensure all terms of select are strings that have a ":" 
    439   for i in range(len(select)): 
    440     if isinstance(select[i], int): select[i] = str(select[i]) 
    441     if select[i] == '-1': select[i] = 'len(x)-1:len(x)' 
    442     elif not select[i].count(':'): 
    443       select[i] += ':' + str(int(select[i])+1) 
    444  
    445   # take only the selected 'id' 
    446   if id != None: 
    447     param = [] 
    448     for j in range(len(params)): 
    449       param.append([p[id] for p in params[j]]) 
    450     params = param[:] 
    451  
    452   # at this point, we should have: 
    453   #bounds = [(60,105),(0,30),(2.1,2.8)] or [(None,None),(None,None),(None,None)] 
    454   #select = ['-1:'] or [':'] or [':1','1:2','2:3','3:'] or similar 
    455   #id = 0 or None 
    456  
    457   # get dataset coords (and values) for selected axes 
    458   if data: 
    459     slope = get_slope(data, xs, cs) 
    460     coords = get_coords(data, xs, cs) 
    461     #print "bounds: %s" % bounds 
    462     #print "slope: %s" % slope 
    463     #print "coords: %s" % coords 
    464  #else: 
    465  #  slope = [] 
    466  #  coords = [] 
    467  
    468   import matplotlib.pyplot as plt 
    469   if not _2D: 
    470     from mpl_toolkits.mplot3d import Axes3D 
    471     from matplotlib.axes import subplot_class_factory 
    472     Subplot3D = subplot_class_factory(Axes3D) 
    473  
    474   plots = len(select) 
    475   if not flatten: 
    476     dim1,dim2 = best_dimensions(plots) 
    477   else: dim1,dim2 = 1,1 
    478  
    479   # use the default bounds where not specified 
    480   bounds = [list(i) for i in bounds] 
    481   for i in range(len(bounds)): 
    482     if bounds[i][0] is None: bounds[i][0] = 0 
    483     if bounds[i][1] is None: bounds[i][1] = 1 
    484  
    485   # swap the axes appropriately when plotting cones (when replacing an axis) 
    486   if _2D and xs == 0: 
    487     if data: 
    488       slope[0],slope[1] = slope[1],slope[0] 
    489       coords = [list(reversed(pt[:2]))+pt[2:] for pt in coords] 
    490    #bounds[0],bounds[1] = bounds[1],bounds[0] 
    491    #label[0],label[1] = label[1],label[0] 
    492     axis = xs #None 
    493   elif not _2D and vertical_cones and xs in range(len(bounds)): 
    494     # adjust slope, bounds, and data so cone axis is last  
    495     if data: 
    496       slope = swap(slope, xs)  
    497       coords = [swap(pt,xs) for pt in coords] 
    498     bounds = swap(bounds, xs)  
    499     label = swap(label, xs) 
    500     axis = None 
    501   else: 
    502     axis = xs 
    503  
    504   # correctly bound the first plot.  there must be at least one plot 
    505   fig = plt.figure()  
    506   if _2D: 
    507     ax1 = fig.add_subplot(dim1,dim2,1) 
    508     ax1.plot([bounds[0][0]],[bounds[1][0]]) 
    509     ax1.plot([bounds[0][1]],[bounds[1][1]]) 
    510   else: 
    511     ax1 = Subplot3D(fig, dim1,dim2,1) 
    512     ax1.plot([bounds[0][0]],[bounds[1][0]],[bounds[2][0]]) 
    513     ax1.plot([bounds[0][1]],[bounds[1][1]],[bounds[2][1]]) 
    514   if not flatten: 
    515     exec "plt.title('iterations[%s]')" % select[0] 
    516   else: 
    517     exec "plt.title('iterations[*]')" 
    518   if cones and data and xs in range(len(bounds)): 
    519     if _2D: 
    520       plot_bowtie(ax1,coords,slope,bounds,axis=axis,tol=gap) 
    521     else: 
    522       plot_cones(ax1,coords,slope,bounds,axis=axis,strict=strict,tol=gap) 
    523   if legacy and data: 
    524     plot_data(ax1,coords,bounds,strict=strict) 
    525  #clip_axes(ax1,bounds) 
    526   label_axes(ax1,label) 
    527   a = [ax1] 
    528  
    529   # set up additional plots 
    530   if not flatten: 
    531     for i in range(2, plots + 1): 
    532       if _2D: 
    533         exec "ax%d = fig.add_subplot(dim1,dim2,%d)" % (i,i) 
    534         exec "ax%d.plot([bounds[0][0]],[bounds[1][0]])" % i 
    535         exec "ax%d.plot([bounds[0][1]],[bounds[1][1]])" % i 
    536       else: 
    537         exec "ax%d = Subplot3D(fig, dim1,dim2,%d)" % (i,i) 
    538         exec "ax%d.plot([bounds[0][0]],[bounds[1][0]],[bounds[2][0]])" % i 
    539         exec "ax%d.plot([bounds[0][1]],[bounds[1][1]],[bounds[2][1]])" % i 
    540       exec "plt.title('iterations[%s]')" % select[i - 1] 
    541       if cones and data and xs in range(len(bounds)): 
    542         if _2D: 
    543           exec "plot_bowtie(ax%d,coords,slope,bounds,axis=axis,tol=gap)" % i 
    544         else: 
    545           exec "plot_cones(ax%d,coords,slope,bounds,axis=axis,strict=strict,tol=gap)" % i 
    546       if legacy and data: 
    547         exec "plot_data(ax%d,coords,bounds,strict=strict)" % i 
    548      #exec "clip_axes(ax%d,bounds)" % i 
    549       exec "label_axes(ax%d,label)" % i 
    550       exec "a.append(ax%d)" % i 
    551  
    552   # turn each "n:m" in select to a list 
    553   _select = [] 
    554   for sel in select: 
    555     if sel[0] == ':': _select.append("0"+sel) 
    556     else: _select.append(sel) 
    557   for i in range(len(_select)): 
    558     if _select[i][-1] == ':': select[i] = _select[i]+str(len(x)) 
    559     else: select[i] = _select[i] 
    560   for i in range(len(select)): 
    561     p = select[i].split(":") 
    562     if p[0][0] == '-': p[0] = "len(x)"+p[0] 
    563     if p[1][0] == '-': p[1] = "len(x)"+p[1] 
    564     select[i] = p[0]+":"+p[1] 
    565   steps = [eval("range(%s)" % sel.replace(":",",")) for sel in select] 
    566  
    567   # at this point, we should have: 
    568   #steps = [[0,1],[1,2],[2,3],[3,4,5,6,7,8]] or similar 
    569   if flatten: 
    570     from mystic.tools import flatten 
    571     steps = [list(flatten(steps))] 
    572  
    573   # plot all the scenario "data" 
    574   from numpy import inf, e 
    575   scale = e**(scale - 1.0) 
    576   for v in range(len(steps)): 
    577     if len(steps[v]) > 1: qp = float(max(steps[v])) 
    578     else: qp = inf 
    579     for s in steps[v]: 
    580       par = eval("[params[q][%s][0] for q in range(len(params))]" % s) 
    581       pm = scenario() 
    582       pm.load(par, npts) 
    583       d = dataset() 
    584       d.load(pm.coords, pm.values) 
    585       # dot color determined by number of simultaneous iterations 
    586       t = str((s/qp)**scale) 
    587       # get and plot dataset coords for selected axes       
    588       _coords = get_coords(d, xs, cs) 
    589       # check if we are replacing an axis 
    590       if _2D and xs == 0: 
    591         if data: # adjust data so cone axis is last  
    592           _coords = [list(reversed(pt[:2]))+pt[2:] for pt in _coords] 
    593       elif not _2D and vertical_cones and xs in range(len(bounds)): 
    594         if data: # adjust data so cone axis is last  
    595           _coords = [swap(pt,xs) for pt in _coords] 
    596       plot_data(a[v], _coords, bounds, color=t, strict=strict) 
    597  
    598   if _2D: # strict = True 
    599     plt.xlim((bounds[0][0],bounds[0][1])) 
    600     plt.ylim((bounds[1][0],bounds[1][1])) 
    601   plt.show()  
    602  
    60319# EOF 
Note: See TracChangeset for help on using the changeset viewer.