In [1]:
import scotch
import matplotlib.pyplot as plt
%matplotlib inline
import numpy as np

Things to do

set json files for each version of the model -- just change the parameter values

get actors and states dicts

figure out how to translate them into the right file patterns


In [2]:
model = scotch.model('HCVLinTrackSS.json')


scotch.scotch.model with 9 states, 13 parameters, and 25 events.

States :
T, initial condition = 10000
E1, initial condition = 0
E1X, initial condition = 0
E2, initial condition = 0
E2X, initial condition = 0
I1, initial condition = 10
I2, initial condition = 1
DT, initial condition = 0
DI, initial condition = 0

Parameters :
kappa = 0
phi_T = 0.14
beta_V = 0.000000005
nu_I = 1/7
phi_I = 0.14 * 0.8
mu = 0.00001
khi = 8.18 * 100000000000 / 10000* 4.1825 / (22.3*1500)
s = 1.01
eta = 0.01
nu_T = 0.014
beta_I = 4.1825 * 0.00001 / (60*24)
alpha = 1
alpha_X = .013

Events :
Number 0: 	probability (1 - eta) * beta_I * I1 * T * (1 - mu)
Number 1: 	probability (1 - eta) * beta_V * khi * I1 * T * (1 - mu)
Number 2: 	probability eta * beta_I * I1 * T * (1 - mu)
Number 3: 	probability eta * beta_V * I1 * khi * T * (1 - mu)
Number 4: 	probability mu * (1 - eta) * beta_I * I1 * T
Number 5: 	probability mu * (1 - eta) * beta_V * khi * I1 * T
Number 6: 	probability (1 - eta) * s * beta_I * I2 * T
Number 7: 	probability (1 - eta) * s * beta_V * khi * I2 * T
Number 8: 	probability mu * eta * beta_I * I1 * T
Number 9: 	probability mu * eta * beta_V * khi * I1 * T
Number 10: 	probability eta * s * beta_I * I2 * T
Number 11: 	probability eta * s * beta_V * khi * I2 * T
Number 12: 	probability nu_T * T
Number 13: 	probability alpha * E1
Number 14: 	probability alpha_X * E1X
Number 15: 	probability alpha * E2
Number 16: 	probability alpha_X * E2
Number 17: 	probability nu_I * I1
Number 18: 	probability nu_I * I2
Number 19: 	probability nu_T * E1
Number 20: 	probability nu_T * E2
Number 21: 	probability nu_T * E1X
Number 22: 	probability nu_T * E2X
Number 23: 	probability phi_T * DT
Number 24: 	probability phi_I * DI


In [3]:
scotch.helpers.add_actors_wizard(model)


This wizard will help assign states as actors for each events
The states to choose from are ['T', 'E1', 'E1X', 'E2', 'E2X', 'I1', 'I2', 'DT', 'DI']
Enter actor for the following event (enter '[]' without quotes if no actor)
(1 - eta) * beta_I * I1 * T * (1 - mu)
I1
Enter actor for the following event (enter '[]' without quotes if no actor)
(1 - eta) * beta_V * khi * I1 * T * (1 - mu)
I1
Enter actor for the following event (enter '[]' without quotes if no actor)
eta * beta_I * I1 * T * (1 - mu)
I1
Enter actor for the following event (enter '[]' without quotes if no actor)
eta * beta_V * I1 * khi * T * (1 - mu)
I1
Enter actor for the following event (enter '[]' without quotes if no actor)
mu * (1 - eta) * beta_I * I1 * T
I1
Enter actor for the following event (enter '[]' without quotes if no actor)
mu * (1 - eta) * beta_V * khi * I1 * T
I1
Enter actor for the following event (enter '[]' without quotes if no actor)
(1 - eta) * s * beta_I * I2 * T
I2
Enter actor for the following event (enter '[]' without quotes if no actor)
(1 - eta) * s * beta_V * khi * I2 * T
I2
Enter actor for the following event (enter '[]' without quotes if no actor)
mu * eta * beta_I * I1 * T
I1
Enter actor for the following event (enter '[]' without quotes if no actor)
mu * eta * beta_V * khi * I1 * T
I1
Enter actor for the following event (enter '[]' without quotes if no actor)
eta * s * beta_I * I2 * T
I2
Enter actor for the following event (enter '[]' without quotes if no actor)
eta * s * beta_V * khi * I2 * T
I2
Enter actor for the following event (enter '[]' without quotes if no actor)
nu_T * T
[]
Enter actor for the following event (enter '[]' without quotes if no actor)
alpha * E1
[]
Enter actor for the following event (enter '[]' without quotes if no actor)
alpha_X * E1X
[]
Enter actor for the following event (enter '[]' without quotes if no actor)
alpha * E2
[]
Enter actor for the following event (enter '[]' without quotes if no actor)
alpha_X * E2
[]
Enter actor for the following event (enter '[]' without quotes if no actor)
nu_I * I1
[]
Enter actor for the following event (enter '[]' without quotes if no actor)
nu_I * I2
[]
Enter actor for the following event (enter '[]' without quotes if no actor)
nu_T * E1
[]
Enter actor for the following event (enter '[]' without quotes if no actor)
nu_T * E2
[]
Enter actor for the following event (enter '[]' without quotes if no actor)
nu_T * E1X
[]
Enter actor for the following event (enter '[]' without quotes if no actor)
nu_T * E2X
[]
Enter actor for the following event (enter '[]' without quotes if no actor)
phi_T * DT
[]
Enter actor for the following event (enter '[]' without quotes if no actor)
phi_I * DI
[]

In [4]:
model.save('HCVLinTrackSS.json')

In [3]:
[t,trace, tracked_events] =scotch.simulate.tauLeap(model,60,1, True)


[                         ]

In [4]:
plt.plot(t,trace)
plt.legend(model.states)


Out[4]:
<matplotlib.legend.Legend at 0x105dcb450>

In [5]:
trace[:,0]


Out[5]:
array([ 10000.,   9869.,   9741.,   9611.,   9516.,   9461.,   9401.,
         9337.,   9267.,   9205.,   9135.,   9073.,   9035.,   9011.,
         8947.,   8884.,   8818.,   8712.,   8626.,   8531.,   8394.,
         8217.,   7993.,   7771.,   7535.,   7239.,   6897.,   6539.,
         6156.,   5739.,   5298.,   4837.,   4429.,   4041.,   3736.,
         3416.,   3152.,   2976.,   2787.,   2638.,   2526.,   2514.,
         2524.,   2457.,   2456.,   2421.,   2432.,   2454.,   2425.,
         2463.,   2513.,   2559.,   2630.,   2688.,   2753.,   2770.,
         2789.,   2834.,   2886.,   2907.])

In [6]:
[statesDict, actorsDict] = scotch.helpers.trackIndividuals(model,tracked_events,t, True,True)

In [17]:
for x in statesDict['I1']:
    x.sort()
    (key=lambda tup: tup[1])

In [31]:
actorsDict['(1 - eta) * beta_V * khi * I1 * T * (1 - mu)'][0].sort(key = lambda tup: (tup[0],tup[1]))

In [7]:
actorsDict['(1 - eta) * beta_V * khi * I1 * T * (1 - mu)'][0]


Out[7]:
[(0, 6108),
 (1, 8423),
 (2, 21),
 (2, 8808),
 (4, 2133),
 (4, 5684),
 (8, 4367),
 (8, 6691)]

In [26]:
a = sorted(actorsDict['(1 - eta) * beta_V * khi * I1 * T * (1 - mu)'][0], key = lambda tup: tup[0])
x = [b for c in a for b in c]

In [30]:
print(x)
print(a)


[0, 4944, 3, 6943, 3, 1720, 6, 2705, 8, 8047]
[(0, 4944), (3, 6943), (3, 1720), (6, 2705), (8, 8047)]

In [20]:
trace[:,5]


Out[20]:
array([   10.,     9.,    11.,    13.,    17.,    22.,    33.,    34.,
          46.,    53.,    63.,    73.,   101.,   120.,   147.,   179.,
         222.,   258.,   315.,   365.,   432.,   530.,   614.,   763.,
         857.,  1025.,  1185.,  1375.,  1544.,  1815.,  2062.,  2333.,
        2517.,  2722.,  2902.,  3004.,  3105.,  3181.,  3142.,  3152.,
        3078.,  3049.,  3006.,  2888.,  2860.,  2836.,  2774.,  2707.,
        2696.,  2661.,  2662.,  2654.,  2640.,  2634.,  2625.,  2588.,
        2587.,  2584.,  2588.,  2613.])

In [10]:
[x[50] for x in actorsDict]


---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-10-20f82ef4d436> in <module>()
----> 1 [x[50] for x in actorsDict]

IndexError: string index out of range

In [21]:
[x[0] for x in [actorsDict[x[0]] for x in model.events]]


Out[21]:
[[],
 [(0, 6108),
  (1, 8423),
  (2, 21),
  (2, 8808),
  (4, 2133),
  (4, 5684),
  (8, 4367),
  (8, 6691)],
 [],
 [],
 [],
 [],
 [],
 [],
 [],
 [],
 [],
 [],
 [],
 [],
 [],
 [],
 [],
 [],
 [],
 [],
 [],
 [],
 [],
 [],
 []]

In [ ]: