Source code for pynucastro.reduction.generate_data
#!/usr/bin/env python3
import numpy as np
from pynucastro import Composition, Nucleus
from pynucastro.reduction.load_network import load_network
[docs]
def dataset(network, n=10, permute=True, b_rho=None, b_T=None, b_Z=None):
"""
Generate a dataset with *n* datapoints. Will either be returned as a sequence of tuples,
each with order (composition, density, temperature) if *permute* is *True* (default),
or the transpose of that if *permute* is *False*. The parameters *b_rho*, *b_T*, and
*b_Z* are tuples giving the bounds on density, temperature, and metallicity respectively.
"""
if isinstance(n, int):
n = np.ones(3, dtype=np.int32) * n
else:
n = np.array(list(n), dtype=np.int32)
# Bounds on each variable
if b_rho is None:
b_rho = (1e2, 1e6) # density (g/cm^3)
if b_T is None:
b_T = (8e6, 1.5e9) # temperature (K)
if b_Z is None:
b_Z = (0.02, 0.2) # metallicity
rho = np.geomspace(*b_rho, num=n[1])
T = np.geomspace(*b_T, num=n[2])
comp_list = []
for Z in np.linspace(*b_Z, num=n[0]):
comp = Composition(network.get_nuclei())
# 75/25 H1/He4 ratio
comp.X[Nucleus("p")] = (1. - Z) * 0.75
comp.X[Nucleus("he4")] = (1. - Z) * 0.25
# 50% of remaining stuff in carbon-12
comp.X[Nucleus("c12")] = Z * 0.5
# Split the rest across all other nuclei
omit = {Nucleus("p"), Nucleus("he4"), Nucleus("c12")}
rem = 0.5 * Z / (len(comp.X) - len(omit))
for nuc in comp.X:
if nuc not in omit:
comp.X[nuc] = rem
comp.normalize()
comp_list.append(comp)
if not permute:
yield comp_list
yield rho
yield T
return
for rho_i in rho:
for T_i in T:
for comp in comp_list:
yield (comp, rho_i, T_i)
[docs]
def main():
network = load_network(Nucleus("ni56"))
conds_list = list(dataset(network))
rho = sorted({conds[1] for conds in conds_list})
T = sorted({conds[2] for conds in conds_list})
comp = {conds[0] for conds in conds_list}
print("ρ")
print(rho)
print()
print("T")
print(T)
print()
def comp_converter(comp):
X = comp.X[Nucleus("p")]
Y = comp.X[Nucleus("he4")]
ZC12 = comp.X[Nucleus("c12")]
Zmin = min(comp.X.values())
return (X, Y, ZC12, Zmin)
print("X")
print(np.array(list(map(comp_converter, comp))))
if __name__ == "__main__":
main()