Merge pull request #55 from esahon/rjt

Add rooted junction tree models. Note: the documentation must still be updated before a release

Author: solliolli

Project.toml

@@ -11,6 +11,7 @@ JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
+OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 
 [compat]
 DataFrames = "1.3"

docs/src/api.md

@@ -58,6 +58,7 @@ UtilityMatrix(::InfluenceDiagram, ::Name)
 add_utilities!
 generate_arcs!
 generate_diagram!
+RJT
 indices
 I_j_indices
 indices_in_vector
@@ -93,10 +94,18 @@ conditional_value_at_risk(::Model, ::InfluenceDiagram, ::PathCompatibilityVariab
 
 ### Decision Strategy from Variables
 ```@docs
-LocalDecisionStrategy(::Node, ::Vector{VariableRef})
+LocalDecisionStrategy(::Node, ::Array{VariableRef})
 DecisionStrategy(::InfluenceDiagram, ::OrderedDict{Name, DecisionProgramming.DecisionVariable})
 ```
 
+### RJT model
+```@docs
+RJTVariables
+expected_value(::Model, ::InfluenceDiagram, ::DecisionProgramming.RJTVariables)
+conditional_value_at_risk(::Model, ::InfluenceDiagram, ::DecisionProgramming.RJTVariables, ::Float64)
+generate_model
+```
+
 ## `heuristics.jl`
 ### Single policy update
 ```@docs

docs/src/decision-programming/RJT-model.md

@@ -0,0 +1,108 @@
+# [RJT model](@id RJT-model)
+## Introduction
+Influence diagrams can be represented as directed rooted trees composed of clusters, which can be transformed into gradual rooted junction trees (RJTs) by imposing additional constraints. These can then be used to formulate an optimization model. Solving for optimal decision strategies using these formulations can be done with significantly less computing time than for full path based formulations. Using RJT based formulations is thus generally preferable.
+
+The explanations for RJT construction and RJT based model formulation largely follow that of Herrala et al. (2024). [^1]
+
+## Converting influence diagrams to RJTs
+
+An influence diagram $G = (N, A)$ can be represented as a directed rooted tree $\mathscr{G} = (\mathscr{V}, \mathscr{A})$ composed of *clusters* $C \isin V$, which are subsets of the nodes of the ID, that is, $C \subset \mathscr{V}$. Both $G$ and $\mathscr{G}$ are directed acyclic graphs whose vertices are connected with directed arcs in $A$ and $\mathscr{A}$ , respectively. The main difference between these diagrams lies in the nature of the vertices. In an influence diagram, the set of nodes $N$ consists of individual chance events, decisions and consequences, while the clusters in $\mathscr{V}$ comprise multiple nodes, hence the notational distinction between $N$ and $\mathscr{V}$.
+
+In order to reformulate this tree into a MIP model, additional constraints need to be imposed, making $\mathscr{G}$ a *gradual rooted junction tree*. A directed rooted tree $\mathscr{G} = (\mathscr{V}, \mathscr{A})$ consisting of clusters $C \in \mathscr{V}$ of nodes $j \in N$ is a gradual rooted junction tree corresponding to the influence diagram $G$ if:
+
+ 1. Given two clusters $C_1$ and $C_2$ in the junction tree, any cluster $C$ on the unique undirected path between $C_1$ and $C_2$ satisfies $C_1 \cap C_2 \subset C$. 
+ 2. Each cluster $C \in \mathscr{V}$ is the root cluster of exactly one node $j \in N$ (that is, the root of the subgraph induced by the clusters with node $j$) and all nodes $j \in N$ appear in at least one of the clusters.
+ 3. For each cluster, $I(j) \in C_j$, where $C_j$ is the _root cluster_ of $j \in N$.
+
+A rooted tree satisfying condition (1) is said to satisfy the *running intersection property*. This condition is sufficient for making $\mathscr{G}$ a rooted junction tree (RJT). In addition, as a consequence of condition (2), we see that a gradual RJT has as many clusters as the original influence diagram has nodes, and each node $j \in N$ can be thought as corresponding to one of the clusters $C \in \mathscr{V}$. Because of this, we refer to clusters using the corresponding nodes $j \in N$ in the influence diagram as the *root cluster* of node $j \in N$, which is denoted as $C_j \in \mathscr{V}$.
+
+An example of influence diagram (upper figure) conversion to RJT (lower figure) for pig breeding problem with $N=4$ is shown below.
+
+<img src="figures/pig_breeding_N=4.jpg" width="400"> <br>
+
+<img src="figures/pig_breeding_rjt_N=4.jpg" width="330">
+
+## Formulating an optimization problem based on gradual RJT
+
+Formulating an optimization model based on the gradual RJT representation starts by introducing a vector of moments $\mu_{C_j}$ for each cluster $C_j, \ j \in N$. Parmentier et al. (2020) [^2] show that for RJTs, we can impose constraints so that these become moments of a distribution $\mu_N$ that factorizes according to $G(N,A)$. The joint distribution $\mathbb{P}$ is said to factorize [^3] according to $G$ if $\mathbb{P}$ can be expressed as 
+
+$$\mathbb{P}(X_N = s_N) = \prod_{j \in N}\mathbb{P}(X_j=s_j \mid X_{I(j)}=s_{I(j)}).\tag{1}$$
+
+In the formulation, $\mu_{C_j}(s_{C_j})$ represents the probability of the nodes within the cluster $C_j$ being in states $s_{C_j}$ and condition (3) of definition above ensures that $\mathbb{P}(X_j=s_j \mid X_{I(j)}=s_{I(j)})$ can thus be obtained from $\mu_{C_j}(s_{C_j})$ for each $j \in N$.
+
+Variable definitions are given in the influence diagram section.
+
+## MIP model formulation
+
+### Objective function
+
+Based on the observations above, the MIP model can be formulated. The objective function becomes:
+
+$${\text{maximize}} \sum_{j \in N^V} \sum_{s_{C_j} \in S_{C_j}} \mu_{C_j}(s_{C_j}) u_{C_j}(s_{C_j})\tag{2}$$
+
+This represents an expected utility maximization problem where $u_{C_j}$ represent the utility values associated with different realizations of the nodes within the cluster $C_j$.
+
+### Constraints
+#### Decision variables
+
+Only one positive decision can be taken given the same information:
+
+$$\sum_{s_j \in S_{I(j)}}\delta (s_j \mid s_{I(j)})=1, \ \forall j \in N^D, s_{I(j)} \in S_{I(j)}\tag{3}$$
+
+#### μ-variables
+
+μ-variables need to represent valid probability distributions with probabilities summing up to 1. Thus, we get constraint:
+
+$$\sum_{s_{C_j} \in S_{C_j}} \mu_{C_j}(s_{C_j}) = 1, \ \forall j \in N\\\tag{3}$$
+
+We use moments $\mu_{\overline{C}_j}$ (μ_bar in the code) to ease the notation. These are defined as
+
+$$\mu_{\overline{C}_j}(s_{\overline{C}_j}) = \sum_{s_j \in S_j} \mu_{C_j}(s_{C_j})\tag{4}$$
+
+
+representing the marginal distribution for cluster $C_j$ with the node $j$ marginalized out. $\overline{C}_j = C_j \setminus j$ is used for notational brevity.
+
+
+#### Local consistency
+
+The following constraint enforces local consistency between adjacent clusters, meaning that for a pair $C_i, C_j$ of adjacent clusters, the marginal distribution for the nodes in both $C_i$ and $C_j$ (that is, $C_i \cap C_j$) must be the same when obtained from either $C_i$ or $C_j$. This is formulated as:
+
+$$\sum_{\substack{s_{C_i} \in S_{C_i}, \\ s_{C_i \cap C_j} = s^*_{C_i \cap C_j}}} \mu_{C_i}(s_{C_i})= \sum_{\substack{s_{C_j} \in S_{C_j}, \\ s_{C_i \cap C_j} = s^*_{C_i \cap C_j}}} \mu_{C_j}(s_{C_j}),  \nonumber\\ \tag{5}$$
+
+$$\qquad\qquad\qquad\qquad  \ \forall (C_i,C_j) \in \mathscr{A}, s^*_{C_i \cap C_j} \in S_{C_i \cap C_j}\\$$
+
+#### Conditional probabilities and decision strategies
+
+Moments $\mu_{\overline{C}_j}$ defined above are used here. We get constraints:
+
+$$\mu_{C_j}(s_{C_j}) = \mu_{\overline{C}_j}(s_{\overline{C}_j}) p(s_j \mid s_{I(j)}), \ \forall j \in N^C \cup N^V, s_{C_j} \in S_{C_j} \tag{6}\\$$
+
+$$\mu_{C_j}(s_{C_j}) \le \delta(s_j \mid s_{I(j)}), \ \forall j \in N^D, s_{C_j} \in S_{C_j} \tag{7}\\$$
+
+The value $p(s_j \mid s_{I(j)})$ is the conditional probability of a state $s_j$ given the information state $s_{I(j)}$ and $\delta(s_j \mid s_{I(j)})$ the decision strategy in node $j \in N^D$.
+
+#### Non-negativity and integer constraints
+
+All probability and decision variables have to be non-negative. Decision variables are binary.
+
+$$\mu_{C_j}(s_{C_j}) \ge 0, \ \forall j \in N, s_{C_j} \in S_{C_j}\tag{8}$$
+
+$$\delta(s_j \mid s_{I(j)}) \in \{0,1\}, \ \forall j \in N^D, s_j \in S_j, s_{I(j)} \in S_{I(j)}\tag{9}$$
+
+## Limitations
+
+Currently, the RJT formulation commands in the package do not support forbidden path or fixed path features.
+
+## Computational considerations
+
+Around 2-3 magnitudes faster solving times are expected using RJT formulations. [^1] In problems with small treewidths, such as the pig breeding problem, the solution times hardly even change when increasing the number of nodes.
+
+## References
+
+[^1]: Herrala, O., Terho, T., Oliveira, F., 2024. Risk-averse decision strategies for influence diagrams using rooted junction trees. Retrieved from [https://arxiv.org/abs/2401.03734]
+
+[^2]: Parmentier, A., Cohen, V., Leclere, V., Obozinski, G., Salmon, J., 2020. `
+Integer programming on the junction tree polytope for influence diagrams. INFORMS Journal on Optimization 2, 209–228.
+
+[^3]: Koller, D., Friedman, N., 2009. Probabilistic graphical models: principles
+and techniques. MIT press

docs/src/decision-programming/decision-model.md

@@ -74,6 +74,10 @@ $$\mathcal{U}^-(𝐬) = \mathcal{U}(𝐬) - \max_{𝐬∈𝐒} \mathcal{U}(𝐬)
 ## Conditional Value-at-Risk
 The section [Measuring Risk](@ref) explains and visualizes the relationships between the formulation of expected value, value-at-risk and conditional value-at-risk for discrete probability distribution.
 
+In this section, CVaR models are defined for both path-based and RJT models.
+
+### Path-based model
+
 Given decision strategy $Z,$ we define the cumulative distribution of compatible paths' probabilities as
 
 $$F_Z(t) = ∑_{𝐬∈𝐒∣\mathcal{U}(𝐬)≤t} x(𝐬) p(𝐬).$$
@@ -138,10 +142,44 @@ We can express the conditional value-at-risk objective as
 
 $$\operatorname{CVaR}_α(Z)=\frac{1}{α}∑_{𝐬∈𝐒}\bar{ρ}(𝐬) \mathcal{U}(𝐬)\tag{25}.$$
 
+### RJT model
+
+CVaR formulation for the RJT model is close to that of path-based model. A diagram can have only a single value node, when using RJT-based CVaR. Trying to call the RJT-based CVaR function using a diagram with more than one value node results in an error.
+
+We denote the possible utility values with $u ∈ U$ and suppose we can define the probability $p(u)$ of attaining a given utility value. In the presence of a single value node, we define $p(u) = ∑_{s_{C_v}∈ \text{\{} S_{C_v} \vert U(s_{C_v})=u \text{\}} }µ(s_{C_v})$. We can then pose the constraints
+
+$$η-u≤M λ(u),\quad ∀u∈U \tag{26}$$
+
+$$η-u≥(M+ϵ) λ(u) - M,\quad ∀u∈U \tag{27}$$
+
+$$η-u≤(M+ϵ) \bar{λ}(u) - ϵ,\quad ∀u∈U \tag{28}$$
+
+$$η-u≥M (\bar{λ}(u) - 1),\quad ∀u∈U \tag{29}$$
+
+$$\bar{ρ}(u) ≤ \bar{λ}(u),\quad ∀u∈U \tag{30}$$
+
+$$p(u) - (1 - λ(u)) ≤ ρ(u) ≤ λ(u),\quad ∀u∈U \tag{31}$$
+
+$$ρ(u) ≤ \bar{ρ}(u) ≤ p(u),\quad ∀u∈U \tag{32}$$
+
+$$∑_{u∈U}\bar{ρ}(u) = α \tag{33}$$
+
+$$\bar{λ}(u), λ(u)∈\{0, 1\},\quad ∀u∈U \tag{34}$$
+
+$$\bar{ρ}(u),ρ(u)∈[0, 1],\quad ∀u∈U \tag{35}$$
+
+$$η∈\mathbb{R} \tag{36}$$
+
+where where α is the probability level in VaR<sub>α</sub>.
+
+$CVaR_α$ can be obtained as $1/α ∑_{u∈U} \bar{ρ}(u)u$.
+
+More details, including explanations of variables and constraints, can be found from Herrala et al. (2024)[^4].
+
 ## Convex Combination
 We can combine expected value and conditional value-at-risk using a convex combination at a fixed probability level $α∈(0, 1]$ as follows
 
-$$w \operatorname{E}(Z) + (1-w) \operatorname{CVaR}_α(Z), \tag{26}$$
+$$w \operatorname{E}(Z) + (1-w) \operatorname{CVaR}_α(Z), \tag{37}$$
 
 where the parameter $w∈[0, 1]$ expresses the decision maker's **risk tolerance**.
 
@@ -152,3 +190,5 @@ where the parameter $w∈[0, 1]$ expresses the decision maker's **risk tolerance
 [^2]: Hölsä, O. (2020). Decision Programming Framework for Evaluating Testing Costs of Disease-Prone Pigs. Retrieved from [http://urn.fi/URN:NBN:fi:aalto-202009295618](http://urn.fi/URN:NBN:fi:aalto-202009295618)
 
 [^3]: Hankimaa, H., Herrala, O., Oliveira, F., Tollander de Balsch, J. (2023). DecisionProgramming.jl -- A framework for modelling decision problems using mathematical programming. Retrieved from [https://arxiv.org/abs/2307.13299](https://arxiv.org/abs/2307.13299)
+
+[^4]: Herrala, O., Terho, T., Oliveira, F., 2024. Risk-averse decision strategies for influence diagrams using rooted junction trees. Retrieved from [https://arxiv.org/abs/2401.03734]

docs/src/decision-programming/figures/pig_breeding_N=4.jpg

@@ -189,8 +189,19 @@ The expected utility is used as the objective and the problem is solved using Gu
 ```julia
 EV = expected_value(model, diagram, x_s)
 @objective(model, Max, EV)
+```
 
+Alternatively, RJT formulation can be used by replacing commands on path compatibility variables and objective function creation with commands
 
+```julia
+μ_s = RJTVariables(model, diagram, z)
+EV = expected_value(model, diagram, μ_s)
+@objective(model, Max, EV)
+```
+
+and then solving using the solver. Significantly faster solving times are expected using RJT formulation.
+
+```julia
 optimizer = optimizer_with_attributes(
     () -> HiGHS.Optimizer()
 )

docs/src/decision-programming/figures/pig_breeding_rjt_N=4.jpg

@@ -185,6 +185,16 @@ EV = expected_value(model, diagram, x_s)
 
 and set up the solver.
 
+Alternatively, RJT formulation can be used by replacing commands on path compatibility variables and objective function creation with commands
+
+```julia
+μ_s = RJTVariables(model, diagram, z)
+EV = expected_value(model, diagram, μ_s)
+@objective(model, Max, EV)
+```
+
+and then solving using the solver. Significantly faster solving times are expected using RJT formulation.
+
 ```julia
 optimizer = optimizer_with_attributes(
     () -> HiGHS.Optimizer()
@@ -198,6 +208,15 @@ spu = singlePolicyUpdate(diagram, model, z, x_s)
 optimize!(model)
 ```
 
+<!-- Onko tämä hyvä, voisi tehdä kunnon esimerkin myös CVaRista, mutta onko nyt tarpeen? -->
+
+CVaR model can be created by adding the following constraint to the model. The model has to be built so that there is only one value node. The constraint with this specific numerical value here is tested and meaningful for N = 6.
+
+```
+α = 0.05
+CVaR = conditional_value_at_risk(model, diagram, μ_s, α; probability_scale_factor = 1.0)
+@constraint(model, CVaR>=300.0)
+```
 
 ## Analyzing results
 

docs/src/examples/n-monitoring.md

@@ -139,6 +139,16 @@ EV = expected_value(model, diagram, x_s)
 @objective(model, Max, EV)
 ```
 
+Alternatively, RJT formulation can be used by replacing commands on path compatibility variables and objective function creation with commands
+
+```julia
+μ_s = RJTVariables(model, diagram, z)
+EV = expected_value(model, diagram, μ_s)
+@objective(model, Max, EV)
+```
+
+and then solving using the solver. Significantly faster solving times are expected using RJT formulation.
+
 We can perform the optimization using an optimizer such as HiGHS.
 
 ```julia

docs/src/examples/pig-breeding.md

@@ -171,21 +171,12 @@ Y_HB["no CHD", "treatment"] = 7.64528451705134
 Y_HB["no CHD", "no treatment"] = 7.70088349200034
 add_utilities!(diagram, "HB", Y_HB)
 
-generate_diagram!(diagram)
-
-
-@info("Creating the decision model.")
-model = Model()
-z = DecisionVariables(model, diagram)
-
 # Defining forbidden paths to include all those where a test is repeated twice
 forbidden_tests = ForbiddenPath(diagram, ["T1","T2"], [("TRS", "TRS"),("GRS", "GRS"),("no test", "TRS"), ("no test", "GRS")])
 fixed_R0 = FixedPath(diagram, Dict("R0" => chosen_risk_level))
-scale_factor = 10000.0
-x_s = PathCompatibilityVariables(model, diagram, z; fixed = fixed_R0, forbidden_paths = [forbidden_tests], probability_cut=false)
 
-EV = expected_value(model, diagram, x_s)
-@objective(model, Max, EV)
+@info("Creating the decision model.")
+model, z, x_s = generate_model(diagram, model_type="DP", forbidden_paths=[forbidden_tests], fixed=fixed_R0)
 
 @info("Starting the optimization process.")
 optimizer = optimizer_with_attributes(