<jupyter_start><jupyter_text>function as an arg (solve_ivp)investigate default valuesuse of kwargs (large option set)function as a return value (dense_output)function as an arg 2 (events). + Higher-order : optional level function factory (modèle à compartiment)<jupyter_code>from numpy import *
from numpy.linalg import *
from matplotlib.pyplot import *
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
WEEK = 7
YEAR = 365
N = 100
beta = 1 / (WEEK)
gamma = 1 / (2 * WEEK)
omega = 1 / YEAR
def dSIR(t, SIR):
    S, I, R = SIR
    dS = omega * R - beta * I * S / N
    dI = beta * I * S / N - gamma * I
    dR = gamma * I - omega * R  
    return (dS, dI, dR)
S0, I0 = 99.0, 1.0
R0 = N - S0 - I0
t_span = [0.0, 5*YEAR]
results = solve_ivp(dSIR, t_span=t_span, y0=(S0, I0, R0))
list(results.keys())
def plot(t, y):
    pass
# plot(**results)
def plot(t, y, **_):
    print(list(_.keys()))

plot(**results)
def plot(t, y, **_):
    plt.plot(t, y.T, "-+", label=["S", "I", "R"])
    plt.xlabel("Temps (années)")
    plt.xticks(arange(0, 5*YEAR+1, YEAR), [str(i) for i in range(5+1)])
    plt.ylabel("Populations (%)")
    plt.yticks(arange(0, 101, 10))
    plt.axis([0, 5*YEAR, 0, 100])
    plt.legend()
    plt.grid(True)
plot(**results)
plt.savefig("images/SIR-discrete.svg")
options = eval(open("options.py").read())
options
results = solve_ivp(dSIR, t_span=t_span, y0=(S0, I0, R0), **options)
plot(**results)
plt.savefig("images/SIR-options.svg")
plot(**results)
plt.savefig("images/SIR-options.svg")<jupyter_output><empty_output><jupyter_text>-----<jupyter_code>t_eval = arange(0.0, 5*YEAR, 30)
results = solve_ivp(dSIR, t_span=t_span, t_eval=t_eval, y0=(S0, I0, R0))
t, y = results.t, results.y
plt.plot(t, y.T, "-+", label=["S", "I", "R"])
plt.xlabel("Temps en jours")
plt.ylabel("Populations")
plt.legend()
pass
t = arange(0.0, 5*YEAR, 1)
dense_output = True
t
results = solve_ivp(dSIR, t_span=t_span, dense_output=dense_output, y0=(S0, I0, R0))
tc, yc, sol = results["t"], results["y"], results["sol"]
tc
sol
callable(sol)
def plot(t_, sol, **_):
    plt.plot(t_, sol(t_).T, "-", label=["S", "I", "R"])
    plt.xlabel("Temps (années)")
    plt.xticks(arange(0, 5*YEAR+1, YEAR), [str(i) for i in range(5+1)])
    plt.ylabel("Populations (%)")
    plt.yticks(arange(0, 101, 10))
    plt.axis([0, 5*YEAR, 0, 100])
    plt.legend()
    plt.grid(True)
plot(arange(0, 5*YEAR), **results)
plt.savefig("images/SIR-dense.svg")
def event(t, SIR):
    S, I, R = SIR
    return I - 10.0
results = solve_ivp(dSIR, y0=(S0, I0, R0), t_span=t_span, dense_output=True, events=[event])
def plot(t_, sol, t_events, **_):
    plt.plot(t_, sol(t_).T, "-", label=["S", "I", "R"])
    t_event = t_events[0]
    plt.plot(t_event, sol(t_event)[1], "k.")
    plt.xlabel("Temps (années)")
    plt.xticks(arange(0, 5*YEAR+1, YEAR), [str(i) for i in range(5+1)])
    plt.ylabel("Populations (%)")
    plt.yticks(arange(0, 101, 10))
    plt.axis([0, 5*YEAR, 0, 100])
    plt.legend()
    plt.grid(True)
plot(t, **results)
savefig("images/SIR-events.svg")
I_max = 2.0

def event(t, SIR):
    S, I, R = SIR
    return I - I_max
results = solve_ivp(dSIR, y0=(S0, I0, R0), t_span=t_span, dense_output=True, events=[event])
plot(t, **results)
savefig("images/SIR-events-2.svg")
def I_crosses(I_threshold):
    def event(t, SIR):
        S, I, R = SIR
        return I - I_threshold
    return event
results = solve_ivp(dSIR, y0=(S0, I0, R0), t_span=t_span, dense_output=True, events=[I_crosses(2.0)])
plot(t, **results)
#savefig("images/SIR-events-2bis.svg")
t_events
plt.plot(t, sol(t).T[:,1], "-", color="orange", label="I")
plt.plot(tc, yc.T[:, 1], "+", color="orange")
plt.plot([t[0], t[-1]], [2.0, 2.0], "k", alpha=0.25)
plt.plot(t_events[0], 2.0 * ones_like(t_events[0]), "k+")
plt.xlabel("Temps en jours")
plt.ylabel("Population")
plt.grid(True)
plt.legend()
pass
SIR_dynamics = {
 ("S", "I"): "beta * I * S / N",
 ("I", "R"): "gamma * I",
 ("R", "S"): "omega * R"
}
def get_vars(dynamics):
    vars = []
    for pair in dynamics.keys():
        vars.extend(pair)
    # remove duplicated vars
    return list(set(vars))

def make_fun(dynamics):
    vars = get_vars(dynamics)
    n = len(vars)
    
    def fun(t, state):
        ns = globals().copy()
        for var, value in zip(vars, state):
            ns[var] = value
        dstate = []
        for i in range(n):
            d = 0
            var = vars[i]
            for (edge, expr) in dynamics.items():
                source, target = edge
                if source == var:
                    d -= eval(expr, ns)
                if target == var:
                    d += eval(expr, ns)
            dstate.append(d)
        return dstate
        
    return fun 
dSIR2 = make_fun(SIR_dynamics)
dSIR(None, (1/3, 1/3, 1/3))
dSIR2(None, (1/3, 1/3, 1/3))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.19648760330578513, 'fraction_numerical': 0.032231404958677684, 'mean_word_length': 3.333930170098478, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '261309', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 2143, 'n_tokens_neox': 2024, 'n_words': 464, 'path': 'tps/data/tmp.ipynb', 'repo_name': 'boisgera/python-fr', 'size': '4840'}