Initial Conditions for estimating parameters of the SIR Model.

Hi Prof. @lululxvi , I was trying to infer the parameters of the SIR model based on my implementation which is similar to Systems biology informed deep learning for inferring parameters and hidden dynamics. The SIR ODE equations are generally represented as follows: SIR Model Equations. I did not apply any input feature transformation function as I want to process my raw input of time as it is and thus also did My code is as follows:

import numpy as np
from scipy.integrate import odeint
import random
import deepxde as dde
from deepxde.backend import tf
random.seed(5)
# Total population, N.
N = 7900010
# Initial number of infected and recovered individuals, I0 and R0.
I0, R0 = 10, 0
# Everyone else, S0, is susceptible to infection initially.
S0 = N - I0 - R0
# Contact rate, beta, and mean recovery rate, gamma, (in 1/days).
beta, gamma = 0.5, 0.33  
# A grid of time points (in days)
t = np.linspace(1, 150, 150)

# The SIR model differential equations.
def SIR_model(y, t, N, beta, gamma):
    S, I, R = y
    dSdt = -beta * S * I/N
    dIdt = beta * S * I / N - gamma * I
    dRdt = gamma * I
    return dSdt, dIdt, dRdt

# Initial conditions vector
y0 = S0, I0, R0
# Integrate the SIR equations over the time grid, t.
#t = np.arange(1, 365,1)[:, None]
ret = odeint(SIR_model, y0, t, args=(N, beta, gamma))
S, I, R = ret.T
t = t[:, None]
y=ret
def pinn(data_t, data_y):
  beta=tf.math.softplus(tf.Variable(0, trainable=True, dtype=tf.float32))
  gamma=tf.math.softplus(tf.Variable(0, trainable=True, dtype=tf.float32))
  var_list=[beta,gamma]
  def ODE(t,y):
    return [tf.gradients(y[:, 0:1], t)[0]-(-beta * y[0] * y[1]/N),
            tf.gradients(y[:, 1:2], t)[0]-((beta * y[0] * y[1] / N) - gamma * y[1]),
            tf.gradients(y[:, 2:3], t)[0]-(gamma*y[1])]
  geom = dde.geometry.TimeDomain(data_t[0, 0], data_t[-1, 0])
  def boundary(x, _):
    return np.isclose(x[0], data_t[-1, 0])
  y1 = data_y[-1]
  bc0 = dde.DirichletBC(geom, lambda X: y1[0], boundary, component=0)
  bc1 = dde.DirichletBC(geom, lambda X: y1[1], boundary, component=1)
  bc2 = dde.DirichletBC(geom, lambda X: y1[2], boundary, component=2)
  # Observes
  n = len(data_t)
  idx = np.append(
      np.random.choice(np.arange(1, n - 1), size=n // 4, replace=False), [0, n - 1])
  ptset = dde.utils.external.PointSet(data_t[idx])
  inside = lambda x, _: ptset.inside(x)
  observe_I=dde.DirichletBC(geom, ptset.values_to_func(data_y[idx, 0:1]), inside, component=1)
  np.savetxt("SIR_input.dat", np.hstack((data_t[idx], data_y[idx, 0:1])))
  data = dde.data.PDE(geom,ODE,[bc0,bc1,observe_I],anchors=data_t)
  net = dde.maps.FNN([1] + [128] * 3 + [2], "swish", "Glorot normal")
  # No Feature Transformation
  # No Output Transformation
  model = dde.Model(data, net)
  checkpointer = dde.callbacks.ModelCheckpoint("./model/model.ckpt", verbose=1, save_better_only=True, period=1000)
  variable = dde.callbacks.VariableValue(var_list, period=1000, filename="variables.dat", precision=3,)
  callbacks = [checkpointer, variable]
  bc_weights = [1, 1]
  data_weights = [1e3]
  model.compile("adam", lr=1e-3,loss_weights=[0] * 3 + bc_weights + data_weights)
  model.train(epochs=10000, display_every=1000)
  losshistory, train_state = model.train(epochs=90000,display_every=1000,callbacks=callbacks,disregard_previous_best=True,)
  #dde.saveplot(losshistory, train_state, issave=True, isplot=True)
  #var_list1 = [model.sess.run(v) for v in var_list]
  #print(var_list1)
  ode_weights = [1e-3, 1e-3, 1e-3]
  model.compile("adam", lr=1e-3, loss_weights=ode_weights + bc_weights + data_weights)
  losshistory, train_state = model.train(epochs=90000,display_every=1000,callbacks=callbacks,disregard_previous_best=True,)
  dde.saveplot(losshistory, train_state, issave=True, isplot=True)
  var_list = [model.sess.run(v) for v in var_list]
  return var_list
var_list = pinn(t, y)

However, after the number of iterations, the inferred parameters [beta, gamma] are only [1.1191349e-16, 0.0049484214], which is much less than the original parameters I used to create the synthetic data with. Is there some sort of mistake made by me while modeling (wrt to the boundary conditions etc) or should I try to increase the number of epochs (which I did in fact, to no better results)? I will be extremely grateful to you if you guide me towards the correct direction here.

Best, Akash Choudhuri.

lululxvi / deepxde

Initial Conditions for estimating parameters of the SIR Model. #515