diff --git a/docs/SMF_example 2.png b/docs/SMF_example 2.png new file mode 100644 index 000000000..c983c8561 Binary files /dev/null and b/docs/SMF_example 2.png differ diff --git a/docs/SMF_example.png b/docs/SMF_example.png new file mode 100644 index 000000000..c983c8561 Binary files /dev/null and b/docs/SMF_example.png differ diff --git a/docs/active_set.rst b/docs/active_set.rst new file mode 100644 index 000000000..c210bd9a6 --- /dev/null +++ b/docs/active_set.rst @@ -0,0 +1,81 @@ +Solver: Active-set Method (ACTIVE-SET) +============================================================================================= + +Description: +------------ +The solver is dedicated for problems with linear constraints in the form of :math:`C_e x = de, C_i x \leq di`. +It employs the active-set method mentioned in Jorge Norcedal's Numerical Optimization. + +Modifications & Implementation: +------------------------------- + +**find_feasible_initial**: Find an initial feasible solution if one is not provided. + +**line_search**: A back-tracking line search method. + +**compute_search_direction**: Compute a search direction by solving a direction-finding quadratic subproblem at solution x. +It takes the index set of the active constraints (`W`), objective gradient at solution x (`g(x)`), and the concatenated constraint +coefficient matrix (`C`) as inputs and returns the optimal search direction along with the Lagrange multipliers of the active constraints. + +.. math:: + + \begin{align} + \min && (1/2)|| d ||^2+ g(x)^T d \\\\ + \text{s.t.} & C_k^T d = 0, \quad \forall k \in W \\ + \end{align} + +**_feasible**: Check whether a solution x is feasible to the problem. + +**finite_diff**: Approximate objective gradient using the finite difference method. + +Scope: +------ +* objective_type: single + +* constraint_type: box, deterministic (linear) + +* variable_type: continuous + +* gradient_observations: not available + +Solver Factors: +--------------- +* crn_across_solns: Use CRN across solutions? + + * Default: True + +* r: number of replications taken at each solution + + * Default: 30 + +* alpha: Tolerance for sufficient decrease condition. + + * Default: 0.2 + +* beta: Step size reduction factor in line search. + + * Default: 0.9 + +* alpha_max: Maximum step size + + * Default: 10.0 + +* lambda: Magnifying factor for r inside the finite difference function + + * Default: 2 + +* tol: Floating point tolerance for checking tightness of constraints + + * Default: 1e-7 + +* tol2: Floating point tolerance for checking closeness of dot product to zero + + * Default: 1e-7 + +* finite_diff_step: Step size for finite difference + + * Default: 1e-5 + +References: +=========== +Nocedal, J., & Wright, S. J. (2006). Numerical optimization (2nd ed.). Springer. \ No newline at end of file diff --git a/docs/covid.rst b/docs/covid.rst new file mode 100644 index 000000000..dd149832c --- /dev/null +++ b/docs/covid.rst @@ -0,0 +1,461 @@ +Model: COVID-19 Disease Progression with Testing and Vaccination (COVID) +================================================================= + +Description: +------------ +COVID-19 is a contagious respiratory disease with a high trasmission rate. A college campus implements +regular survelliance testing to identify, isolate, and reduce disease spread. However, the model can also +be generalized and simulated the interaction among any given groups of people. +The initial prevalance level of the disease is :math:`init_infect_percent`. The population is divided +into different groups with intra-group interaction matrix :math:`inter_rate`. There is a probability of :math:`p_trans` (or :math:`p_trans_vac`) +transmissions per interaction. The transmission rate per individual can be calculated by multiplying +:math:`inter_rate` by :math:`p_trans` (or :math:`p_trans_vac`). The disease progression for each individual follows the following semi-Markov process: + +.. image:: covid_test_compartments.png + :width: 400 + + +Clarifications: + +1. All recovered people will not be infected again; + +2. In the exposed state, people cannot infect others; during the infectious/asymptomatic state, people can infect other; + +3. Remove all infected individuals (except those in symptomatic state) to isolation if they are tested and have result of positive at each day; + +4. All individuals in symptomatic state will be isolated; + +5. For individuals in isolation state, assume they follow the same disease progression; + +6. Given being infected, the transmission probabilities for vaccinated and non-vaccinated individuals are different; + +7. Given being infected, receiving vaccination does influence the probability of being asymptomatic but not influence process of illness; + +8. Assume all the N vaccinations are received at the beginning of the cycle, and we do not consider multiple doses here; + +9. Assume all initial infectious people are unvaccinated; + +10. Assume all symptomatic people are isolated; + +11. Assume all testing happens at the end of the day. (i.e. people will be isolated at the end of the day); + +12. At day :math:`t`, for individuals who will transit to next the state in day :math:`t+1`, assume they will not be tested at day :math:`t` but at day :math:`t+1`; + +13. At the end of isolation, individuals will return to recovered state. + + +The simulation is generated as follows: + +1. For each day in :math:`n` and for each group :math:`g`, generate newly exposed individuals. + +2. At day :math:`t`, for each state, each group :math:`g`, consider the newly coming individuals and the remaining scheduled people: + + (a) Schedule for newly coming individuals: generate :math:`day_exp` number of days remaining in exposed and :math:`gen_exp_num` number of people who are going to move to infectious state at date :math:`day_exp + t`. + + (b) Generate isolation starting schedule for the scheduled exposed people: generate day :math:`day_move_iso_exp` and :math:`num_iso_exp` number of people who are going to move to isolation at date :math:`day_move_iso_exp + t`. + + (c) Remove the scheduled isolated people from the scheduled exposed group and the total exposed group. + + (d) At the end of :math:`day_exp + t`, the remaining scheduled exposed people are going to move to infectious. + + (e) Generate the schedule for exposed isolated people: generate :math:`day_exp_iso` number of days remaining in exposed state and :math:`gen_exp_num_iso` number of people who are going to move to infectious isolation state at date :math:`day_exp_iso + t`. + + (f) At the end of :math:`day_exp_iso + t`, the :math:`gen_exp_num_iso` number of isolated exposed individuals are going to move to infectious isolation state. + + (g) Schedule for newly coming individuals to infectious state: generate :math:`day_inf` number of days remaining in infectious and :math:`gen_inf_num` number of people who are going to move to symptomatic /asymptomatic state at date :math:`day_inf + t`. + + (h) Generate isolation starting schedule for the scheduled infectious people: generate day :math:`day_move_iso_inf` and :math:`num_iso_inf` number of people who are going to move to isolation at date :math:`day_move_iso_inf + t`. + + (i) Remove the scheduled isolated people from the scheduled infectious group and the total infectious group. + + (j) At the end of :math:`day_inf + t`, generate the number of asymptomatic people from the remaining scheduled infectious people, move them to asymptomatic state. Symptomatic people are all going to symptomatic isolation. + + (k) Generate the schedule for infectious isolated people: generate :math:`day_inf_iso` number of days remaining in infectious state and :math:`gen_inf_num_iso` number of people who are going to move to symptomatic /asymptomatic isolation state at date :math:`day_inf_iso + t`. + + (l) At the end of :math:`day_inf_iso + t`, generate the number of symptomatic and asymptomatic people from the :math:`gen_inf_num_iso` number of isolated infectious people, move them to symptomatic/asymptomatic isolation state. + + (m) Generate the schedule for symptomatic individuals: generate :math:`day_sym` number of days remaining in symptomatic state and :math:`gen_sym_num` number of people who are going to move to recover at date :math:`day_sym + t`. + + (n) Generate the schedule for asymptomatic individuals: generate :math:`day_asym ` number of days remaining in asymptomatic state and :math:`gen_asym_num` number of people who are going to move to recover at date :math:`day_asym + t`. + + (o) Generate isolation starting schedule for the scheduled asymptomatic people: generate day :math:`day_move_iso_asym` and :math:`num_iso_asym` number of people who are going to move to asymptomatic isolation state at date :math:`day_move_iso_asym + t`. + + (p) Remove the scheduled isolated asymptomatic people from the scheduled asymptomatic group and the total asymptomatic group. + + (q) At the end of :math:`day_asym + t`, the remaining scheduled asymptomatic people are going to move to recover state. + + (r) Generate the schedule for asymptomatic isolated individuals: generate :math:`day_asym_iso` number of days remaining in asymptomatic isolation state and :math:`gen_asym_num_iso` number of poeple who are going to recover at date :math:`day_asym_iso + t`. + + (s) At the end of :math:`day_asym_iso + t`, the remaining isolated asymptomatic people are going to move to recover state. + + +Sources of Randomness: +---------------------- +There are six sources of randomness. + +1. The number of newly exposed (non)vaccinated individuals on each day follows a Poisson distribution with mean equals transmission +rate times number of infected (infectious + symptomatic + asymptomatic) individuals times fraction of susceptible (non)vaccinated individuals. + + mean of new exposed people: + For vaccinated group: + :math:`new_exp_vac = p_trans_vac * iter_rate * number of infected and asymptomatic individuals * vac_sus/ (sum(group_size)- # isolated individuals)` + For non-vaccinated group: + :math:`new_exp_nonvac = p_trans * iter_rate * number of infected and asymptomatic individuals * non_vac_sus / (sum(group_size)- # isolated individuals)` + +2. For each day, assume individuals can only stay in exposed state for no longer than :math:`max_exp` days, the number of individuals move to infectious state follows multinomial distribution with :math:`n=#new_exp_(non)vac` +and :math:`p=[p_1, p_2, ..., p_max_exp]` with :math:`p_i` is the pmf of possion distribution with :math:`x=i` and normalize the sum of :math:`p` to be 1. + +3. For each day, assume individuals can only stay in infectious state for no longer than :math:`max_inf` days, the number of individuals move to (a)symptomatic state follows multinomial distribution with :math:`n=#new_inf_(non)vac` +and :math:`p=[p_1, p_2, ..., p_max_inf]` with :math:`p_i` is the pmf of possion distribution with :math:`x=i` and normalize the sum of :math:`p` to be 1. + +4. For each day, assume individuals can only stay in symptomatic state for no longer than :math:`max_symp` days, the number of individuals move to recover state follows multinomial distribution with :math:`n=#new_sym_(non)vac` +and :math:`p=[p_1, p_2, ..., p_max_sym]` with :math:`p_i` is the pmf of possion distribution with :math:`x=i` and normalize the sum of :math:`p` to be 1. + +5. For each day, assume individuals can only stay in asymptomatic state for no longer than :math:`max_asymp` days, the number of individuals move to recover state follows multinomial distribution with :math:`n=#new_asym_(non)vac` +and :math:`p=[p_1, p_2, ..., p_max_asym]` with :math:`p_i` is the pmf of possion distribution with :math:`x=i` and normalize the sum of :math:`p` to be 1. + +6. An exposed/infectious/asymptomatic individual in group :math:`g` has a probability :math:`freq_g` of being tested and moved to the isolated states. + +7. At day :math:`t`, for state exposed/infectious/asymptomatic, assume :math:`n_day` number of individuals are going to leave to next infected states at date :math:`t + day`, the number of individuals move to isolation state follows multinomial distribution with :math:`n=n_day` +and :math:`p=[p_0, p_1, ..., p_{day-1}]` with :math:`p_i = (freq_g)^i` for :math:`0 0` for each :math:`i`, +but there is an associated cost. Unlike the original san problem, we now want to minimize :math:`ET(\theta)`, +where :math:`T(\theta)` is the (random) duration of the longest path from :math:`a` +to :math:`i`. + +The objective function is convex in :math:`\theta`. An IPA estimator of the gradient +is also given in the code. + +Constraints: +------------ +We require that :math:`\theta_i > 0` for each :math:`i`. +Additionaly, we include another constraint to impose a lower bound to the sum of arc_means. +which is :math:`\sum_i \theta_i \geq sum\_lb` + +Problem Factors: +---------------- +* budget: Max # of replications for a solver to take. + + * Default: 10000 + +* arc_costs: Cost associated to each arc. + + * Default: (1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1) + +* sum_lb: The lower bound for sum of arc_means. + + * Default: 30.0 + +Fixed Model Factors: +-------------------- +* N/A + +Starting Solution: +------------------ +* initial_solution: (8,) * 13 + +Random Solutions: +----------------- +Use acceptance-rejection to sample each arc_mean uniformly from a lognormal distribution with +2.5- and 97.5-percentiles at 0.1 and 10 respectively such that the arc_mean remains on the feasible region. + +Optimal Solution: +----------------- +Unknown + Optimal Objective Function Value: --------------------------------- Unknown \ No newline at end of file diff --git a/docs/smf.rst b/docs/smf.rst new file mode 100644 index 000000000..7c23370c1 --- /dev/null +++ b/docs/smf.rst @@ -0,0 +1,116 @@ +Model: Stochastic Max Flow (smf) +======================================== +Description: +------------ +Consider a max-flow problem with source node :math:`s` +and sink node :math:`t` and each arc :math:`i` is associated +with a capacity :math:`C_i` + +An example of a network with 10 total nodes and 20 arcs is shown below + +.. image:: SMF_example.png + :alt: The SMF diagram has failed to display + :width: 500 + +Sources of Randomness: +---------------------- +Noise in reduction of arc capacity is multivariate normally distributed +with mean mean_noise and covariance matrix of cov_noise. +The capacity of arc i to the max-flow solver model would be :math:`1000*max(assigned\_capacities[i]-noise[i],0)`, it was scaled up 1000 times since the maxflow solver would only solve for integer flow + + +Model Factors: +--------------- +* num_nodes: Total number of nodes. + + * Default: 10 + +* arcs: List of arcs + + * Default: 20 arcs, [(S,1),(S,2),(S,3),(1,2),(1,4),(2,4),(4,2),(3,2),(2,5),(4,5),(3,6),(3,7),(6,5),(6,7),(5,8),(6,8),(6,t),(7,t),(8,t)] + +* source_index": Source node index + + * default: 0 + +* sink_index: Sink node index + + * default: 9 + +* assigned_capacities: Assigned capacity of each arc + + * default: [5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5] + +* mean_noise: The mean noise in reduction of arc capacities + + * default: [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] + +* cov_noise: Covariance matrix of noise + + * default: Diagonal of 4, the rest 0 + +Responses: +---------- +* max_flow: the maximum flow from the source node S to the sink node t (scaled back down 1000 times) + + :math:`\sum_{(s,j) \in arcs} \nobreakspace \phi_si \nobreakspace / \nobreakspace 1000` + +References: +----------- +No reference + +Optimization Problem: Maximize system flow with random noice in capacity (SMF-1) +================================================================================ + +Decision Variables: +-------------------- +* assigned_capacities :math:`x_i` + +Objectives: +------------ +We want to maximize the expected max flow from the source node s to the sink node t. + +The maximum flow of the system, which is the objective function can then be represented +by the total flow out from the source node, which is + +:math:`f(\phi) = \sum_{(s,i) \in arcs}\nobreakspace \ \phi_si` + +Constraints: +------------ + +* We constrain the total capacity to be allocated over the arcs to be :math:`b` + + * :math:`\sum_{(i,j) \in arcs} \nobreakspace x_ij \nobreakspace\leq \nobreakspace b` + + * :math:`0 \nobreakspace \leq \nobreakspace x_ij\ \ \ \forall (s,i) \in arcs` + + +Problem Factors: +---------------- +* budget: Max # of replications for a solver to take. + + * Default: 10000 + +* b: Total set-capacity that can be allocated to arcs + + * Default: 100 + +Fixed Model Factors: +-------------------- +* N/A + +Starting Solution: +------------------ +* initial_solution: (1,) * 20 + +Random Solutions: +----------------- +Sample uniformly from the simplex of :math:`\sum_{(i,j) \in arcs} \nobreakspace x_ij \nobreakspace \leq \nobreakspace b` and :math:`0 \nobreakspace \leq \nobreakspace x_ij\ \forall (s,i) \in arcs` + +Optimal Solution: +----------------- +Unknown + +Optimal Objective Function Value: +--------------------------------- +Unknown \ No newline at end of file diff --git a/simopt/models/covid.py b/simopt/models/covid.py new file mode 100644 index 000000000..de02dcf5e --- /dev/null +++ b/simopt/models/covid.py @@ -0,0 +1,2247 @@ +""" +Summary +------- +Under some testing and vaccination policy, simulate the spread of COVID-19 +over a period of time by multinomial distribution. +A detailed description of the model/problem can be found +`here `_. +""" +import numpy as np + +from ..base import Model, Problem + + +class COVID_vac(Model): + """ + A model that simulates the spread of COVID disease under isolation policy and vaccination policy + with a multinomial distribution. + Returns the total number of infected people and the total number of patients who have observable symptoms + during the whole period. + + Attributes + ---------- + name : string + name of model + n_rngs : int + number of random-number generators used to run a simulation replication + n_responses : int + number of responses (performance measures) + factors : dict + changeable factors of the simulation model + specifications : dict + details of each factor (for GUI, data validation, and defaults) + check_factor_list : dict + switch case for checking factor simulatability + + Arguments + --------- + fixed_factors : dict + fixed_factors of the simulation model + + See also + -------- + base.Model + """ + def __init__(self, fixed_factors={}): + self.name = "COVID_vac" + self.n_rngs = 7 + self.n_responses = 6 + self.factors = fixed_factors + self.specifications = { + "num_groups": { + "description": "Number of groups.", + "datatype": int, + "default": 3 + # "default": 7 + }, + "p_trans": { + "description": "Probability of transmission per interaction.", + "datatype": float, + "default": 0.018 + }, + "p_trans_vac": { + "description": "Probability of transmission per interaction if being vaccinated", + "datatype": float, + "default": 0.01 + }, + "inter_rate": { + "description": "Interaction rates between two groups per day", + "datatype": tuple, + "default": (10.58, 5, 2, 4, 6.37, 3, 6.9, 4, 2) + # "default": (10.57, 4, 1, 1, 2, 1, 1, 6, 12.8, 2, 2, 3.5, 1, 1, 2, 2.5, 8.4, 1.5, 1.5, 1, 2, 1.5, 1.5, 1, + # 11.61, 3.31, 3.5, 1, 1, 4.29, 4.91, 2.8, 2.9, 2.22, 1, 2.3, 1, 1, 1, 1, 2.2, 8.1, 2.75, 1, 1, 2.7, 1, 2.99, 1) # For 7 groups + }, + "group_size": { + "description": "Size of each group.", + "datatype": tuple, + "default": (8123, 4921, 3598) + # "default": (8123, 3645, 4921, 3598, 3598, 1907, 4478) #total 30270, for 7 gruops + }, + "lamb_exp_inf": { + "description": "Mean number of days from exposed to infectious.", + "datatype": float, + "default": 2.0 + }, + "lamb_inf_sym": { + "description": "Mean number of days from infectious to symptomatic.", + "datatype": float, + "default": 3.0 + }, + "lamb_sym": { + "description": "Mean number of days from symptomatic/asymptomatic to recovered.", + "datatype": float, + "default": 12.0 + }, + "n": { + "description": "Number of days to simulate.", + "datatype": int, + "default": 200 + }, + "init_infect_percent": { + "description": "Initial prevalance level.", + "datatype": tuple, + "default": (0.00200, 0.00121, 0.0008) + # "default": (0.00200, 0.00121, 0.0008, 0.002, 0.00121, 0.008, 0.002) # For 7 gruops + }, + "asymp_rate": { + "description": "Probability of being asymptomatic.", + "datatype": float, + "default": 0.35 + }, + "asymp_rate_vac": { + "description": "Probability of being asymptomatic for vaccinated individuals.", + "datatype": float, + "default": 0.5 + }, + "vac": { + "description": "The fraction of people being vaccinated in each group.", + "datatype": tuple, + "default": (0.8, 0.3, 0) + # "default": (1/7,1/7,1/7,1/7,1/7,1/7,1/7) # For 7 groups + }, + "total_vac": { + "description": "The total number of vaccines.", + "datatype": int, + "default": 8000 + }, + "freq": { + "description": "Testing frequency of each group.", + "datatype": tuple, + "default": (0/7, 0/7, 0/7) + # "default": (0/7,0/7,0/7,0/7,0/7,0/7,0/7) # For 7 groups + }, + "freq_vac": { + "description": "Testing frequency of each group.", + "datatype": tuple, + "default": (0/7, 0/7, 0/7, 0.8, 0.3, 0) + # "default": (0/7,0/7,0/7,0/7,0/7,0/7,0/7,0,0,0,0,0,0,0) # For 7 groups + }, + "false_neg": { + "description": "False negative rate.", + "datatype": float, + "default": 0.12 + }, + "w": { + "description": "Protection weight for each group.", + "datatype": tuple, + "default": (1, 1, 1) + } + } + + self.check_factor_list = { + "num_groups": self.check_num_groups, + "p_trans": self.check_p_trans, + "p_trans_vac": self.check_p_trans_vac, + "inter_rate": self.check_inter_rate, + "group_size": self.check_group_size, + "lamb_exp_inf": self.check_lamb_exp_inf, + "lamb_inf_sym": self.check_lamb_inf_sym, + "lamb_sym": self.check_lamb_sym, + "n": self.check_n, + "init_infect_percent": self.check_init_infect_percent, + "asymp_rate": self.check_asymp_rate, + "asymp_rate_vac": self.check_asymp_rate_vac, + "vac": self.check_vac, + "total_vac": self.check_total_vac, + "freq": self.check_freq, + "freq_vac": self.check_freq_vac, + "false_neg": self.check_false_neg, + "w": self.check_w + } + # Set factors of the simulation model. + super().__init__(fixed_factors) + + def check_num_groups(self): + return self.factors["num_groups"] > 0 + + def check_p_trans(self): + return self.factors["p_trans"] > 0 + + def check_p_trans_vac(self): + return self.factors["p_trans_vac"] > 0 + + def check_inter_rate(self): + return all(np.array(self.factors["inter_rate"]) >= 0) & (len(self.factors["inter_rate"]) == self.factors["num_groups"]**2) + + def check_group_size(self): + return all(np.array(self.factors["group_size"]) >= 0) & (len(self.factors["group_size"]) == self.factors["num_groups"]) + + def check_lamb_exp_inf(self): + return self.factors["lamb_exp_inf"] > 0 + + def check_lamb_inf_sym(self): + return self.factors["lamb_inf_sym"] > 0 + + def check_lamb_sym(self): + return self.factors["lamb_sym"] > 0 + + def check_n(self): + return self.factors["n"] > 0 + + def check_init_infect_percent(self): + return all(np.array(self.factors["init_infect_percent"]) >= 0) & (len(self.factors["init_infect_percent"]) == self.factors["num_groups"]) + + def check_asymp_rate(self): + return (self.factors["asymp_rate"] > 0) & (self.factors["asymp_rate"] < 1) + + def check_asymp_rate_vac(self): + return (self.factors["asymp_rate_vac"] > 0) & (self.factors["asymp_rate_vac"] < 1) + + def check_vac(self): + return all(np.array(self.factors["vac"]) >= 0) & (len(self.factors["init_infect_percent"]) == self.factors["num_groups"]) + + def check_total_vac(self): + return self.factors["total_vac"] > 0 + + def check_freq(self): + return all(np.array(self.factors["freq"]) >= 0) & (len(self.factors["freq"]) == self.factors["num_groups"]) + + def check_freq_vac(self): + return all(np.array(self.factors["freq_vac"]) >= 0) & (len(self.factors["freq_vac"]) == self.factors["num_groups"] * 2) + + def check_false_neg(self): + return (self.factors["false_neg"] > 0) & (self.factors["false_neg"] < 1) + + def check_w(self): + return (self.factors["w"] >= 0) + + def replicate(self, rng_list): + """ + Simulate a single replication for the current model factors. + + Arguments + --------- + rng_list : list of rng.MRG32k3a objects + rngs for model to use when simulating a replication + + Returns + ------- + responses : dict + performance measures of interest + "num_infected" = Number of infected individuals per day + "num_susceptible" = Number of susceptible individuals per day + "num_exposed" = Number of exposed individuals per day + "num_recovered" = Number of recovered individuals per day + "num_symptomatic" = Number of symptomatic individuals per day + "total_symp" = Total number of symptomatic individuals + """ + # Designate random number generator for generating Poisson random variables. + poisson_numexp_rng = rng_list[0] + binomial_asymp_rng = rng_list[1] + multinomial_exp_rng = rng_list[2] + multinomial_inf_rng = rng_list[3] + multinomial_symp_rng = rng_list[4] + multinomial_asymp_rng = rng_list[5] + multinomial_iso_rng = rng_list[6] + + # Reshape the transmission rate + inter_rate = np.reshape(np.array(self.factors["inter_rate"]), (self.factors["num_groups"], self.factors["num_groups"])) + # Calculate the transmission rate + t_rate = inter_rate * self.factors["p_trans"] + t_rate = np.sum(t_rate, axis=1) + t_rate_vac = inter_rate * self.factors["p_trans_vac"] + t_rate_vac = np.sum(t_rate_vac, axis=1) + + # Initialize states, each row is one day, each column is one group + susceptible = np.zeros((self.factors["n"], self.factors["num_groups"])) + sus_vac = np.zeros((self.factors["n"], self.factors["num_groups"])) + sus_non_vac = np.zeros((self.factors["n"], self.factors["num_groups"])) + + exposed_vac = np.zeros((self.factors["n"], self.factors["num_groups"])) + exposed_non_vac = np.zeros((self.factors["n"], self.factors["num_groups"])) + + infectious_vac = np.zeros((self.factors["n"], self.factors["num_groups"])) + infectious_non_vac = np.zeros((self.factors["n"], self.factors["num_groups"])) + + asymptomatic_vac = np.zeros((self.factors["n"], self.factors["num_groups"])) + asymptomatic_non_vac = np.zeros((self.factors["n"], self.factors["num_groups"])) + + symptomatic_vac = np.zeros((self.factors["n"], self.factors["num_groups"])) + symptomatic_non_vac = np.zeros((self.factors["n"], self.factors["num_groups"])) + + isolation_exp_non_vac = np.zeros((self.factors["n"], self.factors["num_groups"])) + isolation_inf_non_vac = np.zeros((self.factors["n"], self.factors["num_groups"])) + isolation_symp_non_vac = np.zeros((self.factors["n"], self.factors["num_groups"])) + isolation_asymp_non_vac = np.zeros((self.factors["n"], self.factors["num_groups"])) + + isolation_exp_vac = np.zeros((self.factors["n"], self.factors["num_groups"])) + isolation_inf_vac = np.zeros((self.factors["n"], self.factors["num_groups"])) + isolation_symp_vac = np.zeros((self.factors["n"], self.factors["num_groups"])) + isolation_asymp_vac = np.zeros((self.factors["n"], self.factors["num_groups"])) + + recovered = np.zeros((self.factors["n"], self.factors["num_groups"])) + + # Initialize the performance measures of interest + num_infected = np.zeros((self.factors["n"], self.factors["num_groups"])) + num_exposed = np.zeros((self.factors["n"], self.factors["num_groups"])) + num_recovered = np.zeros((self.factors["n"], self.factors["num_groups"])) + num_susceptible = np.zeros((self.factors["n"], self.factors["num_groups"])) + last_susceptible = np.zeros(self.factors["num_groups"]) + num_symptomatic = np.zeros((self.factors["n"], self.factors["num_groups"])) + tot_num_syptomatic = 0 + tot_num_isolated = 0 + num_isolation = np.zeros((self.factors["n"], self.factors["num_groups"])) + + # Number of newly coming people for each day, each state + inf_arr = np.zeros((self.factors["n"], self.factors["num_groups"])) + sym_arr = np.zeros((self.factors["n"], self.factors["num_groups"])) + asym_arr = np.zeros((self.factors["n"], self.factors["num_groups"])) + inf_arr_vac = np.zeros((self.factors["n"], self.factors["num_groups"])) + sym_arr_vac = np.zeros((self.factors["n"], self.factors["num_groups"])) + asym_arr_vac = np.zeros((self.factors["n"], self.factors["num_groups"])) + + # Number of newly isolated people for each day, each state + exp_arr_iso = np.zeros((self.factors["n"], self.factors["num_groups"])) + inf_arr_iso = np.zeros((self.factors["n"], self.factors["num_groups"])) + asym_arr_iso = np.zeros((self.factors["n"], self.factors["num_groups"])) + exp_arr_vac_iso = np.zeros((self.factors["n"], self.factors["num_groups"])) + inf_arr_vac_iso = np.zeros((self.factors["n"], self.factors["num_groups"])) + asym_arr_vac_iso = np.zeros((self.factors["n"], self.factors["num_groups"])) + + # Get the current vaccine policy + vac = self.factors["freq_vac"][self.factors["num_groups"]:] + freq = self.factors["freq_vac"][:self.factors["num_groups"]] + + # Add day 0 num infections + infectious_non_vac[0, :] = np.ceil(np.multiply(list(self.factors["group_size"]), list(self.factors["init_infect_percent"]))) + infectious_vac[0, :] = np.multiply(np.multiply(list(self.factors["group_size"]), list(self.factors["init_infect_percent"])), 0) + susceptible[0, :] = np.subtract(list(self.factors["group_size"]), infectious_non_vac[0, :]) + sus_vac[0, :] = np.ceil(np.multiply(list(self.factors["group_size"]), list(vac))) + sus_non_vac[0, :] = np.subtract(susceptible[0, :], sus_vac[0, :]) + + for g in range(self.factors["num_groups"]): + infect = infectious_vac[0, g] + infectious_non_vac[0, g] + num_infected[0][g] = np.sum(infect) + num_susceptible[0][g] = sus_vac[0, g] + sus_non_vac[0, g] + + # Initialization -- calculate the probability distribution for alias + # For exposed: + p_exp = [] # Min=1, max=7 (by original setting) + exp_max = 7 # Min=1, max=7 (by original setting) + lamb_exp = self.factors["lamb_exp_inf"] + p0_e = lamb_exp ** 0 * np.exp(-lamb_exp) / 1 + for i in range(1, 1 + exp_max): + p0_e = p0_e * lamb_exp / i + p_exp.append(p0_e) + p_exp = p_exp / np.sum(p_exp) # Normalize to sum equal to 1 + + dict_exp = {} + for i in range(1, 1 + exp_max): + dict_exp[i] = p_exp[i - 1] + + exp_table, exp_alias, exp_value = multinomial_exp_rng.alias_init(dict_exp) # New alias: exp_value + + # For infected: + p_inf = [] + inf_max = 8 + lamb_inf = self.factors["lamb_inf_sym"] + p0_i = lamb_inf ** 0 * np.exp(-lamb_inf) / 1 + for i in range(1, 1 + inf_max): + p0_i = p0_i * lamb_inf / i + p_inf.append(p0_i) + p_inf = p_inf / np.sum(p_inf) + + dict_inf = {} + for i in range(1, 1 + inf_max): + dict_inf[i] = p_inf[i - 1] + + inf_table, inf_alias, inf_value = multinomial_inf_rng.alias_init(dict_inf) + + # For symptomatic: + p_sym = [] + asym_sym_max = 20 + lamb_asy_sym = self.factors["lamb_sym"] + p0_as_s = lamb_asy_sym ** 0 * np.exp(-lamb_asy_sym) / 1 + for i in range(1, 1 + asym_sym_max): + p0_as_s = p0_as_s * lamb_asy_sym / i + p_sym.append(p0_as_s) + p_sym = p_sym / np.sum(p_sym) + + dict_sym = {} + for i in range(1, 1 + asym_sym_max): + dict_sym[i] = p_sym[i - 1] + + sym_table, sym_alias, sym_value = multinomial_symp_rng.alias_init(dict_sym) + asym_table, asym_alias, asym_value = multinomial_asymp_rng.alias_init(dict_sym) + + # Isolation schedule generator + def calculate_iso_mul(day_max, pb, n): # Note: day_max is the number of days before change, n is the number of people + p = [] + p0 = pb # pb = freq*(1-self.factors["false_neg"]) + p.append(p0) + for i in range(1, day_max): + p0 = p0 * p0 + p.append(p0) + p.append(1 - np.sum(p)) # p is a list with length day_max+1 + + dict_p = {} + for i in range(len(p)): + dict_p[i] = p[i] + + table, alias, value = multinomial_iso_rng.alias_init(dict_p) + + # For base case that we do not do testing + if pb == 0: + day_iso = np.zeros(len(p)).astype(int) + gen_iso_num = np.zeros(len(p)).astype(int) + else: + gen_iso = [] + for i in range(n): + gen_iso.append(multinomial_iso_rng.alias(table, alias, value)) + day_iso, gen_iso_num = np.unique(gen_iso, return_counts=True) # Note, day_iso start from today. last entry represents no isolation. + + return day_iso, gen_iso_num + + # Disease schedule generator + def calculate_mul(n, inf_table, inf_alias, inf_value): + gen_inf = [] + for i in range(n): + gen_inf.append(multinomial_inf_rng.alias(inf_table, inf_alias, inf_value)) + day_inf, gen_inf_num = np.unique(gen_inf, return_counts=True) + + return day_inf, gen_inf_num + + # Loop through day 1 - day n-1 + for day in range(0, self.factors['n']): + + if day == 0: + num_exp_non_vac = np.array(infectious_non_vac[0, :], int) + num_exp_vac = np.array(infectious_vac[0, :], int) + + else: + # Update the states from the day before + sus_vac[day, :] += sus_vac[day - 1, :] + sus_non_vac[day, :] += sus_non_vac[day - 1, :] + exposed_vac[day, :] += exposed_vac[day - 1, :] + exposed_non_vac[day, :] += exposed_non_vac[day - 1, :] + infectious_vac[day, :] += infectious_vac[day - 1, :] + infectious_non_vac[day, :] += infectious_non_vac[day - 1, :] + asymptomatic_vac[day, :] += asymptomatic_vac[day - 1, :] + asymptomatic_non_vac[day, :] += asymptomatic_non_vac[day - 1, :] + symptomatic_vac[day, :] += symptomatic_vac[day - 1, :] + symptomatic_non_vac[day, :] += symptomatic_non_vac[day - 1, :] + isolation_exp_non_vac[day, :] += isolation_exp_non_vac[day - 1, :] + isolation_exp_vac[day, :] += isolation_exp_vac[day - 1, :] + isolation_inf_non_vac[day, :] += isolation_inf_non_vac[day - 1, :] + isolation_inf_vac[day, :] += isolation_inf_vac[day - 1, :] + isolation_symp_non_vac[day, :] += isolation_symp_non_vac[day - 1, :] + isolation_symp_vac[day, :] += isolation_symp_vac[day - 1, :] + isolation_asymp_non_vac[day, :] += isolation_asymp_non_vac[day - 1, :] + isolation_asymp_vac[day, :] += isolation_asymp_vac[day - 1, :] + recovered[day, :] += recovered[day - 1, :] + + # Potential new exposed people (include those exposed but healthy and exposed will become infectious later) + new_exp_vac = np.multiply(np.multiply(t_rate_vac, (infectious_vac[day, :] + infectious_non_vac[day, :] + asymptomatic_vac[day, :] + asymptomatic_non_vac[day, :])), (sus_vac[day, :] / (sus_vac[day, :] + sus_non_vac[day, :] + exposed_vac[day, :] + exposed_non_vac[day, :] + infectious_vac[day, :] + infectious_non_vac[day, :] + asymptomatic_vac[day, :] + asymptomatic_non_vac[day, :] + recovered[day, :]))) + new_exp_non_vac = np.multiply(np.multiply(t_rate, (infectious_vac[day, :] + infectious_non_vac[day, :] + asymptomatic_vac[day, :] + asymptomatic_non_vac[day, :])), (sus_non_vac[day, :] / (sus_vac[day, :] + sus_non_vac[day, :] + exposed_vac[day, :] + exposed_non_vac[day, :] + infectious_vac[day, :] + infectious_non_vac[day, :] + asymptomatic_vac[day, :] + asymptomatic_non_vac[day, :] + recovered[day, :]))) + + num_exp_vac = [poisson_numexp_rng.poissonvariate(new_exp_vac[i]) for i in range(self.factors["num_groups"])] # Number of newly exposed people (will have disease) this day + num_exp_non_vac = [poisson_numexp_rng.poissonvariate(new_exp_non_vac[i]) for i in range(self.factors["num_groups"])] + + exposed_vac[day, :] = np.add(exposed_vac[day, :], num_exp_vac) + exposed_non_vac[day, :] = np.add(exposed_non_vac[day, :], num_exp_non_vac) + + sus_vac[day, :] = np.subtract(sus_vac[day, :], num_exp_vac) # Update the suspectible people (healthy) + sus_non_vac[day, :] = np.subtract(sus_non_vac[day, :], num_exp_non_vac) + + # Get vector of testing frequency from decision variable + freq = self.factors["freq_vac"][:self.factors["num_groups"]] + + # For exposed_non_vac: + for g in range(self.factors["num_groups"]): + + if day > 0: + # Exp_num: unscheduled people in state exposed + exp_num = num_exp_non_vac[g] + + day_exp, gen_exp_num = calculate_mul(exp_num, exp_table, exp_alias, exp_value) + + for i in range(len(gen_exp_num)): + n = gen_exp_num[i] # Subgroup of new comings according to their schedule + day_max = day_exp[i] + day_move_iso_exp, num_iso_exp = calculate_iso_mul(day_max, freq[g] * (1 - self.factors["false_neg"]), n) # Calculate isolation schedule for new comings + # The last entry in num_iso_exp represents the number of people not isolated. + for j in range(len(num_iso_exp) - 1): + if day + day_move_iso_exp[j] < self.factors["n"]: + exposed_non_vac[day + day_move_iso_exp[j], g] -= num_iso_exp[j] # Remove from exposed state + isolation_exp_non_vac[day + day_move_iso_exp[j], g] += num_iso_exp[j] # Move to isolatioin exposed + exp_arr_iso[day + day_move_iso_exp[j], g] += num_iso_exp[j] # Track the new comings for isolation exposed(to do disease schedule) + gen_exp_num[i] -= num_iso_exp[j] # Remove them from the free newcoming + tot_num_isolated += num_iso_exp[j] + + # Update for free new exposed people: + if day > 0: + for i in range(len(gen_exp_num)): + if day + day_exp[i] < self.factors["n"]: + infectious_non_vac[day + day_exp[i], g] += gen_exp_num[i] # Move to infectious state + exposed_non_vac[day + day_exp[i], g] -= gen_exp_num[i] # Remove from exposed state + inf_arr[day + day_exp[i], g] += gen_exp_num[i] # Track the new comings for infectious state + + # Schedule for isolated_exposed: + exp_iso_num = int(exp_arr_iso[day, g]) + day_exp_iso, gen_exp_num_iso = calculate_mul(exp_iso_num, exp_table, exp_alias, exp_value) # Disease schedule for isolation_exposed + + if day > 0: + for i in range(len(gen_exp_num_iso)): + if day + day_exp_iso[i] < self.factors["n"]: + isolation_inf_non_vac[day + day_exp_iso[i], g] += gen_exp_num_iso[i] # Move to isolation_infectious + isolation_exp_non_vac[day + day_exp_iso[i], g] -= gen_exp_num_iso[i] # Remove from isolation_exposed + inf_arr_iso[day + day_exp_iso[i], g] += gen_exp_num_iso[i] # Track the new comings for isolation_infectious + + # Get the number of new comings in infectious/symp/asymp states: + if day == 0: + inf_num = int(infectious_non_vac[day, g]) + asym_num = 0 + sym_num = 0 + # We do not do test at day0 + + else: + inf_num = int(inf_arr[day, g]) # The number of new comings for infectious state + asym_num = int(asym_arr[day, g]) + sym_num = int(sym_arr[day, g]) + + # Do schedule for the free unscheduled people in infectious state + day_inf, gen_inf_num = calculate_mul(inf_num, inf_table, inf_alias, inf_value) + + # Do test and isolation for new coming people: + for i in range(len(gen_inf_num)): + # The number of new coming individuals that will become infectious in day_max days + n = gen_inf_num[i] + day_max = day_inf[i] + # Assume the day of leaving state do not do test at this state + day_move_iso_inf, num_iso_inf = calculate_iso_mul(day_max, freq[g] * (1 - self.factors["false_neg"]), n) # Store the isolation schedule for this group of new comings + for j in range(len(num_iso_inf) - 1): # Move them to isolation at corrosponding days + if day + day_move_iso_inf[j] < self.factors["n"]: + infectious_non_vac[day + day_move_iso_inf[j], g] -= num_iso_inf[j] + isolation_inf_non_vac[day + day_move_iso_inf[j], g] += num_iso_inf[j] + inf_arr_iso[day + day_move_iso_inf[j], g] += num_iso_inf[j] # Track the number of new isolated + gen_inf_num[i] -= num_iso_inf[j] # Remove them from the free newcoming of day i subgroup + tot_num_isolated += num_iso_inf[j] + + # Move the number of people to the symp/asymp at correponding days + for i in range(len(gen_inf_num)): + if day + day_inf[i] < self.factors["n"]: + if gen_inf_num[i] != 0: + # Generate the number of asymptomatic people + asym_num_inf = binomial_asymp_rng.binomialvariate(int(gen_inf_num[i]), self.factors["asymp_rate"]) + sym_num_inf = gen_inf_num[i] - asym_num_inf + symptomatic_non_vac[day + day_inf[i], g] += sym_num_inf # Move to symptomatic state + asymptomatic_non_vac[day + day_inf[i], g] += asym_num_inf # Move to asymptomatic state + sym_arr[day + day_inf[i], g] += sym_num_inf + asym_arr[day + day_inf[i], g] += asym_num_inf + tot_num_syptomatic += sym_num_inf # Count the total symptomatic number + infectious_non_vac[day + day_inf[i], g] -= gen_inf_num[i] # Remove from infectious state + + # Schedule for infectious isolated + if day > 0: + inf_iso_num = int(inf_arr_iso[day, g]) + day_inf_iso, gen_inf_num_iso = calculate_mul(inf_iso_num, inf_table, inf_alias, inf_value) + + # If we only update the next day condition and then check the next + for i in range(len(gen_inf_num_iso)): + if day + day_inf_iso[i] < self.factors["n"]: + if gen_inf_num_iso[i] != 0: + asym_num_inf = binomial_asymp_rng.binomialvariate(int(gen_inf_num_iso[i]), self.factors["asymp_rate"]) # Number of people will go to isolation_asymp + sym_num_inf = gen_inf_num_iso[i] - asym_num_inf + symptomatic_non_vac[day + day_inf_iso[i], g] += sym_num_inf # All symptomatic people are isolated + isolation_asymp_non_vac[day + day_inf_iso[i], g] += asym_num_inf + sym_arr[day + day_inf_iso[i], g] += sym_num_inf # All symptomatic people are isolated + asym_arr_iso[day + day_inf_iso[i], g] += asym_num_inf + tot_num_syptomatic += sym_num_inf + isolation_inf_non_vac[day + day_inf_iso[i], g] -= gen_inf_num_iso[i] + + # For symptomatic: + day_sym, gen_sym_num = calculate_mul(sym_num, sym_table, sym_alias, sym_value) + + for t in range(len(gen_sym_num)): + if day + day_sym[t] < self.factors["n"]: + symptomatic_non_vac[day + day_sym[t], g] -= gen_sym_num[t] + recovered[day + day_sym[t], g] += gen_sym_num[t] + + # For asymptomatic: + day_asym, gen_asym_num = calculate_mul(asym_num, asym_table, asym_alias, asym_value) + + for i in range(len(gen_asym_num)): + n = gen_asym_num[i] + day_max = day_asym[i] + day_move_iso_asym, num_iso_asym = calculate_iso_mul(day_max, freq[g] * (1 - self.factors["false_neg"]), n) + for j in range(len(num_iso_asym) - 1): + if day + day_move_iso_asym[j] < self.factors["n"]: + asymptomatic_non_vac[day + day_move_iso_asym[j], g] -= num_iso_asym[j] + isolation_asymp_non_vac[day + day_move_iso_asym[j], g] += num_iso_asym[j] + asym_arr_iso[day + day_move_iso_asym[j], g] += num_iso_asym[j] + gen_asym_num[i] -= num_iso_asym[j] + tot_num_isolated += num_iso_asym[j] + + for t in range(len(gen_asym_num)): + if day + day_asym[t] < self.factors["n"]: + asymptomatic_non_vac[day + day_asym[t], g] -= gen_asym_num[t] + recovered[day + day_asym[t], g] += gen_asym_num[t] + + # For asymptomatic isolation: + if day > 0: + asym_num_iso = int(asym_arr_iso[day, g]) + day_asym_iso, gen_asym_num_iso = calculate_mul(asym_num_iso, asym_table, asym_alias, asym_value) + + for t in range(len(gen_asym_num_iso)): + if day + day_asym_iso[t] < self.factors["n"]: + isolation_asymp_non_vac[day + day_asym_iso[t], g] -= gen_asym_num_iso[t] + recovered[day + day_asym_iso[t], g] += gen_asym_num_iso[t] + + # For vaccinated group: + # For exposed_vac: + for g in range(self.factors["num_groups"]): + + if day > 0: + # Exp_num: unscheduled people in state exposed + exp_num = num_exp_vac[g] + day_exp, gen_exp_num = calculate_mul(exp_num, exp_table, exp_alias, exp_value) + + for i in range(len(gen_exp_num)): + n = gen_exp_num[i] # Subgroup of new comings according to their schedule + day_max = day_exp[i] + day_move_iso_exp, num_iso_exp = calculate_iso_mul(day_max, freq[g] * (1 - self.factors["false_neg"]), n) # Calculate isolation schedule for new comings + # The last entry in num_iso_exp represents the number of people not isolated. + for j in range(len(num_iso_exp) - 1): + if day + day_move_iso_exp[j] < self.factors["n"]: + exposed_vac[day + day_move_iso_exp[j], g] -= num_iso_exp[j] #Remove from exposed state + isolation_exp_vac[day + day_move_iso_exp[j], g] += num_iso_exp[j] # Move to isolatioin exposed + exp_arr_vac_iso[day + day_move_iso_exp[j], g] += num_iso_exp[j] + gen_exp_num[i] -= num_iso_exp[j] # Remove them from the free newcoming + tot_num_isolated += num_iso_exp[j] + + # Update for free new exposed people: + if day > 0: + for i in range(len(gen_exp_num)): + if day + day_exp[i] < self.factors["n"]: + infectious_vac[day + day_exp[i], g] += gen_exp_num[i] + exposed_vac[day + day_exp[i], g] -= gen_exp_num[i] + inf_arr_vac[day + day_exp[i], g] += gen_exp_num[i] + + # Schedule for isolated_exposed: + exp_iso_num = int(exp_arr_vac_iso[day, g]) + day_exp_iso, gen_exp_num_iso = calculate_mul(exp_iso_num, exp_table, exp_alias, exp_value) + + if day > 0: + for i in range(len(gen_exp_num_iso)): + if day + day_exp_iso[i] < self.factors["n"]: + isolation_inf_vac[day + day_exp_iso[i], g] += gen_exp_num_iso[i] + isolation_exp_vac[day + day_exp_iso[i], g] -= gen_exp_num_iso[i] + inf_arr_vac_iso[day + day_exp_iso[i], g] += gen_exp_num_iso[i] + + # For infectious: + if day == 0: + inf_num = int(infectious_vac[day, g]) + asym_num = 0 + sym_num = 0 + # We do not do test at day 0 + + else: + inf_num = int(inf_arr_vac[day, g]) + asym_num = int(asym_arr_vac[day, g]) + sym_num = int(sym_arr_vac[day, g]) + + # Do schedule for the free unscheduled people in infectious state + day_inf, gen_inf_num = calculate_mul(inf_num, inf_table, inf_alias, inf_value) + + # Do test and isolation for new coming people: + for i in range(len(gen_inf_num)): + # The number of new coming individuals that will become infectious in day_max days + n = gen_inf_num[i] + day_max = day_inf[i] + # Assume the day of leaving state do not do test at this state + day_move_iso_inf, num_iso_inf = calculate_iso_mul(day_max, freq[g] * (1 - self.factors["false_neg"]), n) # Store the isolation schedule for this group of new comings + for j in range(len(num_iso_inf) - 1): # Move them to isolation at corrosponding days + if day + day_move_iso_inf[j] < self.factors["n"]: + infectious_vac[day + day_move_iso_inf[j], g] -= num_iso_inf[j] + isolation_inf_vac[day + day_move_iso_inf[j], g] += num_iso_inf[j] + inf_arr_vac_iso[day + day_move_iso_inf[j], g] += num_iso_inf[j] # Track the number of new isolated + gen_inf_num[i] -= num_iso_inf[j] # Remove them from the free newcoming of day i subgroup + tot_num_isolated += num_iso_inf[j] + + # Move the number of people to the symp/asymp at correponding days + for i in range(len(gen_inf_num)): + if day + day_inf[i] < self.factors["n"]: + if gen_inf_num[i] != 0: + # Generate the number of asymptomatic people + asym_num_inf = binomial_asymp_rng.binomialvariate(int(gen_inf_num[i]), self.factors["asymp_rate_vac"]) + sym_num_inf = gen_inf_num[i] - asym_num_inf + symptomatic_vac[day + day_inf[i], g] += sym_num_inf + asymptomatic_vac[day + day_inf[i], g] += asym_num_inf + sym_arr_vac[day + day_inf[i], g] += sym_num_inf + asym_arr_vac[day + day_inf[i], g] += asym_num_inf + tot_num_syptomatic += sym_num_inf + infectious_vac[day + day_inf[i], g] -= gen_inf_num[i] + + # Schedule for infectious isolated + if day > 0: + inf_iso_num = int(inf_arr_vac_iso[day, g]) + day_inf_iso, gen_inf_num_iso = calculate_mul(inf_iso_num, inf_table, inf_alias, inf_value) + + # If we only update the next day condition and then check the next + for i in range(len(gen_inf_num_iso)): + if day + day_inf_iso[i] < self.factors["n"]: + if gen_inf_num_iso[i] != 0: + asym_num_inf = binomial_asymp_rng.binomialvariate(int(gen_inf_num_iso[i]), self.factors["asymp_rate_vac"]) + sym_num_inf = gen_inf_num_iso[i] - asym_num_inf + symptomatic_vac[day + day_inf_iso[i], g] += sym_num_inf # All symptomatic people are isolated + isolation_asymp_vac[day + day_inf_iso[i], g] += asym_num_inf + sym_arr_vac[day + day_inf_iso[i], g] += sym_num_inf # All symptomatic people are isolated + asym_arr_vac_iso[day + day_inf_iso[i], g] += asym_num_inf + tot_num_syptomatic += sym_num_inf + isolation_inf_vac[day + day_inf_iso[i], g] -= gen_inf_num_iso[i] + + # For symptomatic: + day_sym, gen_sym_num = calculate_mul(sym_num, sym_table, sym_alias, sym_value) + + for t in range(len(gen_sym_num)): + if day + day_sym[t] < self.factors["n"]: + symptomatic_vac[day + day_sym[t], g] -= gen_sym_num[t] + recovered[day + day_sym[t], g] += gen_sym_num[t] + + # For asymptomatic: + day_asym, gen_asym_num = calculate_mul(asym_num, asym_table, asym_alias, asym_value) + + for i in range(len(gen_asym_num)): + n = gen_asym_num[i] + day_max = day_asym[i] + day_move_iso_asym, num_iso_asym = calculate_iso_mul(day_max, freq[g] * (1 - self.factors["false_neg"]), n) + for j in range(len(num_iso_asym) - 1): + if day + day_move_iso_asym[j] < self.factors["n"]: + asymptomatic_vac[day + day_move_iso_asym[j], g] -= num_iso_asym[j] + isolation_asymp_vac[day + day_move_iso_asym[j], g] += num_iso_asym[j] + asym_arr_vac_iso[day + day_move_iso_asym[j], g] += num_iso_asym[j] + gen_asym_num[i] -= num_iso_asym[j] + tot_num_isolated += num_iso_asym[j] + + for t in range(len(gen_asym_num)): + if day + day_asym[t] < self.factors["n"]: + asymptomatic_vac[day + day_asym[t], g] -= gen_asym_num[t] + recovered[day + day_asym[t], g] += gen_asym_num[t] + + # For asymptomatic isolation: + if day > 0: + asym_num_iso = int(asym_arr_vac_iso[day, g]) + day_asym_iso, gen_asym_num_iso = calculate_mul(asym_num_iso, asym_table, asym_alias, asym_value) + + for t in range(len(gen_asym_num_iso)): + if day + day_asym_iso[t] < self.factors["n"]: + isolation_asymp_vac[day + day_asym_iso[t], g] -= gen_asym_num_iso[t] + recovered[day + day_asym_iso[t], g] += gen_asym_num_iso[t] + + # Update performance measures + for g in range(self.factors["num_groups"]): + num_exposed[day][g] = np.sum(exposed_non_vac[day, g] + exposed_vac[day, g] + isolation_exp_non_vac[day, g] + isolation_exp_vac[day, g]) + num_susceptible[day][g] = np.sum(sus_vac[day, g] + sus_non_vac[day, g]) + num_symptomatic[day][g] = np.sum(symptomatic_vac[day, g] + symptomatic_non_vac[day, g]) # All symp people are isolated. isolated_symp = symp + num_recovered[day][g] = np.sum(recovered[day, g]) + num_infected[day][g] = np.sum(infectious_non_vac[day, g] + infectious_vac[day, g] + symptomatic_non_vac[day, g] + symptomatic_vac[day, g] + asymptomatic_non_vac[day, g] + asymptomatic_vac[day, g] + + isolation_inf_non_vac[day, g] + isolation_inf_vac[day, g] + isolation_asymp_non_vac[day, g] + isolation_asymp_vac[day, g]) + num_isolation[day][g] = np.sum(isolation_exp_non_vac[day, g] + isolation_inf_non_vac[day, g] + symptomatic_non_vac[day, g] + isolation_asymp_non_vac[day, g] + + isolation_exp_vac[day, g] + isolation_inf_vac[day, g] + symptomatic_vac[day, g] + isolation_asymp_vac[day, g]) + + # Compose responses and gradients. + last_susceptible = num_susceptible[self.factors["n"] - 1] # Number of suspectible(unaffected) people at the end of the period + inf = self.factors["group_size"] - last_susceptible # Number of total infected people during the period + w = np.array(self.factors["w"]) # Protection-rate for each group, here we heavy scale for the group3 as example w=[1,1,5]; if no scaling, set w=np.ones + tot_infected = np.dot(w, inf) # After adjusting for the protection rate, the weighted objective value (adjusted infected) + + responses = {"num_infected": num_infected, + "num_exposed": num_exposed, + "num_susceptible": num_susceptible, + "num_recovered": num_recovered, + "num_symptomatic": num_symptomatic, + "free_population": last_susceptible, + "num_isolation": num_isolation, + "total_symp": tot_num_syptomatic, + "total_isolated": tot_num_isolated + tot_num_syptomatic, + "total_infected": tot_infected} + gradients = {response_key: + {factor_key: np.nan for factor_key in self.specifications} + for response_key in responses + } + return responses, gradients + + +""" +Summary +------- +Minimize the average number of daily infected people. +""" + + +class CovidMinInfectVac(Problem): + """ + Base class to implement simulation-optimization problems. + + Attributes + ---------- + name : string + name of problem + dim : int + number of decision variables + n_objectives : int + number of objectives + n_stochastic_constraints : int + number of stochastic constraints + minmax : tuple of int (+/- 1) + indicator of maximization (+1) or minimization (-1) for each objective + constraint_type : string + description of constraints types: + "unconstrained", "box", "deterministic", "stochastic" + variable_type : string + description of variable types: + "discrete", "continuous", "mixed" + lower_bounds : tuple + lower bound for each decision variable + upper_bounds : tuple + upper bound for each decision variable + Ci : ndarray (or None) + Coefficient matrix for linear inequality constraints of the form Ci@x <= di + Ce : ndarray (or None) + Coefficient matrix for linear equality constraints of the form Ce@x = de + di : ndarray (or None) + Constraint vector for linear inequality constraints of the form Ci@x <= di + de : ndarray (or None) + Constraint vector for linear equality constraints of the form Ce@x = de + gradient_available : bool + indicates if gradient of objective function is available + optimal_value : float + optimal objective function value + optimal_solution : tuple + optimal solution + model : Model object + associated simulation model that generates replications + model_default_factors : dict + default values for overriding model-level default factors + model_fixed_factors : dict + combination of overriden model-level factors and defaults + model_decision_factors : set of str + set of keys for factors that are decision variables + rng_list : list of rng.MRG32k3a objects + list of RNGs used to generate a random initial solution + or a random problem instance + factors : dict + changeable factors of the problem + initial_solution : tuple + default initial solution from which solvers start + budget : int > 0 + max number of replications (fn evals) for a solver to take + specifications : dict + details of each factor (for GUI, data validation, and defaults) + + Arguments + --------- + name : str + user-specified name for problem + fixed_factors : dict + dictionary of user-specified problem factors + model_fixed factors : dict + subset of user-specified non-decision factors to pass through to the model + + See also + -------- + base.Problem + """ + def __init__(self, name="COVID-vac-test-1", fixed_factors={}, model_fixed_factors={}): + self.name = name + self.n_objectives = 1 + self.n_stochastic_constraints = 0 + self.minmax = (-1,) + self.constraint_type = "box" + self.variable_type = "continuous" + self.gradient_available = True + self.optimal_value = None + self.optimal_solution = None + self.model_default_factors = {} + self.model_fixed_factors = {} + self.model_decision_factors = {"freq_vac"} + self.factors = fixed_factors + self.specifications = { + "initial_solution": { + "description": "Initial solution from which solvers start.", + "datatype": tuple, + # "default": (0/7, 0/7, 0/7, 0/7, 0/7, 0/7, 0/7, 0, 0, 0, 0, 0, 0, 0) #for 7 groups + "default": (0/7, 0/7, 0/7, 0, 0, 0) + }, + "budget": { + "description": "Max # of replications for a solver to take.", + "datatype": int, + "default": 600 + }, + "vaccine_cap": { + "description": "Maximum number of vaccines.", + "datatype": int, + "default": 8000 + }, + "testing_cap": { + "description": "Maximum testing capacity per day.", + "datatype": int, + "default": 5000 + }, + "pen_coef": { + "description": "Penalty coefficient.", + "datatype": int, + "default": 100000 + }, + } + self.check_factor_list = { + "initial_solution": self.check_initial_solution, + "budget": self.check_budget, + "testing_cap": self.check_testing_cap, + "vaccine_cap": self.check_vaccine_cap, + "pen_coef": self.check_pen_coef + } + super().__init__(fixed_factors, model_fixed_factors) + # Instantiate model with fixed factors and overwritten defaults. + self.model = COVID_vac(self.model_fixed_factors) + self.dim = len(self.model.factors["group_size"]) + self.lower_bounds = (0,) * self.dim + self.upper_bounds = (1,) * self.dim + self.Ci = None + self.Ce = None + self.di = None + self.de = None + + def check_vaccine_cap(self): + return self.factors["vaccine_cap"] > 0 + + def check_testing_cap(self): + return self.factors["testing_cap"] > 0 + + def check_pen_coef(self): + return self.factors["pen_coef"] > 0 + + def vector_to_factor_dict(self, vector): + """ + Convert a vector of variables to a dictionary with factor keys + + Arguments + --------- + vector : tuple + vector of values associated with decision variables + + Returns + ------- + factor_dict : dictionary + dictionary with factor keys and associated values + """ + factor_dict = { + "freq_vac": vector + } + return factor_dict + + def factor_dict_to_vector(self, factor_dict): + """ + Convert a dictionary with factor keys to a vector + of variables. + + Arguments + --------- + factor_dict : dictionary + dictionary with factor keys and associated values + + Returns + ------- + vector : tuple + vector of values associated with decision variables + """ + vector = (factor_dict["freq_vac"], ) + return vector + + def response_dict_to_objectives(self, response_dict): + """ + Convert a dictionary with response keys to a vector + of objectives. + + Arguments + --------- + response_dict : dictionary + dictionary with response keys and associated values + + Returns + ------- + objectives : tuple + vector of objectives + """ + objectives = (response_dict["total_symp"], ) + return objectives + + def response_dict_to_stoch_constraints(self, response_dict): + """ + Convert a dictionary with response keys to a vector + of left-hand sides of stochastic constraints: E[Y] >= 0 + + Arguments + --------- + response_dict : dictionary + dictionary with response keys and associated values + + Returns + ------- + stoch_constraints : tuple + vector of LHSs of stochastic constraint + """ + stoch_constraints = None + return stoch_constraints + + def deterministic_objectives_and_gradients(self, x): + """ + Compute deterministic components of objectives for a solution `x`. + + Arguments + --------- + x : tuple + vector of decision variables + + Returns + ------- + det_objectives : tuple + vector of deterministic components of objectives + det_objectives_gradients : tuple + vector of gradients of deterministic components of objectives + """ + # Checking: + # print(self.factors['pen_coef']) + # for g in range(self.dim): + # print(np.dot(self.model.factors["group_size"][g],x[g])) + # print(np.sum(np.dot(self.model.factors["group_size"][g],x[g]) for g in range(self.dim))) + # print(max(np.sum(np.dot(self.model.factors["group_size"][g],x[g]) for g in range(self.dim)) - self.factors["vaccine_cap"],0)) + + x_v = x[self.dim: self.dim * 2] + x_t = x[: self.dim] + det_objectives_test = (self.factors["pen_coef"] * (max(np.sum(np.ceil(np.dot(self.model.factors["group_size"][g], x_t[g])) for g in range(self.dim)) - self.factors["testing_cap"], 0)),) + det_objectives_vac = (self.factors["pen_coef"] * (max(np.sum(np.ceil(np.dot(self.model.factors["group_size"][g], x_v[g])) for g in range(self.dim)) - self.factors["vaccine_cap"], 0)),) + det_objectives = max(det_objectives_test, det_objectives_vac) + det_objectives_gradients = ((0,),) + return det_objectives, det_objectives_gradients + + def deterministic_stochastic_constraints_and_gradients(self, x): + """ + Compute deterministic components of stochastic constraints + for a solution `x`. + + Arguments + --------- + x : tuple + vector of decision variables + + Returns + ------- + det_stoch_constraints : tuple + vector of deterministic components of stochastic constraints + det_stoch_constraints_gradients : tuple + vector of gradients of deterministic components of + stochastic constraints + """ + det_stoch_constraints = None + det_stoch_constraints_gradients = None + return det_stoch_constraints, det_stoch_constraints_gradients + + def check_deterministic_constraints(self, x): + """ + Check if a solution `x` satisfies the problem's deterministic + constraints. + + Arguments + --------- + x : tuple + vector of decision variables + + Returns + ------- + satisfies : bool + indicates if solution `x` satisfies the deterministic constraints. + """ + return np.all(x >= 0) + + def get_random_solution(self, rand_sol_rng): + """ + Generate a random solution for starting or restarting solvers. + + Arguments + --------- + rand_sol_rng : rng.MRG32k3a object + random-number generator used to sample a new random solution + + Returns + ------- + x : tuple + vector of decision variables + """ + # If uniformly find freq & vac: + # x = tuple([rand_sol_rng.uniform(0,1) for _ in range(self.dim)]) + # f = tuple([rand_sol_rng.uniform(0,1) for _ in range(self.dim)]) + # dx = f + x + + # By redistribution method: (reallocate the resource(vac/freq) if some group have more resource than their population size) + # Linear constraint + # For vaccine: + x = tuple(rand_sol_rng.integer_random_vector_from_simplex(self.model.factors["num_groups"], self.factors["vaccine_cap"])) + v = np.divide(x, self.model.factors["group_size"]) + + x = np.array(x) # For the temporal calculation + v = np.array(v) + + n = len(self.model.factors["group_size"]) + while True: + xn = x - np.array(self.model.factors["group_size"]) + idx = np.argmax(xn) + if v[idx] > 1: + remain = x[idx] - self.model.factors["group_size"][idx] + x[idx] -= remain + x += int(np.floor(remain / n)) + short_idx = np.argmax(np.array(self.model.factors["group_size"]) - x) + rr = remain - int(np.floor(remain / 3) * n) + if rr > 0: + x[short_idx] += rr + v = np.divide(x, self.model.factors["group_size"]) + else: + break + + x = tuple(x) + v = tuple(v) + + # If modify testing freq by redistribution, for testing: + xf = tuple(rand_sol_rng.integer_random_vector_from_simplex(self.model.factors["num_groups"], self.factors["testing_cap"])) + f = np.divide(xf, self.model.factors["group_size"]) + + xf = np.array(xf) # For the temporal calculation + f = np.array(f) + + n = len(self.model.factors["group_size"]) + while True: + xn = xf - np.array(self.model.factors["group_size"]) + idx = np.argmax(xn) + if f[idx] > 1: + remain = xf[idx] - self.model.factors["group_size"][idx] + xf[idx] -= remain + xf += int(np.floor(remain / n)) + short_idx = np.argmax(np.array(self.model.factors["group_size"]) - xf) + rr = remain - int(np.floor(remain / 3) * n) + if rr > 0: + xf[short_idx] += rr + f = np.divide(xf, self.model.factors["group_size"]) + else: + break + + xf = tuple(xf) + f = tuple(f) + + print('initial f, v: ', f, v) + + # Final decision variable freq_vac: + dx = f + v + + return dx + + +# Other problem classes: +# Sub-problem 2: problem with more groups. +class CovidMinInfectVac2(Problem): + """ + Base class to implement simulation-optimization problems. + In this problem classes, consider there are 7 groups in total. + Thus, to study this problem class, one need to modify the model part default values to those corresponding to 7 groups. + + Attributes + ---------- + name : string + name of problem + dim : int + number of decision variables + n_objectives : int + number of objectives + n_stochastic_constraints : int + number of stochastic constraints + minmax : tuple of int (+/- 1) + indicator of maximization (+1) or minimization (-1) for each objective + constraint_type : string + description of constraints types: + "unconstrained", "box", "deterministic", "stochastic" + variable_type : string + description of variable types: + "discrete", "continuous", "mixed" + lower_bounds : tuple + lower bound for each decision variable + upper_bounds : tuple + upper bound for each decision variable + Ci : ndarray (or None) + Coefficient matrix for linear inequality constraints of the form Ci@x <= di + Ce : ndarray (or None) + Coefficient matrix for linear equality constraints of the form Ce@x = de + di : ndarray (or None) + Constraint vector for linear inequality constraints of the form Ci@x <= di + de : ndarray (or None) + Constraint vector for linear equality constraints of the form Ce@x = de + gradient_available : bool + indicates if gradient of objective function is available + optimal_value : float + optimal objective function value + optimal_solution : tuple + optimal solution + model : Model object + associated simulation model that generates replications + model_default_factors : dict + default values for overriding model-level default factors + model_fixed_factors : dict + combination of overriden model-level factors and defaults + model_decision_factors : set of str + set of keys for factors that are decision variables + rng_list : list of rng.MRG32k3a objects + list of RNGs used to generate a random initial solution + or a random problem instance + factors : dict + changeable factors of the problem + initial_solution : tuple + default initial solution from which solvers start + budget : int > 0 + max number of replications (fn evals) for a solver to take + specifications : dict + details of each factor (for GUI, data validation, and defaults) + + Arguments + --------- + name : str + user-specified name for problem + fixed_factors : dict + dictionary of user-specified problem factors + model_fixed factors : dict + subset of user-specified non-decision factors to pass through to the model + + See also + -------- + base.Problem + """ + def __init__(self, name="COVID-vac-test-2", fixed_factors={}, model_fixed_factors={}): + self.name = name + self.n_objectives = 1 + self.n_stochastic_constraints = 0 + self.minmax = (-1,) + self.constraint_type = "box" + self.variable_type = "continuous" + self.gradient_available = True + self.optimal_value = None + self.optimal_solution = None + self.model_default_factors = {} + self.model_fixed_factors = {} + self.model_decision_factors = {"freq_vac"} + self.factors = fixed_factors + self.specifications = { + "initial_solution": { + "description": "Initial solution from which solvers start.", + "datatype": tuple, + "default": (0/7, 0/7, 0/7, 0/7, 0/7, 0/7, 0/7, 0, 0, 0, 0, 0, 0, 0) #for 7 groups + }, + "budget": { + "description": "Max # of replications for a solver to take.", + "datatype": int, + "default": 600 + }, + "vaccine_cap": { + "description": "Maximum number of vaccines.", + "datatype": int, + "default": 8000 + }, + "testing_cap": { + "description": "Maximum testing capacity per day.", + "datatype": int, + "default": 5000 + }, + "pen_coef": { + "description": "Penalty coefficient.", + "datatype": int, + "default": 100000 + }, + } + self.check_factor_list = { + "initial_solution": self.check_initial_solution, + "budget": self.check_budget, + "testing_cap": self.check_testing_cap, + "vaccine_cap": self.check_vaccine_cap, + "pen_coef": self.check_pen_coef + } + super().__init__(fixed_factors, model_fixed_factors) + # Instantiate model with fixed factors and overwritten defaults. + self.model = COVID_vac(self.model_fixed_factors) + self.dim = len(self.model.factors["group_size"]) + self.lower_bounds = (0,) * self.dim + self.upper_bounds = (1,) * self.dim + self.Ci = None + self.Ce = None + self.di = None + self.de = None + + def check_vaccine_cap(self): + return self.factors["vaccine_cap"] > 0 + + def check_testing_cap(self): + return self.factors["testing_cap"] > 0 + + def check_pen_coef(self): + return self.factors["pen_coef"] > 0 + + def vector_to_factor_dict(self, vector): + """ + Convert a vector of variables to a dictionary with factor keys + + Arguments + --------- + vector : tuple + vector of values associated with decision variables + + Returns + ------- + factor_dict : dictionary + dictionary with factor keys and associated values + """ + factor_dict = { + "freq_vac": vector + } + return factor_dict + + def factor_dict_to_vector(self, factor_dict): + """ + Convert a dictionary with factor keys to a vector + of variables. + + Arguments + --------- + factor_dict : dictionary + dictionary with factor keys and associated values + + Returns + ------- + vector : tuple + vector of values associated with decision variables + """ + vector = (factor_dict["freq_vac"], ) + return vector + + def response_dict_to_objectives(self, response_dict): + """ + Convert a dictionary with response keys to a vector + of objectives. + + Arguments + --------- + response_dict : dictionary + dictionary with response keys and associated values + + Returns + ------- + objectives : tuple + vector of objectives + """ + objectives = (response_dict["total_infected"], ) + return objectives + + def response_dict_to_stoch_constraints(self, response_dict): + """ + Convert a dictionary with response keys to a vector + of left-hand sides of stochastic constraints: E[Y] >= 0 + + Arguments + --------- + response_dict : dictionary + dictionary with response keys and associated values + + Returns + ------- + stoch_constraints : tuple + vector of LHSs of stochastic constraint + """ + stoch_constraints = None + return stoch_constraints + + def deterministic_objectives_and_gradients(self, x): + """ + Compute deterministic components of objectives for a solution `x`. + + Arguments + --------- + x : tuple + vector of decision variables + + Returns + ------- + det_objectives : tuple + vector of deterministic components of objectives + det_objectives_gradients : tuple + vector of gradients of deterministic components of objectives + """ + # Checking: + # print(self.factors['pen_coef']) + # for g in range(self.dim): + # print(np.dot(self.model.factors["group_size"][g],x[g])) + # print(np.sum(np.dot(self.model.factors["group_size"][g],x[g]) for g in range(self.dim))) + # print(max(np.sum(np.dot(self.model.factors["group_size"][g],x[g]) for g in range(self.dim)) - self.factors["vaccine_cap"],0)) + + x_v = x[self.dim: self.dim * 2] + x_t = x[: self.dim] + det_objectives_test = (self.factors["pen_coef"] * (max(np.sum(np.ceil(np.dot(self.model.factors["group_size"][g], x_t[g])) for g in range(self.dim)) - self.factors["testing_cap"], 0)),) + det_objectives_vac = (self.factors["pen_coef"] * (max(np.sum(np.ceil(np.dot(self.model.factors["group_size"][g], x_v[g])) for g in range(self.dim)) - self.factors["vaccine_cap"], 0)),) + det_objectives = max(det_objectives_test, det_objectives_vac) + det_objectives_gradients = ((0,),) + return det_objectives, det_objectives_gradients + + def deterministic_stochastic_constraints_and_gradients(self, x): + """ + Compute deterministic components of stochastic constraints + for a solution `x`. + + Arguments + --------- + x : tuple + vector of decision variables + + Returns + ------- + det_stoch_constraints : tuple + vector of deterministic components of stochastic constraints + det_stoch_constraints_gradients : tuple + vector of gradients of deterministic components of + stochastic constraints + """ + det_stoch_constraints = None + det_stoch_constraints_gradients = None + return det_stoch_constraints, det_stoch_constraints_gradients + + def check_deterministic_constraints(self, x): + """ + Check if a solution `x` satisfies the problem's deterministic + constraints. + + Arguments + --------- + x : tuple + vector of decision variables + + Returns + ------- + satisfies : bool + indicates if solution `x` satisfies the deterministic constraints. + """ + return np.all(x >= 0) + + def get_random_solution(self, rand_sol_rng): + """ + Generate a random solution for starting or restarting solvers. + + Arguments + --------- + rand_sol_rng : rng.MRG32k3a object + random-number generator used to sample a new random solution + + Returns + ------- + x : tuple + vector of decision variables + """ + # If uniformly find freq & vac: + # x = tuple([rand_sol_rng.uniform(0,1) for _ in range(self.dim)]) + # f = tuple([rand_sol_rng.uniform(0,1) for _ in range(self.dim)]) + # dx = f + x + + # By redistribution method: (reallocate the resource(vac/freq) if some group have more resource than their population size) + # Linear constraint + # For vaccine: + x = tuple(rand_sol_rng.integer_random_vector_from_simplex(self.model.factors["num_groups"], self.factors["vaccine_cap"])) + v = np.divide(x, self.model.factors["group_size"]) + + x = np.array(x) # For the temporal calculation + v = np.array(v) + + n = len(self.model.factors["group_size"]) + while True: + xn = x - np.array(self.model.factors["group_size"]) + idx = np.argmax(xn) + if v[idx] > 1: + remain = x[idx] - self.model.factors["group_size"][idx] + x[idx] -= remain + x += int(np.floor(remain / n)) + short_idx = np.argmax(np.array(self.model.factors["group_size"]) - x) + rr = remain - int(np.floor(remain / 3) * n) + if rr > 0: + x[short_idx] += rr + v = np.divide(x, self.model.factors["group_size"]) + else: + break + + x = tuple(x) + v = tuple(v) + + # If modify testing freq by redistribution, for testing: + xf = tuple(rand_sol_rng.integer_random_vector_from_simplex(self.model.factors["num_groups"], self.factors["testing_cap"])) + f = np.divide(xf, self.model.factors["group_size"]) + + xf = np.array(xf) # For the temporal calculation + f = np.array(f) + + n = len(self.model.factors["group_size"]) + while True: + xn = xf - np.array(self.model.factors["group_size"]) + idx = np.argmax(xn) + if f[idx] > 1: + remain = xf[idx] - self.model.factors["group_size"][idx] + xf[idx] -= remain + xf += int(np.floor(remain / n)) + short_idx = np.argmax(np.array(self.model.factors["group_size"]) - xf) + rr = remain - int(np.floor(remain / 3) * n) + if rr > 0: + xf[short_idx] += rr + f = np.divide(xf, self.model.factors["group_size"]) + else: + break + + xf = tuple(xf) + f = tuple(f) + + print('initial f, v: ', f, v) + + # Final decision variable freq_vac: + dx = f + v + + return dx + + +# Sub-problem 3: problem that consider different protection weight for different groups. +class CovidMinInfectVac3(Problem): + """ + Base class to implement simulation-optimization problems. + This problem studies the case when different groups are influenced by the disease to various degree. + To study this problem, one need to modify the factor "w" in model classes to assign a protection weight for each group. + + Attributes + ---------- + name : string + name of problem + dim : int + number of decision variables + n_objectives : int + number of objectives + n_stochastic_constraints : int + number of stochastic constraints + minmax : tuple of int (+/- 1) + indicator of maximization (+1) or minimization (-1) for each objective + constraint_type : string + description of constraints types: + "unconstrained", "box", "deterministic", "stochastic" + variable_type : string + description of variable types: + "discrete", "continuous", "mixed" + lower_bounds : tuple + lower bound for each decision variable + upper_bounds : tuple + upper bound for each decision variable + Ci : ndarray (or None) + Coefficient matrix for linear inequality constraints of the form Ci@x <= di + Ce : ndarray (or None) + Coefficient matrix for linear equality constraints of the form Ce@x = de + di : ndarray (or None) + Constraint vector for linear inequality constraints of the form Ci@x <= di + de : ndarray (or None) + Constraint vector for linear equality constraints of the form Ce@x = de + gradient_available : bool + indicates if gradient of objective function is available + optimal_value : float + optimal objective function value + optimal_solution : tuple + optimal solution + model : Model object + associated simulation model that generates replications + model_default_factors : dict + default values for overriding model-level default factors + model_fixed_factors : dict + combination of overriden model-level factors and defaults + model_decision_factors : set of str + set of keys for factors that are decision variables + rng_list : list of rng.MRG32k3a objects + list of RNGs used to generate a random initial solution + or a random problem instance + factors : dict + changeable factors of the problem + initial_solution : tuple + default initial solution from which solvers start + budget : int > 0 + max number of replications (fn evals) for a solver to take + specifications : dict + details of each factor (for GUI, data validation, and defaults) + + Arguments + --------- + name : str + user-specified name for problem + fixed_factors : dict + dictionary of user-specified problem factors + model_fixed factors : dict + subset of user-specified non-decision factors to pass through to the model + + See also + -------- + base.Problem + """ + def __init__(self, name="COVID-vac-test-3", fixed_factors={}, model_fixed_factors={}): + self.name = name + self.n_objectives = 1 + self.n_stochastic_constraints = 0 + self.minmax = (-1,) + self.constraint_type = "box" + self.variable_type = "continuous" + self.gradient_available = True + self.optimal_value = None + self.optimal_solution = None + self.model_default_factors = {} + self.model_fixed_factors = {} + self.model_decision_factors = {"freq_vac"} + self.factors = fixed_factors + self.specifications = { + "initial_solution": { + "description": "Initial solution from which solvers start.", + "datatype": tuple, + "default": (0/7, 0/7, 0/7, 0, 0, 0) + }, + "budget": { + "description": "Max # of replications for a solver to take.", + "datatype": int, + "default": 600 + }, + "vaccine_cap": { + "description": "Maximum number of vaccines.", + "datatype": int, + "default": 8000 + }, + "testing_cap": { + "description": "Maximum testing capacity per day.", + "datatype": int, + "default": 5000 + }, + "pen_coef": { + "description": "Penalty coefficient.", + "datatype": int, + "default": 100000 + }, + } + self.check_factor_list = { + "initial_solution": self.check_initial_solution, + "budget": self.check_budget, + "testing_cap": self.check_testing_cap, + "vaccine_cap": self.check_vaccine_cap, + "pen_coef": self.check_pen_coef + } + super().__init__(fixed_factors, model_fixed_factors) + # Instantiate model with fixed factors and overwritten defaults. + self.model = COVID_vac(self.model_fixed_factors) + self.dim = len(self.model.factors["group_size"]) + self.lower_bounds = (0,) * self.dim + self.upper_bounds = (1,) * self.dim + self.Ci = None + self.Ce = None + self.di = None + self.de = None + + def check_vaccine_cap(self): + return self.factors["vaccine_cap"] > 0 + + def check_testing_cap(self): + return self.factors["testing_cap"] > 0 + + def check_pen_coef(self): + return self.factors["pen_coef"] > 0 + + def vector_to_factor_dict(self, vector): + """ + Convert a vector of variables to a dictionary with factor keys + + Arguments + --------- + vector : tuple + vector of values associated with decision variables + + Returns + ------- + factor_dict : dictionary + dictionary with factor keys and associated values + """ + factor_dict = { + "freq_vac": vector + } + return factor_dict + + def factor_dict_to_vector(self, factor_dict): + """ + Convert a dictionary with factor keys to a vector + of variables. + + Arguments + --------- + factor_dict : dictionary + dictionary with factor keys and associated values + + Returns + ------- + vector : tuple + vector of values associated with decision variables + """ + vector = (factor_dict["freq_vac"], ) + return vector + + def response_dict_to_objectives(self, response_dict): + """ + Convert a dictionary with response keys to a vector + of objectives. + + Arguments + --------- + response_dict : dictionary + dictionary with response keys and associated values + + Returns + ------- + objectives : tuple + vector of objectives + """ + objectives = (response_dict["total_infected"], ) # here, we use the weighted number of infected people as objective; + return objectives + + def response_dict_to_stoch_constraints(self, response_dict): + """ + Convert a dictionary with response keys to a vector + of left-hand sides of stochastic constraints: E[Y] >= 0 + + Arguments + --------- + response_dict : dictionary + dictionary with response keys and associated values + + Returns + ------- + stoch_constraints : tuple + vector of LHSs of stochastic constraint + """ + stoch_constraints = None + return stoch_constraints + + def deterministic_objectives_and_gradients(self, x): + """ + Compute deterministic components of objectives for a solution `x`. + + Arguments + --------- + x : tuple + vector of decision variables + + Returns + ------- + det_objectives : tuple + vector of deterministic components of objectives + det_objectives_gradients : tuple + vector of gradients of deterministic components of objectives + """ + # Checking: + # print(self.factors['pen_coef']) + # for g in range(self.dim): + # print(np.dot(self.model.factors["group_size"][g],x[g])) + # print(np.sum(np.dot(self.model.factors["group_size"][g],x[g]) for g in range(self.dim))) + # print(max(np.sum(np.dot(self.model.factors["group_size"][g],x[g]) for g in range(self.dim)) - self.factors["vaccine_cap"],0)) + + x_v = x[self.dim: self.dim * 2] + x_t = x[: self.dim] + det_objectives_test = (self.factors["pen_coef"] * (max(np.sum(np.ceil(np.dot(self.model.factors["group_size"][g], x_t[g])) for g in range(self.dim)) - self.factors["testing_cap"], 0)),) + det_objectives_vac = (self.factors["pen_coef"] * (max(np.sum(np.ceil(np.dot(self.model.factors["group_size"][g], x_v[g])) for g in range(self.dim)) - self.factors["vaccine_cap"], 0)),) + det_objectives = max(det_objectives_test, det_objectives_vac) + det_objectives_gradients = ((0,),) + return det_objectives, det_objectives_gradients + + def deterministic_stochastic_constraints_and_gradients(self, x): + """ + Compute deterministic components of stochastic constraints + for a solution `x`. + + Arguments + --------- + x : tuple + vector of decision variables + + Returns + ------- + det_stoch_constraints : tuple + vector of deterministic components of stochastic constraints + det_stoch_constraints_gradients : tuple + vector of gradients of deterministic components of + stochastic constraints + """ + det_stoch_constraints = None + det_stoch_constraints_gradients = None + return det_stoch_constraints, det_stoch_constraints_gradients + + def check_deterministic_constraints(self, x): + """ + Check if a solution `x` satisfies the problem's deterministic + constraints. + + Arguments + --------- + x : tuple + vector of decision variables + + Returns + ------- + satisfies : bool + indicates if solution `x` satisfies the deterministic constraints. + """ + return np.all(x >= 0) + + def get_random_solution(self, rand_sol_rng): + """ + Generate a random solution for starting or restarting solvers. + + Arguments + --------- + rand_sol_rng : rng.MRG32k3a object + random-number generator used to sample a new random solution + + Returns + ------- + x : tuple + vector of decision variables + """ + # If uniformly find freq & vac: + # x = tuple([rand_sol_rng.uniform(0,1) for _ in range(self.dim)]) + # f = tuple([rand_sol_rng.uniform(0,1) for _ in range(self.dim)]) + # dx = f + x + + # By redistribution method: (reallocate the resource(vac/freq) if some group have more resource than their population size) + # Linear constraint + # For vaccine: + x = tuple(rand_sol_rng.integer_random_vector_from_simplex(self.model.factors["num_groups"], self.factors["vaccine_cap"])) + v = np.divide(x, self.model.factors["group_size"]) + + x = np.array(x) # For the temporal calculation + v = np.array(v) + + n = len(self.model.factors["group_size"]) + while True: + xn = x - np.array(self.model.factors["group_size"]) + idx = np.argmax(xn) + if v[idx] > 1: + remain = x[idx] - self.model.factors["group_size"][idx] + x[idx] -= remain + x += int(np.floor(remain / n)) + short_idx = np.argmax(np.array(self.model.factors["group_size"]) - x) + rr = remain - int(np.floor(remain / 3) * n) + if rr > 0: + x[short_idx] += rr + v = np.divide(x, self.model.factors["group_size"]) + else: + break + + x = tuple(x) + v = tuple(v) + + # If modify testing freq by redistribution, for testing: + xf = tuple(rand_sol_rng.integer_random_vector_from_simplex(self.model.factors["num_groups"], self.factors["testing_cap"])) + f = np.divide(xf, self.model.factors["group_size"]) + + xf = np.array(xf) # For the temporal calculation + f = np.array(f) + + n = len(self.model.factors["group_size"]) + while True: + xn = xf - np.array(self.model.factors["group_size"]) + idx = np.argmax(xn) + if f[idx] > 1: + remain = xf[idx] - self.model.factors["group_size"][idx] + xf[idx] -= remain + xf += int(np.floor(remain / n)) + short_idx = np.argmax(np.array(self.model.factors["group_size"]) - xf) + rr = remain - int(np.floor(remain / 3) * n) + if rr > 0: + xf[short_idx] += rr + f = np.divide(xf, self.model.factors["group_size"]) + else: + break + + xf = tuple(xf) + f = tuple(f) + + print('initial f, v: ', f, v) + + # Final decision variable freq_vac: + dx = f + v + + return dx + + +# Sub-problem 4: with given testing frequency for each group, find the optimal vaccination policy. +class CovidMinInfectVac4(Problem): + """ + Base class to implement simulation-optimization problems. + + Attributes + ---------- + name : string + name of problem + dim : int + number of decision variables + n_objectives : int + number of objectives + n_stochastic_constraints : int + number of stochastic constraints + minmax : tuple of int (+/- 1) + indicator of maximization (+1) or minimization (-1) for each objective + constraint_type : string + description of constraints types: + "unconstrained", "box", "deterministic", "stochastic" + variable_type : string + description of variable types: + "discrete", "continuous", "mixed" + lower_bounds : tuple + lower bound for each decision variable + upper_bounds : tuple + upper bound for each decision variable + Ci : ndarray (or None) + Coefficient matrix for linear inequality constraints of the form Ci@x <= di + Ce : ndarray (or None) + Coefficient matrix for linear equality constraints of the form Ce@x = de + di : ndarray (or None) + Constraint vector for linear inequality constraints of the form Ci@x <= di + de : ndarray (or None) + Constraint vector for linear equality constraints of the form Ce@x = de + gradient_available : bool + indicates if gradient of objective function is available + optimal_value : float + optimal objective function value + optimal_solution : tuple + optimal solution + model : Model object + associated simulation model that generates replications + model_default_factors : dict + default values for overriding model-level default factors + model_fixed_factors : dict + combination of overriden model-level factors and defaults + model_decision_factors : set of str + set of keys for factors that are decision variables + rng_list : list of rng.MRG32k3a objects + list of RNGs used to generate a random initial solution + or a random problem instance + factors : dict + changeable factors of the problem + initial_solution : tuple + default initial solution from which solvers start + budget : int > 0 + max number of replications (fn evals) for a solver to take + specifications : dict + details of each factor (for GUI, data validation, and defaults) + + Arguments + --------- + name : str + user-specified name for problem + fixed_factors : dict + dictionary of user-specified problem factors + model_fixed factors : dict + subset of user-specified non-decision factors to pass through to the model + + See also + -------- + base.Problem + """ + def __init__(self, name="COVID-vac-test-4", fixed_factors={}, model_fixed_factors={}): + self.name = name + self.n_objectives = 1 + self.n_stochastic_constraints = 0 + self.minmax = (-1,) + self.constraint_type = "box" + self.variable_type = "continuous" + self.gradient_available = True + self.optimal_value = None + self.optimal_solution = None + self.model_default_factors = {} + self.model_fixed_factors = {} + self.model_decision_factors = {"freq_vac"} + self.factors = fixed_factors + self.specifications = { + "initial_solution": { + "description": "Initial solution from which solvers start.", + "datatype": tuple, + "default": (1/7, 1/7, 1/7, 0, 0, 0) + }, + "budget": { + "description": "Max # of replications for a solver to take.", + "datatype": int, + "default": 600 + }, + "vaccine_cap": { + "description": "Maximum number of vaccines.", + "datatype": int, + "default": 8000 + }, + "testing_cap": { + "description": "Maximum testing capacity per day.", + "datatype": int, + "default": 5000 + }, + "testing_freq": { + "description": "Testing frequency for each group each day.", + "datatype": tuple, + "default": (1/7, 1/7, 1/7) + }, + "pen_coef": { + "description": "Penalty coefficient.", + "datatype": int, + "default": 100000 + }, + } + self.check_factor_list = { + "initial_solution": self.check_initial_solution, + "budget": self.check_budget, + "testing_cap": self.check_testing_cap, + "vaccine_cap": self.check_vaccine_cap, + "testing_freq": self.check_testing_freq, + "pen_coef": self.check_pen_coef + } + super().__init__(fixed_factors, model_fixed_factors) + # Instantiate model with fixed factors and overwritten defaults. + self.model = COVID_vac(self.model_fixed_factors) + self.dim = len(self.model.factors["group_size"]) + self.lower_bounds = (0,) * self.dim + self.upper_bounds = (1,) * self.dim + self.Ci = None + self.Ce = None + self.di = None + self.de = None + + def check_vaccine_cap(self): + return self.factors["vaccine_cap"] > 0 + + def check_testing_cap(self): + return self.factors["testing_cap"] > 0 + + def check_testing_freq(self): + return all(np.array(self.factors["testing_freq"]) >= 0) + + def check_pen_coef(self): + return self.factors["pen_coef"] > 0 + + def vector_to_factor_dict(self, vector): + """ + Convert a vector of variables to a dictionary with factor keys + + Arguments + --------- + vector : tuple + vector of values associated with decision variables + + Returns + ------- + factor_dict : dictionary + dictionary with factor keys and associated values + """ + factor_dict = { + "freq_vac": vector + } + return factor_dict + + def factor_dict_to_vector(self, factor_dict): + """ + Convert a dictionary with factor keys to a vector + of variables. + + Arguments + --------- + factor_dict : dictionary + dictionary with factor keys and associated values + + Returns + ------- + vector : tuple + vector of values associated with decision variables + """ + vector = (factor_dict["freq_vac"], ) + return vector + + def response_dict_to_objectives(self, response_dict): + """ + Convert a dictionary with response keys to a vector + of objectives. + + Arguments + --------- + response_dict : dictionary + dictionary with response keys and associated values + + Returns + ------- + objectives : tuple + vector of objectives + """ + objectives = (response_dict["total_infected"], ) + return objectives + + def response_dict_to_stoch_constraints(self, response_dict): + """ + Convert a dictionary with response keys to a vector + of left-hand sides of stochastic constraints: E[Y] >= 0 + + Arguments + --------- + response_dict : dictionary + dictionary with response keys and associated values + + Returns + ------- + stoch_constraints : tuple + vector of LHSs of stochastic constraint + """ + stoch_constraints = None + return stoch_constraints + + def deterministic_objectives_and_gradients(self, x): + """ + Compute deterministic components of objectives for a solution `x`. + + Arguments + --------- + x : tuple + vector of decision variables + + Returns + ------- + det_objectives : tuple + vector of deterministic components of objectives + det_objectives_gradients : tuple + vector of gradients of deterministic components of objectives + """ + # Checking: + # print(self.factors['pen_coef']) + # for g in range(self.dim): + # print(np.dot(self.model.factors["group_size"][g],x[g])) + # print(np.sum(np.dot(self.model.factors["group_size"][g],x[g]) for g in range(self.dim))) + # print(max(np.sum(np.dot(self.model.factors["group_size"][g],x[g]) for g in range(self.dim)) - self.factors["vaccine_cap"],0)) + + x_v = x[self.dim: self.dim * 2] + x_t = x[: self.dim] + det_objectives_test = (self.factors["pen_coef"] * (max(np.sum(np.ceil(np.dot(self.model.factors["group_size"][g], x_t[g])) for g in range(self.dim)) - self.factors["testing_cap"], 0)),) + det_objectives_vac = (self.factors["pen_coef"] * (max(np.sum(np.ceil(np.dot(self.model.factors["group_size"][g], x_v[g])) for g in range(self.dim)) - self.factors["vaccine_cap"], 0)),) + det_objectives = max(det_objectives_test, det_objectives_vac) + det_objectives_gradients = ((0,),) + return det_objectives, det_objectives_gradients + + def deterministic_stochastic_constraints_and_gradients(self, x): + """ + Compute deterministic components of stochastic constraints + for a solution `x`. + + Arguments + --------- + x : tuple + vector of decision variables + + Returns + ------- + det_stoch_constraints : tuple + vector of deterministic components of stochastic constraints + det_stoch_constraints_gradients : tuple + vector of gradients of deterministic components of + stochastic constraints + """ + det_stoch_constraints = None + det_stoch_constraints_gradients = None + return det_stoch_constraints, det_stoch_constraints_gradients + + def check_deterministic_constraints(self, x): + """ + Check if a solution `x` satisfies the problem's deterministic + constraints. + + Arguments + --------- + x : tuple + vector of decision variables + + Returns + ------- + satisfies : bool + indicates if solution `x` satisfies the deterministic constraints. + """ + return np.all(x >= 0) + + def get_random_solution(self, rand_sol_rng): + """ + Generate a random solution for starting or restarting solvers. + + Arguments + --------- + rand_sol_rng : rng.MRG32k3a object + random-number generator used to sample a new random solution + + Returns + ------- + x : tuple + vector of decision variables + """ + # If uniformly find freq & vac: + # x = tuple([rand_sol_rng.uniform(0,1) for _ in range(self.dim)]) + # f = tuple([rand_sol_rng.uniform(0,1) for _ in range(self.dim)]) + # dx = f + x + + # By redistribution method: (reallocate the resource(vac/freq) if some group have more resource than their population size) + # Linear constraint + # For vaccine: + x = tuple(rand_sol_rng.integer_random_vector_from_simplex(self.model.factors["num_groups"], self.factors["vaccine_cap"])) + v = np.divide(x, self.model.factors["group_size"]) + + x = np.array(x) # For the temporal calculation + v = np.array(v) + + n = len(self.model.factors["group_size"]) + while True: + xn = x - np.array(self.model.factors["group_size"]) + idx = np.argmax(xn) + if v[idx] > 1: + remain = x[idx] - self.model.factors["group_size"][idx] + x[idx] -= remain + x += int(np.floor(remain / n)) + short_idx = np.argmax(np.array(self.model.factors["group_size"]) - x) + rr = remain - int(np.floor(remain / 3) * n) + if rr > 0: + x[short_idx] += rr + v = np.divide(x, self.model.factors["group_size"]) + else: + break + + x = tuple(x) + v = tuple(v) + + # with given testing policy for the period: + f = self.factors["testing_freq"] + + xf = tuple(xf) + f = tuple(f) + + print('initial f, v: ', f, v) + + # Final decision variable freq_vac: + dx = f + v + + return dx \ No newline at end of file diff --git a/simopt/models/san.py b/simopt/models/san.py index 8d219cd3f..d83e08ede 100644 --- a/simopt/models/san.py +++ b/simopt/models/san.py @@ -129,7 +129,7 @@ def replicate(self, rng_list): graph_in[a[1]].add(a[0]) graph_out[a[0]].add(a[1]) indegrees = [len(graph_in[n]) for n in range(1, self.factors["num_nodes"] + 1)] - # outdegrees = [len(graph_out[n]) for n in range(1, self.factors["num_nodes"]+1)] + queue = [] topo_order = [] for n in range(self.factors["num_nodes"]): @@ -212,6 +212,14 @@ class SANLongestPath(Problem): lower bound for each decision variable upper_bounds : tuple upper bound for each decision variable + Ci : ndarray (or None) + Coefficient matrix for linear inequality constraints of the form Ci@x <= di + Ce : ndarray (or None) + Coefficient matrix for linear equality constraints of the form Ce@x = de + di : ndarray (or None) + Constraint vector for linear inequality constraints of the form Ci@x <= di + de : ndarray (or None) + Constraint vector for linear equality constraints of the form Ce@x = de gradient_available : bool indicates if gradient of objective function is available optimal_value : tuple @@ -252,6 +260,289 @@ class SANLongestPath(Problem): base.Problem """ def __init__(self, name="SAN-1", fixed_factors=None, model_fixed_factors=None): + if fixed_factors is None: + fixed_factors = {} + if model_fixed_factors is None: + model_fixed_factors = {} + self.name = name + self.n_objectives = 1 + self.n_stochastic_constraints = 0 + self.minmax = (-1,) + self.constraint_type = "box" + self.variable_type = "continuous" + self.gradient_available = True + self.optimal_value = None + self.optimal_solution = None + self.model_default_factors = {} + self.model_decision_factors = {"arc_means"} + self.factors = fixed_factors + self.specifications = { + "initial_solution": { + "description": "initial solution", + "datatype": tuple, + "default": (8,) * 13 + }, + "budget": { + "description": "max # of replications for a solver to take", + "datatype": int, + "default": 2000 + }, + "arc_costs": { + "description": "Cost associated to each arc.", + "datatype": tuple, + "default": (1, ) * 13 + } + } + self.check_factor_list = { + "initial_solution": self.check_initial_solution, + "budget": self.check_budget, + "arc_costs": self.check_arc_costs + } + super().__init__(fixed_factors, model_fixed_factors) + # Instantiate model with fixed factors and over-riden defaults. + self.model = SAN(self.model_fixed_factors) + self.dim = len(self.model.factors["arcs"]) + self.lower_bounds = (1e-2, ) * self.dim + self.upper_bounds = (np.inf, ) * self.dim + self.Ci = None + self.Ce = None + self.di = None + self.de = None + + def check_arc_costs(self): + positive = True + for x in list(self.factors["arc_costs"]): + positive = positive & x > 0 + return (len(self.factors["arc_costs"]) != self.model.factors["num_arcs"]) & positive + + def vector_to_factor_dict(self, vector): + """ + Convert a vector of variables to a dictionary with factor keys + + Arguments + --------- + vector : tuple + vector of values associated with decision variables + + Returns + ------- + factor_dict : dictionary + dictionary with factor keys and associated values + """ + factor_dict = { + "arc_means": vector[:] + } + return factor_dict + + def factor_dict_to_vector(self, factor_dict): + """ + Convert a dictionary with factor keys to a vector + of variables. + + Arguments + --------- + factor_dict : dictionary + dictionary with factor keys and associated values + + Returns + ------- + vector : tuple + vector of values associated with decision variables + """ + vector = tuple(factor_dict["arc_means"]) + return vector + + def response_dict_to_objectives(self, response_dict): + """ + Convert a dictionary with response keys to a vector + of objectives. + + Arguments + --------- + response_dict : dictionary + dictionary with response keys and associated values + + Returns + ------- + objectives : tuple + vector of objectives + """ + objectives = (response_dict["longest_path_length"],) + return objectives + + def response_dict_to_stoch_constraints(self, response_dict): + """ + Convert a dictionary with response keys to a vector + of left-hand sides of stochastic constraints: E[Y] <= 0 + + Arguments + --------- + response_dict : dictionary + dictionary with response keys and associated values + + Returns + ------- + stoch_constraints : tuple + vector of LHSs of stochastic constraint + """ + stoch_constraints = None + return stoch_constraints + + def deterministic_stochastic_constraints_and_gradients(self, x): + """ + Compute deterministic components of stochastic constraints for a solution `x`. + + Arguments + --------- + x : tuple + vector of decision variables + + Returns + ------- + det_stoch_constraints : tuple + vector of deterministic components of stochastic constraints + det_stoch_constraints_gradients : tuple + vector of gradients of deterministic components of stochastic constraints + """ + det_stoch_constraints = None + det_stoch_constraints_gradients = ((0,) * self.dim,) # tuple of tuples – of sizes self.dim by self.dim, full of zeros + return det_stoch_constraints, det_stoch_constraints_gradients + + def deterministic_objectives_and_gradients(self, x): + """ + Compute deterministic components of objectives for a solution `x`. + + Arguments + --------- + x : tuple + vector of decision variables + + Returns + ------- + det_objectives : tuple + vector of deterministic components of objectives + det_objectives_gradients : tuple + vector of gradients of deterministic components of objectives + """ + det_objectives = (np.sum(np.array(self.factors["arc_costs"]) / np.array(x)),) + det_objectives_gradients = (-np.array(self.factors["arc_costs"]) / (np.array(x) ** 2),) + return det_objectives, det_objectives_gradients + + def check_deterministic_constraints(self, x): + """ + Check if a solution `x` satisfies the problem's deterministic constraints. + + Arguments + --------- + x : tuple + vector of decision variables + + Returns + ------- + satisfies : bool + indicates if solution `x` satisfies the deterministic constraints. + """ + return np.all(np.array(x) >= 0) + + def get_random_solution(self, rand_sol_rng): + """ + Generate a random solution for starting or restarting solvers. + + Arguments + --------- + rand_sol_rng : mrg32k3a.mrg32k3a.MRG32k3a object + random-number generator used to sample a new random solution + + Returns + ------- + x : tuple + vector of decision variables + """ + x = tuple([rand_sol_rng.lognormalvariate(lq=0.1, uq=10) for _ in range(self.dim)]) + return x + + +""" +Summary +------- +Minimize the duration of the longest path from a to i subject to a lower bound in sum of arc_means. +""" + + +class SANLongestPathConstr(Problem): + """ + Base class to implement simulation-optimization problems. + + Attributes + ---------- + name : string + name of problem + dim : int + number of decision variables + n_objectives : int + number of objectives + n_stochastic_constraints : int + number of stochastic constraints + minmax : tuple of int (+/- 1) + indicator of maximization (+1) or minimization (-1) for each objective + constraint_type : string + description of constraints types: + "unconstrained", "box", "deterministic", "stochastic" + variable_type : string + description of variable types: + "discrete", "continuous", "mixed" + lower_bounds : tuple + lower bound for each decision variable + upper_bounds : tuple + upper bound for each decision variable + Ci : ndarray (or None) + Coefficient matrix for linear inequality constraints of the form Ci@x <= di + Ce : ndarray (or None) + Coefficient matrix for linear equality constraints of the form Ce@x = de + di : ndarray (or None) + Constraint vector for linear inequality constraints of the form Ci@x <= di + de : ndarray (or None) + Constraint vector for linear equality constraints of the form Ce@x = de + gradient_available : bool + indicates if gradient of objective function is available + optimal_value : tuple + optimal objective function value + optimal_solution : tuple + optimal solution + model : Model object + associated simulation model that generates replications + model_default_factors : dict + default values for overriding model-level default factors + model_fixed_factors : dict + combination of overriden model-level factors and defaults + model_decision_factors : set of str + set of keys for factors that are decision variables + rng_list : list of mrg32k3a.mrg32k3a.MRG32k3a objects + list of RNGs used to generate a random initial solution + or a random problem instance + factors : dict + changeable factors of the problem + initial_solution : list + default initial solution from which solvers start + budget : int > 0 + max number of replications (fn evals) for a solver to take + specifications : dict + details of each factor (for GUI, data validation, and defaults) + + Arguments + --------- + name : str + user-specified name for problem + fixed_factors : dict + dictionary of user-specified problem factors + model_fixed factors : dict + subset of user-specified non-decision factors to pass through to the model + + See also + -------- + base.Problem + """ + def __init__(self, name="SAN-2", fixed_factors=None, model_fixed_factors=None): if fixed_factors is None: fixed_factors = {} if model_fixed_factors is None: @@ -283,6 +574,11 @@ def __init__(self, name="SAN-1", fixed_factors=None, model_fixed_factors=None): "description": "Cost associated to each arc.", "datatype": tuple, "default": (1,) * 13 + }, + "sum_lb": { + "description": "Lower bound for the sum of arc means", + "datatype": float, + "default": 30.0 } } self.check_factor_list = { @@ -296,6 +592,10 @@ def __init__(self, name="SAN-1", fixed_factors=None, model_fixed_factors=None): self.dim = len(self.model.factors["arcs"]) self.lower_bounds = (1e-2,) * self.dim self.upper_bounds = (np.inf,) * self.dim + self.Ci = -1 * np.ones(13) + self.Ce = None + self.di = -1 * np.array([self.factors["sum_lb"]]) + self.de = None def check_arc_costs(self): positive = True @@ -412,8 +712,8 @@ def deterministic_objectives_and_gradients(self, x): det_objectives_gradients : tuple vector of gradients of deterministic components of objectives """ - det_objectives = (np.sum(np.array(self.factors["arc_costs"]) / np.array(x)),) - det_objectives_gradients = (-np.array(self.factors["arc_costs"]) / (np.array(x) ** 2),) + det_objectives = (0,) + det_objectives_gradients = ((0,) * self.dim,) return det_objectives, det_objectives_gradients def check_deterministic_constraints(self, x): @@ -446,5 +746,9 @@ def get_random_solution(self, rand_sol_rng): x : tuple vector of decision variables """ - x = tuple([rand_sol_rng.lognormalvariate(lq=0.1, uq=10) for _ in range(self.dim)]) + while True: + x = [rand_sol_rng.lognormalvariate(lq = 0.1, uq = 10) for _ in range(self.dim)] + if np.sum(x) >= self.factors['sum_lb']: + break + x= tuple(x) return x diff --git a/simopt/models/smf.py b/simopt/models/smf.py new file mode 100644 index 000000000..aa21f853e --- /dev/null +++ b/simopt/models/smf.py @@ -0,0 +1,480 @@ +""" +Summary +------- +Simulate duration of a stochastic Max-Flow network (SMF). +A detailed description of the model/problem can be found +`here `_. +""" + +import numpy as np +from ortools.graph.python import max_flow +from ..base import Model, Problem + + +class SMF(Model): + """ + A model that simulates a stochastic Max-Flow problem with + capacities deducted with multivariate distributed noise distributed durations + + Attributes + ---------- + name : string + name of model + n_rngs : int + number of random-number generators used to run a simulation replication + n_responses : int + number of responses (performance measures) + factors : dict + changeable factors of the simulation model + specifications : dict + details of each factor (for GUI and data validation) + check_factor_list : dict + switch case for checking factor simulatability + + Arguments + --------- + fixed_factors : nested dict + fixed factors of the simulation model + + See also + -------- + base.Model + """ + def __init__(self, fixed_factors=None): + if fixed_factors is None: + fixed_factors = {} + self.name = "SMF" + self.n_rngs = 1 + self.n_responses = 1 + cov_fac = np.zeros((20, 20)) + np.fill_diagonal(cov_fac, 4) + cov_fac = cov_fac.tolist() + self.specifications = { + "num_nodes": { + "description": "number of nodes, 0 being the source, highest being the sink", + "datatype": int, + "default": 10 + }, + "source_index": { + "description": "source node index", + "datatype": int, + "default": 0 + }, + "sink_index": { + "description": "sink node index", + "datatype": int, + "default": 9 + }, + "arcs": { + "description": "list of arcs", + "datatype": list, + "default": [(0, 1), (0, 2), (0, 3), (1, 2), (1, 4), (2, 4), (4, 2), (3, 2), (2, 5), (4, 5), (3, 6), (3, 7), (6, 2), (6, 5), (6, 7), (5, 8), (6, 8), (6, 9), (7, 9), (8, 9)] + }, + "assigned_capacities": { + "description": "Assigned capacity of each arc", + "datatype": list, + "default": [5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5] + }, + "mean_noise": { + "description": "The mean noise in reduction of arc capacities", + "datatype": list, + "default": [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] + }, + "cov_noise": { + "description": "Covariance matrix of noise", + "datatype": list, + "default": cov_fac + } + + } + self.check_factor_list = { + "num_nodes": self.check_num_nodes, + "arcs": self.check_arcs, + "assigned_capacities": self.check_capacities, + "mean_noise": self.check_mean, + "cov_noise": self.check_cov, + "source_index": self.check_s, + "sink_index": self.check_t + } + # Set factors of the simulation model. + super().__init__(fixed_factors) + + def check_num_nodes(self): + return self.factors["num_nodes"] > 0 + + def dfs(self, graph, start, visited=None): + if visited is None: + visited = set() + visited.add(start) + for next in graph[start] - visited: + self.dfs(graph, next, visited) + return visited + + def check_arcs(self): + if len(self.factors["arcs"]) <= 0: + return False + # Check source is connected to the sink. + graph = {node: set() for node in range(0, self.factors["num_nodes"])} + for a in self.factors["arcs"]: + graph[a[0]].add(a[1]) + visited = self.dfs(graph, self.factors["source_index"]) + if self.factors["source_index"] in visited and self.factors["sink_index"] in visited: + return True + return False + + def check_capacities(self): + positive = True + for x in list(self.factors["assigned_capacities"]): + positive = positive & (x > 0) + return (len(self.factors["assigned_capacities"]) == len(self.factors["arcs"])) & positive + + def check_mean(self): + return len(self.factors["mean_noise"]) == len(self.factors["arcs"]) + + def check_cov(self): + return np.array(self.factors["cov_noise"]).shape == (len(self.factors["arcs"]), len(self.factors["arcs"])) + + def check_s(self): + return self.factors["source_index"] >= 0 and self.factors["source_index"] <= self.factors["num_nodes"] + + def check_t(self): + return self.factors["sink_index"] >= 0 and self.factors["sink_index"] <= self.factors["num_nodes"] + + def replicate(self, rng_list): + """ + Simulate a single replication for the current model factors. + + Arguments + --------- + rng_list : list of mrg32k3a.mrg32k3a.MRG32k3a + rngs for model to use when simulating a replication + + Returns + ------- + responses : dict + performance measures of interest + "longest_path_length" = length/duration of longest path + gradients : dict of dicts + gradient estimates for each response + """ + # Designate separate random number generators. + solver = max_flow.SimpleMaxFlow() + exp_rng = rng_list[0] + # From input graph generate start end end nodes. + start_nodes = [] + end_nodes = [] + for i, j in self.factors["arcs"]: + start_nodes.append(i) + end_nodes.append(j) + # Generate actual capacity. + for i in range(len(self.factors["arcs"])): + noise = exp_rng.mvnormalvariate(self.factors["mean_noise"], np.array(self.factors["cov_noise"])) + capacities = [] + for i in range(len(noise)): + capacities.append(max(1000 * (self.factors["assigned_capacities"][i] - noise[i]), 0)) + # Add arcs in bulk. + solver.add_arcs_with_capacity(start_nodes, end_nodes, capacities) + status = solver.solve(self.factors["source_index"], self.factors["sink_index"]) + if status != solver.OPTIMAL: + print('There was an issue with the max flow input.') + print(f'Status: {status}') + exit(1) + + # Construct gradient vector (=1 if has a outflow from min-cut nodes). + gradient = np.zeros(len(self.factors["arcs"])) + grad_arclist = [] + min_cut_nodes = solver.get_source_side_min_cut() + for i in min_cut_nodes: + for j in range(self.factors['num_nodes']): + if j not in min_cut_nodes: + grad_arc = (i, j) + if (i, j) in self.factors['arcs']: + grad_arclist.append(grad_arc) + for arc in grad_arclist: + gradient[self.factors['arcs'].index(arc)] = 1 + + responses = {"Max Flow": solver.optimal_flow() / 1000} + gradients = {response_key: {factor_key: np.nan for factor_key in self.specifications} for response_key in responses} + gradients["Max Flow"]["assigned_capacities"] = gradient + return responses, gradients + + +""" +Summary +------- +Maximize the expected max flow from the source node s to the sink node t. +""" + + +class SMF_Max(Problem): + """ + Base class to implement simulation-optimization problems. + + Attributes + ---------- + name : string + name of problem + dim : int + number of decision variables + n_objectives : int + number of objectives + n_stochastic_constraints : int + number of stochastic constraints + minmax : tuple of int (+/- 1) + indicator of maximization (+1) or minimization (-1) for each objective + constraint_type : string + description of constraints types: + "unconstrained", "box", "deterministic", "stochastic" + variable_type : string + description of variable types: + "discrete", "continuous", "mixed" + lower_bounds : tuple + lower bound for each decision variable + upper_bounds : tuple + upper bound for each decision variable + Ci : ndarray (or None) + Coefficient matrix for linear inequality constraints of the form Ci@x <= di + Ce : ndarray (or None) + Coefficient matrix for linear equality constraints of the form Ce@x = de + di : ndarray (or None) + Constraint vector for linear inequality constraints of the form Ci@x <= di + de : ndarray (or None) + Constraint vector for linear equality constraints of the form Ce@x = de + gradient_available : bool + indicates if gradient of objective function is available + optimal_value : tuple + optimal objective function value + optimal_solution : tuple + optimal solution + model : Model object + associated simulation model that generates replications + model_default_factors : dict + default values for overriding model-level default factors + model_fixed_factors : dict + combination of overriden model-level factors and defaults + model_decision_factors : set of str + set of keys for factors that are decision variables + rng_list : list of mrg32k3a.mrg32k3a.MRG32k3a objects + list of RNGs used to generate a random initial solution + or a random problem instance + factors : dict + changeable factors of the problem + initial_solution : list + default initial solution from which solvers start + budget : int > 0 + max number of replications (fn evals) for a solver to take + specifications : dict + details of each factor (for GUI, data validation, and defaults) + + Arguments + --------- + name : str + user-specified name for problem + fixed_factors : dict + dictionary of user-specified problem factors + model_fixed factors : dict + subset of user-specified non-decision factors to pass through to the model + + See also + -------- + base.Problem + """ + def __init__(self, name="SMF-1", fixed_factors=None, model_fixed_factors=None): + if fixed_factors is None: + fixed_factors = {} + if model_fixed_factors is None: + model_fixed_factors = {} + self.name = name + self.n_objectives = 1 + self.n_stochastic_constraints = 0 + self.minmax = (1, ) + self.constraint_type = "deterministic" + self.variable_type = "continuous" + self.gradient_available = True + self.optimal_value = None + self.optimal_solution = None + self.model_default_factors = {} + self.model_decision_factors = {"assigned_capacities"} + self.factors = fixed_factors + self.specifications = { + "initial_solution": { + "description": "initial solution", + "datatype": tuple, + "default": (1, ) * 20 + }, + "budget": { + "description": "max # of replications for a solver to take", + "datatype": int, + "default": 10000 + }, + "cap": { + "description": "total set-capacity to be allocated to arcs.", + "datatype": int, + "default": 100 + } + } + self.check_factor_list = { + "initial_solution": self.check_initial_solution, + "budget": self.check_budget, + "cap": self.check_cap + } + super().__init__(fixed_factors, model_fixed_factors) + # Instantiate model with fixed factors and over-riden defaults. + self.model = SMF(self.model_fixed_factors) + self.dim = len(self.model.factors["arcs"]) + self.lower_bounds = (0, ) * self.dim + self.upper_bounds = (np.inf, ) * self.dim + self.Ci = np.ones(20) + self.Ce = None + self.di = np.array([self.factors["cap"]]) + self.de = None + + def check_cap(self): + return self.factors["cap"] >= 0 + + def vector_to_factor_dict(self, vector): + """ + Convert a vector of variables to a dictionary with factor keys + + Arguments + --------- + vector : tuple + vector of values associated with decision variables + + Returns + ------- + factor_dict : dictionary + dictionary with factor keys and associated values + """ + factor_dict = { + "assigned_capacities": vector[:] + } + return factor_dict + + def factor_dict_to_vector(self, factor_dict): + """ + Convert a dictionary with factor keys to a vector + of variables. + + Arguments + --------- + factor_dict : dictionary + dictionary with factor keys and associated values + + Returns + ------- + vector : tuple + vector of values associated with decision variables + """ + vector = tuple(factor_dict["assigned_capacities"]) + return vector + + def response_dict_to_objectives(self, response_dict): + """ + Convert a dictionary with response keys to a vector + of objectives. + + Arguments + --------- + response_dict : dictionary + dictionary with response keys and associated values + + Returns + ------- + objectives : tuple + vector of objectives + """ + objectives = (response_dict["Max Flow"], ) + return objectives + + def response_dict_to_stoch_constraints(self, response_dict): + """ + Convert a dictionary with response keys to a vector + of left-hand sides of stochastic constraints: E[Y] <= 0 + + Arguments + --------- + response_dict : dictionary + dictionary with response keys and associated values + + Returns + ------- + stoch_constraints : tuple + vector of LHSs of stochastic constraint + """ + stoch_constraints = None + return stoch_constraints + + def deterministic_stochastic_constraints_and_gradients(self, x): + """ + Compute deterministic components of stochastic constraints for a solution `x`. + + Arguments + --------- + x : tuple + vector of decision variables + + Returns + ------- + det_stoch_constraints : tuple + vector of deterministic components of stochastic constraints + det_stoch_constraints_gradients : tuple + vector of gradients of deterministic components of stochastic constraints + """ + det_stoch_constraints = None + det_stoch_constraints_gradients = None + return det_stoch_constraints, det_stoch_constraints_gradients + + def deterministic_objectives_and_gradients(self, x): + """ + Compute deterministic components of objectives for a solution `x`. + + Arguments + --------- + x : tuple + vector of decision variables + + Returns + ------- + det_objectives : tuple + vector of deterministic components of objectives + det_objectives_gradients : tuple + vector of gradients of deterministic components of objectives + """ + det_objectives = (0, ) + det_objectives_gradients = ((0, ) * self.dim,) + return det_objectives, det_objectives_gradients + + def check_deterministic_constraints(self, x): + """ + Check if a solution `x` satisfies the problem's deterministic constraints. + + Arguments + --------- + x : tuple + vector of decision variables + + Returns + ------- + satisfies : bool + indicates if solution `x` satisfies the deterministic constraints. + """ + + return sum(self.factors["assigned_capacities"]) <= self.factors["cap"] + + def get_random_solution(self, rand_sol_rng): + """ + Generate a random solution for starting or restarting solvers. + + Arguments + --------- + rand_sol_rng : mrg32k3a.mrg32k3a.MRG32k3a object + random-number generator used to sample a new random solution + + Returns + ------- + x : tuple + vector of decision variables + """ + x = rand_sol_rng.continuous_random_vector_from_simplex(len(self.model.factors["arcs"]), self.factors["cap"], False) + return x diff --git a/simopt/solvers/active_set.py b/simopt/solvers/active_set.py new file mode 100644 index 000000000..48488ebad --- /dev/null +++ b/simopt/solvers/active_set.py @@ -0,0 +1,596 @@ +""" +Summary +------- +ACTIVESET: An active set algorithm for problems with linear constraints i.e., Ce@x = de, Ci@x <= di. +A detailed description of the solver can be found `here `_. +""" +import numpy as np +import cvxpy as cp +import warnings +warnings.filterwarnings("ignore") + +from ..base import Solver + + +class ACTIVESET(Solver): + """ + The Active Set solver. + + Attributes + ---------- + name : string + name of solver + objective_type : string + description of objective types: + "single" or "multi" + constraint_type : string + description of constraints types: + "unconstrained", "box", "deterministic", "stochastic" + variable_type : string + description of variable types: + "discrete", "continuous", "mixed" + gradient_needed : bool + indicates if gradient of objective function is needed + factors : dict + changeable factors (i.e., parameters) of the solver + specifications : dict + details of each factor (for GUI, data validation, and defaults) + rng_list : list of rng.MRG32k3a objects + list of RNGs used for the solver's internal purposes + + Arguments + --------- + name : str + user-specified name for solver + fixed_factors : dict + fixed_factors of the solver + + See also + -------- + base.Solver + """ + def __init__(self, name="ACTIVESET", fixed_factors={}): + self.name = name + self.objective_type = "single" + self.constraint_type = "deterministic" + self.variable_type = "continuous" + self.gradient_needed = False + self.specifications = { + "crn_across_solns": { + "description": "use CRN across solutions?", + "datatype": bool, + "default": True + }, + "r": { + "description": "number of replications taken at each solution", + "datatype": int, + "default": 30 + }, + "alpha": { + "description": "tolerance for sufficient decrease condition.", + "datatype": float, + "default": 0.2 + }, + "beta": { + "description": "step size reduction factor in line search.", + "datatype": float, + "default": 0.9 + }, + "alpha_max": { + "description": "maximum step size.", + "datatype": float, + "default": 10.0 + }, + "lambda": { + "description": "magnifying factor for r inside the finite difference function", + "datatype": int, + "default": 2 + }, + "tol": { + "description": "floating point tolerance for checking tightness of constraints", + "datatype": float, + "default": 1e-7 + }, + "tol2": { + "description": "floating point tolerance for checking closeness of dot product to zero", + "datatype": float, + "default": 1e-7 + }, + "finite_diff_step": { + "description": "step size for finite difference", + "datatype": float, + "default": 1e-5 + } + } + self.check_factor_list = { + "crn_across_solns": self.check_crn_across_solns, + "r": self.check_r, + "alpha": self.check_alpha, + "beta": self.check_beta, + "alpha_max": self.check_alpha_max, + "lambda": self.check_lambda, + "tol": self.check_tol, + "tol2": self.check_tol2, + "finite_diff_step": self.check_finite_diff_step + } + super().__init__(fixed_factors) + + def check_r(self): + return self.factors["r"] > 0 + + def check_alpha(self): + return self.factors["alpha"] > 0 + + def check_beta(self): + return self.factors["beta"] > 0 & self.factors["beta"] < 1 + + def check_alpha_max(self): + return self.factors["alpha_max"] > 0 + + def check_lambda(self): + return self.factors["lambda"] > 0 + + def check_tol(self): + return self.factors["tol"] > 0 + + def check_tol2(self): + return self.factors["tol2"] > 0 + + def check_finite_diff_step(self): + return self.factors["finite_diff_step"] > 0 + + def solve(self, problem): + """ + Run a single macroreplication of a solver on a problem. + + Arguments + --------- + problem : Problem object + simulation-optimization problem to solve + crn_across_solns : bool + indicates if CRN are used when simulating different solutions + + Returns + ------- + recommended_solns : list of Solution objects + list of solutions recommended throughout the budget + intermediate_budgets : list of ints + list of intermediate budgets when recommended solutions changes + """ + recommended_solns = [] + intermediate_budgets = [] + expended_budget = 0 + + # Default values. + r = self.factors["r"] + alpha = self.factors["alpha"] + beta = self.factors["beta"] + tol = self.factors["tol"] + tol2 = self.factors["tol2"] + max_step = self.factors["alpha_max"] # Maximum step size + + # Upper bound and lower bound. + lower_bound = np.array(problem.lower_bounds) + upper_bound = np.array(problem.upper_bounds) + + # Input inequality and equlaity constraint matrix and vector. + # Cix <= di + # Cex = de + Ci = problem.Ci + di = problem.di + Ce = problem.Ce + de = problem.de + + # Remove redundant upper/lower bounds. + ub_inf_idx = np.where(~np.isinf(upper_bound))[0] + lb_inf_idx = np.where(~np.isinf(lower_bound))[0] + + # Form a constraint coefficient matrix where all the equality constraints are put on top and + # all the bound constraints in the bottom and a constraint coefficient vector. + if (Ce is not None) and (de is not None) and (Ci is not None) and (di is not None): + C = np.vstack((Ce, Ci)) + d = np.vstack((de.T, di.T)) + elif (Ce is not None) and (de is not None): + C = Ce + d = de.T + elif (Ci is not None) and (di is not None): + C = Ci + d = di.T + else: + C = np.empty([1, problem.dim]) + d = np.empty([1, 1]) + + if len(ub_inf_idx) > 0: + C = np.vstack((C, np.identity(upper_bound.shape[0]))) + d = np.vstack((d, upper_bound[np.newaxis].T)) + if len(lb_inf_idx) > 0: + C = np.vstack((C, -np.identity(lower_bound.shape[0]))) + d = np.vstack((d, -lower_bound[np.newaxis].T)) + + # Number of equality constraints. + if (Ce is not None) and (de is not None): + neq = len(de) + else: + neq = 0 + + # Checker for whether the problem is unconstrained. + unconstr_flag = (Ce is None) & (Ci is None) & (di is None) & (de is None) & (all(np.isinf(lower_bound))) & (all(np.isinf(upper_bound))) + + # Start with the initial solution. + new_solution = self.create_new_solution(problem.factors["initial_solution"], problem) + new_x = new_solution.x + + # If the initial solution is not feasible, generate one using phase one simplex. + if (not unconstr_flag) & (not self._feasible(new_x, problem, tol)): + new_x = self.find_feasible_initial(problem, Ce, Ci, de, di, tol) + new_solution = self.create_new_solution(tuple(new_x), problem) + + # Use r simulated observations to estimate the objective value. + problem.simulate(new_solution, r) + expended_budget += r + best_solution = new_solution + recommended_solns.append(new_solution) + intermediate_budgets.append(expended_budget) + + # Active constraint index vector. + acidx = [] + if not unconstr_flag: + # Initialize the active set to be the set of indices of the tight constraints. + cx = np.dot(C, new_x) + for j in range(cx.shape[0]): + if j < neq or np.isclose(cx[j], d[j], rtol=0, atol= tol): + acidx.append(j) + + while expended_budget < problem.factors["budget"]: + new_x = new_solution.x + # Check variable bounds. + forward = np.isclose(new_x, lower_bound, atol = tol).astype(int) + backward = np.isclose(new_x, upper_bound, atol = tol).astype(int) + # BdsCheck: 1 stands for forward, -1 stands for backward, 0 means central diff. + BdsCheck = np.subtract(forward, backward) + + if problem.gradient_available: + # Use IPA gradient if available. + grad = -1 * problem.minmax[0] * new_solution.objectives_gradients_mean[0] + else: + # Use finite difference to estimate gradient if IPA gradient is not available. + grad = self.finite_diff(new_solution, BdsCheck, problem, r, self.factors["finite_diff_step"]) + expended_budget += (2 * problem.dim - np.sum(BdsCheck != 0)) * r + # A while loop to prevent zero gradient + while np.all((grad == 0)): + if expended_budget > problem.factors["budget"]: + break + grad = self.finite_diff(new_solution, BdsCheck, problem, alpha, r) + expended_budget += (2 * problem.dim - np.sum(BdsCheck != 0)) * r + # Update r after each iteration. + r = int(self.factors["lambda"] * r) + + # If the active set is empty, search on negative gradient. + if len(acidx) == 0: + dir = -grad + else: + # Find the search direction and Lagrange multipliers of the direction-finding problem. + dir, lmbd, = self.compute_search_direction(acidx, grad, problem, C) + # If the optimal search direction is 0 + if (np.isclose(np.dot(grad, dir), 0, rtol=0, atol=tol2)): + # Terminate if Lagrange multipliers of the inequality constraints in the active set are all nonnegative. + if unconstr_flag or np.all(lmbd[neq:] >= 0): + break + # Otherwise, drop the inequality constraint in the active set with the most negative Lagrange multiplier. + else: + q = acidx[neq + np.argmin(lmbd[neq:][lmbd[neq:] < 0])] + acidx.remove(q) + else: + if not unconstr_flag: + idx = list(set(np.arange(C.shape[0])) - set(acidx)) # Constraints that are not in the active set. + # If all constraints are feasible. + if unconstr_flag or np.all(C[idx,:] @ dir <= 0): + # Line search to determine a step_size. + new_solution, step_size, expended_budget = self.line_search(problem, expended_budget, r, grad, new_solution, max_step, dir, alpha, beta) + # Update maximum step size for the next iteration. + max_step = step_size + + # Ratio test to determine the maximum step size possible + else: + # Get all indices not in the active set such that Ai^Td>0 + r_idx = list(set(idx).intersection(set((C @ dir > 0).nonzero()[0]))) + # Compute the ratio test + ra = d[r_idx,:].flatten() - C[r_idx, :] @ new_x + ra_d = C[r_idx, :] @ dir + # Initialize maximum step size. + s_star = np.inf + # Initialize blocking constraint index. + q = -1 + # Perform ratio test. + for i in range(len(ra)): + if ra_d[i] - tol > 0: + s = ra[i]/ra_d[i] + if s < s_star: + s_star = s + q = r_idx[i] + # If there is no blocking constraint (i.e., s_star >= 1) + if s_star >= 1: + # Line search to determine a step_size. + new_solution, step_size, expended_budget = self.line_search(problem, expended_budget, r, grad, new_solution, s_star, dir, alpha, beta) + # If there is a blocking constraint (i.e., s_star < 1) + else: + # Add blocking constraint to the active set. + if q not in acidx: + acidx.append(q) + # No need to do line search if s_star is 0. + if s_star > 0: + # Line search to determine a step_size. + new_solution, step_size, expended_budget = self.line_search(problem, expended_budget, r, grad, new_solution, s_star, dir, alpha, beta) + + # Append new solution. + if (problem.minmax[0] * new_solution.objectives_mean > problem.minmax[0] * best_solution.objectives_mean): + best_solution = new_solution + recommended_solns.append(new_solution) + intermediate_budgets.append(expended_budget) + + return recommended_solns, intermediate_budgets + + + def compute_search_direction(self, acidx, grad, problem, C): + ''' + Compute a search direction by solving a direction-finding quadratic subproblem at solution x. + + Arguments + --------- + acidx: list + list of indices of active constraints + grad : ndarray + the estimated objective gradient at new_solution + problem : Problem object + simulation-optimization problem to solve + C : ndarray + constraint matrix + + Returns + ------- + d : ndarray + search direction + lmbd : ndarray + Lagrange multipliers for this LP + ''' + # Define variables. + d = cp.Variable(problem.dim) + + # Define constraints. + constraints = [C[acidx, :] @ d == 0] + + # Define objective. + obj = cp.Minimize(grad @ d + 0.5 * cp.quad_form(d, np.identity(problem.dim))) + prob = cp.Problem(obj, constraints) + prob.solve() + # Get Lagrange multipliers + lmbd = prob.constraints[0].dual_value + + dir = np.array(d.value) + dir[np.abs(dir) < self.factors["tol"]] = 0 + + return dir, lmbd + + + def finite_diff(self, new_solution, BdsCheck, problem, r, stepsize = 1e-5): + ''' + Finite difference for approximating objective gradient at new_solution. + + Arguments + --------- + new_solution : Solution object + a solution to the problem + BdsCheck : ndarray + an array that checks for lower/upper bounds at each dimension + problem : Problem object + simulation-optimization problem to solve + r : int + number of replications taken at each solution + stepsize: float + step size for finite differences + + Returns + ------- + grad : ndarray + the estimated objective gradient at new_solution + ''' + lower_bound = problem.lower_bounds + upper_bound = problem.upper_bounds + fn = -1 * problem.minmax[0] * new_solution.objectives_mean + new_x = new_solution.x + # Store values for each dimension. + FnPlusMinus = np.zeros((problem.dim, 3)) + grad = np.zeros(problem.dim) + + for i in range(problem.dim): + # Initialization. + x1 = list(new_x) + x2 = list(new_x) + # Forward stepsize. + steph1 = stepsize + # Backward stepsize. + steph2 = stepsize + + # Check variable bounds. + if x1[i] + steph1 > upper_bound[i]: + steph1 = np.abs(upper_bound[i] - x1[i]) + if x2[i] - steph2 < lower_bound[i]: + steph2 = np.abs(x2[i] - lower_bound[i]) + + # Decide stepsize. + # Central diff. + if BdsCheck[i] == 0: + FnPlusMinus[i, 2] = min(steph1, steph2) + x1[i] = x1[i] + FnPlusMinus[i, 2] + x2[i] = x2[i] - FnPlusMinus[i, 2] + # Forward diff. + elif BdsCheck[i] == 1: + FnPlusMinus[i, 2] = steph1 + x1[i] = x1[i] + FnPlusMinus[i, 2] + # Backward diff. + else: + FnPlusMinus[i, 2] = steph2 + x2[i] = x2[i] - FnPlusMinus[i, 2] + x1_solution = self.create_new_solution(tuple(x1), problem) + if BdsCheck[i] != -1: + problem.simulate_up_to([x1_solution], r) + fn1 = -1 * problem.minmax[0] * x1_solution.objectives_mean + # First column is f(x+h,y). + FnPlusMinus[i, 0] = fn1 + x2_solution = self.create_new_solution(tuple(x2), problem) + if BdsCheck[i] != 1: + problem.simulate_up_to([x2_solution], r) + fn2 = -1 * problem.minmax[0] * x2_solution.objectives_mean + # Second column is f(x-h,y). + FnPlusMinus[i, 1] = fn2 + # Calculate gradient. + if BdsCheck[i] == 0: + grad[i] = (fn1 - fn2) / (2 * FnPlusMinus[i, 2]) + elif BdsCheck[i] == 1: + grad[i] = (fn1 - fn) / FnPlusMinus[i, 2] + elif BdsCheck[i] == -1: + grad[i] = (fn - fn2) / FnPlusMinus[i, 2] + + return grad + + def line_search(self, problem, expended_budget, r, grad, cur_sol, alpha_0, d, alpha, beta): + """ + A backtracking line-search along [x, x + rd] assuming all solution on the line are feasible. + + Arguments + --------- + problem : Problem object + simulation-optimization problem to solve + expended_budget: int + current expended budget + r : int + number of replications taken at each solution + grad : ndarray + objective gradient of cur_sol + cur_sol : Solution object + current solution + alpha_0 : float + maximum step size allowed + d : ndarray + search direction + alpha: float + tolerance for sufficient decrease condition + beta: float + step size reduction factor + + Returns + ------- + x_new_solution : Solution + a solution obtained by line search + step_size : float + computed step size + expended_budget : int + updated expended budget + """ + x = cur_sol.x + fx = -1 * problem.minmax[0] * cur_sol.objectives_mean + step_size = alpha_0 + count = 0 + x_new_solution = self.create_new_solution(tuple(x), problem) + while True: + if expended_budget > problem.factors["budget"]: + break + x_new = x + step_size * d + # Create a solution object for x_new. + x_new_solution = self.create_new_solution(tuple(x_new), problem) + # Use r simulated observations to estimate the objective value. + problem.simulate(x_new_solution, r) + expended_budget += r + # Check the sufficient decrease condition. + f_new = -1 * problem.minmax[0] * x_new_solution.objectives_mean + if f_new < fx + alpha * step_size * np.dot(grad, d): + break + step_size *= beta + count += 1 + return x_new_solution, step_size, expended_budget + + def find_feasible_initial(self, problem, Ae, Ai, be, bi, tol): + ''' + Find an initial feasible solution (if not user-provided) + by solving phase one simplex. + + Arguments + --------- + problem : Problem object + simulation-optimization problem to solve + C: ndarray + constraint coefficient matrix + d: ndarray + constraint coefficient vector + + Returns + ------- + x0 : ndarray + an initial feasible solution + tol: float + Floating point comparison tolerance + ''' + upper_bound = np.array(problem.upper_bounds) + lower_bound = np.array(problem.lower_bounds) + + # Define decision variables. + x = cp.Variable(problem.dim) + + # Define constraints. + constraints = [] + + if (Ae is not None) and (be is not None): + constraints.append(Ae @ x == be.ravel()) + if (Ai is not None) and (bi is not None): + constraints.append(Ai @ x <= bi.ravel()) + + # Removing redundant bound constraints. + ub_inf_idx = np.where(~np.isinf(upper_bound))[0] + if len(ub_inf_idx) > 0: + for i in ub_inf_idx: + constraints.append(x[i] <= upper_bound[i]) + lb_inf_idx = np.where(~np.isinf(lower_bound)) + if len(lb_inf_idx) > 0: + for i in lb_inf_idx: + constraints.append(x[i] >= lower_bound[i]) + + # Define objective function. + obj = cp.Minimize(0) + + # Create problem. + model = cp.Problem(obj, constraints) + + # Solve problem. + model.solve(solver = cp.SCIPY) + + # Check for optimality. + if model.status not in [cp.OPTIMAL, cp.OPTIMAL_INACCURATE] : + raise ValueError("Could not find feasible x0") + x0 = x.value + if not self._feasible(x0, problem, tol): + raise ValueError("Could not find feasible x0") + + return x0 + + def _feasible(self, x, problem, tol): + """ + Check whether a solution x is feasible to the problem. + + Arguments + --------- + x : tuple + a solution vector + problem : Problem object + simulation-optimization problem to solve + tol: float + Floating point comparison tolerance + """ + x = np.asarray(x) + lb = np.asarray(problem.lower_bounds) + ub = np.asarray(problem.upper_bounds) + res = True + if (problem.Ci is not None) and (problem.di is not None): + res = res & np.all(problem.Ci @ x <= problem.di + tol) + if (problem.Ce is not None) and (problem.de is not None): + res = res & (np.allclose(np.dot(problem.Ce, x), problem.de, rtol=0, atol=tol)) + return res & (np.all(x >= lb)) & (np.all(x <= ub)) \ No newline at end of file diff --git a/simopt/solvers/pgd.py b/simopt/solvers/pgd.py new file mode 100644 index 000000000..2725ec122 --- /dev/null +++ b/simopt/solvers/pgd.py @@ -0,0 +1,538 @@ +""" +Summary +------- +PGD: A projected gradient descent algorithm for problems with linear constraints, i.e., Ce@x = de, Ci@x <= di. +A detailed description of the solver can be found `here `_. +""" +import numpy as np +import cvxpy as cp +import warnings +warnings.filterwarnings("ignore") + +from ..base import Solver + + +class PGD(Solver): + """ + The PGD solver. + + Attributes + ---------- + name : string + name of solver + objective_type : string + description of objective types: + "single" or "multi" + constraint_type : string + description of constraints types: + "unconstrained", "box", "deterministic", "stochastic" + variable_type : string + description of variable types: + "discrete", "continuous", "mixed" + gradient_needed : bool + indicates if gradient of objective function is needed + factors : dict + changeable factors (i.e., parameters) of the solver + specifications : dict + details of each factor (for GUI, data validation, and defaults) + rng_list : list of rng.MRG32k3a objects + list of RNGs used for the solver's internal purposes + + Arguments + --------- + name : str + user-specified name for solver + fixed_factors : dict + fixed_factors of the solver + + See also + -------- + base.Solver + """ + def __init__(self, name="PGD", fixed_factors={}): + self.name = name + self.objective_type = "single" + self.constraint_type = "deterministic" + self.variable_type = "continuous" + self.gradient_needed = False + self.specifications = { + "crn_across_solns": { + "description": "use CRN across solutions?", + "datatype": bool, + "default": True + }, + "r": { + "description": "number of replications taken at each solution", + "datatype": int, + "default": 30 + }, + "alpha": { + "description": "tolerance for sufficient decrease condition", + "datatype": float, + "default": 0.2 + }, + "beta": { + "description": "step size reduction factor in line search", + "datatype": float, + "default": 0.9 + }, + "alpha_max": { + "description": "maximum step size", + "datatype": float, + "default": 10.0 + }, + "lambda": { + "description": "magnifying factor for r inside the finite difference function", + "datatype": int, + "default": 2 + }, + "tol": { + "description": "floating point comparison tolerance", + "datatype": float, + "default": 1e-7 + }, + "finite_diff_step": { + "description": "step size for finite difference", + "datatype": float, + "default": 1e-5 + } + } + self.check_factor_list = { + "crn_across_solns": self.check_crn_across_solns, + "r": self.check_r, + "alpha": self.check_alpha, + "beta": self.check_beta, + "alpha_max": self.check_alpha_max, + "lambda": self.check_lambda, + "tol": self.check_tol, + "finite_diff_step": self.check_finite_diff_step + } + super().__init__(fixed_factors) + + def check_r(self): + return self.factors["r"] > 0 + + def check_alpha(self): + return self.factors["alpha"] > 0 + + def check_beta(self): + return self.factors["beta"] > 0 & self.factors["beta"] < 1 + + def check_alpha_max(self): + return self.factors["alpha_max"] > 0 + + def check_lambda(self): + return self.factors["lambda"] > 0 + + def check_tol(self): + return self.factors["tol"] > 0 + + def check_finite_diff_step(self): + return self.factors["finite_diff_step"] > 0 + + def solve(self, problem): + """ + Run a single macroreplication of a solver on a problem. + + Arguments + --------- + problem : Problem object + simulation-optimization problem to solve + crn_across_solns : bool + indicates if CRN are used when simulating different solutions + + Returns + ------- + recommended_solns : list of Solution objects + list of solutions recommended throughout the budget + intermediate_budgets : list of ints + list of intermediate budgets when recommended solutions changes + """ + recommended_solns = [] + intermediate_budgets = [] + expended_budget = 0 + + # Default values. + r = self.factors["r"] + alpha = self.factors["alpha"] + beta = self.factors["beta"] + tol = self.factors["tol"] + max_step = self.factors["alpha_max"] + + # Upper bound and lower bound. + lower_bound = np.array(problem.lower_bounds) + upper_bound = np.array(problem.upper_bounds) + + # Input inequality and equlaity constraint matrix and vector. + # Cix <= di + # Cex = de + Ci = problem.Ci + di = problem.di + Ce = problem.Ce + de = problem.de + + # Checker for whether the problem is unconstrained. + unconstr_flag = (Ce is None) & (Ci is None) & (di is None) & (de is None) & (all(np.isinf(lower_bound))) & (all(np.isinf(upper_bound))) + + # Start with the initial solution. + new_solution = self.create_new_solution(problem.factors["initial_solution"], problem) + new_x = new_solution.x + + # If the initial solution is not feasible, generate one using phase one simplex. + if (not unconstr_flag) & (not self._feasible(new_x, problem, tol)): + new_x = self.find_feasible_initial(problem, Ce, Ci, de, di, tol) + new_solution = self.create_new_solution(tuple(new_x), problem) + + # Use r simulated observations to estimate the objective value. + problem.simulate(new_solution, r) + expended_budget += r + best_solution = new_solution + recommended_solns.append(new_solution) + intermediate_budgets.append(expended_budget) + + while expended_budget < problem.factors["budget"]: + new_x = new_solution.x + # Check variable bounds. + forward = np.isclose(new_x, lower_bound, atol = tol).astype(int) + backward = np.isclose(new_x, upper_bound, atol = tol).astype(int) + # BdsCheck: 1 stands for forward, -1 stands for backward, 0 means central diff. + BdsCheck = np.subtract(forward, backward) + + if problem.gradient_available: + # Use IPA gradient if available. + grad = -1 * problem.minmax[0] * new_solution.objectives_gradients_mean[0] + else: + # Use finite difference to estimate gradient if IPA gradient is not available. + grad = self.finite_diff(new_solution, BdsCheck, problem, r, stepsize = self.factors["finite_diff_step"]) + expended_budget += (2 * problem.dim - np.sum(BdsCheck != 0)) * r + # A while loop to prevent zero gradient. + while np.all((grad == 0)): + if expended_budget > problem.factors["budget"]: + break + grad = self.finite_diff(new_solution, BdsCheck, problem, r) + expended_budget += (2 * problem.dim - np.sum(BdsCheck != 0)) * r + # Update r after each iteration. + r = int(self.factors["lambda"] * r) + + # Get search direction by taking negative normalized gradient. + dir = -grad / np.linalg.norm(grad) + + # Get a temp solution. + temp_x = new_x + max_step * dir + + # Check feasibility of temp_x. + if unconstr_flag or self._feasible(temp_x, problem, tol): + # Perform line search. + new_solution, step_size, expended_budget = self.line_search(problem, expended_budget, r, grad, new_solution, max_step, dir, alpha, beta) + # Update maximum step size for the next iteration. + max_step = step_size + else: + # If not feasible, project temp_x back to the feasible set. + proj_x = self.project_grad(problem, temp_x, Ce, Ci, de, di) + # Get new search direction based on projection. + dir = proj_x - new_x + # Perform line search. + new_solution, step_size, expended_budget = self.line_search(problem, expended_budget, r, grad, new_solution, 1, dir, alpha, beta) + # Update maximum step size for the next iteration. + max_step = step_size + + # Append new solution. + if (problem.minmax[0] * new_solution.objectives_mean > problem.minmax[0] * best_solution.objectives_mean): + best_solution = new_solution + recommended_solns.append(new_solution) + intermediate_budgets.append(expended_budget) + + return recommended_solns, intermediate_budgets + + + def finite_diff(self, new_solution, BdsCheck, problem, r, stepsize = 1e-5): + ''' + Finite difference for approximating objective gradient at new_solution. + + Arguments + --------- + new_solution : Solution object + a solution to the problem + BdsCheck : ndarray + an array that checks for lower/upper bounds at each dimension + problem : Problem object + simulation-optimization problem to solve + r : int + number of replications taken at each solution + stepsize: float + step size for finite differences + + Returns + ------- + grad : ndarray + the estimated objective gradient at new_solution + ''' + lower_bound = problem.lower_bounds + upper_bound = problem.upper_bounds + fn = -1 * problem.minmax[0] * new_solution.objectives_mean + new_x = new_solution.x + # Store values for each dimension. + FnPlusMinus = np.zeros((problem.dim, 3)) + grad = np.zeros(problem.dim) + + for i in range(problem.dim): + # Initialization. + x1 = list(new_x) + x2 = list(new_x) + # Forward stepsize. + steph1 = stepsize + # Backward stepsize. + steph2 = stepsize + + # Check variable bounds. + if x1[i] + steph1 > upper_bound[i]: + steph1 = np.abs(upper_bound[i] - x1[i]) + if x2[i] - steph2 < lower_bound[i]: + steph2 = np.abs(x2[i] - lower_bound[i]) + + # Decide stepsize. + # Central diff. + if BdsCheck[i] == 0: + FnPlusMinus[i, 2] = min(steph1, steph2) + x1[i] = x1[i] + FnPlusMinus[i, 2] + x2[i] = x2[i] - FnPlusMinus[i, 2] + # Forward diff. + elif BdsCheck[i] == 1: + FnPlusMinus[i, 2] = steph1 + x1[i] = x1[i] + FnPlusMinus[i, 2] + # Backward diff. + else: + FnPlusMinus[i, 2] = steph2 + x2[i] = x2[i] - FnPlusMinus[i, 2] + x1_solution = self.create_new_solution(tuple(x1), problem) + if BdsCheck[i] != -1: + problem.simulate_up_to([x1_solution], r) + fn1 = -1 * problem.minmax[0] * x1_solution.objectives_mean + # First column is f(x+h,y). + FnPlusMinus[i, 0] = fn1 + x2_solution = self.create_new_solution(tuple(x2), problem) + if BdsCheck[i] != 1: + problem.simulate_up_to([x2_solution], r) + fn2 = -1 * problem.minmax[0] * x2_solution.objectives_mean + # Second column is f(x-h,y). + FnPlusMinus[i, 1] = fn2 + # Calculate gradient. + if BdsCheck[i] == 0: + grad[i] = (fn1 - fn2) / (2 * FnPlusMinus[i, 2]) + elif BdsCheck[i] == 1: + grad[i] = (fn1 - fn) / FnPlusMinus[i, 2] + elif BdsCheck[i] == -1: + grad[i] = (fn - fn2) / FnPlusMinus[i, 2] + + return grad + + def _feasible(self, x, problem, tol): + """ + Check whether a solution x is feasible to the problem. + + Arguments + --------- + x : tuple + a solution vector + problem : Problem object + simulation-optimization problem to solve + tol: float + Floating point comparison tolerance + """ + x = np.asarray(x) + lb = np.asarray(problem.lower_bounds) + ub = np.asarray(problem.upper_bounds) + res = True + if (problem.Ci is not None) and (problem.di is not None): + res = res & np.all(problem.Ci @ x <= problem.di + tol) + if (problem.Ce is not None) and (problem.de is not None): + res = res & (np.allclose(np.dot(problem.Ce, x), problem.de, rtol=0, atol=tol)) + return res & (np.all(x >= lb)) & (np.all(x <= ub)) + + def project_grad(self, problem, x, Ae, Ai, be, bi): + """ + Project the vector x onto the hyperplane H: Ae x = be, Ai x <= bi by solving a quadratic projection problem: + + min d^Td + s.t. Ae(x + d) = be + Ai(x + d) <= bi + (x + d) >= lb + (x + d) <= ub + + Arguments + --------- + problem : Problem object + simulation-optimization problem to solve + x : ndarray + vector to be projected + Ae: ndarray + equality constraint coefficient matrix + be: ndarray + equality constraint coefficient vector + Ai: ndarray + inequality constraint coefficient matrix + bi: ndarray + inequality constraint coefficient vector + Returns + ------- + x_new : ndarray + the projected vector + """ + # Define variables. + d = cp.Variable(problem.dim) + + # Define objective. + obj = cp.Minimize(cp.quad_form(d, np.identity(problem.dim))) + + # Define constraints. + constraints = [] + if (Ae is not None) and (be is not None): + constraints.append(Ae @ (x + d) == be.ravel()) + if (Ai is not None) and (bi is not None): + constraints.append(Ai @ (x + d) <= bi.ravel()) + + upper_bound = np.array(problem.upper_bounds) + lower_bound = np.array(problem.lower_bounds) + + # Removing redundant bound constraints. + ub_inf_idx = np.where(~np.isinf(upper_bound))[0] + if len(ub_inf_idx) > 0: + for i in ub_inf_idx: + constraints.append((x + d)[i] <= upper_bound[i]) + lb_inf_idx = np.where(~np.isinf(lower_bound))[0] + if len(lb_inf_idx) > 0: + for i in lb_inf_idx: + constraints.append((x + d)[i] >= lower_bound[i]) + + # Form and solve problem. + prob = cp.Problem(obj, constraints) + prob.solve() + + dir = d.value + + # Get the projected vector. + x_new = x + dir + + # Avoid floating point error + x_new[np.abs(x_new) < self.factors["tol"]] = 0 + + return x_new + + def line_search(self, problem, expended_budget, r, grad, cur_sol, alpha_0, d, alpha, beta): + """ + A backtracking line-search along [x, x + rd] assuming all solution on the line are feasible. + + Arguments + --------- + problem : Problem object + simulation-optimization problem to solve + expended_budget: int + current expended budget + r : int + number of replications taken at each solution + grad : ndarray + objective gradient of cur_sol + cur_sol : Solution object + current solution + alpha_0 : float + maximum step size allowed + d : ndarray + search direction + alpha: float + tolerance for sufficient decrease condition + beta: float + step size reduction factor + + Returns + ------- + x_new_solution : Solution + a solution obtained by line search + step_size : float + computed step size + expended_budget : int + updated expended budget + """ + x = cur_sol.x + fx = -1 * problem.minmax[0] * cur_sol.objectives_mean + step_size = alpha_0 + count = 0 + x_new_solution = cur_sol + while True: + if expended_budget > problem.factors["budget"]: + break + x_new = x + step_size * d + # Create a solution object for x_new. + x_new_solution = self.create_new_solution(tuple(x_new), problem) + # Use r simulated observations to estimate the objective value. + problem.simulate(x_new_solution, r) + expended_budget += r + # Check the sufficient decrease condition. + f_new = -1 * problem.minmax[0] * x_new_solution.objectives_mean + if f_new < fx + alpha * step_size * np.dot(grad, d): + break + step_size *= beta + count += 1 + return x_new_solution, step_size, expended_budget, + + def find_feasible_initial(self, problem, Ae, Ai, be, bi, tol): + ''' + Find an initial feasible solution (if not user-provided) + by solving phase one simplex. + + Arguments + --------- + problem : Problem object + simulation-optimization problem to solve + C: ndarray + constraint coefficient matrix + d: ndarray + constraint coefficient vector + + Returns + ------- + x0 : ndarray + an initial feasible solution + tol: float + Floating point comparison tolerance + ''' + upper_bound = np.array(problem.upper_bounds) + lower_bound = np.array(problem.lower_bounds) + + # Define decision variables. + x = cp.Variable(problem.dim) + + # Define constraints. + constraints = [] + + if (Ae is not None) and (be is not None): + constraints.append(Ae @ x == be.ravel()) + if (Ai is not None) and (bi is not None): + constraints.append(Ai @ x <= bi.ravel()) + + # Removing redundant bound constraints. + ub_inf_idx = np.where(~np.isinf(upper_bound))[0] + if len(ub_inf_idx) > 0: + for i in ub_inf_idx: + constraints.append(x[i] <= upper_bound[i]) + lb_inf_idx = np.where(~np.isinf(lower_bound))[0] + if len(lb_inf_idx) > 0: + for i in lb_inf_idx: + constraints.append(x[i] >= lower_bound[i]) + + # Define objective function. + obj = cp.Minimize(0) + + # Create problem. + model = cp.Problem(obj, constraints) + + # Solve problem. + model.solve(solver = cp.SCIPY) + + # Check for optimality. + if model.status not in [cp.OPTIMAL, cp.OPTIMAL_INACCURATE] : + raise ValueError("Could not find feasible x0") + x0 = x.value + if not self._feasible(x0, problem, tol): + raise ValueError("Could not find feasible x0") + + return x0 diff --git a/simopt/solvers/pgdss.py b/simopt/solvers/pgdss.py new file mode 100644 index 000000000..6d73fa34a --- /dev/null +++ b/simopt/solvers/pgdss.py @@ -0,0 +1,512 @@ +""" +Summary +------- +PGD-SS: A projected gradient descent algorithm with adaptive step search +for problems with linear constraints, i.e., Ce@x = de, Ci@x <= di. +A detailed description of the solver can be found `here `_. +""" +import numpy as np +import cvxpy as cp +import warnings +warnings.filterwarnings("ignore") + +from ..base import Solver + + +class PGDSS(Solver): + """ + The PGD solver with adaptive step search. + + Attributes + ---------- + name : string + name of solver + objective_type : string + description of objective types: + "single" or "multi" + constraint_type : string + description of constraints types: + "unconstrained", "box", "deterministic", "stochastic" + variable_type : string + description of variable types: + "discrete", "continuous", "mixed" + gradient_needed : bool + indicates if gradient of objective function is needed + factors : dict + changeable factors (i.e., parameters) of the solver + specifications : dict + details of each factor (for GUI, data validation, and defaults) + rng_list : list of rng.MRG32k3a objects + list of RNGs used for the solver's internal purposes + + Arguments + --------- + name : str + user-specified name for solver + fixed_factors : dict + fixed_factors of the solver + + See also + -------- + base.Solver + """ + def __init__(self, name="PGD-SS", fixed_factors={}): + self.name = name + self.objective_type = "single" + self.constraint_type = "deterministic" + self.variable_type = "continuous" + self.gradient_needed = False + self.specifications = { + "crn_across_solns": { + "description": "use CRN across solutions?", + "datatype": bool, + "default": True + }, + "r": { + "description": "number of replications taken at each solution", + "datatype": int, + "default": 30 + }, + "theta": { + "description": "constant in the Armijo condition", + "datatype": int, + "default": 0.2 + }, + "gamma": { + "description": "constant for shrinking the step size", + "datatype": int, + "default": 0.8 + }, + "alpha_max": { + "description": "maximum step size", + "datatype": int, + "default": 10 + }, + "alpha_0": { + "description": "initial step size", + "datatype": int, + "default": 1 + }, + "epsilon_f": { + "description": "additive constant in the Armijo condition", + "datatype": int, + "default": 1e-3 # In the paper, this value is estimated for every epoch but a value > 0 is justified in practice. + }, + "lambda": { + "description": "magnifying factor for r inside the finite difference function", + "datatype": int, + "default": 2 + }, + "tol": { + "description": "floating point comparison tolerance", + "datatype": float, + "default": 1e-7 + }, + "finite_diff_step": { + "description": "step size for finite difference", + "datatype": float, + "default": 1e-5 + } + + } + self.check_factor_list = { + "crn_across_solns": self.check_crn_across_solns, + "r": self.check_r, + "theta": self.check_theta, + "gamma": self.check_gamma, + "alpha_max": self.check_alpha_max, + "alpha_0": self.check_alpha_0, + "epsilon_f": self.check_epsilon_f, + "lambda": self.check_lambda, + "tol": self.check_tol, + "finite_diff_step": self.check_finite_diff_step + } + super().__init__(fixed_factors) + + def check_r(self): + return self.factors["r"] > 0 + + def check_theta(self): + return self.factors["theta"] > 0 & self.factors["theta"] < 1 + + def check_gamma(self): + return self.factors["gamma"] > 0 & self.factors["gamma"] < 1 + + def check_alpha_max(self): + return self.factors["alpha_max"] > 0 + + def check_alpha_0(self): + return self.factors["alpha_0"] > 0 + + def check_epsilon_f(self): + return self.factors["epsilon_f"] > 0 + + def check_tol(self): + return self.factors["tol"] > 0 + + def check_lambda(self): + return self.factors["lambda"] > 0 + + def check_finite_diff_step(self): + return self.factors["finite_diff_step"] > 0 + + def solve(self, problem): + """ + Run a single macroreplication of a solver on a problem. + + Arguments + --------- + problem : Problem object + simulation-optimization problem to solve + crn_across_solns : bool + indicates if CRN are used when simulating different solutions + + Returns + ------- + recommended_solns : list of Solution objects + list of solutions recommended throughout the budget + intermediate_budgets : list of ints + list of intermediate budgets when recommended solutions changes + """ + recommended_solns = [] + intermediate_budgets = [] + expended_budget = 0 + + # Default values. + r = self.factors["r"] + tol = self.factors["tol"] + theta = self.factors["theta"] + gamma = self.factors["gamma"] + alpha_max = self.factors["alpha_max"] + alpha_0 = self.factors["alpha_0"] + epsilon_f = self.factors["epsilon_f"] + + # Upper bound and lower bound. + lower_bound = np.array(problem.lower_bounds) + upper_bound = np.array(problem.upper_bounds) + + # Initialize stepsize. + alpha = alpha_0 + + # Input inequality and equlaity constraint matrix and vector. + # Cix <= di + # Cex = de + Ci = problem.Ci + di = problem.di + Ce = problem.Ce + de = problem.de + + # Checker for whether the problem is unconstrained. + unconstr_flag = (Ce is None) & (Ci is None) & (di is None) & (de is None) & (all(np.isinf(lower_bound))) & (all(np.isinf(upper_bound))) + + # Start with the initial solution. + new_solution = self.create_new_solution(problem.factors["initial_solution"], problem) + new_x = new_solution.x + + # If the initial solution is not feasible, generate one using phase one simplex. + if (not unconstr_flag) & (not self._feasible(new_x, problem, tol)): + new_x = self.find_feasible_initial(problem, Ce, Ci, de, di, tol) + new_solution = self.create_new_solution(tuple(new_x), problem) + + # Use r simulated observations to estimate the objective value. + problem.simulate(new_solution, r) + expended_budget += r + best_solution = new_solution + recommended_solns.append(new_solution) + intermediate_budgets.append(expended_budget) + + while expended_budget < problem.factors["budget"]: + new_x = new_solution.x + # Check variable bounds. + forward = np.isclose(new_x, lower_bound, atol = tol).astype(int) + backward = np.isclose(new_x, upper_bound, atol = tol).astype(int) + # BdsCheck: 1 stands for forward, -1 stands for backward, 0 means central diff. + BdsCheck = np.subtract(forward, backward) + + if problem.gradient_available: + # Use IPA gradient if available. + grad = -1 * problem.minmax[0] * new_solution.objectives_gradients_mean[0] + else: + # Use finite difference to estimate gradient if IPA gradient is not available. + grad = self.finite_diff(new_solution, BdsCheck, problem, r, stepsize = self.factors["finite_diff_step"]) + expended_budget += (2 * problem.dim - np.sum(BdsCheck != 0)) * r + # A while loop to prevent zero gradient. + while np.all((grad == 0)): + if expended_budget > problem.factors["budget"]: + break + grad = self.finite_diff(new_solution, BdsCheck, problem, r) + expended_budget += (2 * problem.dim - np.sum(BdsCheck != 0)) * r + # Update r after each iteration. + r = int(self.factors["lambda"] * r) + + # Get search direction by taking negative normalized gradient. + dir = -grad / np.linalg.norm(grad) + + # Get a temp solution. + temp_x = new_x + alpha * dir + + if unconstr_flag or self._feasible(temp_x, problem, tol): + candidate_solution = self.create_new_solution(tuple(temp_x), problem) + else: + # If not feasible, project temp_x back to the feasible set. + proj_x = self.project_grad(problem, temp_x, Ce, Ci, de, di) + candidate_solution = self.create_new_solution(tuple(proj_x), problem) + # Get new search direction based on projection. + dir = proj_x - new_x + + # Use r simulated observations to estimate the objective value. + problem.simulate(candidate_solution, r) + expended_budget += r + + # Check the modified Armijo condition for sufficient decrease. + if (-1 * problem.minmax[0] * candidate_solution.objectives_mean) <= ( + -1 * problem.minmax[0] * new_solution.objectives_mean + alpha * theta * np.dot(grad, dir) + 2 * epsilon_f): + # Successful step + new_solution = candidate_solution + # Enlarge step size. + alpha = min(alpha_max, alpha / gamma) + else: + # Unsuccessful step - reduce step size. + alpha = gamma * alpha + + # Append new solution. + if (problem.minmax[0] * new_solution.objectives_mean > problem.minmax[0] * best_solution.objectives_mean): + best_solution = new_solution + recommended_solns.append(new_solution) + intermediate_budgets.append(expended_budget) + + return recommended_solns, intermediate_budgets + + + def finite_diff(self, new_solution, BdsCheck, problem, r, stepsize = 1e-5): + ''' + Finite difference for approximating objective gradient at new_solution. + + Arguments + --------- + new_solution : Solution object + a solution to the problem + BdsCheck : ndarray + an array that checks for lower/upper bounds at each dimension + problem : Problem object + simulation-optimization problem to solve + r : int + number of replications taken at each solution + stepsize: float + step size for finite differences + + Returns + ------- + grad : ndarray + the estimated objective gradient at new_solution + ''' + lower_bound = problem.lower_bounds + upper_bound = problem.upper_bounds + fn = -1 * problem.minmax[0] * new_solution.objectives_mean + new_x = new_solution.x + # Store values for each dimension. + FnPlusMinus = np.zeros((problem.dim, 3)) + grad = np.zeros(problem.dim) + + for i in range(problem.dim): + # Initialization. + x1 = list(new_x) + x2 = list(new_x) + # Forward stepsize. + steph1 = stepsize + # Backward stepsize. + steph2 = stepsize + + # Check variable bounds. + if x1[i] + steph1 > upper_bound[i]: + steph1 = np.abs(upper_bound[i] - x1[i]) + if x2[i] - steph2 < lower_bound[i]: + steph2 = np.abs(x2[i] - lower_bound[i]) + + # Decide stepsize. + # Central diff. + if BdsCheck[i] == 0: + FnPlusMinus[i, 2] = min(steph1, steph2) + x1[i] = x1[i] + FnPlusMinus[i, 2] + x2[i] = x2[i] - FnPlusMinus[i, 2] + # Forward diff. + elif BdsCheck[i] == 1: + FnPlusMinus[i, 2] = steph1 + x1[i] = x1[i] + FnPlusMinus[i, 2] + # Backward diff. + else: + FnPlusMinus[i, 2] = steph2 + x2[i] = x2[i] - FnPlusMinus[i, 2] + x1_solution = self.create_new_solution(tuple(x1), problem) + if BdsCheck[i] != -1: + problem.simulate_up_to([x1_solution], r) + fn1 = -1 * problem.minmax[0] * x1_solution.objectives_mean + # First column is f(x+h,y). + FnPlusMinus[i, 0] = fn1 + x2_solution = self.create_new_solution(tuple(x2), problem) + if BdsCheck[i] != 1: + problem.simulate_up_to([x2_solution], r) + fn2 = -1 * problem.minmax[0] * x2_solution.objectives_mean + # Second column is f(x-h,y). + FnPlusMinus[i, 1] = fn2 + # Calculate gradient. + if BdsCheck[i] == 0: + grad[i] = (fn1 - fn2) / (2 * FnPlusMinus[i, 2]) + elif BdsCheck[i] == 1: + grad[i] = (fn1 - fn) / FnPlusMinus[i, 2] + elif BdsCheck[i] == -1: + grad[i] = (fn - fn2) / FnPlusMinus[i, 2] + + return grad + + def _feasible(self, x, problem, tol): + """ + Check whether a solution x is feasible to the problem. + + Arguments + --------- + x : tuple + a solution vector + problem : Problem object + simulation-optimization problem to solve + tol: float + Floating point comparison tolerance + """ + x = np.asarray(x) + lb = np.asarray(problem.lower_bounds) + ub = np.asarray(problem.upper_bounds) + res = True + if (problem.Ci is not None) and (problem.di is not None): + res = res & np.all(problem.Ci @ x <= problem.di + tol) + if (problem.Ce is not None) and (problem.de is not None): + res = res & (np.allclose(np.dot(problem.Ce, x), problem.de, rtol=0, atol=tol)) + return res & (np.all(x >= lb)) & (np.all(x <= ub)) + + def project_grad(self, problem, x, Ae, Ai, be, bi): + """ + Project the vector x onto the hyperplane H: Ae x = be, Ai x <= bi by solving a quadratic projection problem: + + min d^Td + s.t. Ae(x + d) = be + Ai(x + d) <= bi + (x + d) >= lb + (x + d) <= ub + + Arguments + --------- + problem : Problem object + simulation-optimization problem to solve + x : ndarray + vector to be projected + Ae: ndarray + equality constraint coefficient matrix + be: ndarray + equality constraint coefficient vector + Ai: ndarray + inequality constraint coefficient matrix + bi: ndarray + inequality constraint coefficient vector + Returns + ------- + x_new : ndarray + the projected vector + """ + # Define variables. + d = cp.Variable(problem.dim) + + # Define objective. + obj = cp.Minimize(cp.quad_form(d, np.identity(problem.dim))) + + # Define constraints. + constraints = [] + if (Ae is not None) and (be is not None): + constraints.append(Ae @ (x + d) == be.ravel()) + if (Ai is not None) and (bi is not None): + constraints.append(Ai @ (x + d) <= bi.ravel()) + + upper_bound = np.array(problem.upper_bounds) + lower_bound = np.array(problem.lower_bounds) + # Removing redundant bound constraints. + ub_inf_idx = np.where(~np.isinf(upper_bound))[0] + if len(ub_inf_idx) > 0: + for i in ub_inf_idx: + constraints.append((x + d)[i] <= upper_bound[i]) + lb_inf_idx = np.where(~np.isinf(lower_bound))[0] + if len(lb_inf_idx) > 0: + for i in lb_inf_idx: + constraints.append((x + d)[i] >= lower_bound[i]) + + # Form and solve problem. + prob = cp.Problem(obj, constraints) + prob.solve() + + # Get the projected vector. + x_new = x + d.value + + # Avoid floating point error + x_new[np.abs(x_new) < self.factors["tol"]] = 0 + + return x_new + + def find_feasible_initial(self, problem, Ae, Ai, be, bi, tol): + ''' + Find an initial feasible solution (if not user-provided) + by solving phase one simplex. + + Arguments + --------- + problem : Problem object + simulation-optimization problem to solve + C: ndarray + constraint coefficient matrix + d: ndarray + constraint coefficient vector + + Returns + ------- + x0 : ndarray + an initial feasible solution + tol: float + Floating point comparison tolerance + ''' + upper_bound = np.array(problem.upper_bounds) + lower_bound = np.array(problem.lower_bounds) + + # Define decision variables. + x = cp.Variable(problem.dim) + + # Define constraints. + constraints = [] + + if (Ae is not None) and (be is not None): + constraints.append(Ae @ x == be.ravel()) + if (Ai is not None) and (bi is not None): + constraints.append(Ai @ x <= bi.ravel()) + + # Removing redundant bound constraints. + ub_inf_idx = np.where(~np.isinf(upper_bound))[0] + if len(ub_inf_idx) > 0: + for i in ub_inf_idx: + constraints.append(x[i] <= upper_bound[i]) + lb_inf_idx = np.where(~np.isinf(lower_bound))[0] + if len(lb_inf_idx) > 0: + for i in lb_inf_idx: + constraints.append(x[i] >= lower_bound[i]) + + # Define objective function. + obj = cp.Minimize(0) + + # Create problem. + model = cp.Problem(obj, constraints) + + # Solve problem. + model.solve(solver = cp.SCIPY) + + # Check for optimality. + if model.status not in [cp.OPTIMAL, cp.OPTIMAL_INACCURATE] : + raise ValueError("Could not find feasible x0") + x0 = x.value + if not self._feasible(x0, problem, tol): + raise ValueError("Could not find feasible x0") + + return x0