Plotting Examples#

There are a variety of methods that can be used to visualize a network, which we demonstrate below. For most plots, NetworkX is used to create the structure of the network and matplotlib manages the drawing.

Tip

You can interactively zoom and pan a network’s visualization within Jupyter by installing the ipympl package and then adding:

%matplotlib widget

to the notebook (in a cell that is executed).

import pynucastro as pyna

Standard view#

Let’s build a quick network that represents He, C, and O burning

library = pyna.ReacLibLibrary()
sub = library.linking_nuclei(["p", "n", "he4", "c12", "n13", "n14", "o16",
                              "ne20", "ne21", "na23", "mg23", "mg24", "al27",
                              "si27", "si28"])
rc = pyna.RateCollection(libraries=sub)
/opt/hostedtoolcache/Python/3.11.12/x64/lib/python3.11/site-packages/pynucastro/networks/rate_collection.py:689: UserWarning: ReacLib neutron decay rate (<n_to_p_weak_wc12>) does not account for degeneracy at high densities. Consider using tabular rate from Langanke.
  warnings.warn(msg)

The standard network plot of this is:

fig = rc.plot()
_images/fd753319954cdadc0317d25a077b6853938e9bc8e0b3c40f424383bb28fd4db7.png

Rotated view#

Since everything lies close to the diagonal, it is hard to see the structure. We can instead look at a rotated view of the network.

While we’re at it, we’ll also color the rates by reaction rate strength

rho = 1.e4
T = 1.e8
comp = pyna.Composition(rc.get_nuclei())
comp.set_solar_like()
fig = rc.plot(rho, T, comp, rotated=True)
_images/203441ece45cadb36c29bafaad79ac1d286823f98af9101927f15449ed51262f.png

Curved edges#

We can also plot the edges curved, to help distinguish the forward and reverse rates

fig = rc.plot(rho, T, comp, rotated=True, hide_xalpha=True,
              curved_edges=True)
_images/fb2b0b380fc69d8aa2c8f4975562f4850204790f4126dcccd0e6dcacf680e5db.png

Adding a legend#

We can also add a “compass rose” style legend that shows the directions of key captures. We specify the location in terms of (Z, N)

fig = rc.plot(rho, T, comp, rotated=True, hide_xalpha=True,
              curved_edges=True, legend_coord=(3, 2))
_images/45f3781573a06830294f8c32dc59a3f0045341231a161c8832233c3b643b576a.png

Labeling edges#

Any edge can be labelled.

edge_labels = {(pyna.Nucleus("he4"), pyna.Nucleus("c12")):
               r"$\alpha(\alpha\alpha,\gamma){}^{12}\mathrm{C}$"}
fig = rc.plot(rho, T, comp, rotated=True, hide_xalpha=True,
              edge_labels=edge_labels, size=(1200, 700),
              curved_edges=True)
_images/0aba8f74a7e603ee326994c0beec4e6dc690201a4540ed126609c39d9cd56af2.png

Filters#

We can use a rate filter to show only certain links. For example, to show only links involving proton capture, we can do:

fig = rc.plot(rho, T, comp, rotated=True, hide_xalpha=True,
              rate_filter_function=lambda r: pyna.Nucleus("p") in r.reactants)
_images/f2a8caf1f3bc673bdb971844b470f4f896693ecf5ebfc46ee4d0948017df5255.png

We can also highlight rates using the highlight_filter_function mechanism. This takes a function that operates on a rate and returns True if we wish to highlight the rate (i.e., as if we used a highlighter marker over the line).

Here’s an example where we highlight the rates that have \(Q > 0\)

fig = rc.plot(rotated=True, hide_xalpha=True, size=(1500, 600),
              highlight_filter_function=lambda r: r.Q > 0,
              curved_edges=True)
_images/690dcd941134eaf66402114d6d8afc462e94d00bb47ac775b2071a53bd4d6a7a.png

Plotting the Jacobian#

For implicit integration, we solve:

\[\dot{\bf Y} = {\bf f}({\bf Y})\]

and the Jacobian, \({\bf J} = \partial {\bf f} / \partial {\bf Y}\) plays an important role in solving for the update.

We can visualize the Jacobian of our network easily. By default, we will show all elements where

\[\frac{\mathrm{max}_{i,j} (|J_{ij}|)}{|J_{ij}|} > 10^{10}\]

but this ratio can be set via the rate_scaling keyword argument.

rho = 1.e4
T = 2.e9
comp.set_equal()
fig = rc.plot_jacobian(rho, T, comp, rate_scaling=1.e10)
_images/2cecfb8dbad893f66791983792e6a5f081fe4051f2fbc695474cbeaf9b20392b.png

Note

If you set rate_scaling too large, then roundoff error in mapping the Jacobian elements to colors can cause elements that are zero to appear colored. It is best to stay below 1.e15 to minimize roundoff.

Plot nuclides on a grid#

Nuclides in a network may also be visualized as cells on a grid of Z vs. N, colored by some quantity. This can be more interpretable for large networks.

Calling gridplot without any arguments will just plot the grid – to see anything interesting we need to supply some conditions. Here is a plot of nuclide mass fraction on a log scale, with a 36 square inch figure:

rate_names = ["c12(p,g)n13", 
              "c13(p,g)n14",
              "n13(,)c13",
              "n13(p,g)o14",
              "n14(p,g)o15",
              "n15(p,a)c12",
              "o14(,)n14",
              "o15(,)n15"]
rates = library.get_rate_by_name(rate_names)
rc = pyna.RateCollection(rates=rates)
rho = 1.e4
T = 1.e8
comp = pyna.Composition(rc.get_nuclei())
comp.set_solar_like()
fig = rc.gridplot(comp=comp, color_field="X", scale="log", area=36)
_images/90421b02571a8246e652ca04ac3c2377c15c2f3a3e6ac4f5dcfc835289e63c44.png

The plot is configurable through a large number of keyword arguments. Here we want to look at the rates at which nuclides are being created or destroyed, so we color by \(\dot{Y}\), the rate of change of molar abundance. Density and temperature need to be supplied to evaluate the rates. A full list of valid keyword arguments can be found in the API documentation.

fig = rc.gridplot(comp=comp, rho=rho, T=T, color_field="ydot", area=36,
                  cmap="RdGy", cbar_bounds=(-0.1, 0.1))
_images/8256f1bd36f79178c8e706e70f7c0f9585416cd5e21d9028a06fbac3985667e9.png

Unlike the network plot, this won’t omit hydrogen and helium by default. To just look at the heavier nuclides, we can define a function to filter by proton number. Here we also plot activity, which is the sum of all rates contributions to Ydot ignoring the sign.

ff = lambda nuc: nuc.Z > 2
fig = rc.gridplot(comp=comp, rho=rho, T=T, color_field="activity", scale="log",
                  filter_function=ff, area=20, cmap="Blues")
_images/cb6c1c17e1986ed7eaa9b2f8a817ff8e029c40faf9c0d4195b7296536f0a6185.png

Network chart#

A network chart shows how each rate affects each nucleus

fig = rc.plot_network_chart(rho, T, comp)
_images/23ffd338a3c8d776c4df2faeeb2b1de11779e9c796c2e91924caa2bc96e6de75.png