Profiling and debugging

Profiling

One of the main reasons people don’t use ABMs is because they can be very slow. While “vanilla Starsim” is quite fast (10,000 agents running for 100 timesteps should take about a second), custom modules, if not properly written, can be quite slow.

The first step of fixing a slow module is to identify the problem. To do this, Starsim includes some built-in profiling tools.

Let’s look at a simple simulation:

import sciris as sc
import starsim as ss
ss.options(jupyter=True)
sc.options(jupyter=True)

pars = dict(
    start = '2000-01-01',
    stop = '2020-01-01',
    diseases = 'sis',
    networks = 'random'
)

# Profile sim
sim = ss.Sim(pars)
prof = sim.profile()

Initializing sim with 10000 agents

Profiling 15 function(s):

 <bound method Sim.run of Sim(n=10000; 2000.01.01—2020.01.01; networks=randomnet; diseases=sis)>

<bound method Sim.start_step of Sim(n=10000; 2000.01.01—2020.01.01; networks=randomnet; diseases=sis [...]

<function Module.start_step at 0x7f8d79434e00>

<function Module.start_step at 0x7f8d79434e00>

<bound method SIS.step_state of sis(pars=[init_prev, beta, dur_inf, waning, imm_boost, _n_initial_ca [...]

<bound method DynamicNetwork.step of randomnet(n_edges=50000; pars=[n_contacts, dur, beta]; states=[ [...]

<bound method Infection.step of sis(pars=[init_prev, beta, dur_inf, waning, imm_boost, _n_initial_ca [...]

<bound method People.step_die of People(n=10000; age=30.2±17.4)>

<bound method People.update_results of People(n=10000; age=30.2±17.4)>

<bound method Network.update_results of randomnet(n_edges=50000; pars=[n_contacts, dur, beta]; state [...]

<function SIS.update_results at 0x7f8d78e1e0c0>

<function Module.finish_step at 0x7f8d79434fe0>

<function Module.finish_step at 0x7f8d79434fe0>

<bound method People.finish_step of People(n=10000; age=30.2±17.4)>

<bound method Sim.finish_step of Sim(n=10000; 2000.01.01—2020.01.01; networks=randomnet; diseases=si [...] 




  Running 2000.01.01 ( 0/21) (0.00 s)  ———————————————————— 5%


  Running 2010.01.01 (10/21) (0.14 s)  ••••••••••—————————— 52%


  Running 2020.01.01 (20/21) (0.21 s)  •••••••••••••••••••• 100%



Elapsed time: 0.240 s





————————————————————————————————————————————————————————————————————

Profile of networks.Network.update_results: 0.000532634 s (0.22175%)

————————————————————————————————————————————————————————————————————



Total time: 0.000532634 s

File: /home/runner/work/starsim/starsim/starsim/networks.py

Function: Network.update_results at line 278



Line #      Hits         Time  Per Hit   % Time  Line Contents

==============================================================

   278                                               def update_results(self):

   279                                                   """ Store the number of edges in the network """

   280        21     108492.0   5166.3     20.4          super().update_results()

   281        21     415045.0  19764.0     77.9          self.results['n_edges'][self.ti] = len(self)

   282        21       9097.0    433.2      1.7          return







——————————————————————————————————————————————————————————————

Profile of people.People.finish_step: 0.00127495 s (0.530796%)

——————————————————————————————————————————————————————————————



Total time: 0.00127495 s

File: /home/runner/work/starsim/starsim/starsim/people.py

Function: People.finish_step at line 508



Line #      Hits         Time  Per Hit   % Time  Line Contents

==============================================================

   508                                               def finish_step(self):

   509                                                   """ Remove dead agents and run post-step updates """

   510        21     948139.0  45149.5     74.4          self.remove_dead()

   511        21     318114.0  15148.3     25.0          self.update_post()

   512        21       8698.0    414.2      0.7          return







——————————————————————————————————————————————————————————

Profile of people.People.step_die: 0.00249719 s (1.03965%)

——————————————————————————————————————————————————————————



Total time: 0.00249719 s

File: /home/runner/work/starsim/starsim/starsim/people.py

Function: People.step_die at line 468



Line #      Hits         Time  Per Hit   % Time  Line Contents

==============================================================

   468                                               def step_die(self):

   469                                                   """ Carry out any deaths or removals that took place this timestep """

   470        21    2153418.0 102543.7     86.2          death_uids = ((self.ti_dead <= self.sim.ti) | (self.ti_removed <= self.sim.ti)).uids

   471        21      91310.0   4348.1      3.7          self.alive[death_uids] = False  # Whilst not dead, removed agents should not be included in alive totals

   472                                           

   473                                                   # Execute deaths that took place this timestep (i.e., changing the `alive` state of the agents). This is executed

   474                                                   # before analyzers have run so that analyzers are able to inspect and record outcomes for agents that died this timestep

   475        42     187330.0   4460.2      7.5          for disease in self.sim.diseases():

   476        21      22060.0   1050.5      0.9              if isinstance(disease, ss.Disease):

   477        21      33645.0   1602.1      1.3                  disease.step_die(death_uids)

   478                                           

   479        21       9430.0    449.0      0.4          return death_uids







———————————————————————————————————————————————————————————————

Profile of diseases.SIS.update_results: 0.00373117 s (1.55339%)

———————————————————————————————————————————————————————————————



Total time: 0.00373117 s

File: /home/runner/work/starsim/starsim/starsim/diseases.py

Function: SIS.update_results at line 803



Line #      Hits         Time  Per Hit   % Time  Line Contents

==============================================================

   803                                               @ss.required()

   804                                               def update_results(self):

   805                                                   """ Store the population immunity (susceptibility) """

   806        21    2687481.0 127975.3     72.0          super().update_results()

   807        21    1034232.0  49249.1     27.7          self.results['rel_sus'][self.ti] = self.rel_sus.mean()

   808        21       9457.0    450.3      0.3          return







————————————————————————————————————————————————————————————————

Profile of people.People.update_results: 0.00464806 s (1.93511%)

————————————————————————————————————————————————————————————————



Total time: 0.00464806 s

File: /home/runner/work/starsim/starsim/starsim/people.py

Function: People.update_results at line 497



Line #      Hits         Time  Per Hit   % Time  Line Contents

==============================================================

   497                                               def update_results(self):

   498                                                   """ Record per-timestep population counts into the simulation results """

   499        21      32212.0   1533.9      0.7          ti = self.sim.ti

   500        21      10309.0    490.9      0.2          res = self.sim.results

   501        63     478074.0   7588.5     10.3          for state in self.auto_state_list: # Count each auto-generated BoolState result, e.g. n_alive, n_female

   502        42    1047650.0  24944.0     22.5              res[f'n_{state.name}'][ti] = np.count_nonzero(getattr(self, state.name))

   503        21    1143791.0  54466.2     24.6          res.new_deaths[ti] = np.count_nonzero(self.ti_dead == ti)

   504        21     991708.0  47224.2     21.3          res.new_emigrants[ti] = np.count_nonzero(self.ti_removed == ti)

   505        21     934400.0  44495.2     20.1          res.cum_deaths[ti] = np.sum(res.new_deaths[:ti]) # TODO: inefficient to compute the cumulative sum on every timestep!

   506        21       9919.0    472.3      0.2          return







——————————————————————————————————————————————————————

Profile of sim.Sim.start_step: 0.00550509 s (2.29192%)

——————————————————————————————————————————————————————



Total time: 0.00550509 s

File: /home/runner/work/starsim/starsim/starsim/sim.py

Function: Sim.start_step at line 451



Line #      Hits         Time  Per Hit   % Time  Line Contents

==============================================================

   451                                               def start_step(self):

   452                                                   """ Start the step -- only print progress; all actual changes happen in the modules """

   453                                           

   454                                                   # Set the time and if we have reached the end of the simulation, then do nothing

   455        21      10431.0    496.7      0.2          if self.complete:

   456                                                       errormsg = 'Simulation already complete (call sim.init() to re-run)'

   457                                                       raise AlreadyRunError(errormsg)

   458                                           

   459                                                   # Print out progress if needed

   460        21    3773842.0 179706.8     68.6          self.elapsed = self.timer.toc(output=True)

   461        21      13474.0    641.6      0.2          if self.verbose: # Print progress

   462        21      12025.0    572.6      0.2              t = self.t

   463        21     287019.0  13667.6      5.2              simlabel = f'"{self.label}": ' if self.label else ''

   464        21     916864.0  43660.2     16.7              string = f'  Running {simlabel}{t.now("str")} ({t.ti:2.0f}/{t.npts}) ({self.elapsed:0.2f} s) '

   465        21      18173.0    865.4      0.3              if self.verbose >= 1:

   466                                                           sc.heading(string)

   467        21      11900.0    566.7      0.2              elif self.verbose > 0:

   468        21      24488.0   1166.1      0.4                  if not (t.ti % int(1.0 / self.verbose)):

   469         3     426197.0 142065.7      7.7                      sc.progressbar(t.ti + 1, t.npts, label=string, length=20, newline=True)

   470        21      10672.0    508.2      0.2          return







—————————————————————————————————————————————————————————

Profile of modules.Module.start_step: 0.0126 s (5.24572%)

—————————————————————————————————————————————————————————



Total time: 0.0126 s

File: /home/runner/work/starsim/starsim/starsim/modules.py

Function: Module.start_step at line 723



Line #      Hits         Time  Per Hit   % Time  Line Contents

==============================================================

   723                                               @required()

   724                                               def start_step(self):

   725                                                   """ Tasks to perform at the beginning of the step """

   726        42      35111.0    836.0      0.3          if self.finalized:

   727                                                       errormsg = f'The module {self._debug_name} has already been run. Did you mean to copy it before running it?'

   728                                                       raise RuntimeError(errormsg)

   729        42      21599.0    514.3      0.2          if self.dists is not None: # Will be None if no distributions are defined

   730        42   12521547.0 298132.1     99.4              self.dists.jump_dt() # Advance random number generators forward for calls on this step

   731        42      21785.0    518.7      0.2          return







—————————————————————————————————————————————————————————

Profile of diseases.SIS.step_state: 0.017993 s (7.49097%)

—————————————————————————————————————————————————————————



Total time: 0.017993 s

File: /home/runner/work/starsim/starsim/starsim/diseases.py

Function: SIS.step_state at line 763



Line #      Hits         Time  Per Hit   % Time  Line Contents

==============================================================

   763                                               def step_state(self):

   764                                                   """ Progress infectious -> recovered """

   765        21    1823312.0  86824.4     10.1          recovered = (self.infected & (self.ti_recovered <= self.ti)).uids

   766        21     109612.0   5219.6      0.6          self.infected[recovered] = False

   767        21      67135.0   3196.9      0.4          self.susceptible[recovered] = True

   768        21   15983191.0 761104.3     88.8          self.update_immunity()

   769        21       9781.0    465.8      0.1          return







———————————————————————————————————————————————————————————————

Profile of networks.DynamicNetwork.step: 0.0552302 s (22.9938%)

———————————————————————————————————————————————————————————————



Total time: 0.0552302 s

File: /home/runner/work/starsim/starsim/starsim/networks.py

Function: DynamicNetwork.step at line 437



Line #      Hits         Time  Per Hit   % Time  Line Contents

==============================================================

   437                                               def step(self):

   438                                                   """ Remove expired partnerships and add new ones """

   439        21   30140525.0 1.44e+06     54.6          self.end_pairs()

   440        21   25078103.0 1.19e+06     45.4          self.add_pairs()

   441        21      11543.0    549.7      0.0          return







—————————————————————————————————————————————————————————

Profile of diseases.Infection.step: 0.111573 s (46.4509%)

—————————————————————————————————————————————————————————



Total time: 0.111573 s

File: /home/runner/work/starsim/starsim/starsim/diseases.py

Function: Infection.step at line 216



Line #      Hits         Time  Per Hit   % Time  Line Contents

==============================================================

   216                                               def step(self):

   217                                                   """

   218                                                   Perform key infection updates, including infection and setting prognoses

   219                                                   """

   220                                                   # Create new cases

   221        21   92493086.0  4.4e+06     82.9          new_cases, sources, networks = self.infect() # TODO: store outputs in self or use objdict rather than 3 returns

   222                                           

   223                                                   # Set prognoses

   224        21      17003.0    809.7      0.0          if len(new_cases):

   225        21   19051008.0 907190.9     17.1              self.set_outcomes(new_cases, sources)

   226                                           

   227        21      11743.0    559.2      0.0          return new_cases, sources, networks



Figure(672x480)

This graph (which is a shortcut to sim.loop.plot_cpu()) shows us how much time each step in the integration loop takes. We can get line-by-line detail of where each function is taking time, though:

prof.disp(maxentries=5)

——————————————————————————————————————————————————————

Profile of sim.Sim.start_step: 0.00550509 s (2.29192%)

——————————————————————————————————————————————————————



Total time: 0.00550509 s

File: /home/runner/work/starsim/starsim/starsim/sim.py

Function: Sim.start_step at line 451



Line #      Hits         Time  Per Hit   % Time  Line Contents

==============================================================

   451                                               def start_step(self):

   452                                                   """ Start the step -- only print progress; all actual changes happen in the modules """

   453                                           

   454                                                   # Set the time and if we have reached the end of the simulation, then do nothing

   455        21      10431.0    496.7      0.2          if self.complete:

   456                                                       errormsg = 'Simulation already complete (call sim.init() to re-run)'

   457                                                       raise AlreadyRunError(errormsg)

   458                                           

   459                                                   # Print out progress if needed

   460        21    3773842.0 179706.8     68.6          self.elapsed = self.timer.toc(output=True)

   461        21      13474.0    641.6      0.2          if self.verbose: # Print progress

   462        21      12025.0    572.6      0.2              t = self.t

   463        21     287019.0  13667.6      5.2              simlabel = f'"{self.label}": ' if self.label else ''

   464        21     916864.0  43660.2     16.7              string = f'  Running {simlabel}{t.now("str")} ({t.ti:2.0f}/{t.npts}) ({self.elapsed:0.2f} s) '

   465        21      18173.0    865.4      0.3              if self.verbose >= 1:

   466                                                           sc.heading(string)

   467        21      11900.0    566.7      0.2              elif self.verbose > 0:

   468        21      24488.0   1166.1      0.4                  if not (t.ti % int(1.0 / self.verbose)):

   469         3     426197.0 142065.7      7.7                      sc.progressbar(t.ti + 1, t.npts, label=string, length=20, newline=True)

   470        21      10672.0    508.2      0.2          return







—————————————————————————————————————————————————————————

Profile of modules.Module.start_step: 0.0126 s (5.24572%)

—————————————————————————————————————————————————————————



Total time: 0.0126 s

File: /home/runner/work/starsim/starsim/starsim/modules.py

Function: Module.start_step at line 723



Line #      Hits         Time  Per Hit   % Time  Line Contents

==============================================================

   723                                               @required()

   724                                               def start_step(self):

   725                                                   """ Tasks to perform at the beginning of the step """

   726        42      35111.0    836.0      0.3          if self.finalized:

   727                                                       errormsg = f'The module {self._debug_name} has already been run. Did you mean to copy it before running it?'

   728                                                       raise RuntimeError(errormsg)

   729        42      21599.0    514.3      0.2          if self.dists is not None: # Will be None if no distributions are defined

   730        42   12521547.0 298132.1     99.4              self.dists.jump_dt() # Advance random number generators forward for calls on this step

   731        42      21785.0    518.7      0.2          return







—————————————————————————————————————————————————————————

Profile of diseases.SIS.step_state: 0.017993 s (7.49097%)

—————————————————————————————————————————————————————————



Total time: 0.017993 s

File: /home/runner/work/starsim/starsim/starsim/diseases.py

Function: SIS.step_state at line 763



Line #      Hits         Time  Per Hit   % Time  Line Contents

==============================================================

   763                                               def step_state(self):

   764                                                   """ Progress infectious -> recovered """

   765        21    1823312.0  86824.4     10.1          recovered = (self.infected & (self.ti_recovered <= self.ti)).uids

   766        21     109612.0   5219.6      0.6          self.infected[recovered] = False

   767        21      67135.0   3196.9      0.4          self.susceptible[recovered] = True

   768        21   15983191.0 761104.3     88.8          self.update_immunity()

   769        21       9781.0    465.8      0.1          return







———————————————————————————————————————————————————————————————

Profile of networks.DynamicNetwork.step: 0.0552302 s (22.9938%)

———————————————————————————————————————————————————————————————



Total time: 0.0552302 s

File: /home/runner/work/starsim/starsim/starsim/networks.py

Function: DynamicNetwork.step at line 437



Line #      Hits         Time  Per Hit   % Time  Line Contents

==============================================================

   437                                               def step(self):

   438                                                   """ Remove expired partnerships and add new ones """

   439        21   30140525.0 1.44e+06     54.6          self.end_pairs()

   440        21   25078103.0 1.19e+06     45.4          self.add_pairs()

   441        21      11543.0    549.7      0.0          return







—————————————————————————————————————————————————————————

Profile of diseases.Infection.step: 0.111573 s (46.4509%)

—————————————————————————————————————————————————————————



Total time: 0.111573 s

File: /home/runner/work/starsim/starsim/starsim/diseases.py

Function: Infection.step at line 216



Line #      Hits         Time  Per Hit   % Time  Line Contents

==============================================================

   216                                               def step(self):

   217                                                   """

   218                                                   Perform key infection updates, including infection and setting prognoses

   219                                                   """

   220                                                   # Create new cases

   221        21   92493086.0  4.4e+06     82.9          new_cases, sources, networks = self.infect() # TODO: store outputs in self or use objdict rather than 3 returns

   222                                           

   223                                                   # Set prognoses

   224        21      17003.0    809.7      0.0          if len(new_cases):

   225        21   19051008.0 907190.9     17.1              self.set_outcomes(new_cases, sources)

   226                                           

   227        21      11743.0    559.2      0.0          return new_cases, sources, networks

(Note that the names of the functions here refer to the actual functions called, which may not match the graph above. That’s because, for example, ss.SIS does not define its own step() method, but instead inherits step() from Infection. In the graph, this is shown as sis.step(), but is listed in the table as Infection.step(). This is because it’s referring to the actual code being run, so refers to where those lines of code exist in the codebase; there is no code corresponding to SIS.step() since it’s just inherited from Infection.step().)

If you want more detail, you can also define custom functions to follow. For example, we can see that ss.SIS.infect() takes the most time in ss.SIS.step(), so let’s profile that:

prof = sim.profile(follow=ss.SIS.infect, plot=False)
prof.disp()

Initializing sim with 10000 agents

Profiling 1 function(s):

 <function Infection.infect at 0x7f8d78e1ca40> 




  Running 2000.01.01 ( 0/21) (0.00 s)  ———————————————————— 5%


  Running 2010.01.01 (10/21) (0.07 s)  ••••••••••—————————— 52%


  Running 2020.01.01 (20/21) (0.14 s)  •••••••••••••••••••• 100%



Elapsed time: 0.166 s





————————————————————————————————————————————————————————————

Profile of diseases.Infection.infect: 0.0647676 s (38.9776%)

————————————————————————————————————————————————————————————



Total time: 0.0647676 s

File: /home/runner/work/starsim/starsim/starsim/diseases.py

Function: Infection.infect at line 257



Line #      Hits         Time  Per Hit   % Time  Line Contents

==============================================================

   257                                               def infect(self):

   258                                                   """ Determine who gets infected on this timestep via transmission on the network """

   259        21      19657.0    936.0      0.0          new_cases = []

   260        21       8498.0    404.7      0.0          sources = []

   261        21       8409.0    400.4      0.0          networks = []

   262        21     591735.0  28177.9      0.9          betamap = self.validate_beta()

   263                                           

   264                                                   # Compute effective transmissibility and susceptibility directly on the raw

   265                                                   # (full-length) arrays. This avoids the gather/scatter, full-length astype copy,

   266                                                   # and extra wrapper allocations of the Arr math operators; edges only ever index

   267                                                   # living agents, so stale raw values for inactive agents are never used.

   268        21     432553.0  20597.8      0.7          rel_trans = self.rel_trans.asnew(self.infectious.raw * self.rel_trans.raw, copy=False)

   269        21     264078.0  12575.1      0.4          rel_sus   = self.rel_sus.asnew(self.susceptible.raw * self.rel_sus.raw, copy=False)

   270                                           

   271        42     151366.0   3604.0      0.2          for i, (nkey,route) in enumerate(self.sim.networks.items()):

   272        21      34617.0   1648.4      0.1              nk = ss.standardize_netkey(nkey)

   273                                           

   274                                                       # Main use case: networks

   275        21      14637.0    697.0      0.0              if isinstance(route, ss.Network):

   276        21     254847.0  12135.6      0.4                  if len(route): # Skip networks with no edges

   277        21       9786.0    466.0      0.0                      edges = route.edges

   278        21     357757.0  17036.0      0.6                      p1_to_p2 = [edges.p1, edges.p2, betamap[nk][0]]  # p1→p2 direction, beta 0

   279        21     343433.0  16354.0      0.5                      p2_to_p1 = [edges.p2, edges.p1, betamap[nk][1]]  # p2→p1 direction, beta 1

   280        63      42450.0    673.8      0.1                      for src, trg, beta in [p1_to_p2, p2_to_p1]:

   281        42      72054.0   1715.6      0.1                          if beta: # Skip networks with no transmission

   282        42    1157897.0  27569.0      1.8                              disease_beta = beta.to_prob(self.t.dt) if isinstance(beta, ss.Rate) else beta

   283        42    1907900.0  45426.2      2.9                              beta_per_dt = route.net_beta(disease_beta=disease_beta, disease=self) # Compute beta for this network and timestep

   284        42   46388479.0  1.1e+06     71.6                              randvals = self.trans_rng.rvs(src, trg) # Generate a new random number based on the two other random numbers

   285        42      37014.0    881.3      0.1                              args = (src, trg, rel_trans, rel_sus, beta_per_dt, randvals) # Set up the arguments to calculate transmission

   286        42    9961974.0 237189.9     15.4                              target_uids, source_uids = self.compute_transmission(*args) # Actually calculate it

   287        42      27183.0    647.2      0.0                              new_cases.append(target_uids)

   288        42      17189.0    409.3      0.0                              sources.append(source_uids)

   289        42     425045.0  10120.1      0.7                              networks.append(np.full(len(target_uids), dtype=ss_int, fill_value=i))

   290                                           

   291                                                       # Handle everything else: mixing pools, environmental transmission, etc.

   292                                                       elif isinstance(route, ss.Route):

   293                                                           # Mixing pools are unidirectional, only use the first beta value

   294                                                           disease_beta = betamap[nk][0].to_prob(self.t.dt) if isinstance(betamap[nk][0], ss.Rate) else betamap[nk][0]

   295                                                           target_uids = route.compute_transmission(rel_sus, rel_trans, disease_beta, disease=self)

   296                                                           new_cases.append(target_uids)

   297                                                           sources.append(np.full(len(target_uids), dtype=ss_int, fill_value=ss.dtypes.int_nan))

   298                                                           networks.append(np.full(len(target_uids), dtype=ss_int, fill_value=i))

   299                                                       else:

   300                                                           errormsg = f'Cannot compute transmission via route {type(route)}; please subclass ss.Route and define a compute_transmission() method'

   301                                                           raise TypeError(errormsg)

   302                                           

   303                                                   # Finalize

   304        21      15325.0    729.8      0.0          if len(new_cases) and len(sources):

   305        21     357042.0  17002.0      0.6              new_cases = ss.uids.concatenate(new_cases)

   306        21    1431352.0  68159.6      2.2              new_cases, inds = new_cases.unique(return_index=True)

   307        21     319734.0  15225.4      0.5              sources = ss.uids.concatenate(sources)[inds]

   308        21     103734.0   4939.7      0.2              networks = np.concatenate(networks)[inds]

   309                                                   else:

   310                                                       new_cases = ss.uids()

   311                                                       sources = ss.uids()

   312                                                       networks = np.empty(0, dtype=ss_int)

   313                                           

   314        21      11902.0    566.8      0.0          return new_cases, sources, networks







————————————————————————————————————————————————————————————

Profile of diseases.Infection.infect: 0.0647676 s (38.9776%)

————————————————————————————————————————————————————————————



Total time: 0.0647676 s

File: /home/runner/work/starsim/starsim/starsim/diseases.py

Function: Infection.infect at line 257



Line #      Hits         Time  Per Hit   % Time  Line Contents

==============================================================

   257                                               def infect(self):

   258                                                   """ Determine who gets infected on this timestep via transmission on the network """

   259        21      19657.0    936.0      0.0          new_cases = []

   260        21       8498.0    404.7      0.0          sources = []

   261        21       8409.0    400.4      0.0          networks = []

   262        21     591735.0  28177.9      0.9          betamap = self.validate_beta()

   263                                           

   264                                                   # Compute effective transmissibility and susceptibility directly on the raw

   265                                                   # (full-length) arrays. This avoids the gather/scatter, full-length astype copy,

   266                                                   # and extra wrapper allocations of the Arr math operators; edges only ever index

   267                                                   # living agents, so stale raw values for inactive agents are never used.

   268        21     432553.0  20597.8      0.7          rel_trans = self.rel_trans.asnew(self.infectious.raw * self.rel_trans.raw, copy=False)

   269        21     264078.0  12575.1      0.4          rel_sus   = self.rel_sus.asnew(self.susceptible.raw * self.rel_sus.raw, copy=False)

   270                                           

   271        42     151366.0   3604.0      0.2          for i, (nkey,route) in enumerate(self.sim.networks.items()):

   272        21      34617.0   1648.4      0.1              nk = ss.standardize_netkey(nkey)

   273                                           

   274                                                       # Main use case: networks

   275        21      14637.0    697.0      0.0              if isinstance(route, ss.Network):

   276        21     254847.0  12135.6      0.4                  if len(route): # Skip networks with no edges

   277        21       9786.0    466.0      0.0                      edges = route.edges

   278        21     357757.0  17036.0      0.6                      p1_to_p2 = [edges.p1, edges.p2, betamap[nk][0]]  # p1→p2 direction, beta 0

   279        21     343433.0  16354.0      0.5                      p2_to_p1 = [edges.p2, edges.p1, betamap[nk][1]]  # p2→p1 direction, beta 1

   280        63      42450.0    673.8      0.1                      for src, trg, beta in [p1_to_p2, p2_to_p1]:

   281        42      72054.0   1715.6      0.1                          if beta: # Skip networks with no transmission

   282        42    1157897.0  27569.0      1.8                              disease_beta = beta.to_prob(self.t.dt) if isinstance(beta, ss.Rate) else beta

   283        42    1907900.0  45426.2      2.9                              beta_per_dt = route.net_beta(disease_beta=disease_beta, disease=self) # Compute beta for this network and timestep

   284        42   46388479.0  1.1e+06     71.6                              randvals = self.trans_rng.rvs(src, trg) # Generate a new random number based on the two other random numbers

   285        42      37014.0    881.3      0.1                              args = (src, trg, rel_trans, rel_sus, beta_per_dt, randvals) # Set up the arguments to calculate transmission

   286        42    9961974.0 237189.9     15.4                              target_uids, source_uids = self.compute_transmission(*args) # Actually calculate it

   287        42      27183.0    647.2      0.0                              new_cases.append(target_uids)

   288        42      17189.0    409.3      0.0                              sources.append(source_uids)

   289        42     425045.0  10120.1      0.7                              networks.append(np.full(len(target_uids), dtype=ss_int, fill_value=i))

   290                                           

   291                                                       # Handle everything else: mixing pools, environmental transmission, etc.

   292                                                       elif isinstance(route, ss.Route):

   293                                                           # Mixing pools are unidirectional, only use the first beta value

   294                                                           disease_beta = betamap[nk][0].to_prob(self.t.dt) if isinstance(betamap[nk][0], ss.Rate) else betamap[nk][0]

   295                                                           target_uids = route.compute_transmission(rel_sus, rel_trans, disease_beta, disease=self)

   296                                                           new_cases.append(target_uids)

   297                                                           sources.append(np.full(len(target_uids), dtype=ss_int, fill_value=ss.dtypes.int_nan))

   298                                                           networks.append(np.full(len(target_uids), dtype=ss_int, fill_value=i))

   299                                                       else:

   300                                                           errormsg = f'Cannot compute transmission via route {type(route)}; please subclass ss.Route and define a compute_transmission() method'

   301                                                           raise TypeError(errormsg)

   302                                           

   303                                                   # Finalize

   304        21      15325.0    729.8      0.0          if len(new_cases) and len(sources):

   305        21     357042.0  17002.0      0.6              new_cases = ss.uids.concatenate(new_cases)

   306        21    1431352.0  68159.6      2.2              new_cases, inds = new_cases.unique(return_index=True)

   307        21     319734.0  15225.4      0.5              sources = ss.uids.concatenate(sources)[inds]

   308        21     103734.0   4939.7      0.2              networks = np.concatenate(networks)[inds]

   309                                                   else:

   310                                                       new_cases = ss.uids()

   311                                                       sources = ss.uids()

   312                                                       networks = np.empty(0, dtype=ss_int)

   313                                           

   314        21      11902.0    566.8      0.0          return new_cases, sources, networks

(Note: you can only follow functions that are called as part of sim.run() this way. To follow other functions, such as those run by sim.init(), you can use sc.profile() directly.)

Debugging

When figuring out what your sim is doing – whether it’s doing something it shouldn’t be, or not doing something it should – sim.loop is your friend. It shows everything that will happen in the sim, and in what order:

import starsim as ss

sim = ss.Sim(
    start = 2000,
    stop = 2002,
    diseases = 'sis',
    networks = 'random',
    verbose = 0,
)
sim.run()
sim.loop.df.disp()
# %%

    time  ti  func_order                     label     module       func_name    cpu_time
0   2000   0           0            sim.start_step        sim      start_step  9.1792e-05
1   2000   0           1      randomnet.start_step  randomnet      start_step  1.3736e-04
2   2000   0           2            sis.start_step        sis      start_step  2.3152e-04
3   2000   0           3            sis.step_state        sis      step_state  1.1375e-04
4   2000   0           4            randomnet.step  randomnet            step  1.0088e-03
5   2000   0           5                  sis.step        sis            step  3.2489e-03
6   2000   0           6           people.step_die     people        step_die  5.6526e-05
7   2000   0           7     people.update_results     people  update_results  8.9737e-05
8   2000   0           8  randomnet.update_results  randomnet  update_results  9.2780e-06
9   2000   0           9        sis.update_results        sis  update_results  9.0309e-05
10  2000   0          10     randomnet.finish_step  randomnet     finish_step  3.6370e-06
11  2000   0          11           sis.finish_step        sis     finish_step  1.1120e-06
12  2000   0          12        people.finish_step     people     finish_step  3.1929e-05
13  2000   0          13           sim.finish_step        sim     finish_step  1.8840e-06
14  2001   1           0            sim.start_step        sim      start_step  6.6625e-05
15  2001   1           1      randomnet.start_step  randomnet      start_step  1.5379e-04
16  2001   1           2            sis.start_step        sis      start_step  2.0448e-04
17  2001   1           3            sis.step_state        sis      step_state  9.7101e-05
18  2001   1           4            randomnet.step  randomnet            step  9.9110e-04
19  2001   1           5                  sis.step        sis            step  2.5808e-03
20  2001   1           6           people.step_die     people        step_die  5.0926e-05
21  2001   1           7     people.update_results     people  update_results  1.0437e-04
22  2001   1           8  randomnet.update_results  randomnet  update_results  7.4140e-06
23  2001   1           9        sis.update_results        sis  update_results  7.6052e-05
24  2001   1          10     randomnet.finish_step  randomnet     finish_step  2.7750e-06
25  2001   1          11           sis.finish_step        sis     finish_step  1.2420e-06
26  2001   1          12        people.finish_step     people     finish_step  2.7462e-05
27  2001   1          13           sim.finish_step        sim     finish_step  1.0520e-06
28  2002   2           0            sim.start_step        sim      start_step  6.5041e-05
29  2002   2           1      randomnet.start_step  randomnet      start_step  1.3303e-04
30  2002   2           2            sis.start_step        sis      start_step  1.9686e-04
31  2002   2           3            sis.step_state        sis      step_state  9.6109e-05
32  2002   2           4            randomnet.step  randomnet            step  9.7613e-04
33  2002   2           5                  sis.step        sis            step  2.8875e-03
34  2002   2           6           people.step_die     people        step_die  5.1045e-05
35  2002   2           7     people.update_results     people  update_results  7.7325e-05
36  2002   2           8  randomnet.update_results  randomnet  update_results  6.8730e-06
37  2002   2           9        sis.update_results        sis  update_results  7.3598e-05
38  2002   2          10     randomnet.finish_step  randomnet     finish_step  2.3540e-06
39  2002   2          11           sis.finish_step        sis     finish_step  1.0920e-06
40  2002   2          12        people.finish_step     people     finish_step  2.7321e-05
41  2002   2          13           sim.finish_step        sim     finish_step  7.8200e-07

As you can see, it’s a lot – this is only three timesteps and two modules, and it’s already 41 steps.

The typical way to do debugging is to insert breakpoints or print statements into your modules for custom debugging (e.g., to print a value), or to use analyzers for heavier-lift debugging. Starsim also lets you manually modify the loop by inserting “probes” or other arbitrary functions. For example, if you wanted to check the population size after each time the People object is updated:

def check_pop_size(sim):
    print(f'Population size is {len(sim.people)}')

sim = ss.Sim(diseases='sir', networks='random', demographics=True, dur=10)
sim.init()
sim.loop.insert(check_pop_size, label='people.finish_step')
sim.run()

Initializing sim with 10000 agents

  Running 2000 ( 0/11) (0.00 s)  •——————————————————— 9%
Population size is 10074
Population size is 10191
Population size is 10281
Population size is 10376
Population size is 10454
Population size is 10556
Population size is 10652
Population size is 10742
Population size is 10853
Population size is 10946

  Running 2010 (10/11) (0.07 s)  •••••••••••••••••••• 100%

Population size is 11041

Sim(n=10000; 2000—2010.0; demographics=births, deaths; networks=randomnet; diseases=sir)

In this case, you get the same output as using an analyzer:

def check_pop_size(sim):
    print(f'Population size is {len(sim.people)}')

sim = ss.Sim(diseases='sir', networks='random', demographics=True, dur=10, analyzers=check_pop_size)
sim.run()

Initializing sim with 10000 agents

  Running 2000 ( 0/11) (0.00 s)  •——————————————————— 9%
Population size is 10179
Population size is 10287
Population size is 10379
Population size is 10485
Population size is 10559
Population size is 10665
Population size is 10750
Population size is 10857
Population size is 10951
Population size is 11066

  Running 2010 (10/11) (0.07 s)  •••••••••••••••••••• 100%

Population size is 11163

Sim(n=10000; 2000—2010.0; demographics=births, deaths; networks=randomnet; diseases=sir; analyzers=check_pop_size)

However, inserting functions directly in the loop gives you more control over their exact placement, whereas analyzers are always executed last in the timestep.

The loop also has methods for visualizing itself. You can get a simple representation of the loop with loop.plot():

sim.loop.plot()

Figure(672x480)

Or a slightly more detailed one with loop.plot_step_order():

sim.loop.plot_step_order()

Figure(672x480)

This is especially useful if your simulation has modules with different timesteps, e.g.:

sis = ss.SIS(dt=0.1)
net = ss.RandomNet(dt=0.5)
births = ss.Births(dt=1)
sim = ss.Sim(dt=0.1, dur=5, diseases=sis, networks=net, demographics=births)
sim.init()
sim.loop.plot_step_order()

Initializing sim with 10000 agents
Figure(672x480)

(Note: this is a 3D plot, so it helps if you can plot it in a separate window interactively to be able to move it around, rather than just in a notebook.)