# Code Review
###### tags: `others`

## Simulation

### Original

Maintain $I\cap g$ for all contact group $g$

For each $v\in S\cap g$, the infected prob is

$$
\begin{split}
&1-\prod_{u\in I\cap g}(1-c_{age(u),age(v)}\cdot p_{age(v)})\\
=&1-\prod_{i\in age}(1-c_{i,age(v)}\cdot p_{age(v)})^{n_i}=f(age(v))
\end{split}
$$

where $n_i$ is the number of infectious nodes of age $i$.

Thus the number of nodes infected of age $j$ is $\sim \mbox{Binomial}(m_j,f(j))$ where $m_j$ is the number of susceptible nodes of age $j$.

inline std::vector<uint> infection(char src, char des, uint period, const std::vector<double>& ptrans) {
        std::vector<uint> re;
#pragma omp parallel num_threads(16)
        {
            std::vector<uint> tre;
#pragma omp for schedule(dynamic, 128) nowait
            for (uint i = 0; i < N_gp; ++i) {
                if (!icnt_all[i] || cgpp[i].period != period) continue;
                for (uint v = 0; v < N_ag; ++v) {
                    auto& aggp = cgpp[i].nds[v];
                    if (aggp.size() == 0) continue;

                    double m = 1;
                    for (uint u = 0; u < N_ag; ++u) {
                        m *= std::pow(1 - cgpp[i].cm.getRate(u, v) * ptrans[v], icnt[i][u]);
                    }
                    int k = Random::bino_dis(aggp.size(), 1 - m);
                    for (auto c : Random::choose(aggp.size(), k)) {
                        if (__sync_bool_compare_and_swap(&ndp[aggp[c]].stateID, src, des))
                            tre.push_back(aggp[c]);
                    }
                }
            }
            
#pragma omp critical
            {
                re.insert(re.end(), tre.begin(), tre.end());
            }
        }
        return re;
    }

std::vector<uint> s2e = infection('S', 'E', ts.getPeriod(), prob_transmission); // chged
std::vector<uint> f2e = infection('F', 'E', ts.getPeriod(), prob_transmission); // chged

### New

#### New Feature

- two dose of vaccination $S, V_1, V_2,F$
- vaccination efficiency wanning $\epsilon_{d}$
- $I_{pre}, I_{asym}, I_{sym}$
- relative infectiuousness $\tau$

#### Infection Algorithm
- For each $u\in I$
- For each $g\ni u$
- $k=\lceil|g|\cdot \max_{i,j}c_{i,j}\rceil$
- choose $k$ nodes among $g$
- for each v in these $k$ nodes, infect it if
- it is of state $S$ ot $V_1$ or $V_2$
- with prob $p(u,v)\cdot|g|/k$
- where $p(u,v)=(1-\epsilon_{v})\cdot c_{age(u),age(v)}\cdot p_{age(v)}\cdot\tau_{state(u)}$

inline void infection(ExpiringState& ext, double tau, const Time::TimeStep& ts, Nodes& s2e, Nodes& v2e, Nodes& w2e, Nodes& f2e) {
        for (auto u : ext) {
            for (auto& cgp : ndp[u].gp[ts.getPeriod()]) {
                /*
                std::vector<uint> re_s2e, re_v2e, re_w2e;
                */

                uint n = cgp.size();
                double k = ceil(n * cgp.getContactMatrix().getPmax());
                for (auto idx : Random::choose(n, k)) {
                    uint v = cgp.at(idx);
                    uint dt = ts - ndp[v].ts;
                    double p = epsilon_bar(ndp[v].stateID, dt) * cgp.getContactMatrix().getRate(ndp[u].age, ndp[v].age) * tau;
                    if (Random::trial(p * n / k)) {
                        //
                        switch (ndp[v].stateID) {
                            case 'S':
                                s2e.push_back(v);
                                ndp[v].stateID = 'E';
                                break;
                            case 'V':
                                v2e.push_back(v);
                                ndp[v].stateID = 'E';
                                break;
                            case 'W':
                                w2e.push_back(v);
                                ndp[v].stateID = 'E';
                                break;
                            case 'F':
                                f2e.push_back(v);
                                ndp[v].stateID = 'F';
                                break;
                        }
                        //

                        /*
                        if (__sync_bool_compare_and_swap(&ndp[v].stateID, 'S', 'E'))
                            tre_s2e.push_back(aggp[c]);
                        else if (__sync_bool_compare_and_swap(&ndp[v].stateID, 'V', 'E'))
                            tre_v2e.push_back(aggp[c]);
                        else if (__sync_bool_compare_and_swap(&ndp[aggp[c]].stateID, 'W', 'E'))
                            tre_w2e.push_back(aggp[c]);
                        else if (__sync_bool_compare_and_swap(&ndp[aggp[c]].stateID, 'F', 'E'))
                            tre_f2e.push_back(aggp[c]);
                        */
                    }
                }

                /*
                {
                    s2e.insert(s2e.end(), tre_s2e.begin(), tre_s2e.end());
                    v2e.insert(v2e.end(), tre_v2e.begin(), tre_v2e.end());
                    w2e.insert(w2e.end(), tre_w2e.begin(), tre_w2e.end());
                }
                */
            }
        }
    }

Nodes s2e, v2e, w2e, f2e;
infection(I_pre, tau_I_pre, ts, s2e, v2e, w2e, f2e)
infection(I_asym, tau_I_asym, ts, s2e, v2e, w2e, f2e)
infection(I_sym, tau_I_sym, ts, s2e, v2e, w2e, f2e)