SloppyCell.ReactionNetworks.Plotting

1 import sets 2 3 import scipy 4 5 import SloppyCell.Plotting 6 from SloppyCell.Plotting import * 7 import Network_mod 8 9 import logging 10 logger = logging.getLogger('ReactionNetworks.Plotting') 11

12 -def PlotEigenvectors(eigVects, net = None, title = None):

13 nEv = 3 14 nOv = len(eigVects[:,0]) 15 for jj in range(nEv): 16 subplot(nEv, 1, jj+1) 17 if jj == 0 and title is not None: 18 title(title) 19 20 bar(range(nOv), eigVects[:,jj]/scipy.linalg.norm(eigVects[:,jj])) 21 axis([-1, nOv] + axis()[2:]) 22 23 if net is not None: 24 mags = zip(abs(eigVects[:,jj]), range(nOv), eigVects[:,jj]) 25 mags.sort() 26 mags.reverse() 27 for mag, index, val in mags[:5]: 28 name = net.optimizableVars[index].name 29 if name is None: 30 name = net.optimizableVars[index].id 31 text(index, val + scipy.sign(val)*0.05, 32 name, 33 horizontalalignment='center', 34 verticalalignment='center') 35 36 a = list(axis()) 37 a[0:2] = [-.03*nOv, nOv*1.03] 38 a[2] -= 0.1 39 a[3] += 0.1 40 axis(a)

41

42 -def plotStateSpaceTrajectoriesForVariables(traj, id1, id2, thresholds = None):

43 xx = traj.getVariableTrajectory(id1) 44 yy = traj.getVariableTrajectory(id2) 45 plot(xx, yy) 46 plot([xx[0]], [yy[0]], 'or') 47 xlabel(id1) 48 ylabel(id2) 49 if thresholds is not None: 50 a = axis() 51 vlines([thresholds[0]], a[2], a[3]) 52 hlines([thresholds[1]], a[0], a[1])

53

54 -def plotTrajectoriesForVariables(traj, ids = None, showLegend = True):

55 if ids is None: 56 ids = traj.net.variables.keys() 57 58 cW = ColorWheel() 59 60 lines = [] 61 legend = [] 62 for id in ids: 63 line = plot(traj.timepoints, traj.getVariableTrajectory(id), 64 cW.next()[::2]) 65 lines.append(line) 66 legend.append(id) 67 68 if showLegend: 69 legend(lines, legend)

70 71

72 -def PlotTrajectoriesForExperiments(model, experiments, params = None, with_data=True, 73 plotPts = 100, overlap = .1, skip = 1, showLegend=True):

74 75 # First find the maximum time in our data 76 maxTime = 0 77 chemsNeededByCalc = {} 78 for exptName in experiments: 79 dataByCalc = model.exptColl[exptName].GetData() 80 for calc in dataByCalc: 81 chemsNeededByCalc.setdefault(calc, []) 82 for chem in dataByCalc[calc].keys(): 83 chemsNeededByCalc[calc].append(chem) 84 thisMaxTime = max(dataByCalc[calc][chem].keys()) 85 if thisMaxTime > maxTime: 86 maxTime = thisMaxTime 87 88 lines = [] 89 legend = [] 90 times = scipy.linspace(0, maxTime*(1 + overlap), plotPts) 91 varsByCalc = {} 92 for calc in chemsNeededByCalc: 93 varsByCalc[calc] = {} 94 for chem in chemsNeededByCalc[calc]: 95 varsByCalc[calc][chem] = times 96 97 model.GetCalculationCollection().Calculate(varsByCalc, params) 98 calcVals = model.GetCalculationCollection().GetResults(varsByCalc) 99 cW = ColorWheel() 100 for exptName in experiments: 101 expt = exptColl[exptName] 102 dataByCalc = expt.GetData() 103 for calc in dataByCalc: 104 for chem in dataByCalc[calc]: 105 color, sym, dash = cW.next() 106 107 if with_data: 108 for time, (data, error) in dataByCalc[calc][chem].items()[::skip]: 109 errorbar(time, data, yerr=error, color=color, 110 mfc = color, marker=sym, 111 ecolor=color, capsize=6) 112 113 predicted = scipy.array(calcVals[calc][chem].items()) 114 order = scipy.argsort(predicted[:,0]) 115 predicted = scipy.take(predicted, order) 116 predicted[:,1] = predicted[:,1] *\ 117 model.GetScaleFactors()[exptName][chem] 118 lines.append(plot(predicted[:,0], predicted[:,1], 119 color=color, linestyle=dash, 120 linewidth = 3)) 121 legend.append(chem + ' in ' + str(calc))# + ' for ' + str(exptName)) 122 123 124 if showLegend: 125 legend(lines, legend, loc=4)

126

127 -def PlotDataForExperiments(model, experiments, skip = 1):

128 exptColl = model.GetExperimentCollection() 129 130 cW = ColorWheel() 131 for exptName in experiments: 132 expt = exptColl[exptName] 133 dataByCalc = expt.GetData() 134 for calc in dataByCalc: 135 for chem in dataByCalc[calc]: 136 color, sym, dash = cW.next() 137 d = scipy.zeros((len(dataByCalc[calc][chem].values()[::skip]), 138 3), scipy.Float) 139 for ii, (time, (data, error))\ 140 in enumerate(dataByCalc[calc][chem].items()[::skip]): 141 d[ii] = [time, data, error] 142 143 errorbar(d[:,0], d[:,1], yerr=d[:,2], color=color, 144 mfc=color, marker=sym, 145 ecolor=color, capsize=6)

146

147 -def plot_model_data(model, expts = None, style = 'errorbars', 148 show_legend = True, loc = 'upper left'):

149 """ 150 Plot the data in the given experiments for the given model. 151 152 Note: You may need to run a Plotting.show() to display the plot. 153 154 Inputs: 155 model: Model whose experiments to plot 156 expts: List of experiment IDs to plot 157 style: Style of plot. Currently supported options are: 158 'errorbars': Plots points and bars for each data point 159 'lines': Plots a continuous line for the data 160 show_legend: Boolean that control whether or not to show the legend 161 loc: Location of the legend. See help(Plotting.legend) for options. 162 """ 163 return plot_model_results(model, expts, style, show_legend, loc, 164 plot_trajectories = False)

165

166 -def plot_model_results(model, expts = None, style='errorbars', 167 show_legend = True, loc = 'upper left', 168 plot_data = True, plot_trajectories = True, 169 data_to_plot=None):

170 """ 171 Plot the fits to the given experiments for the last cost evalution of the 172 model. 173 174 Note: You may need to run a Plotting.show() to display the plot. 175 176 Inputs: 177 model: Model whose results to plot 178 expts: List of experiment IDs to plot, if None is specified, all 179 experiments are plotted 180 style: Style of plot. Currently supported options are: 181 'errorbars': Plots points and bars for each data point 182 'lines': Plots a continuous line for the data 183 show_legend: Boolean that control whether or not to show the legend 184 loc: Location of the legend. See help(Plotting.legend) for options. 185 plot_data: Boolean that controls whether the data is plotted 186 plot_trajectories: Boolean that controls whether the trajectories are 187 plotted 188 data_to_plot: If None, all data variables will be plotted. Otherwise, 189 pass a list of id's and only variables in that list 190 will be plotted. 191 """ 192 exptColl = model.get_expts() 193 calcColl = model.get_calcs() 194 195 lines, labels = [], [] 196 cW = ColorWheel() 197 198 if expts is None: 199 expts = exptColl.keys() 200 201 for exptId in expts: 202 expt = exptColl[exptId] 203 dataByCalc = expt.GetData() 204 for ds in expt.scaled_extrema_data: 205 dataByCalc.setdefault(ds['calcKey'], {}) 206 dataByCalc[ds['calcKey']].setdefault(ds['var'], {}) 207 # We sort the calculation names for easier comparison across plots 208 sortedCalcIds = dataByCalc.keys() 209 sortedCalcIds.sort() 210 for calcId in sortedCalcIds: 211 # Pull the trajectory from that calculation, defaulting to None 212 # if it doesn't exist. 213 net = calcColl.get(calcId) 214 traj = getattr(net, 'trajectory', None) 215 for dataId, dataDict in dataByCalc[calcId].items(): 216 # Skip this variable if it's not amongst our data_to_plot 217 # list. 218 if (data_to_plot is not None) and (dataId not in data_to_plot): 219 continue 220 color, sym, dash = cW.next() 221 222 if plot_trajectories: 223 if traj is None: 224 print 'No trajectory in calculation %s!' % calcId 225 print 'The cost must be evaluated before the results', 226 print 'can be plotted.' 227 return 228 229 scaleFactor = model.GetScaleFactors()[exptId][dataId] 230 result = scaleFactor*traj.getVariableTrajectory(dataId) 231 l = plot(traj.timepoints, result, color=color, 232 linestyle=dash,linewidth=3) 233 234 # We superimpose a dotted black line to distinguish 235 # theory from data in this case 236 if style is 'lines': 237 plot(traj.timepoints, result, 'k--', linewidth=3, 238 zorder = 10) 239 240 241 if plot_data and dataDict: 242 # Pull the data out of the dictionary and into an array 243 d = scipy.array([[t, v, e] for (t, (v, e)) 244 in dataDict.items()]) 245 if style is 'errorbars': 246 l = errorbar(d[:,0], d[:,1], yerr=d[:,2], color=color, 247 mfc = color, linestyle='', marker=sym, 248 ecolor='k', capsize=6)[0] 249 elif style is 'lines': 250 # Make sure we order the data before plotting 251 order = scipy.argsort(d[:,0], 0) 252 d = scipy.take(d, order, 0) 253 l = plot(d[:,0], d[:,1], color=color, 254 linestyle=dash) 255 lines.append(l) 256 257 # Plot the extra data points. 258 if plot_data: 259 for res in model.residuals: 260 if isinstance(res, Residuals.ScaledExtremum)\ 261 and res.exptKey == exptId and res.calcKey == calcId: 262 t = res.last_time_result 263 val = res.yMeas 264 sigma = res.ySigma 265 errorbar([t], [val], [sigma], color=color, 266 linestyle='', marker=sym, 267 ecolor='k', capsize=6, mfc='w', 268 mec= color, mew=2, ms=10)[0] 269 270 # Let's print the pretty name for our variable if we can. 271 lines.append(l) 272 name = net.get_component_name(dataId) 273 labels.append('%s in %s for %s' % (name, calcId, exptId)) 274 275 if show_legend: 276 legend(lines, labels, loc=loc) 277 278 return lines, labels

279

280 -def plot_ensemble_results(model, ensemble, expts = None, 281 style='errorbars', 282 show_legend = True, loc = 'upper left', 283 plot_data = True, plot_trajectories = True):

284 """ 285 Plot the fits to the given experiments over an ensemble. 286 287 Note that this recalculates the cost for every member of the ensemble, so 288 it may be very slow. Filtering correlated members from the ensemble is 289 strongly recommended. 290 291 Inputs: 292 model: Model whose results to plot 293 ensemble: Parameter ensemble 294 expts: List of experiment IDs to plot, if None is specified, all 295 experiments are plotted 296 style: Style of plot. Currently supported options are: 297 'errorbars': Plots points and bars for each data point 298 'lines': Plots a continuous line for the data 299 show_legend: Boolean that control whether or not to show the legend 300 loc: Location of the legend. See help(Plotting.legend) for options. 301 plot_data: Boolean that controls whether the data is plotted 302 plot_trajectories: Boolean that controls whether the trajectories are 303 plotted 304 """ 305 exptColl = model.get_expts() 306 nets = model.get_calcs() 307 308 if expts is None: 309 expts = exptColl.keys() 310 311 lines, labels = [], [] 312 cW = ColorWheel() 313 314 Network_mod.Network.pretty_plotting() 315 model.cost(ensemble[0]) 316 timepoints = {} 317 for netId, net in nets.items(): 318 traj = getattr(net, 'trajectory', None) 319 if traj is not None: 320 net.times_to_add = scipy.linspace(traj.timepoints[0], 321 traj.timepoints[-1], 1000) 322 323 Network_mod.Network.full_speed() 324 325 results = {} 326 for params in ensemble: 327 model.cost(params) 328 for exptId in expts: 329 expt = exptColl[exptId] 330 results.setdefault(exptId, {}) 331 dataByCalc = expt.GetData() 332 for netId in dataByCalc.keys(): 333 results[exptId].setdefault(netId, {}) 334 # Pull the trajectory from that calculation, defaulting to None 335 # if it doesn't exist. 336 net = nets.get(netId) 337 traj = net.trajectory 338 for dataId in dataByCalc[netId].keys(): 339 results[exptId][netId].setdefault(dataId, []) 340 341 scaleFactor = model.GetScaleFactors()[exptId][dataId] 342 result = scaleFactor*traj.get_var_traj(dataId) 343 results[exptId][netId][dataId].append(result) 344 345 for exptId in expts: 346 expt = exptColl[exptId] 347 dataByCalc = expt.GetData() 348 # We sort the calculation names for easier comparison across plots 349 sortedCalcIds = dataByCalc.keys() 350 sortedCalcIds.sort() 351 for netId in sortedCalcIds: 352 for dataId, dataDict in dataByCalc[netId].items(): 353 color, sym, dash = cW.next() 354 355 if plot_data: 356 # Pull the data out of the dictionary and into an array 357 d = scipy.array([[t, v, e] for (t, (v, e)) 358 in dataDict.items()]) 359 if style is 'errorbars': 360 l = errorbar(d[:,0], d[:,1], yerr=d[:,2], fmt='o', 361 color=color, markerfacecolor=color, 362 marker=sym, ecolor='k', capsize=6)[0] 363 elif style is 'lines': 364 # Make sure we order the data before plotting 365 order = scipy.argsort(d[:,0], 0) 366 d = scipy.take(d, order, 0) 367 l = plot(d[:,0], d[:,1], color=color, 368 linestyle=dash) 369 lines.append(l) 370 371 if plot_trajectories: 372 times = model.get_calcs().get(netId).trajectory.get_times() 373 mean_vals = scipy.mean(results[exptId][netId][dataId], 0) 374 std_vals = scipy.std(results[exptId][netId][dataId], 0) 375 376 lower_vals = mean_vals - std_vals 377 upper_vals = mean_vals + std_vals 378 379 # Plot the polygon 380 xpts = scipy.concatenate((times, times[::-1])) 381 ypts = scipy.concatenate((lower_vals, upper_vals[::-1])) 382 fill(xpts, ypts, fc=color, alpha=0.4) 383 384 # Let's print the pretty name for our variable if we can. 385 name = net.get_component_name(dataId) 386 labels.append('%s in %s for %s' % (name, netId, exptId)) 387 388 for netId, net in nets.items(): 389 del net.times_to_add 390 391 if show_legend: 392 legend(lines, labels, loc=loc) 393 394 for net in nets.values(): 395 net.times_to_add = None 396 397 return lines, labels

398 399

400 -def plot_trajectory(traj, vars = None, 401 show_legend = True, loc = 'upper left', 402 logx = False, logy = False):

403 if vars is None: 404 vars = traj.dynamicVarKeys 405 406 plot_funcs_dict = {(False, False): plot, 407 (True, False): semilogx, 408 (False, True): semilogy, 409 (True, True): loglog} 410 plot_func = plot_funcs_dict[(logx, logy)] 411 412 cW = ColorWheel(symbols=None) 413 lines, labels = [], [] 414 for id in vars: 415 color, sym, dash = cW.next() 416 label = str(id) 417 line = plot_func(traj.timepoints, traj.getVariableTrajectory(id), 418 color = color, mfc = color, marker=sym, linestyle=dash, 419 linewidth=3, label=label) 420 lines.append(line) 421 labels.append(label) 422 423 if show_legend: 424 legend(loc=loc) 425 426 return (lines, labels)

427

428 -def plot_deriv_trajectory(traj, vars = None, 429 show_legend = True, loc = 'upper left', 430 logx = False, logy = False):

431 if vars is None: 432 vars = traj.dynamicVarKeys 433 434 plot_funcs_dict = {(False, False): plot, 435 (True, False): semilogx, 436 (False, True): semilogy, 437 (True, True): loglog} 438 plot_func = plot_funcs_dict[(logx, logy)] 439 440 cW = ColorWheel(symbols=None) 441 lines, labels = [], [] 442 for id in vars: 443 try: 444 color, sym, dash = cW.next() 445 label = str(id) 446 line = plot_func(traj.timepoints, traj.get_var_traj((id, 'time')), 447 color = color, mfc = color, marker=sym, linestyle=dash, 448 linewidth=3, label=label) 449 lines.append(line) 450 labels.append(label) 451 except ValueError: 452 logger.warn('Derivative of variable %s not present in trajectory. Make sure '\ 453 'return_derivs is set to True in integrate function.' % (id)) 454 455 456 if show_legend: 457 legend(loc=loc) 458 459 return (lines, labels)

460

461 -def plot_ensemble_trajs(best_traj=None, mean_traj=None, 462 std_traj=None, std_devs=1.0, 463 vars=None, 464 show_legend = True, loc = 'upper left'):

465 """ 466 Plot the results of a Network ensemble. 467 468 Inputs: 469 best_traj -- Best-fit trajectory 470 mean_traj -- Mean trajectory 471 std_traj -- Trajectory of standard deviations 472 std_devs -- Number of standard deviations to draw bounds at 473 vars -- List of variable ids to plot the bounds for 474 475 show_legend -- Boolean to show a legend or not 476 loc -- Location code for legend 477 """ 478 cW = ColorWheel() 479 lines = [] 480 labels = [] 481 for var in vars: 482 color, sym, dash = cW.next() 483 if best_traj is not None: 484 l = plot(best_traj.timepoints, best_traj.getVariableTrajectory(var), 485 linestyle='-', color=color, linewidth=2) 486 if mean_traj is not None: 487 times = mean_traj.timepoints 488 mean_vals = mean_traj.getVariableTrajectory(var) 489 plot(times, mean_vals, linestyle = '--', linewidth=2, 490 color=color) 491 if mean_traj is not None and std_traj is not None and std_devs > 0: 492 # Calculate the bounds 493 std_vals = std_traj.getVariableTrajectory(var) 494 lower_vals = mean_vals - std_devs*std_vals 495 upper_vals = mean_vals + std_devs*std_vals 496 497 # Plot the polygon 498 xpts = scipy.concatenate((times, times[::-1])) 499 ypts = scipy.concatenate((lower_vals, upper_vals[::-1])) 500 fill(xpts, ypts, facecolor=color, alpha=0.4) 501 502 lines.append(l) 503 labels.append(var) 504 505 if show_legend: 506 legend(lines, labels, loc=loc)

507

Source Code for Module SloppyCell.ReactionNetworks.Plotting