notes.txt

Runtime modes:

main --+---+--> one --+--> player-type
       |   |
       |   +--> grp --+--> player-counts
       |   |
       |   +--> pop --+--> player-dist
       |
       +--> env [one/grp]



Notes on Planner:

    Last_week.count("W")

    28, 20-27

    1.  There are two important concepts for any person of the Planner population:
        a. how is the work history (X) calculated
        b. what is the distribution of P[work|X]

    2.  For work history calculations to be reasonable, we have to attach exponentially decreasing
        weights to past days. (with a const exp factor: h)
            => X_it = W_i{t-1} + h.X_i{t-1}     for one person
            => X_t = eta_w_{t-1} + h.X_{t-1}    for a population

    3.  To get the macro results (eta_w_t in terms of some aggregate(s) of X_it), compatible with
        micro analysis (P[work]_it for each X_it), we need a very specific distribution of P[work|X]

    4.  If we plan to use X = avg(X_i) as our aggregate, P[work|X] must be linear in X (in [0, 1/1-h))

    5.  Fixing the upper and lower bounds of P[work|X] as 0, 1 doesn't work. This is because of -
        a.  with fixed bounds, all information must come from h, which makes it time-variant. With
            a time-variant h, the formula of 2 doesn't hold. (Imagine h was 0.9 upto t-1, and guy
            keeping track of X in a time variant h requires computing h from the entire history,
            removing any benefits of aggregation, and becoming impossible for part 3)

        b.  with fixed bounds of 0, 1: we have P[work] = 1 - (1 - h)X, we know for symptomatic guys
            P[work] = 0 for all X, but no h can capture this.

    6.  Thus we must have: eta_w = mX + c where m and c are time variant, m <= 0, and the bounds of
        eta_w, ie. [c + m/1-h, c] must be a subset of [0, 1]. Instead of working with m, c, we work
        with the (time variant) bounds themselves, eta_w_min, eta_w_max, with
            eta_w_t = (eta_w_min - eta_w_max) * (1 - h) * X_t + eta_w_max
                where, eta_w_min/max are a function of params and current environment
            also, X_t = eta_w_{t-1} + h.X_{t-1}
        Tracking these two (X and eta_w)_t gives the entire solution.

    7.  Since we are storing the information in the bounds, adding parametric information to h is,
        at best redundant and at worst incorrect. Much better to keep it a constant at data.json.

    8.  The last question that remains is what should the bounds be? The probability that simple
        takes must lie between the bounds. The following seems reasonable -

        Simple   Planner(Max)   Planner(Min)
          95%       97%             93%
          90%       94%             86%
          50%       70%             30%
          10%       14%              6%
           5%        7%              3%


    *.  An important implementation detail is how should X vary in execute_infection()?
        X is U[0 , 1/1-h] -> P is U[a, b] (0 => b, 1/1-h => a)
        X_mean  =  y
        if a < 1/2(1-h), then assume U[0, 2 * a] else assume U[ 2 * a - 1/(1-h), 1/1-h]
        c = (p3/p1)*a      b = (1 + p3/p1)*a

        <X>_S,t-1 known = a     P_S,t-1 = p1
        <X>_S, t = ? (b)        P_S, t  = p2
        <X>_I, t = ? (c)        P_I, t  = p3

        <X>_W , <X>_H
        eta_iw , eta_ih

        <X>_I = eta_iw * <X>_W + eta_ih * <X>_H

        <X>_S = 1-eta_iw * <X>_W + 1-eta_ih * <X>_H


        p1 = p2 + p3, all three are known

        a*p1 = b*p2 + c*p3







        # assuming self.eta_w stores previous day's working ratio
        # Since planner type tries to go to work for a number of target days in a week
        #
        # some way to know avg people who went to work in susceptible population?
        #
        # 0: n0, 1: n1, ...
        # p0*n0 + p1*n1 + ... + p7*n7 (p4 = raghav) p0 > p1 > ... > p7
        #         t     t+1
        # p0=1    n0 -> n0=0
        # p1=1    n1 -> n1=n0
        # p2=1    n2 -> n2=n1
        # p3=0.8  n3 -> n3=n2 + 0.2*n3
        # p4=0.2  n4 -> n4=n3
        # p5=0    n5 ->
        # p6=0    n6 ->
        # p7=0    n7 ->
        #
        # p0..7, n0..7 -> n0..7
        #
        # P[worked exactly 7 days ago] = x
        #
        # n0 = n0*(1-p0) + n1*(1-p1)*x
        # n1 = n0*p0 + n1*(1-p1)*(1 - x) + n1*p1*x + n2*(1-p2)*x
        # n2 = n1*p1 + n2*(1-p2)*(1 - x) + n2*p2*x + n3*(1-p3)*x
        # ...
        # n6 = n5*p5 + n6*(1-p6)*(1 - x) + n6*p6*x + n7*(1-p7)*x
        # n7 = n6*p6 + n7*(1-p7)*(1 - x) + n7*p7*x
        # ________________________________________________________________
        #
        # p = f(x)    -> x is a measure of 'work history'
        # x = [1 if worked yesterday] + [6/7 if worked day before] + ...
        # x(@t+1) = x - [1/7 * n_days worked last week] + [1 if worked yesterday]
        #
        # x = [1 if worked yesterday] + [h if worked day before] + [h^2 ...] + ...
        #
        #
        # x = hx + [1 if worked yesterday] : for one person
        #
        # W.W.H.H.W.W.H
        # x=? (p=0.5)
        # 0|1|1.5|0.75|0.375|1.1875|1.59375|0.796875
        #
        # 1 / (1-h)
        #
        # 1, 2, ..., N=10^10
        # x1, x2, ..., xN
        #
        # E[no. of people working | t=t] = \sum f(xi)
        # E[no. of people working | t=t+1] = \sum f(xi*h + [1 if i worked yesterday])
        #         IF f is linear, f(x) = 1 - (1-h)*x
        #                                  = \sum f(xi)*h + \sum f([1 if i worked yesterday])
        #
        # eta_w = \sum f(xi)*h + \sum f([1 if i worked yesterday else 0])
        # eta_w = eta_w*h + eta_w*f(1) + (1-eta_w)*f(0)
        #       = eta_w*h + eta_w*h + (1-eta_w)
        #       = 1 - eta_w*(1 - 2h)
        #       this works if h < 1/2
        #       otherwise eta_w > 1
        #
        #
        # where p = f(x) in [0,1]
        #
        # wt = 1 - n/7
        # eta_w = N(1-h) - (1-2h)eta_w





        ##########################################################################################################
        #
        # ratio_over_threshold_w = 1 - ((threshold_sw - (1 -(0)*num / unaware_days -  self.fluctuation / 2)) / self.fluctuation)
        # ratio_over_threshold_h = 1 - ((threshold_sh - (1 - self.fluctuation / 2)) / self.fluctuation)
        #
        # w.append(last_w[0] * ratio_over_threshold_w + (1 - last_w[0]) * ratio_over_threshold_h)
        #
        #
        # for i in range(total_days - 1):
        #     if i < unaware_days:
        #         ratio_over_threshold_w = 1 - ((threshold_sw - (1 - (i + 1) / unaware_days - self.fluctuation / 2)) / self.fluctuation)
        #         ratio_over_threshold_w = min(1.0, max(ratio_over_threshold_w, 0))
        #         ratio_over_threshold_h = 1 - (
        #                 (threshold_sh - (1 - (i + 1) / unaware_days - self.fluctuation / 2)) / self.fluctuation)
        #         ratio_over_threshold_h = min(1.0, max(ratio_over_threshold_h, 0))
        #     else:
        #         ratio_over_threshold_w = max((self.fluctuation / 2 - threshold_sw) / self.fluctuation, 0)
        #         ratio_over_threshold_h = max((self.fluctuation / 2 - threshold_sh) / self.fluctuation, 0)
        #     w.append(last_w[i] * ratio_over_threshold_w + (1 - last_w[i]) * ratio_over_threshold_h)

        # ratio_over_threshold_w = 1 - ((threshold_sw - (1 - self.fluctuation / 2)) / self.fluctuation)
        # ratio_over_threshold_h = 1 - ((threshold_sh - (1 - self.fluctuation / 2)) / self.fluctuation)

        # oye @CC ye direct 1 nhi hai??
        # w.append(last_w[total_days] * ratio_over_threshold_w + (1 - last_w[total_days]) * ratio_over_threshold_h)
        # w.append(1)

t = 0

for all types: ratio(W)

t = 1

for all types: ratio(W) ??


TODO BUGS:
    1. sick people don't work in groups, so why raise the chance of infection?

































# info on data.json
    > danger is sqrt death_risk (adjust constant in base player)





# type_shr: periodic?
if len(self.last_week_actions) >= 7:
    self.last_week_actions.pop(0)

assert len(self.action_plan) == 0

cash_work = self.u_economic_w
virus_utility = -self._params["danger"]
death_risk = self.death_risk
active_infection_risk = self.work_infection_risk
passive_infection_risk = self.home_infection_risk
death_utility = self.u_death
health_belief = self.p_healthy
caution_multiplier = 100

# Strategy
# If he hasn't gone to work 4 days in last week, must go to work
# If has gone to work, then compares utility of work and home
# Since he works minimum 4 days a week, is extra cautious about other 3

work = cash_work + virus_utility * active_infection_risk * caution_multiplier + virus_utility * (
        1 - health_belief)

if self.last_week_actions.count("W") < 4:
    action = "W"
    logger.debug("Hasn't gone to work for 4 days in last week, so "
                 "choosing {0}".format(action))
elif work > 0:
    action = "W"
    logger.debug(
            "Working as perceived economic payoff is: {0:.3f}".format(
                    work)
    )
else:
    action = "H"

self.action_plan.append(action)
self.last_week_actions.append(action)


# type ri: linear
class TypeRi(_Person):

    def plan(self):
        if len(self.action_plan) != 0:
            return

        plan_days = 15

        virus_contact_prob_h = 1 - np.power(1 - self.home_infection_risk,
                                            np.arange(plan_days))
        logger.debug("virus contact prob_h = " + str(virus_contact_prob_h))
        virus_contact_prob_w = 1 - np.power(
                1 - (self.home_infection_risk + self.work_infection_risk),
                np.arange(plan_days))
        logger.debug("virus contact prob_w = " + str(virus_contact_prob_w))
        virus_contact_prob_delta = virus_contact_prob_w - virus_contact_prob_h
        logger.debug(
            "virus contact prob_delta = " + str(virus_contact_prob_delta))
        virus_util = virus_contact_prob_delta * self.death_risk * self.u_death
        logger.debug("virus util = " + str(virus_util))
        economic_util = self.u_economic_w * np.arange(plan_days)
        logger.debug("economic util = " + str(economic_util))
        net_util = virus_util + economic_util
        logger.debug(("net util = " + str(net_util)))
        work_days = int(np.argmax(net_util))
        logger.debug("work days = " + str(work_days))

        for i in range(work_days):
            self.action_plan.append("W")
        for i in range(plan_days - work_days):
            self.action_plan.append("H")

## estimating_h.py
import numpy as np

target_days = 4


u_economic_w = 1.2
u_death = 20000
death_prob = 1 - 0.95644

coeff_i = 0.05
coeff_wi =  0.2

total_infectious = 0.04991064123654763
working_infectious = 0.005820863793475713

p_h_max = 23/24

eta_ih = total_infectious * coeff_i
eta_iw = np.clip(1 - (1 - total_infectious * coeff_i) * (1 - working_infectious * coeff_wi), eta_ih, 1)

cutoff_p_h = np.sqrt(u_death * death_prob * (eta_iw - eta_ih) / u_economic_w)


_eta_w = 0.16055732
assert -0.01 < (p_h_max - np.sqrt(u_death * death_prob * (eta_iw - eta_ih) / u_economic_w)) / 0.25 - _eta_w < 0.01

h = 1 - 1/(2*((p_h_max - np.sqrt(u_death * death_prob * (eta_iw - eta_ih) / u_economic_w)) / 0.25))


weights = np.logspace(0,6,num=7,base = h)

loops = 10000
x = 0

for _ in range(loops):
    perm = np.arange(7)
    perm = np.random.permutation(perm)
    perm = perm[:target_days]

    arr = np.zeros(7)
    for idx in perm:
        arr[idx] = 1

    x += np.sum(arr*weights)

x /= loops
p = 1 - (1-h)*x
print(p)

#
# WHWHWHWHWHWHWHWH(50%W).
# p = f(h) for a certain work history X p = 1/2(1-h)
# p = eta_w
# eta_w = g(params)
# h = f-1(g(params)) = 1 - 1/(2*eta_w[Raghav])
#

########################
import argparse
import logging
from typing import List, Dict

import matplotlib.pyplot as plt
import numpy as np

from constants import sections, s_pops, max_utility, job_risk, survival, env_params, player_data, player_types, T

logger: logging.Logger


def _init():
    """ all file related initialization code """
    global logger
    logger = logging.getLogger("Log")


def _add_args(parser: argparse.ArgumentParser):
    parser.add_argument("--eta-types", nargs=4, metavar=('nC', 'nP', 'nS', 'nU'),
                        help="Ratios of the population types (will be scaled to ensure the sum to 1)"
                             " in case of a population simulation",
                        type=float, default=[0.33, 0.33, 0.34, 0.00])
    pass


def _parse_args(args: argparse.Namespace):
    type_pops = np.asarray(args.eta_types)
    type_pops /= np.sum(type_pops)
    i_ratio = args.i_start / args.n_start
    return Population(type_pops=type_pops, i_ratio=i_ratio, t_max=args.t_max)


class Population:
    # @formatter:off

    # size range over which a persons belief of his own health (type) varies
    p_h_delta: float

    # constants used in various archetypes
    coward_data: Dict[str, float]
    player_data: Dict[str, float]

    # player division across ( section x archetype x i_stage )
    n_sections: int
    n_stages: int  # number of stages
    n_types: int = 4  # number of archetypes
    people: np.ndarray  # ratio of population for each cell
    deaths: np.ndarray  # ( section x archetype )
    eta_w: np.ndarray  # ratio of people who worked / will work in that cell
    n_age_groups: int
    n_job_types: int
    # for planner
    X: np.ndarray  # work measure, (section x i_stage)
    h: float = player_types['planner']['h']  # memory over past days
    influence_cap: float = player_types['planner']['cap']  # max influence on P[W/H] due to X

    # utility measures
    net_utility: np.ndarray

    # prototype enum
    C: int = 0  # archetype index for type coward
    P: int = 1  # archetype index for type planner
    S: int = 2  # archetype index for type simple
    U: int = 3  # archetype index for type unaware

    eta_iw: np.ndarray
    eta_ih: float

    T_max: int

    # For graphs
    t_deaths_ot: List[float]  # total deaths over time
    f_deaths_ot: List[float]  # fresh deaths over time
    t_deaths_ot_s: List[float]  # total deaths over time for each section
    t_utility_ot: List[float]  # total utility over
    d_utility_ot: List[float]  # daily utility over
    t_deaths_ot_a: List[float]  # total deaths over time for each age group
    eta_w_ot: List[float]  # Population going to work

    timeline: np.ndarray  # time scale
    infection_rate: List[float]  # infection rate over time
    archetype_dict: Dict[int, str]  # Contains Labels for n_types
    age_group_dict: Dict[int, str]  # Contains Labels for n_age_group
    sections_dict: Dict[int, str]  # Contains Labels for n_job_types

    # @formatter:on

    class ForcedExit(Exception):
        pass

    @staticmethod
    def _safe_divide(a: np.ndarray, b: np.ndarray):
        """ returns a/b if b != 0 else 0 """
        return np.divide(a, b, out=np.zeros_like(a), where=b != 0)

    def __init__(self, type_pops: np.ndarray, i_ratio: float, t_max: int):
        self.coward_data = player_types["coward"]
        self.player_data = player_data
        self.p_h_delta = self.player_data['p-healthy-fluctuation']

        self.n_sections = len(sections)
        self.n_stages = env_params['t-removal'] + 1
        self.n_age_groups = 3
        self.n_job_types = 5
        self.people = np.zeros((self.n_sections, self.n_types, self.n_stages))

        s_type_pops = np.expand_dims(np.array(s_pops), axis=1) * np.expand_dims(type_pops, axis=0)

        self.people[:, :, 1] = s_type_pops * i_ratio
        self.people[:, :, 0] = s_type_pops - self.people[:, :, 1]

        self.eta_w = np.zeros((self.n_sections, self.n_types, self.n_stages))
        self.eta_w[:, :, 0] = 1
        self.eta_w[:, :, -1] = 1

        self.net_utility = np.zeros(self.n_types)
        self.deaths = np.zeros((self.n_sections, self.n_types))

        self.eta_iw = np.zeros(self.n_sections)
        self.eta_ih = 0

        self.X = np.zeros((self.n_sections, self.n_stages))
        h = max_utility / (max_utility + (1 - survival) * self.player_data["u-death"])

        self.X[:, 0] = 1 / (1 - h)

        self.T_max = t_max

        self.t_deaths_ot = []
        self.f_deaths_ot = []
        self.t_deaths_ot_s = []
        self.t_deaths_ot_a = []
        self.t_utility_ot = []
        self.d_utility_ot = []

        self.eta_w_ot = []

        self.timeline = np.arange(self.T_max)
        self.infection_rate = []
        self.archetype_dict = {0: 'Gullible', 1: 'Planner', 2: 'Simple', 3:'Unaware'}
        self.sections_dict = {0: 'Primary', 1: 'Secondary', 2: 'Tertiary', 3: 'Essential', 4: 'Retired'}
        self.age_group_dict = {0: 'Young', 1: 'Middle', 2: 'Old'}
        # TODO: set linear instead of step
        self.threshold_sw = np.where(job_risk < self.coward_data['job-risk-threshold'],
                                     self.coward_data['low-risk-w-threshold'],
                                     self.coward_data['w-threshold'])
        self.threshold_sh = np.where(job_risk < self.coward_data['job-risk-threshold'],
                                     self.coward_data['low-risk-h-threshold'],
                                     self.coward_data['h-threshold'])

    def _execute_infection(self):
        """
        Adjusts self.people, self.eta_w and self.X forward by one day wrt the pandemic.
        """

        # 1. Find the probability of infection for S work and S home, ie. eta_iw, eta_ih respectively
        infectious_mask = np.zeros(self.n_stages)
        infectious_mask[env_params['t-infectious']:-1] = 1
        s_working_infectious = np.einsum("ijk,ijk,k -> j", self.people, self.eta_w, infectious_mask)

        working_infectious = np.sum(s_working_infectious)  # ratio of (working, infected) / total population
        infected = np.sum(self.people, axis=(0, 1))
        assert infected.shape == (self.n_stages,)
        total_infectious = np.sum(infected[env_params['t-infectious']:-1])
        total_infected = np.sum(infected[1:-1])
        self.infection_rate.append(float(total_infectious))
        coeff_i = 0.05
        coeff_wi = 1 * job_risk

        self.eta_ih = total_infectious * coeff_i
        self.eta_iw = np.clip(1 - (1 - total_infectious * coeff_i) * (1 - working_infectious * coeff_wi), self.eta_ih,
                              1)
        logger.info("total-infectious: {0}, working-infectious : {1}".format(total_infectious, working_infectious))
        logger.info("eta_iw: {0}, eta_ih: {1:.2f}".format(str(self.eta_iw), self.eta_ih))

        # 2. find out the ratio of susceptible that are getting infected
        eta_w = self.eta_w[:, :, 0]
        eta_iw = np.expand_dims(self.eta_iw, axis=1)
        eta_i = eta_w * eta_iw + (1 - eta_w) * self.eta_ih

        assert eta_i.shape == (self.n_sections, self.n_types)

        if total_infected == 0:
            raise self.ForcedExit("Infection Eradicated")

        # 3. execute the population transitions
        fresh_infected = eta_i * self.people[:, :, 0]
        fresh_recovered = self.people[:, :, -2] * np.expand_dims(survival, axis=1)
        fresh_deaths = self.people[:, :, -2] - fresh_recovered
        old_recovered = self.people[:, :, -1]

        self.people[:, :, 2:-1] = self.people[:, :, 1:-2]
        self.people[:, :, -1] += fresh_recovered
        self.people[:, :, 0] -= fresh_infected
        self.people[:, :, 1] = fresh_infected
        self.deaths += fresh_deaths

        # @Prajjwal - I think temp need not be multiplied in line 200, but just in case
        temp = np.asarray(s_pops).reshape(self.n_job_types, self.n_age_groups)
        age_group_deaths_s = np.sum(
            np.asarray(fresh_deaths[:, self.S]).reshape(self.n_job_types, self.n_age_groups), axis=0)
        if len(self.t_deaths_ot_s) != 0:
            self.t_deaths_ot_s.append(self.t_deaths_ot_s[-1] + fresh_deaths[:, self.S])
            self.t_deaths_ot_a.append(self.t_deaths_ot_a[-1] + age_group_deaths_s)
        else:
            self.t_deaths_ot_s.append(fresh_deaths[:, self.S])
            self.t_deaths_ot_a.append(age_group_deaths_s)
        self.eta_w_ot.append(np.sum(np.asarray(eta_w) * np.expand_dims(np.asarray(s_pops), axis=1), axis=0))
        self.t_deaths_ot.append(np.sum(self.deaths, axis=0))
        self.f_deaths_ot.append(np.sum(fresh_deaths, axis=0))

        logger.info("Fresh Deaths: " + str(np.sum(fresh_deaths, axis=0)))
        logger.info("Total Deaths: " + str(np.sum(self.deaths, axis=0)))

        # 4. execute the eta_w transitions (note: eta_w will be wrt the new population distribution)
        eta_wi = self._safe_divide(eta_iw * eta_w, eta_i)  # P[W given they get I]
        eta_ws = self._safe_divide((1 - eta_iw) * eta_w, (1 - eta_i))  # P[W given they remain S]
        self.eta_w[:, :, -1] = 1 - self._safe_divide(fresh_recovered, self.people[:, :, -1])
        self.eta_w[:, :, 2:-1] = self.eta_w[:, :, 1:-2]
        self.eta_w[:, :, 1] = eta_wi
        self.eta_w[:, :, 0] = eta_ws

        # 5. execute X transitions (note: this is only for planner archetype)
        self.X[:, -1] = self._safe_divide(
            self.X[:, -1] * old_recovered[:, self.P] + self.X[:, -2] * fresh_recovered[:, self.P],
            self.people[:, self.P, -1]
        )
        self.X[:, 2:-1] = self.X[:, 1:-2]
        pop_infected_today = self._safe_divide(self.eta_w[:, self.P, 1], self.eta_w[:, self.P, 0])
        self.X[:, 1] = pop_infected_today * self.X[:, 0]
        self.X[:, 0] = (1 + pop_infected_today) * self.X[:, 0]
        self.X = np.clip(self.X, 0, 1 / 1 - self.h)

    def _update_utility(self):
        """ Takes the actions supplied, increments the utility and returns the risks of infection """
        i_loss = np.clip(1 - np.arange(self.n_stages) / env_params['t-infectious'], 0, 1)
        i_loss[-1] = 1
        i_loss = i_loss ** 2
        s_utility = np.einsum("ijk,ijk,i,k -> j", self.people, self.eta_w, max_utility, i_loss)
        logger.info("Fresh utility: " + str(s_utility))
        self.net_utility += s_utility
        logger.info("Total utility: " + str(self.net_utility))

        if self.t_utility_ot:
            self.t_utility_ot.append(s_utility + self.t_utility_ot[-1])
        else:
            self.t_utility_ot.append(s_utility)

        self.d_utility_ot.append(s_utility)

    def _get_action_coward(self):
        """
        Forwards eta_w[all_archetypes, C, all_stages] by one day. That is, using self parameters (in
        particular self.eta_w, the ratio of people in each cell of (section x archetype x i_stage)
        who had chosen to work yesterday (they might have been in a different cell per the last axis
        then)) set self.eta_w to the ratio of people in each cell who will choose to work today.
        """
        w = np.zeros((self.n_sections, self.n_stages,))

        ts = env_params["t-symptoms"]
        last_w = self.eta_w[:, self.C, :]

        ratio_over_threshold_w = 1 - ((np.expand_dims(self.threshold_sw, axis=1) - (
                1 - np.expand_dims(np.arange(ts), axis=0) / ts - self.p_h_delta / 2)) / self.p_h_delta)
        ratio_over_threshold_h = 1 - ((np.expand_dims(self.threshold_sh, axis=1) - (
                1 - np.expand_dims(np.arange(ts), axis=0) / ts - self.p_h_delta / 2)) / self.p_h_delta)
        # ratio_over_threshold_h = 1 - ((threshold_sh - (1 - np.arange(ts) / ts - self.p_h_delta / 2)) / self.p_h_delta)
        ratio_over_threshold_w = np.clip(ratio_over_threshold_w, 0, 1)
        ratio_over_threshold_h = np.clip(ratio_over_threshold_h, 0, 1)
        w[:, :ts] = last_w[:, :ts] * ratio_over_threshold_w + (1 - last_w[:, :ts]) * ratio_over_threshold_h

        w[:, ts:-1] = last_w[:, ts:-1] * np.maximum((self.p_h_delta / 2 - self.threshold_sw) / self.p_h_delta, 0) \
                      + (1 - last_w[:, ts:-1]) * np.maximum((self.p_h_delta / 2 - self.threshold_sh) / self.p_h_delta,
                                                            0)
        w[:, -1] = 1

        self.eta_w[:, self.C, :] = np.asarray(w)

    def _get_action_simple(self):
        u_iw = (self.eta_iw - self.eta_ih) * (1 - survival) * self.player_data["u-death"]
        p_max = 1 + self.p_h_delta / 2 - np.arange(self.n_stages) / env_params["t-symptoms"]
        cutoff = np.expand_dims(np.sqrt(u_iw / max_utility), axis=1)
        eta_w = np.where(cutoff > 1, 0, (np.expand_dims(p_max, axis=0) - cutoff) / self.p_h_delta)
        eta_w[:, -1] = 1
        self.eta_w[:, self.S, :] = np.clip(eta_w, 0, 1)

    def _get_action_planner(self):
        self.X = self.eta_w[:, self.P, :] + self.h * self.X

        eta_w_s = self.eta_w[:, self.S, :]
        eta_w_del = np.minimum(1 - eta_w_s, eta_w_s) * 0.4
        eta_w_min = eta_w_s - eta_w_del
        eta_w_max = eta_w_s + eta_w_del

        self.eta_w[:, self.P, :] = (eta_w_min - eta_w_max) * (1 - self.h) * self.X + eta_w_max

    def _get_actions_unaware(self):
        self.eta_w[:, self.U, :] = 1


    def simulate(self):
        try:
            for T[0] in range(self.T_max):
                logger.debug("Population sections:\n" + str(self.people))
                logger.debug("Percentage working:\n" + str(self.eta_w))
                logger.debug(
                    "Final Utility: {}".format(self.net_utility + np.sum(self.deaths, axis=0) * player_data['u-death']))

                self._get_action_coward()
                self._get_action_simple()
                # simple must be called before planner !
                self._get_action_planner()
                self._update_utility()
                self._execute_infection()
        except self.ForcedExit as e:
            logger.info(e)
            print(e)

        logger.info("Population sections:\n" + str(self.people))
        logger.info("Percentage working:\n" + str(self.eta_w))
        logger.info("Final Utility: {}".format(self.net_utility + np.sum(self.deaths, axis=0) * player_data['u-death']))

    def total_death_plot(self):
        fig, ax1 = plt.subplots()
        color = 'tab:red'
        ax1.plot(self.timeline, self.infection_rate,color=color, label="Infection")
        ax1.set_ylabel("Infected Population")
        col = ['tab:blue', 'tab:orange', 'tab:green', 'tab:purple']
        ax2 = ax1.twinx()
        for i in range(self.n_types):
            ax2.plot(self.timeline, np.asarray(self.t_deaths_ot)[:, i],color = col[i], label=self.archetype_dict[i])

        ax2.set_xlabel("Time(days)")
        ax2.set_ylabel("Death percentage")
        plt.grid()
        ax1.legend(loc="center right")
        ax2.legend()
        fig.tight_layout()
        plt.title("Total Death Percentage for 4 archetypes")
        plt.savefig("graphs/total_deaths_with_unaware.jpg", bbox_inches='tight')

    #        plt.show()

    def fresh_death_plot(self):
        fig, ax1 = plt.subplots()
        color = 'tab:brown'
        ax1.plot(self.timeline, self.infection_rate, color=color, label="Infection")
        ax1.set_ylabel("Infected Population")

        ax2 = ax1.twinx()
        for i in range(self.n_types):
            ax2.plot(self.timeline, np.asarray(self.f_deaths_ot)[:, i], label=self.archetype_dict[i])

        ax2.set_xlabel("Time(days)")
        ax2.set_ylabel("Death percentage")
        plt.grid()
        ax1.legend(loc="upper left")
        ax2.legend()
        fig.tight_layout()
        plt.title("Daily Death Percentage for 3 archetypes")
        plt.savefig("graphs/fresh_deaths.jpg", bbox_inches='tight')

    def total_utility_plot(self):
        fig, ax1 = plt.subplots()
        color = 'tab:red'
        ax1.plot(self.timeline, self.infection_rate, color=color, label="Infection")
        ax1.set_ylabel("Infected Population")

        ax2 = ax1.twinx()
        for i in range(self.n_types):
            ax2.plot(self.timeline, np.asarray(self.t_utility_ot)[:, i], label=self.archetype_dict[i])

        ax2.set_xlabel("Time(days)")
        ax2.set_ylabel(" Total Utility")
        plt.grid()
        ax1.legend(loc="upper right")
        ax2.legend()
        fig.tight_layout()
        plt.title("Total Utility for 3 archetypes")
        plt.savefig("graphs/total_utility.jpg", bbox_inches='tight')

    def daily_utility_plot(self):
        fig, ax1 = plt.subplots()
        color = 'tab:red'
        ax1.plot(self.timeline, self.infection_rate, color=color, label="Infection")
        ax1.set_ylabel("Infected Population")

        ax2 = ax1.twinx()
        for i in range(self.n_types):
            ax2.plot(self.timeline, np.asarray(self.d_utility_ot)[:, i], label=self.archetype_dict[i])

        ax2.set_xlabel("Time(days)")
        ax2.set_ylabel("Daily Utility")
        plt.grid()
        ax1.legend(loc="upper right")
        ax2.legend()
        fig.tight_layout()
        plt.title("Daily Utility for 3 archetypes")
        plt.savefig("graphs/daily_utility.jpg", bbox_inches='tight')

    def total_death_plot_sections_simple(self):
        fig, ax1 = plt.subplots()
        for j in range(self.n_sections):
            ax1.plot(self.timeline, np.asarray(self.t_deaths_ot_s)[:, j], label=sections[j])
        plt.tight_layout()
        ax1.set_xlabel("Time(days)")
        ax1.set_ylabel("Death percentage")
        plt.grid()
        ax1.legend()
        plt.title("Death Percentage for all Sections for Archetype Simple")
        plt.savefig("graphs/total_deaths_per_section.jpg", bbox_inches='tight')


    def total_death_plot_age_group_simple(self):
        fig, ax1 = plt.subplots()
        for j in range(self.n_age_groups):
            ax1.plot(self.timeline, np.asarray(self.t_deaths_ot_a)[:, j], label=self.age_group_dict[j])

        ax1.set_xlabel("Time(days)")
        ax1.set_ylabel("Death percentage")
        plt.grid()
        ax1.legend()
        plt.tight_layout()
        plt.title("Death Percentage for all Sections for Archetype Simple")
        plt.savefig("graphs/total_deaths_per_age_group.jpg", bbox_inches='tight')

    def susceptible_eta_w(self):
        fig, ax1 = plt.subplots()
        color = 'tab:red'
        ax1.plot(self.timeline, self.infection_rate, color=color, label="Infection")
        ax1.set_ylabel("Infected Population")
        # TODO : Planner and Simple are literally following each other, tried altering planner ( took 0.8 instead of
        #  0.4, not working), so have only drawn this graph for one type
        ax2 = ax1.twinx()
        # for i in range(self.n_types):
        ax2.plot(self.timeline, np.asarray(self.eta_w_ot)[:, self.P], label=self.archetype_dict[self.P])

        ax2.set_xlabel("Time(days)")
        ax2.set_ylabel("Daily Working Population")
        plt.grid()
        ax1.legend(loc="center right")
        ax2.legend()
        plt.tight_layout()
        # plt.gcf().subplots_adjust(bottom=0.15)
        plt.title("Daily Working Population for Archetype Planner")
        plt.savefig("graphs/daily_working.jpg", bbox_inches='tight')

    def plot_graphs(self):
        self.total_death_plot()
        # self.fresh_death_plot()
        # self.total_utility_plot()
        # self.daily_utility_plot()
        # self.total_death_plot_sections_simple()
        # self.total_death_plot_age_group_simple()
        # self.susceptible_eta_w()