I'm trying to move a Covid-19 model programmed in python to Mathematica and I can't figure out how to translate this segment of the code to Mathematica.

def maxhosp(Ntest, f=1.0):
# Total population, N.
N = 31914989
# Initial number of infected and recovered individuals, I0 and R0.
I0, Q0, R0 = 1000, 0, 0
# Everyone else, S0, is susceptible to infection initially.
S0 = N - I0 - R0 - Q0
# Contact rate, beta, and mean recovery rate, gamma, (in 1/days).
beta, gamma = 1./6.0, 1.0/15.0 
# A grid of time points (in days)
t = np.linspace(0, 365, 365)
# The SIR model differential equations.
hosprate = 0.13
def deriv(y, t, N, beta, gamma):
    S, I, Q, R = y
    dSdt = -beta * S * I / N
    dIdt = beta * S * I / N - Ntest* f*posrate* I/N - gamma * I
    dQdt = Ntest* f*posrate* I/N - gamma * Q
    dRdt = gamma * (I+Q)
    return dSdt, dIdt, dQdt, dRdt
# Initial conditions vector
y0 = S0, I0, Q0, R0
# Integrate the SIR equations over the time grid, t.
ret = odeint(deriv, y0, t, args=(N, beta, gamma))
S, I, Q, R = ret.T
return np.max(hosprate*(I+Q)/N)
ntests = np.logspace(3, np.log10(20e6), 300)
sims1 = [maxhosp(n, f=1.0) for n in ntests]
sims3 = [maxhosp(n, f=3.0) for n in ntests]
sims10 = [maxhosp(n, f=10.0) for n in ntests]
sims30 = [maxhosp(n, f=30.0) for n in ntests]
# plot in frac of population
fig, ax = plt.subplots(figsize=(8,6))
ax.axhline(0.01, color='r', lw=1, linestyle='-', label='#camas en \nel Perú')
ax.plot(ntests/1e6, np.array(sims1)*N/1e6, lw=2, label='$f=1.0$')

Is it possible to create a function like maxhosp that takes ntest and f as inputs and returns the max value hosprate*(I+Q)/N for each given n and be able to graph it?

  • 1
    $\begingroup$ What you are asking is undoubtedly possible, but I don’t read Python, so can you explain what maxhosp should do and how it is applied, ad to what kind of data? It seems to me that you would be carrying out some sort of Map operation, in which you carry out a similar operation on all members of an array, using each item as input to your function in turn. If you want two paired inputs, then look at Fold. If you want all possible combinations of inputs, then probably Outer would help. $\endgroup$ – MarcoB Jun 2 '20 at 3:08
  • 1
    $\begingroup$ It would also be helpful if you could describe your model in the abstract; perhaps an idiomatic Mathematica approach could be very achieved more easily. $\endgroup$ – MarcoB Jun 2 '20 at 3:10
  • $\begingroup$ From your code, it looks like you are solving some ODEs, in WM you can use DSolve[] or NDSolve[] functions to do this, just look at the documentation for there functions. Also, try to search WM code for epidemic dynamics, like demonstrations.wolfram.com/SIREpidemicDynamics looks similar. $\endgroup$ – I.M. Jun 2 '20 at 6:41
  • $\begingroup$ Is there any part of this that you can do yourself? Defining the variables, at least? Make an attempt with NDSolve? For what it's worth, here is an example of how to solve a similar model with NDSolve. $\endgroup$ – C. E. Jun 2 '20 at 6:57
  • $\begingroup$ If you already have the python program, why do you need a Mathematica version? In other words, what is the context of your question? $\endgroup$ – Somos Jun 2 '20 at 12:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.