Commit a72dc4ec authored by Tiago Peixoto's avatar Tiago Peixoto

dynamics: fix behavior of epidemic models with constant_beta=True

parent 7ddc03f0
...@@ -129,7 +129,6 @@ public: ...@@ -129,7 +129,6 @@ public:
for (auto e : edges_range(g)) for (auto e : edges_range(g))
beta[e] = std::log1p(-_beta[e]); beta[e] = std::log1p(-_beta[e]);
_beta = beta; _beta = beta;
}
for (auto v : vertices_range(g)) for (auto v : vertices_range(g))
{ {
...@@ -142,6 +141,7 @@ public: ...@@ -142,6 +141,7 @@ public:
_m_temp[v] = _m[v]; _m_temp[v] = _m[v];
} }
} }
}
}; };
template <class Graph> template <class Graph>
...@@ -182,6 +182,9 @@ public: ...@@ -182,6 +182,9 @@ public:
} }
else else
{ {
if constexpr (!constant_beta)
return;
if constexpr (sync) if constexpr (sync)
{ {
for (auto e : out_edges_range(v, g)) for (auto e : out_edges_range(v, g))
...@@ -233,13 +236,26 @@ public: ...@@ -233,13 +236,26 @@ public:
return 1; return 1;
} }
auto m = _m[v];
double prob = 0; double prob = 0;
if constexpr (!weighted || constant_beta)
{
auto m = _m[v];
if constexpr (!weighted) if constexpr (!weighted)
prob = _prob[m]; prob = _prob[m];
else else
prob = 1 - std::exp(m); prob = 1 - std::exp(m);
}
else
{
for (auto e : in_or_out_edges_range(v, g))
{
auto w = source(e, g);
if (_s[w] == State::I)
prob += std::log1p(-_beta[e]);
}
prob = 1 - std::exp(prob);
}
std::bernoulli_distribution minfect(prob); std::bernoulli_distribution minfect(prob);
if (prob > 0 && minfect(rng)) if (prob > 0 && minfect(rng))
......
...@@ -162,8 +162,8 @@ class EpidemicStateBase(DiscreteStateBase): ...@@ -162,8 +162,8 @@ class EpidemicStateBase(DiscreteStateBase):
if weighted: if weighted:
_check_prop_scalar(beta, "beta") _check_prop_scalar(beta, "beta")
if beta.value_type() != "double": if beta.value_type() != "double":
if constant_beta: if not constant_beta:
raise ValueError("if constant_beta == True, the type of beta must be double") raise ValueError("if constant_beta == False, the type of beta must be double")
beta = beta.copy("double") beta = beta.copy("double")
params["beta"] = beta params["beta"] = beta
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment