diff --git a/README.md b/README.md index f938ca4f..c084ad92 100644 --- a/README.md +++ b/README.md @@ -14,16 +14,22 @@ We can create an influence diagram as follows: ```julia using DecisionProgramming -S = States([2, 2, 2, 2]) -C = [ChanceNode(2, [1]), ChanceNode(3, [1])] -D = [DecisionNode(1, Node[]), DecisionNode(4, [2, 3])] -V = [ValueNode(5, [4])] -X = [Probabilities(2, [0.4 0.6; 0.6 0.4]), Probabilities(3, [0.7 0.3; 0.3 0.7])] -Y = [Consequences(5, [1.5, 1.7])] -validate_influence_diagram(S, C, D, V) -sort!.((C, D, V, X, Y), by = x -> x.j) -P = DefaultPathProbability(C, X) -U = DefaultPathUtility(V, Y) + +diagram = InfluenceDiagram() + +add_node!(diagram, DecisionNode("A", [], ["a", "b"])) +add_node!(diagram, DecisionNode("D", ["B", "C"], ["k", "l"])) +add_node!(diagram, ChanceNode("B", ["A"], ["x", "y"])) +add_node!(diagram, ChanceNode("C", ["A"], ["v", "w"])) +add_node!(diagram, ValueNode("V", ["D"])) + +generate_arcs!(diagram) + +add_probabilities!(diagram, "B", [0.4 0.6; 0.6 0.4]) +add_probabilities!(diagram, "C", [0.7 0.3; 0.3 0.7]) +add_utilities!(diagram, "V", [1.5, 1.7]) + +generate_diagram!(diagram) ``` Using the influence diagram, we create the decision model as follow: @@ -31,13 +37,13 @@ Using the influence diagram, we create the decision model as follow: ```julia using JuMP model = Model() -z = DecisionVariables(model, S, D) -x_s = PathCompatibilityVariables(model, z, S, P) -EV = expected_value(model, x_s, U, P) +z = DecisionVariables(model, diagram) +x_s = PathCompatibilityVariables(model, diagram, z) +EV = expected_value(model, diagram, x_s) @objective(model, Max, EV) ``` -Finally, we can optimize the model using MILP solver. +We can optimize the model using MILP solver. ```julia using Gurobi diff --git a/docs/src/api.md b/docs/src/api.md index c04f7d19..002d88cf 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -1,25 +1,27 @@ # API Reference `DecisionProgramming.jl` API reference. + ## `influence_diagram.jl` ### Nodes ```@docs Node +Name +AbstractNode ChanceNode DecisionNode ValueNode State States -States(::Vector{Tuple{State, Vector{Node}}}) -validate_influence_diagram ``` ### Paths ```@docs Path -paths(::AbstractVector{State}) -paths(::AbstractVector{State}, ::Dict{Node, State}) ForbiddenPath +FixedPath +paths(::AbstractVector{State}) +paths(::AbstractVector{State}, FixedPath) ``` ### Probabilities @@ -33,9 +35,10 @@ AbstractPathProbability DefaultPathProbability ``` -### Consequences +### Utilities ```@docs -Consequences +Utility +Utilities ``` ### Path Utility @@ -44,6 +47,18 @@ AbstractPathUtility DefaultPathUtility ``` +### InfluenceDiagram +InfluenceDiagram +generate_arcs! +generate_diagram! +add_node! +ProbabilityMatrix +set_probability! +add_probabilities! +UtilityMatrix +set_utility! +add_utilities! + ### Decision Strategy ```@docs LocalDecisionStrategy @@ -61,10 +76,8 @@ lazy_probability_cut(::Model, ::PathCompatibilityVariables, ::AbstractPathProbab ### Objective Functions ```@docs -PositivePathUtility -NegativePathUtility -expected_value(::Model, ::PathCompatibilityVariables, ::AbstractPathUtility, ::AbstractPathProbability; ::Float64) -conditional_value_at_risk(::Model, ::PathCompatibilityVariables{N}, ::AbstractPathUtility, ::AbstractPathProbability, ::Float64; ::Float64) where N +expected_value(::Model, ::InfluenceDiagram, ::PathCompatibilityVariables; ::Float64) +conditional_value_at_risk(::Model, ::InfluenceDiagram, ::PathCompatibilityVariables{N}, ::Float64; ::Float64) where N ``` ### Decision Strategy from Variables @@ -76,13 +89,14 @@ DecisionStrategy(::DecisionVariables) ## `analysis.jl` ```@docs CompatiblePaths +CompatiblePaths(::InfluenceDiagram, ::DecisionStrategy, ::FixedPath) UtilityDistribution -UtilityDistribution(::States, ::AbstractPathProbability, ::AbstractPathUtility, ::DecisionStrategy) +UtilityDistribution(::InfluenceDiagram, ::DecisionStrategy) StateProbabilities -StateProbabilities(::States, ::AbstractPathProbability, ::DecisionStrategy) -StateProbabilities(::States, ::AbstractPathProbability, ::DecisionStrategy, ::Node, ::State, ::StateProbabilities) -value_at_risk(::Vector{Float64}, ::Vector{Float64}, ::Float64) -conditional_value_at_risk(::Vector{Float64}, ::Vector{Float64}, ::Float64) +StateProbabilities(::InfluenceDiagram, ::DecisionStrategy) +StateProbabilities(::InfluenceDiagram, ::DecisionStrategy, ::Name, ::Name, ::StateProbabilities) +value_at_risk(::UtilityDistribution, ::Float64) +conditional_value_at_risk(::UtilityDistribution, ::Float64) ``` ## `printing.jl` diff --git a/docs/src/examples/CHD_preventative_care.md b/docs/src/examples/CHD_preventative_care.md index e6d0660a..dc6d7509 100644 --- a/docs/src/examples/CHD_preventative_care.md +++ b/docs/src/examples/CHD_preventative_care.md @@ -6,19 +6,19 @@ In this example, we will showcase the subproblem, which optimises the decision strategy given a single prior risk level. The chosen risk level in this example is 12%. The solution to the main problem is found in [^1]. -## Influence Diagram +## Influence diagram ![](figures/CHD_preventative_care.svg) The influence diagram representation of the problem is seen above. The chance nodes $R$ represent the patient's risk estimate – the prior risk estimate being $R0$. The risk estimate nodes $R0$, $R1$ and $R2$ have 101 states $R = \{0\%, 1\%, ..., 100\%\}$, which are the discretised risk levels for the risk estimates. -The risk estimate is updated according to the first and second test decisions, which are represented by decision nodes $T1$ and $T2$. These nodes have states $T = \{\text{TRS, GRS, no test}\}$. The health of the patient represented by chance node $H$ also affects the update of the risk estimate. In this model, the health of the patient indicates whether they will have a CHD event in the next ten years or not. Thus, the node has states $H = \{\text{CHD, no CHD}\}$. The treatment decision is represented by node $TD$ and it has states $TD = \{\text{treatment, no treatment}\}$. +The risk estimate is updated according to the first and second testing decisions, which are represented by decision nodes $T1$ and $T2$. These nodes have states $T = \{\text{TRS, GRS, no test}\}$. The health of the patient, represented by chance node $H$, also affects the update of the risk estimate. In this model, the health of the patient indicates whether they will have a CHD event in the next ten years or not. Thus, the node has states $H = \{\text{CHD, no CHD}\}$. The treatment decision is represented by node $TD$ and it has states $TD = \{\text{treatment, no treatment}\}$. The prior risk estimate represented by node $R0$ influences the health node $H$, because in the model we make the assumption that the prior risk estimate accurately describes the probability of having a CHD event. The value nodes in the model are $TC$ and $HB$. Node $TC$ represents the testing costs incurred due to the testing decisions $T1$ and $T2$. Node $HB$ represents the health benefits achieved. The testing costs and health benefits are measured in QALYs. These parameter values were evaluated in the study [^2]. -We begin by declaring the chosen prior risk level and reading the conditional probability data for the tests. The risk level 12% is referred to as 13 because indexing begins from 1 and the first risk level is 0\%. Note also that the sample data in this repository is dummy data due to distribution restrictions on the real data. We also define functions ```update_risk_distribution ```, ```state_probabilities``` and ```analysing_results ```. These functions will be discussed in the following sections. +We begin by declaring the chosen prior risk level and reading the conditional probability data for the tests. Note that the sample data in this repository is dummy data due to distribution restrictions on the real data. We also define functions ```update_risk_distribution ``` and ```state_probabilities```. These functions will be discussed in the following sections. ```julia using Logging @@ -27,7 +27,7 @@ using DecisionProgramming using CSV, DataFrames, PrettyTables -const chosen_risk_level = 13 +const chosen_risk_level = 12% const data = CSV.read("risk_prediction_data.csv", DataFrame) function update_risk_distribution(prior::Int64, t::Int64)... @@ -35,69 +35,66 @@ end function state_probabilities(risk_p::Array{Float64}, t::Int64, h::Int64, prior::Int64)... end - -function analysing_results(Z::DecisionStrategy, sprobs::StateProbabilities)... -end ``` - -We define the decision programming model by first defining the node indices and states: +### Initialise influence diagram +We start defining the Decision Programming model by initialising the influence diagram. ```julia -const R0 = 1 -const H = 2 -const T1 = 3 -const R1 = 4 -const T2 = 5 -const R2 = 6 -const TD = 7 -const TC = 8 -const HB = 9 +diagram = InfluenceDiagram() +``` + +For brevity in the next sections, we define the states of the nodes to be readily available. Notice, that $R_{states}$ is a vector with values $0\%, 1\%,..., 100\%$. +```julia const H_states = ["CHD", "no CHD"] const T_states = ["TRS", "GRS", "no test"] const TD_states = ["treatment", "no treatment"] -const R_states = map( x -> string(x) * "%", [0:1:100;]) -const TC_states = ["TRS", "GRS", "TRS & GRS", "no tests"] -const HB_states = ["CHD & treatment", "CHD & no treatment", "no CHD & treatment", "no CHD & no treatment"] - -@info("Creating the influence diagram.") -S = States([ - (length(R_states), [R0, R1, R2]), - (length(H_states), [H]), - (length(T_states), [T1, T2]), - (length(TD_states), [TD]) -]) - -C = Vector{ChanceNode}() -D = Vector{DecisionNode}() -V = Vector{ValueNode}() -X = Vector{Probabilities}() -Y = Vector{Consequences}() +const R_states = [string(x) * "%" for x in [0:1:100;]] +``` + + We then add the nodes. The chance and decision nodes are identified by their names. When declaring the nodes, they are also given information sets and states. Notice that nodes $R0$ and $H$ are root nodes, meaning that their information sets are empty. In Decision Programming, we add the chance and decision nodes in the follwoing way. + ```julia +add_node!(diagram, ChanceNode("R0", [], R_states)) +add_node!(diagram, ChanceNode("R1", ["R0", "H", "T1"], R_states)) +add_node!(diagram, ChanceNode("R2", ["R1", "H", "T2"], R_states)) +add_node!(diagram, ChanceNode("H", ["R0"], H_states)) + +add_node!(diagram, DecisionNode("T1", ["R0"], T_states)) +add_node!(diagram, DecisionNode("T2", ["R1"], T_states)) +add_node!(diagram, DecisionNode("TD", ["R2"], TD_states)) ``` +The value nodes are added in a similar fashion. However, value nodes do not have states because they map their information states to utility values instead. -Next, we define the nodes with their information sets and corresponding probabilities for chance nodes and consequences for value nodes. +```julia +add_node!(diagram, ValueNode("TC", ["T1", "T2"])) +add_node!(diagram, ValueNode("HB", ["H", "TD"])) +``` + +### Generate arcs +Now that all of the nodes have been added to the influence diagram we generate the arcs between the nodes. This step automatically orders the nodes, gives them indices and reorganises the information into the appropriate form. +```julia +generate_arcs!(diagram) +``` -### Prior risk estimate and health of the patient +### Probabilities of the prior risk estimate and health of the patient -In this subproblem, the prior risk estimate is given and therefore the node $R0$ is in effect a deterministic node. In decision programming a deterministic node is added as a chance node for which the probability of one state is set to one and the probabilities of the rest of the states are set to zero. In this case +In this subproblem, the prior risk estimate is given and therefore the node $R0$ is in effect a deterministic node. In Decision Programming a deterministic node is added as a chance node for which the probability of one state is set to one and the probabilities of the rest of the states are set to zero. In this case $$ℙ(R0 = 12\%)=1$$ and $$ℙ(R0 \neq 12\%)= 0. $$ -Notice also that node $R0$ is the first node in the influence diagram, meaning that its information set $I(R0)$ is empty. In decision programming we add node $R0$ and its state probabilities as follows: +The probability matrix of node $R0$ is added in the following way. Remember that the `ProbabilityMatrix` function initialises the matrix with zeros. ```julia -I_R0 = Vector{Node}() -X_R0 = zeros(S[R0]) -X_R0[chosen_risk_level] = 1 -push!(C, ChanceNode(R0, I_R0)) -push!(X, Probabilities(R0, X_R0)) +X_R0 = ProbabilityMatrix(diagram, "R0") +set_probability!(X_R0, [chosen_risk_level], 1) +add_probabilities!(diagram, "R0", X_R0) ``` -Next we add node $H$ and its state probabilities. For modeling purposes, we define the information set of node $H$ to include the prior risk node $R0$. We set the probability that the patient experiences a CHD event in the next ten years according to the prior risk level such that +Next we add the state probabilities of node $H$. For modeling purposes, we define the information set of node $H$ to include the prior risk node $R0$. We set the probability that the patient experiences a CHD event in the next ten years according to the prior risk level such that $$ℙ(H = \text{CHD} | R0 = \alpha) = \alpha.$$ @@ -107,140 +104,110 @@ $$ℙ(H = \text{no CHD} | R0 = \alpha) = 1 - \alpha$$ Since node $R0$ is deterministic and the health node $H$ is defined in this way, in our model the patient has a 12% chance of experiencing a CHD event and 88% chance of remaining healthy. -Node $H$ and its probabilities are added in the following way. - +In Decision Programming the probability matrix of node $H$ has dimensions (101, 2) because its information set consisting of node $R0$ has 101 states and node $H$ has 2 states. We first set the column related to the state $CHD$ with values from `data.risk_levels` which are $0.00, 0.01, ..., 0.99, 1.00$ and the other column as its complement event. ```julia -I_H = [R0] -X_H = zeros(S[R0], S[H]) -X_H[:, 1] = data.risk_levels -X_H[:, 2] = 1 .- X_H[:, 1] -push!(C, ChanceNode(H, I_H)) -push!(X, Probabilities(H, X_H)) +X_H = ProbabilityMatrix(diagram, "H") +set_probability!(X_H, [:, "CHD"], data.risk_levels) +set_probability!(X_H, [:, "no CHD"], 1 .- data.risk_levels) +add_probabilities!(diagram, "H", X_H) ``` -### Test decisions and updating the risk estimate - -The node representing the first test decision is added to the model. - -```julia -I_T1 = [R0] -push!(D, DecisionNode(T1, I_T1)) -``` +### Probabilities of the updated the risk estimates -For node $R1%$, the probabilities of the states are calculated by aggregating the updated risk estimates, after a test is performed, into the risk levels. The updated risk estimates are calculated using the function ```update_risk_distribution```, which calculates the posterior probability distribution for a given health state, test and prior risk estimate. +For node $R1%$, the probabilities of the states are calculated by aggregating the updated risk estimates into the risk levels after a test is performed. The updated risk estimates are calculated using the function ```update_risk_distribution```, which calculates the posterior probability distribution for a given health state, test and prior risk estimate. $$\textit{risk estimate} = P(\text{CHD} \mid \text{test result}) = \frac{P(\text{test result} \mid \text{CHD})P(\text{CHD})}{P(\text{test result})}$$ The probabilities $P(\text{test result} \mid \text{CHD})$ are test specific and these are read from the CSV data file. The updated risk estimates are aggregated according to the risk levels. These aggregated probabilities are then the state probabilities of node $R1$. The aggregating is done using function ```state_probabilities```. -The node $R1$ and its probabilities are added in the following way. +In Decision Programming the probability distribution over the states of node $R1$ is defined into a probability matrix with dimensions $(101,2,3,101)$. This is because its information set consists of nodes $R0, H$ and, $T$ which have 101, 2 and 3 states respectively and the node $R1$ itself has 101 states. Here, one must know that in Decision Programming the states of the nodes are mapped to numbers in the back-end. For instance, the health states $\text{CHD}$ and $\text{no CHD}$ are indexed 1 and 2. The testing decision states TRS, GRS and no test are 1, 2 and 3. The order of the states is determined by the order in which they are defined when adding the nodes. Knowing this, we can set the probability values into the probability matrix using a very compact syntax. Notice that we add 101 probability values at a time into the matrix. ```julia -I_R1 = [R0, H, T1] -X_R1 = zeros(S[I_R1]..., S[R1]) +X_R = ProbabilityMatrix(diagram, "R1") for s_R0 = 1:101, s_H = 1:2, s_T1 = 1:3 - X_R1[s_R0, s_H, s_T1, :] = state_probabilities(update_risk_distribution(s_R0, s_T1), s_T1, s_H, s_R0) + X_R[s_R0, s_H, s_T1, :] = state_probabilities(update_risk_distribution(s_R0, s_T1), s_T1, s_H, s_R0) end -push!(C, ChanceNode(R1, I_R1)) -push!(X, Probabilities(R1, X_R1)) +add_probabilities!(diagram, "R1", X_R) ``` -Nodes $T2$ and $R2$ are added in a similar fashion to nodes $T1$ and $R1$ above. +We notice that the probability distrubtion is identical in $R1$ and $R2$ because their information states are identical. Therefore we can simply add the same matrix from above as the probability matrix of node $R2$. ```julia -I_T2 = [R1] -push!(D, DecisionNode(T2, I_T2)) - - -I_R2 = [H, R1, T2] -X_R2 = zeros(S[I_R2]..., S[R2]) -for s_R1 = 1:101, s_H = 1:2, s_T2 = 1:3 - X_R2[s_H, s_R1, s_T2, :] = state_probabilities(update_risk_distribution(s_R1, s_T2), s_T2, s_H, s_R1) -end -push!(C, ChanceNode(R2, I_R2)) -push!(X, Probabilities(R2, X_R2)) +add_probabilities!(diagram, "R2", X_R) ``` -We also add the treatment decision node $TD$. The treatment decision is made based on the risk estimate achieved with the testing process. +### Utilities of testing costs and health benefits -```julia -I_TD = [R2] -push!(D, DecisionNode(TD, I_TD)) -``` - -### Test costs and health benefits - -To add the value node $TC$, which represents testing costs, we need to define the consequences for its different information states. The node and the consequences are added in the following way. +We define a utility matrix for node $TC$, which maps all its information states to testing +costs. The unit in which the testing costs are added is quality-adjusted life-year (QALYs). The utility matrix is defined and added in the following way. ```julia -I_TC = [T1, T2] -Y_TC = zeros(S[I_TC]...) cost_TRS = -0.0034645 cost_GRS = -0.004 -cost_forbidden = 0 #the cost of forbidden test combinations is negligible -Y_TC[1 , 1] = cost_forbidden -Y_TC[1 , 2] = cost_TRS + cost_GRS -Y_TC[1, 3] = cost_TRS -Y_TC[2, 1] = cost_GRS + cost_TRS -Y_TC[2, 2] = cost_forbidden -Y_TC[2, 3] = cost_GRS -Y_TC[3, 1] = cost_TRS -Y_TC[3, 2] = cost_GRS -Y_TC[3, 3] = 0 -push!(V, ValueNode(TC, I_TC)) -push!(Y, Consequences(TC, Y_TC)) +forbidden = 0 #the cost of forbidden test combinations is negligible + +Y_TC = UtilityMatrix(diagram, "TC") +set_utility!(Y_TC, ["TRS", "TRS"], forbidden) +set_utility!(Y_TC, ["TRS", "GRS"], cost_TRS + cost_GRS) +set_utility!(Y_TC, ["TRS", "no test"], cost_TRS) +set_utility!(Y_TC, ["GRS", "TRS"], cost_TRS + cost_GRS) +set_utility!(Y_TC, ["GRS", "GRS"], forbidden) +set_utility!(Y_TC, ["GRS", "no test"], cost_GRS) +set_utility!(Y_TC, ["no test", "TRS"], cost_TRS) +set_utility!(Y_TC, ["no test", "GRS"], cost_GRS) +set_utility!(Y_TC, ["no test", "no test"], 0) +add_utilities!(diagram, "TC", Y_TC) + ``` -The health benefits that are achieved are determined by whether treatment is administered and by the health of the patient. We add the final node to the model. +The health benefits that are achieved are determined by whether treatment is administered and by the health of the patient. We add the final utility matrix to the model. ```julia -I_HB = [H, TD] -Y_HB = zeros(S[I_HB]...) -Y_HB[1 , 1] = 6.89713671259061 # sick & treat -Y_HB[1 , 2] = 6.65436854256236 # sick & don't treat -Y_HB[2, 1] = 7.64528451705134 # healthy & treat -Y_HB[2, 2] = 7.70088349200034 # healthy & don't treat -push!(V, ValueNode(HB, I_HB)) -push!(Y, Consequences(HB, Y_HB)) +Y_HB = UtilityMatrix(diagram, "HB") +set_utility!(Y_HB, ["CHD", "treatment"], 6.89713671259061) +set_utility!(Y_HB, ["CHD", "no treatment"], 6.65436854256236 ) +set_utility!(Y_HB, ["no CHD", "treatment"], 7.64528451705134) +set_utility!(Y_HB, ["no CHD", "no treatment"], 7.70088349200034) +add_utilities!(diagram, "HB", Y_HB) ``` -### Validating the Influence Diagram -Before creating the decision model, we need to validate the influence diagram and sort the nodes, probabilities and consequences in increasing order by the node indices. +### Generate influence diagram +Finally, we generate the full influence diagram before defining the decision model. By default this function uses the default path probabilities and utilities, which are defined as the joint probability of all chance events in the diagram and the sum of utilities in value nodes, respectively. In the [Contingent Portfolio Programming](contingent-portfolio-programming.md) example, we show how to use a user-defined custom path utility function. + + +## Decision Model +We define the JuMP model and declare the decision variables. ```julia -validate_influence_diagram(S, C, D, V) -sort!.((C, D, V, X, Y), by = x -> x.j) +model = Model() +z = DecisionVariables(model, diagram) ``` -We also define the path probability and the path utility. We use the default path utility, which is the sum of the consequences of the path. +In this problem, we want to forbid the model from choosing paths where the same test is repeated twice and where the first testing decision is not to perform a test but the second testing decision is to perform a test. We forbid the paths by declaring these combinations of states as forbidden paths. + ```julia -P = DefaultPathProbability(C, X) -U = DefaultPathUtility(V, Y) +forbidden_tests = ForbiddenPath(diagram, ["T1","T2"], [("TRS", "TRS"),("GRS", "GRS"),("no test", "TRS"), ("no test", "GRS")]) ``` +We fix the state of the deterministic $R0$ node by declaring it as a fixed path. Fixing the state of node $R0$ is not necessary because of how the probabilities were defined. However, the fixed state reduces the need for some computation in the back-end. -## Decision Model -We define our model and declare the decision variables. ```julia -model = Model() -z = DecisionVariables(model, S, D) +fixed_R0 = FixedPath(diagram, Dict("R0" => chosen_risk_level)) ``` -In this problem, we want to forbid the model from choosing paths where the same test is repeated twice and where the first testing decision is not to perform a test but the second testing decision is to perform a test. We forbid the paths by declaring these combinations of states as forbidden paths. - -We also choose a scale factor of 1000, which will be used to scale the path probabilities. The probabilities need to be scaled because in this specific problem they are very small since the $R$ nodes have many states. Scaling the probabilities helps the solver find an optimal solution. +We also choose a scale factor of 10000, which will be used to scale the path probabilities. The probabilities need to be scaled because in this specific problem they are very small since the $R$ nodes have a large number of states. Scaling the probabilities helps the solver find an optimal solution. -We declare the path compatibility variables. We fix the state of the deterministic $R0$ node and forbid the unwanted testing strategies and scale the probabilities by giving them as parameters in the function call. +We then declare the path compatibility variables. We fix the state of the deterministic $R0$ node , forbid the unwanted testing strategies and scale the probabilities by giving them as parameters in the function call. ```julia -forbidden_tests = ForbiddenPath[([T1,T2], Set([(1,1),(2,2),(3,1), (3,2)]))] scale_factor = 10000.0 -x_s = PathCompatibilityVariables(model, z, S, P; fixed = Dict(1 => chosen_risk_level), forbidden_paths = forbidden_tests, probability_cut=false) +x_s = PathCompatibilityVariables(model, diagram, z; fixed = fixed_R0, forbidden_paths = [forbidden_tests], probability_cut=false) + ``` We define the objective function as the expected value. ```julia -EV = expected_value(model, x_s, U, P, probability_scale_factor= scale_factor) +EV = expected_value(model, diagram, x_s, probability_scale_factor = scale_factor) @objective(model, Max, EV) ``` @@ -258,72 +225,85 @@ optimize!(model) -## Analyzing Results - -### Decision Strategy -We obtain the optimal decision strategy from the z variable values. +## Analyzing results +We extract the results in the following way. ```julia Z = DecisionStrategy(z) +S_probabilities = StateProbabilities(diagram, Z) +U_distribution = UtilityDistribution(diagram, Z) + ``` -We use the function ```analysing_results``` to visualise the results in order to inspect the decision strategy. We use this tailor made function merely for convinience. From the printout, we can see that when the prior risk level is 12% the optimal decision strategy is to first perform TRS testing. At the second decision stage, GRS should be conducted if the updated risk estimate is between 16% and 28% and otherwise no further testing should be conducted. Treatment should be provided to those who have a final risk estimate greater than 18%. Notice that the blank spaces in the table are states which have a probability of zero, which means that given this data it is impossible for the patient to have their risk estimate updated to those risk levels. +### Decision strategy +We inspect the decision strategy. From the printout, we can see that when the prior risk level is 12% the optimal decision strategy is to first perform TRS testing. At the second decision stage, GRS should be conducted if the updated risk estimate is between 16% and 28% and otherwise no further testing should be conducted. Treatment should be provided to those who have a final risk estimate greater than 18%. Notice that the incompatible states are not included in the printout. The incompatible states are those that have a state probability of zero, which means that given this data it is impossible for the patient to have their risk estimate updated to those risk levels. ```julia -sprobs = StateProbabilities(S, P, Z) -``` -```julia -julia> println(analysing_results(Z, sprobs)) -┌─────────────────┬────────┬────────┬────────┐ -│ Information_set │ T1 │ T2 │ TD │ -│ String │ String │ String │ String │ -├─────────────────┼────────┼────────┼────────┤ -│ 0% │ │ 3 │ 2 │ -│ 1% │ │ 3 │ 2 │ -│ 2% │ │ │ 2 │ -│ 3% │ │ 3 │ 2 │ -│ 4% │ │ │ │ -│ 5% │ │ │ │ -│ 6% │ │ 3 │ 2 │ -│ 7% │ │ 3 │ 2 │ -│ 8% │ │ │ 2 │ -│ 9% │ │ │ 2 │ -│ 10% │ │ 3 │ 2 │ -│ 11% │ │ 3 │ 2 │ -│ 12% │ 1 │ │ 2 │ -│ 13% │ │ 3 │ 2 │ -│ 14% │ │ 3 │ 2 │ -│ 15% │ │ │ 2 │ -│ 16% │ │ 2 │ 2 │ -│ 17% │ │ 2 │ 2 │ -│ 18% │ │ 2 │ 1 │ -│ 19% │ │ │ 1 │ -│ 20% │ │ │ 1 │ -│ 21% │ │ 2 │ 1 │ -│ 22% │ │ 2 │ 1 │ -│ 23% │ │ 2 │ 1 │ -│ 24% │ │ │ 1 │ -│ 25% │ │ │ 1 │ -│ 26% │ │ │ 1 │ -│ 27% │ │ │ 1 │ -│ 28% │ │ 3 │ 1 │ -│ 29% │ │ 3 │ 1 │ -│ 30% │ │ │ 1 │ -│ ⋮ │ ⋮ │ ⋮ │ ⋮ │ -└─────────────────┴────────┴────────┴────────┘ - 70 rows omitted +julia> print_decision_strategy(diagram, Z, S_probabilities) +┌────────────────┬────────────────┐ +│ State(s) of R0 │ Decision in T1 │ +├────────────────┼────────────────┤ +│ 12% │ TRS │ +└────────────────┴────────────────┘ +┌────────────────┬────────────────┐ +│ State(s) of R1 │ Decision in T2 │ +├────────────────┼────────────────┤ +│ 0% │ no test │ +│ 1% │ no test │ +│ 3% │ no test │ +│ 6% │ no test │ +│ 7% │ no test │ +│ 10% │ no test │ +│ 11% │ no test │ +│ 13% │ no test │ +│ 14% │ no test │ +│ 16% │ GRS │ +│ 17% │ GRS │ +│ 18% │ GRS │ +│ 21% │ GRS │ +│ 22% │ GRS │ +│ 23% │ GRS │ +│ 28% │ no test │ +│ 29% │ no test │ +│ 31% │ no test │ +│ 34% │ no test │ +│ ⋮ │ ⋮ │ +└────────────────┴────────────────┘ + rows omitted + +┌────────────────┬────────────────┐ +│ State(s) of R2 │ Decision in TD │ +├────────────────┼────────────────┤ +│ 0% │ no treatment │ +│ 1% │ no treatment │ +│ 2% │ no treatment │ +│ 3% │ no treatment │ +│ 6% │ no treatment │ +│ 7% │ no treatment │ +│ 8% │ no treatment │ +│ 9% │ no treatment │ +│ 10% │ no treatment │ +│ 11% │ no treatment │ +│ 12% │ no treatment │ +│ 13% │ no treatment │ +│ 14% │ no treatment │ +│ 15% │ no treatment │ +│ 16% │ no treatment │ +│ 17% │ no treatment │ +│ 18% │ treatment │ +│ 19% │ treatment │ +│ 20% │ treatment │ +│ ⋮ │ ⋮ │ +└────────────────┴────────────────┘ + rows omitted ``` -### Utility Distribution +### Utility distribution We can also print the utility distribution for the optimal strategy and some basic statistics for the distribution. ```julia -udist = UtilityDistribution(S, P, U, Z) -``` - -```julia -julia> print_utility_distribution(udist) +julia> print_utility_distribution(U_distribution) ┌──────────┬─────────────┐ │ Utility │ Probability │ │ Float64 │ Float64 │ @@ -339,7 +319,7 @@ julia> print_utility_distribution(udist) └──────────┴─────────────┘ ``` ```julia -julia> print_statistics(udist) +julia> print_statistics(U_distribution) ┌──────────┬────────────┐ │ Name │ Statistics │ │ String │ Float64 │ diff --git a/docs/src/examples/n-monitoring.md b/docs/src/examples/n-monitoring.md index 4ff0145e..86d10bf3 100644 --- a/docs/src/examples/n-monitoring.md +++ b/docs/src/examples/n-monitoring.md @@ -8,7 +8,7 @@ The $N$-monitoring problem is described in [^1], sections 4.1 and 6.1. The influence diagram of the generalized $N$-monitoring problem where $N≥1$ and indices $k=1,...,N.$ The nodes are associated with states as follows. **Load state** $L=\{high, low\}$ denotes the load on a structure, **report states** $R_k=\{high, low\}$ report the load state to the **action states** $A_k=\{yes, no\}$ which represent different decisions to fortify the structure. The **failure state** $F=\{failure, success\}$ represents whether or not the (fortified) structure fails under the load $L$. Finally, the utility at target $T$ depends on the whether $F$ fails and the fortification costs. -We draw the cost of fortification $c_k∼U(0,1)$ from a uniform distribution, and the magnitude of fortification is directly proportional to the cost. Fortification is defined as +We begin by choosing $N$ and defining our fortification cost function. We draw the cost of fortification $c_k∼U(0,1)$ from a uniform distribution, and the magnitude of fortification is directly proportional to the cost. Fortification is defined as $$f(A_k=yes) = c_k$$ @@ -19,53 +19,67 @@ using Logging, Random using JuMP, Gurobi using DecisionProgramming -Random.seed!(13) - const N = 4 -const L = [1] -const R_k = [k + 1 for k in 1:N] -const A_k = [(N + 1) + k for k in 1:N] -const F = [2*N + 2] -const T = [2*N + 3] -const L_states = ["high", "low"] -const R_k_states = ["high", "low"] -const A_k_states = ["yes", "no"] -const F_states = ["failure", "success"] + +Random.seed!(13) const c_k = rand(N) const b = 0.03 fortification(k, a) = [c_k[k], 0][a] +``` + +### Initialising the influence diagram +We initialise the influence diagram before adding nodes to it. -S = States([ - (length(L_states), L), - (length(R_k_states), R_k), - (length(A_k_states), A_k), - (length(F_states), F) -]) -C = Vector{ChanceNode}() -D = Vector{DecisionNode}() -V = Vector{ValueNode}() -X = Vector{Probabilities}() -Y = Vector{Consequences}() +```julia +diagram = InfluenceDiagram ``` -### Load State Probability -The probability that the load is high, $ℙ(L=high)$, is drawn from a uniform distribution. +### Adding nodes +Add node $L$ which represents the load on the structure. This node is the root node and thus, has an empty information set. Its states describe the state of the load, they are $high$ and $low$. -$$ℙ(L=high)∼U(0,1)$$ +```julia +add_node!(diagram, ChanceNode("L", [], ["high", "low"])) +``` + +The report nodes $R_k$ and action nodes $A_k$ are easily added with a for-loop. The report nodes have node $L$ in their information sets and their states are $high$ and $low$. The actions are made based on these reports, which is represented by the action nodes $A_k$ having the report nodes $R_k$ in their information sets. The action nodes have states $yes$ and $no$, which represents decisions whether to fortify the structure or not. ```julia -for j in L - I_j = Vector{Node}() - X_j = zeros(S[I_j]..., S[j]) - X_j[1] = rand() - X_j[2] = 1.0 - X_j[1] - push!(C, ChanceNode(j, I_j)) - push!(X, Probabilities(j, X_j)) +for i in 1:N + add_node!(diagram, ChanceNode("R$i", ["L"], ["high", "low"])) + add_node!(diagram, DecisionNode("A$i", ["R$i"], ["yes", "no"])) end ``` -### Reporting Probability -The probabilities of the report states correspond to the load state. We draw the values $x∼U(0,1)$ and $y∼U(0,1)$ from uniform distribution. +The failure node $F$ has the load node $L$ and all of the action nodes $A_k$ in its information set. The failure node has states $failure$ and $success$. +```julia +add_node!(diagram, ChanceNode("F", ["L", ["A$i" for i in 1:N]...], ["failure", "success"])) +``` + +The value node $T$ is added as follows. +```julia +add_node!(diagram, ValueNode("T", ["F", ["A$i" for i in 1:N]...])) +``` + +### Generating arcs +Now that all of the nodes have been added to the influence diagram we generate the arcs between the nodes. This step automatically orders the nodes, gives them indices and reorganises the information into the appropriate form. +```julia +generate_arcs!(diagram) +``` + + +### Load State Probabilities +After generating the arcs, the probabilities and utilities can be added. The probability that the load is high, $ℙ(L=high)$, is drawn from a uniform distribution. For different syntax options for adding probabilities and utilities, see the [usage page](../usage.md). + +$$ℙ(L=high)∼U(0,1)$$ + +```julia +X_L = [rand(), 0] +X_L[2] = 1.0 - X_L[1] +add_probabilities!(diagram, "L", X_L) +``` + +### Reporting Probabilities +The probabilities of the report states correspond to the load state. We draw the values $x∼U(0,1)$ and $y∼U(0,1)$ from uniform distributions. $$ℙ(R_k=high∣L=high)=\max\{x,1-x\}$$ @@ -73,124 +87,110 @@ $$ℙ(R_k=low∣L=low)=\max\{y,1-y\}$$ The probability of a correct report is thus in the range [0.5,1]. (This reflects the fact that a probability under 50% would not even make sense, since we would notice that if the test suggests a high load, the load is more likely to be low, resulting in that a low report "turns into" a high report and vice versa.) +In Decision Programming we add these probabilities by declaring probabilty matrices for nodes $R_k$. The probability matrix of a report node $R_k$ has dimensions (2,2), where the rows correspond to the states $high$ and $low$ of its predecessor node $L$ and the columns to its own states $high$ and $low$. + ```julia -for j in R_k - I_j = L +for i in 1:N x, y = rand(2) - X_j = zeros(S[I_j]..., S[j]) - X_j[1, 1] = max(x, 1-x) - X_j[1, 2] = 1.0 - X_j[1, 1] - X_j[2, 2] = max(y, 1-y) - X_j[2, 1] = 1.0 - X_j[2, 2] - push!(C, ChanceNode(j, I_j)) - push!(X, Probabilities(j, X_j)) + X_R = ProbabilityMatrix(diagram, "R$i") + set_probability!(X_R, ["high", "high"], max(x, 1-x)) + set_probability!(X_R, ["high", "low"], 1-max(x, 1-x)) + set_probability!(X_R, ["low", "low"], max(y, 1-y)) + set_probability!(X_R, ["low", "high"], 1-max(y, 1-y)) + add_probabilities!(diagram, "R$i", X_R) end ``` -### Decision to Fortify +### Probability of Failure +The probability of failure is decresead by fortification actions. We draw the values $x∼U(0,1)$ and $y∼U(0,1)$ from uniform distribution. -Only the corresponding load report is known when making the fortification decision, thus $I(A_k)=R_k$. +$$ℙ(F=failure∣A_N,...,A_1,L=high)=\frac{\max{\{x, 1-x\}}}{\exp{(b ∑_{k=1,...,N} f(A_k))}}$$ +$$ℙ(F=failure∣A_N,...,A_1,L=low)=\frac{\min{\{y, 1-y\}}}{\exp{(b ∑_{k=1,...,N} f(A_k))}}$$ + +First we initialise the probability matrix for node $F$. ```julia -for (i, j) in zip(R_k, A_k) - I_j = [i] - push!(D, DecisionNode(j, I_j)) -end +X_F = ProbabilityMatrix(diagram, "F") ``` -### Probability of Failure -The probabilities of failure which are decresead by fortifications. We draw the values $x∼U(0,1)$ and $y∼U(0,1)$ from uniform distribution. +This matrix has dimensions $(2, \textcolor{orange}{2, 2, 2, 2}, 2)$ because node $L$ and nodes $A_k$, which form the information set of $F$, all have 2 states and node $F$ itself also has 2 states. The orange colored dimensions correspond to the states of the action nodes $A_k$. -$$ℙ(F=failure∣A_N,...,A_1,L=high)=\frac{\max{\{x, 1-x\}}}{\exp{(b ∑_{k=1,...,N} f(A_k))}}$$ +To set the probabilities we have to iterate over the information states. Here it helps to know that in Decision Programming the states of each node are mapped to numbers in the back-end. For instance, the load states $high$ and $low$ are referred to as 1 and 2. The same applies for the action states $yes$ and $no$, they are states 1 and 2. The `paths` function allows us to iterate over the subpaths of specific nodes. In these paths, the states are refer to by their indices. Using this information, we can easily iterate over the information states using the `paths` function and set the probability values into the probability matrix. -$$ℙ(F=failure∣A_N,...,A_1,L=low)=\frac{\min{\{y, 1-y\}}}{\exp{(b ∑_{k=1,...,N} f(A_k))}}$$ ```julia -for j in F - I_j = L ∪ A_k - x, y = rand(2) - X_j = zeros(S[I_j]..., S[j]) - for s in paths(S[A_k]) - d = exp(b * sum(fortification(k, a) for (k, a) in enumerate(s))) - X_j[1, s..., 1] = max(x, 1-x) / d - X_j[1, s..., 2] = 1.0 - X_j[1, s..., 1] - X_j[2, s..., 1] = min(y, 1-y) / d - X_j[2, s..., 2] = 1.0 - X_j[2, s..., 1] - end - push!(C, ChanceNode(j, I_j)) - push!(X, Probabilities(j, X_j)) +x_F, y_F = rand(2) +for s in paths([State(2) for i in 1:N]) + denominator = exp(b * sum(fortification(k, a) for (k, a) in enumerate(s))) + X_F[1, s..., 1] = max(x_F, 1-x_F) / denominator + X_F[1, s..., 2] = 1.0 - X_F[1, s..., 1] + X_F[2, s..., 1] = min(y_F, 1-y_F) / denominator + X_F[2, s..., 2] = 1.0 - X_F[2, s..., 1] end ``` -### Consequences -Utility from consequences at target $T$ from failure state $F$ +After declaring the probability matrix, we add it to the influence diagram. +```julia +add_probabilities!(diagram, "F", X_F) +``` + +### Utility +The utility from the different scenarios of the failure state at target $T$ are $$g(F=failure) = 0$$ -$$g(F=success) = 100$$ +$$g(F=success) = 100$$. -Utility from consequences at target $T$ from action states $A_k$ is +Utilities from the action states $A_k$ at target $T$ are $$f(A_k=yes) = c_k$$ -$$f(A_k=no) = 0$$ +$$f(A_k=no) = 0.$$ -Total cost +The total cost is thus -$$Y(F, A_N, ..., A_1) = g(F) + (-f(A_N)) + ... + (-f(A_1))$$ +$$Y(F, A_N, ..., A_1) = g(F) + (-f(A_N)) + ... + (-f(A_1)).$$ +We first declare the utility matrix for node $T$. ```julia -for j in T - I_j = A_k ∪ F - Y_j = zeros(S[I_j]...) - for s in paths(S[A_k]) - cost = sum(-fortification(k, a) for (k, a) in enumerate(s)) - Y_j[s..., 1] = cost + 0 - Y_j[s..., 2] = cost + 100 - end - push!(V, ValueNode(j, I_j)) - push!(Y, Consequences(j, Y_j)) -end +Y_T = UtilityMatrix(diagram, "T") ``` - -### Validating the Influence Diagram - -Finally, we need to validate the influence diagram and sort the nodes, probabilities and consequences in increasing order by the node indices. +This matrix has dimensions $(2, \textcolor{orange}{2, 2, 2, 2})$, where the dimensions correspond to the numbers of states the nodes in the information set have. Similarly as before, the first dimension corresponds to the states of node $F$ and the other 4 dimensions (in oragne) correspond to the states of the $A_k$ nodes. The utilities are set and added similarly to how the probabilities were added above. ```julia -validate_influence_diagram(S, C, D, V) -sort!.((C, D, V, X, Y), by = x -> x.j) +for s in paths([State(2) for i in 1:N]) + cost = sum(-fortification(k, a) for (k, a) in enumerate(s)) + Y_T[1, s...] = 0 + cost + Y_T[2, s...] = 100 + cost +end +add_utilities!(diagram, "T", Y_T) ``` -We define the path probability. -```julia -P = DefaultPathProbability(C, X) -``` +### Generate Influence Diagram +The full influence diagram can now be generated. We use the default path probabilities and utilities, which are the default setting in this function. In the [Contingent Portfolio Programming](contingent-portfolio-programming.md) example, we show how to use a user-defined custom path utility function. + +In this particular problem, some of the path utilities are negative. In this case, we choose to use the [positive path utility](../decision-programming/decision-model.md) transformation, which translates the path utilities to positive values. This allows us to exclude the probability cut in the next section. -As the path utility, we use the default, which is the sum of the consequences given the path. ```julia -U = DefaultPathUtility(V, Y) +generate_diagram!(diagram, positive_path_utility = true) ``` - ## Decision Model - -An affine transformation is applied to the path utility, making all utilities positive. See the [section](../decision-programming/decision-model.md) on positive path utilities for details. +We initialise the JuMP model and declare the decision and path compatibility variables. Since we applied an affine transformation to the utility function, the probability cut can be excluded from the model formulation. ```julia -U⁺ = PositivePathUtility(S, U) model = Model() -z = DecisionVariables(model, S, D) -x_s = PathCompatibilityVariables(model, z, S, P, probability_cut = false) +z = DecisionVariables(model, diagram) +x_s = PathCompatibilityVariables(model, diagram, z, probability_cut = false) ``` -` The expected utility is used as the objective and the problem is solved using Gurobi. ```julia -EV = expected_value(model, π_s, U⁺) +EV = expected_value(model, diagram, x_s) @objective(model, Max, EV) + optimizer = optimizer_with_attributes( () -> Gurobi.Optimizer(Gurobi.Env()), "IntFeasTol" => 1e-9, @@ -202,47 +202,46 @@ optimize!(model) ## Analyzing Results -The decision strategy shows us that the optimal strategy is to make all four fortifications regardless of the reports (state 1 in fortification nodes corresponds to the option "yes"). +We obtain the decision strategy, state probabilities and utility distribution from the solution. ```julia Z = DecisionStrategy(z) +U_distribution = UtilityDistribution(diagram, Z) +S_probabilities = StateProbabilities(diagram, Z) ``` +The decision strategy shows us that the optimal strategy is to make all four fortifications regardless of the reports. + ```julia-repl -julia> print_decision_strategy(S, Z) -┌────────┬──────┬───┐ -│ Nodes │ (2,) │ 6 │ -├────────┼──────┼───┤ -│ States │ (1,) │ 1 │ -│ States │ (2,) │ 1 │ -└────────┴──────┴───┘ -┌────────┬──────┬───┐ -│ Nodes │ (3,) │ 7 │ -├────────┼──────┼───┤ -│ States │ (1,) │ 1 │ -│ States │ (2,) │ 1 │ -└────────┴──────┴───┘ -┌────────┬──────┬───┐ -│ Nodes │ (4,) │ 8 │ -├────────┼──────┼───┤ -│ States │ (1,) │ 1 │ -│ States │ (2,) │ 1 │ -└────────┴──────┴───┘ -┌────────┬──────┬───┐ -│ Nodes │ (5,) │ 9 │ -├────────┼──────┼───┤ -│ States │ (1,) │ 1 │ -│ States │ (2,) │ 1 │ -└────────┴──────┴───┘ +julia> print_decision_strategy(diagram, Z, S_probabilities) +┌────────────────┬────────────────┐ +│ State(s) of R1 │ Decision in A1 │ +├────────────────┼────────────────┤ +│ high │ yes │ +│ low │ yes │ +└────────────────┴────────────────┘ +┌────────────────┬────────────────┐ +│ State(s) of R2 │ Decision in A2 │ +├────────────────┼────────────────┤ +│ high │ yes │ +│ low │ yes │ +└────────────────┴────────────────┘ +┌────────────────┬────────────────┐ +│ State(s) of R3 │ Decision in A3 │ +├────────────────┼────────────────┤ +│ high │ yes │ +│ low │ yes │ +└────────────────┴────────────────┘ +┌────────────────┬────────────────┐ +│ State(s) of R4 │ Decision in A4 │ +├────────────────┼────────────────┤ +│ high │ yes │ +│ low │ yes │ +└────────────────┴────────────────┘ ``` -The state probabilities for the strategy $Z$ can also be obtained. These tell the probability of each state in each node, given the strategy $Z$. - - -```julia -sprobs = StateProbabilities(S, P, Z) -``` +The state probabilities for strategy $Z$ are also obtained. These tell the probability of each state in each node, given strategy $Z$. ```julia-repl julia> print_state_probabilities(sprobs, L) @@ -283,12 +282,8 @@ julia> print_state_probabilities(sprobs, F) We can also print the utility distribution for the optimal strategy and some basic statistics for the distribution. -```julia -udist = UtilityDistribution(S, P, U, Z) -``` - ```julia-repl -julia> print_utility_distribution(udist) +julia> print_utility_distribution(U_distribution) ┌───────────┬─────────────┐ │ Utility │ Probability │ │ Float64 │ Float64 │ @@ -299,7 +294,7 @@ julia> print_utility_distribution(udist) ``` ```julia-repl -julia> print_statistics(udist) +julia> print_statistics(U_distribution) ┌──────────┬────────────┐ │ Name │ Statistics │ │ String │ Float64 │ diff --git a/docs/src/examples/pig-breeding.md b/docs/src/examples/pig-breeding.md index 446825ee..29f8282c 100644 --- a/docs/src/examples/pig-breeding.md +++ b/docs/src/examples/pig-breeding.md @@ -2,56 +2,75 @@ ## Description The pig breeding problem as described in [^1]. -"A pig breeder is growing pigs for a period of four months and subsequently selling them. During this period the pig may or may not develop a certain disease. If the pig has the disease at the time it must be sold, the pig must be sold for slaughtering, and its expected market price is then 300 DKK (Danish kroner). If it is disease free, its expected market price as a breeding animal is 1000 DKK +>A pig breeder is growing pigs for a period of four months and subsequently selling them. During this period the pig may or may not develop a certain disease. If the pig has the disease at the time it must be sold, the pig must be sold for slaughtering, and its expected market price is then 300 DKK (Danish kroner). If it is disease free, its expected market price as a breeding animal is 1000 DKK +> +>Once a month, a veterinary doctor sees the pig and makes a test for presence of the disease. If the pig is ill, the test will indicate this with probability 0.80, and if the pig is healthy, the test will indicate this with probability 0.90. At each monthly visit, the doctor may or may not treat the pig for the disease by injecting a certain drug. The cost of an injection is 100 DKK. +> +>A pig has the disease in the first month with probability 0.10. A healthy pig develops the disease in the subsequent month with probability 0.20 without injection, whereas a healthy and treated pig develops the disease with probability 0.10, so the injection has some preventive effect. An untreated pig that is unhealthy will remain so in the subsequent month with probability 0.90, whereas the similar probability is 0.50 for an unhealthy pig that is treated. Thus spontaneous cure is possible, but treatment is beneficial on average. -Once a month, a veterinary doctor sees the pig and makes a test for presence of the disease. If the pig is ill, the test will indicate this with probability 0.80, and if the pig is healthy, the test will indicate this with probability 0.90. At each monthly visit, the doctor may or may not treat the pig for the disease by injecting a certain drug. The cost of an injection is 100 DKK. -A pig has the disease in the first month with probability 0.10. A healthy pig develops the disease in the subsequent month with probability 0.20 without injection, whereas a healthy and treated pig develops the disease with probability 0.10, so the injection has some preventive effect. An untreated pig that is unhealthy will remain so in the subsequent month with probability 0.90, whereas the similar probability is 0.50 for an unhealthy pig that is treated. Thus spontaneous cure is possible, but treatment is beneficial on average." - - -## Influence Diagram +## Influence diagram ![](figures/n-month-pig-breeding.svg) -The influence diagram for the the generalized $N$-month pig breeding. The nodes are associated with the following states. **Health states** $h_k=\{ill,healthy\}$ represent the health of the pig at month $k=1,...,N$. **Test states** $t_k=\{positive,negative\}$ represent the result from testing the pig at month $k=1,...,N-1$. **Treat states** $d_k=\{treat, pass\}$ represent the decision to treat the pig with an injection at month $k=1,...,N-1$. +The influence diagram for the generalized $N$-month pig breeding problem. The nodes are associated with the following states. **Health states** $h_k=\{ill,healthy\}$ represent the health of the pig at month $k=1,...,N$. **Test states** $t_k=\{positive,negative\}$ represent the result from testing the pig at month $k=1,...,N-1$. **Treatment states** $d_k=\{treat, pass\}$ represent the decision to treat the pig with an injection at month $k=1,...,N-1$. > The dashed arcs represent the no-forgetting principle and we can toggle them on and off in the formulation. -In decision programming, we start by defining the node indices and states, as follows: +In this example, we solve the 4 month pig breeding problem and thus, declare $N = 4$. ```julia using JuMP, Gurobi using DecisionProgramming const N = 4 -const health = [3*k - 2 for k in 1:N] -const test = [3*k - 1 for k in 1:(N-1)] -const treat = [3*k for k in 1:(N-1)] -const cost = [(3*N - 2) + k for k in 1:(N-1)] -const price = [(3*N - 2) + N] -const health_states = ["ill", "healthy"] -const test_states = ["positive", "negative"] -const treat_states = ["treat", "pass"] - -S = States([ - (length(health_states), health), - (length(test_states), test), - (length(treat_states), treat), -]) -C = Vector{ChanceNode}() -D = Vector{DecisionNode}() -V = Vector{ValueNode}() -X = Vector{Probabilities}() -Y = Vector{Consequences}() +``` +In Decision Programming, we start by initialising an empty influence diagram. Then we define the nodes with their information sets and states and add them to the influence diagram. +```julia +diagram = InfluenceDiagram() ``` -Next, we define the nodes with their information sets and corresponding probabilities or consequences. -### Health at First Month +### Health at first month -As seen in the influence diagram, the node $h_1$ has no arcs into it, making it a root node. Therefore, the information set $I(h_1)$ is empty. +As seen in the influence diagram, the node $h_1$ has no arcs into it making it a root node. Therefore, the information set $I(h_1)$ is empty. The states of this node are $ill$ and $healthy$. -The probability that pig is ill in the first month is + +```julia +add_node!(diagram, ChanceNode("H1", [], ["ill", "healthy"])) +``` + +### Health, test results and treatment decisions at subsequent months +The chance and decision nodes representing the health, test results, treatment decisions for the following months can be added easily using a for-loop. The value node representing the treatment costs in each month is also added. Each node is given a name, its information set and states. Remember that value nodes do not have states. Notice that we do not assume the no-forgetting principle and thus, the information sets of the treatment decisions only contain the previous test result. + +```julia +for i in 1:N-1 + # Testing result + add_node!(diagram, ChanceNode("T$i", ["H$i"], ["positive", "negative"])) + # Decision to treat + add_node!(diagram, DecisionNode("D$i", ["T$i"], ["treat", "pass"])) + # Cost of treatment + add_node!(diagram, ValueNode("C$i", ["D$i"])) + # Health of next period + add_node!(diagram, ChanceNode("H$(i+1)", ["H$(i)", "D$(i)"], ["ill", "healthy"])) +end +``` + +### Market price +The final value node, representing the market price, is added. It has the final health node $h_N$ as its information set. +```julia +add_node!(diagram, ValueNode("MP", ["H$N"])) +``` + +### Generate arcs +Now that all of the nodes have been added to the influence diagram, we generate the arcs between the nodes. This step automatically orders the nodes, gives them indices and reorganises the information into the correct form. +```julia +generate_arcs!(diagram) +``` + +### Probabilities + +We define probability distributions for all chance nodes. For the first health node, the probability distribution is defined over its two states $ill$ and $healthy$. The probability that pig is ill in the first month is $$ℙ(h_1 = ill)=0.1.$$ @@ -59,159 +78,106 @@ We obtain the complement probabilities for binary states by subtracting from one $$ℙ(h_1 = healthy)=1-ℙ(h_1 = ill).$$ -In decision programming, we add the nodes and probabilities as follows: - +In Decision Programming, we add these probabilities for node $h_1$ as follows. Notice, that the probability vector is ordered according to the order that the states were given in when defining node $h_1$. More information on the syntax of adding probabilities is found on the [usage page](../usage.md). ```julia -for j in health[[1]] - I_j = Vector{Node}() - X_j = zeros(S[I_j]..., S[j]) - X_j[1] = 0.1 - X_j[2] = 1.0 - X_j[1] - push!(C, ChanceNode(j, I_j)) - push!(X, Probabilities(j, X_j)) -end +add_probabilities!(diagram, "H1", [0.1, 0.9]) ``` -### Health at Subsequent Months -The probability that the pig is ill in the subsequent months $k=2,...,N$ depends on the treatment decision and state of health in the previous month $k-1$. The nodes $h_{k-1}$ and $d_{k-1}$ are thus in the information set $I(h_k)$, meaning that the probability distribution of $h_k$ is conditional on these nodes: - -$$ℙ(h_k = ill ∣ d_{k-1} = pass, h_{k-1} = healthy)=0.2,$$ +The probability distributions for the other health nodes are identical. Thus, we define one probability matrix and use it for all the subsequent months' health nodes. The probability that the pig is ill in the subsequent months $k=2,...,N$ depends on the treatment decision and state of health in the previous month $k-1$. The nodes $h_{k-1}$ and $d_{k-1}$ are thus in the information set $I(h_k)$, meaning that the probability distribution of $h_k$ is conditional on these nodes: -$$ℙ(h_k = ill ∣ d_{k-1} = treat, h_{k-1} = healthy)=0.1,$$ +$$ℙ(h_k = ill ∣ h_{k-1} = healthy, \ d_{k-1} = pass)=0.2,$$ -$$ℙ(h_k = ill ∣ d_{k-1} = pass, h_{k-1} = ill)=0.9,$$ +$$ℙ(h_k = ill ∣ h_{k-1} = healthy, \ d_{k-1} = treat)=0.1,$$ -$$ℙ(h_k = ill ∣ d_{k-1} = treat, h_{k-1} = ill)=0.5.$$ +$$ℙ(h_k = ill ∣ h_{k-1} = ill, \ d_{k-1} = pass)=0.9,$$ -In decision programming: +$$ℙ(h_k = ill ∣ h_{k-1} = ill, \ d_{k-1} = treat)=0.5.$$ +In Decision Programming, the probability matrix is define in the following way. Notice, that the ordering of the information state corresponds to the order in which the information set was defined when adding the health nodes. ```julia -for (i, k, j) in zip(health[1:end-1], treat, health[2:end]) - I_j = [i, k] - X_j = zeros(S[I_j]..., S[j]) - X_j[2, 2, 1] = 0.2 - X_j[2, 2, 2] = 1.0 - X_j[2, 2, 1] - X_j[2, 1, 1] = 0.1 - X_j[2, 1, 2] = 1.0 - X_j[2, 1, 1] - X_j[1, 2, 1] = 0.9 - X_j[1, 2, 2] = 1.0 - X_j[1, 2, 1] - X_j[1, 1, 1] = 0.5 - X_j[1, 1, 2] = 1.0 - X_j[1, 1, 1] - push!(C, ChanceNode(j, I_j)) - push!(X, Probabilities(j, X_j)) -end +X_H = ProbabilityMatrix(diagram, "H2") +set_probability!(X_H, ["healthy", "pass", :], [0.2, 0.8]) +set_probability!(X_H, ["healthy", "treat", :], [0.1, 0.9]) +set_probability!(X_H, ["ill", "pass", :], [0.9, 0.1]) +set_probability!(X_H, ["ill", "treat", :], [0.5, 0.5]) ``` -Note that the order of states indexing the probabilities is reversed compared to the mathematical definition. - -### Health Test -For the probabilities that the test indicates a pig's health correctly at month $k=1,...,N-1$, we have +Next we define the probability matrix for the test results. Here again, we note that the probability distributions for all test results are identical, and thus we only define the probability matrix once. For the probabilities that the test indicates a pig's health correctly at month $k=1,...,N-1$, we have $$ℙ(t_k = positive ∣ h_k = ill) = 0.8,$$ $$ℙ(t_k = negative ∣ h_k = healthy) = 0.9.$$ -In decision programming: +In Decision Programming: ```julia -for (i, j) in zip(health, test) - I_j = [i] - X_j = zeros(S[I_j]..., S[j]) - X_j[1, 1] = 0.8 - X_j[1, 2] = 1.0 - X_j[1, 1] - X_j[2, 2] = 0.9 - X_j[2, 1] = 1.0 - X_j[2, 2] - push!(C, ChanceNode(j, I_j)) - push!(X, Probabilities(j, X_j)) -end +X_T = ProbabilityMatrix(diagram, "T1") +set_probability!(X_T, ["ill", "positive"], 0.8) +set_probability!(X_T, ["ill", "negative"], 0.2) +set_probability!(X_T, ["healthy", "negative"], 0.9) +set_probability!(X_T, ["healthy", "positive"], 0.1) ``` -### Decision to Treat -In decision programing, we add the decision nodes for decision to treat the pig as follows: +We add the probability matrices into the influence diagram using a for-loop. ```julia -for (i, j) in zip(test, treat) - I_j = [i] - push!(D, DecisionNode(j, I_j)) +for i in 1:N-1 + add_probabilities!(diagram, "T$i", X_T) + add_probabilities!(diagram, "H$(i+1)", X_H) end ``` -The no-forgetting assumption does not hold, and the information set $I(d_k)$ only comprises the previous test result. -### Cost of Treatment -The cost of treatment decision for the pig at month $k=1,...,N-1$ is defined +### Utilities + +The cost incurred by the treatment decision at month $k=1,...,N-1$ is $$Y(d_k=treat) = -100,$$ $$Y(d_k=pass) = 0.$$ -In decision programming: +In Decision Programming the utility values are added using utility matrices. Notice that the utility values in the matrix are given in the same order as the states of node $h_N$ were defined when node $h_N$ was added. ```julia -for (i, j) in zip(treat, cost) - I_j = [i] - Y_j = zeros(S[I_j]...) - Y_j[1] = -100 - Y_j[2] = 0 - push!(V, ValueNode(j, I_j)) - push!(Y, Consequences(j, Y_j)) +for i in 1:N-1 + add_utilities!(diagram, "C$i", [-100.0, 0.0]) end ``` - -### Selling Price -The price of given the pig health at month $N$ is defined +The market price of the pig given its health at month $N$ is $$Y(h_N=ill) = 300,$$ $$Y(h_N=healthy) = 1000.$$ -In decision programming: +In Decision Programming: ```julia -for (i, j) in zip(health[end], price) - I_j = [i] - Y_j = zeros(S[I_j]...) - Y_j[1] = 300 - Y_j[2] = 1000 - push!(V, ValueNode(j, I_j)) - push!(Y, Consequences(j, Y_j)) -end +add_utilities!(diagram, "MP", [300.0, 1000.0]) ``` -### Validating Influence Diagram -Finally, we need to validate the influence diagram and sort the nodes, probabilities and consequences in increasing order by the node indices. +### Generate influence diagram +After adding nodes, generating arcs and defining probability and utility values, we generate the full influence diagram. By default this function uses the default path probabilities and utilities, which are defined as the joint probability of all chance events in the diagram and the sum of utilities in value nodes, respectively. In the [Contingent Portfolio Programming](contingent-portfolio-programming.md) example, we show how to use a user-defined custom path utility function. -```julia -validate_influence_diagram(S, C, D, V) -sort!.((C, D, V, X, Y), by = x -> x.j) -``` +In the pig breeding problem, when the $N$ is large some of the path utilities become negative. In this case, we choose to use the [positive path utility](../decision-programming/decision_model.md) transformation, which allows us to exclude the probability cut in the next section. -We define the path probability. ```julia -P = DefaultPathProbability(C, X) +generate_diagram!(diagram, positive_path_utility = true) ``` -As the path utility, we use the default, which is the sum of the consequences given the path. -```julia -U = DefaultPathUtility(V, Y) -``` - - -## Decision Model +## Decision model -We apply an affine transformation to the utility function, making all path utilities positive. Now that all path utilities are positive, the probability cut can be excluded from the model. The purpose of this is discussed in the [theoretical section](../decision-programming/decision-model.md) of this documentation. +Next we initialise the JuMP model and add the decision variables. Then we add the path compatibility variables. Since we applied an affine transformation to the utility function, making all path utilities positive, the probability cut can be excluded from the model. The purpose of this is discussed in the [theoretical section](../decision-programming/decision-model.md) of this documentation. ```julia -U⁺ = PositivePathUtility(S, U) model = Model() -z = DecisionVariables(model, S, D) -x_s = PathCompatibilityVariables(model, z, S, P, probability_cut = false) +z = DecisionVariables(model, diagram) +x_s = PathCompatibilityVariables(model, diagram, z, probability_cut = false) ``` We create the objective function ```julia -EV = expected_value(model, x_s, U⁺, P) +EV = expected_value(model, diagram, x_s) @objective(model, Max, EV) ``` @@ -227,88 +193,92 @@ optimize!(model) ``` -## Analyzing Results -### Decision Strategy +## Analyzing results -We obtain the optimal decision strategy: +Once the model is solved, we extract the results. The results are the decision strategy, state probabilities and utility distribution. ```julia Z = DecisionStrategy(z) +S_probabilities = StateProbabilities(diagram, Z) +U_distribution = UtilityDistribution(diagram, Z) ``` +### Decision strategy + +The optimal decision strategy is: + ```julia-repl -julia> print_decision_strategy(S, Z) -┌────────┬──────┬───┐ -│ Nodes │ (2,) │ 3 │ -├────────┼──────┼───┤ -│ States │ (1,) │ 2 │ -│ States │ (2,) │ 2 │ -└────────┴──────┴───┘ -┌────────┬──────┬───┐ -│ Nodes │ (5,) │ 6 │ -├────────┼──────┼───┤ -│ States │ (1,) │ 1 │ -│ States │ (2,) │ 2 │ -└────────┴──────┴───┘ -┌────────┬──────┬───┐ -│ Nodes │ (8,) │ 9 │ -├────────┼──────┼───┤ -│ States │ (1,) │ 1 │ -│ States │ (2,) │ 2 │ -└────────┴──────┴───┘ +julia> print_decision_strategy(diagram, Z, S_probabilities) +┌────────────────┬────────────────┐ +│ State(s) of T1 │ Decision in D1 │ +├────────────────┼────────────────┤ +│ positive │ pass │ +│ negative │ pass │ +└────────────────┴────────────────┘ +┌────────────────┬────────────────┐ +│ State(s) of T2 │ Decision in D2 │ +├────────────────┼────────────────┤ +│ positive │ treat │ +│ negative │ pass │ +└────────────────┴────────────────┘ +┌────────────────┬────────────────┐ +│ State(s) of T3 │ Decision in D3 │ +├────────────────┼────────────────┤ +│ positive │ treat │ +│ negative │ pass │ +└────────────────┴────────────────┘ ``` -The optimal strategy is as follows. In the first period, state 2 (no treatment) is chosen in node 3 ($d_1$) regardless of the state of node 2 ($t_1$). In other words, the pig is not treated in the first month. In the two subsequent months, state 1 (treat) is chosen if the corresponding test result is 1 (positive). +The optimal strategy is to not treat the pig in the first month regardless of if it is sick or not. In the two subsequent months, the pig should be treated if the test result is positive. -### State Probabilities +### State probabilities -The state probabilities for the strategy $Z$ can also be obtained. These tell the probability of each state in each node, given the strategy $Z$. +The state probabilities for strategy $Z$ tell the probability of each state in each node, given strategy $Z$. -```julia -sprobs = StateProbabilities(S, P, Z) -``` ```julia-repl -julia> print_state_probabilities(sprobs, health) -┌───────┬──────────┬──────────┬─────────────┐ -│ Node │ State 1 │ State 2 │ Fixed state │ -│ Int64 │ Float64 │ Float64 │ String │ -├───────┼──────────┼──────────┼─────────────┤ -│ 1 │ 0.100000 │ 0.900000 │ │ -│ 4 │ 0.270000 │ 0.730000 │ │ -│ 7 │ 0.295300 │ 0.704700 │ │ -│ 10 │ 0.305167 │ 0.694833 │ │ -└───────┴──────────┴──────────┴─────────────┘ -julia> print_state_probabilities(sprobs, test) -┌───────┬──────────┬──────────┬─────────────┐ -│ Node │ State 1 │ State 2 │ Fixed state │ -│ Int64 │ Float64 │ Float64 │ String │ -├───────┼──────────┼──────────┼─────────────┤ -│ 2 │ 0.170000 │ 0.830000 │ │ -│ 5 │ 0.289000 │ 0.711000 │ │ -│ 8 │ 0.306710 │ 0.693290 │ │ -└───────┴──────────┴──────────┴─────────────┘ -julia> print_state_probabilities(sprobs, treat) -┌───────┬──────────┬──────────┬─────────────┐ -│ Node │ State 1 │ State 2 │ Fixed state │ -│ Int64 │ Float64 │ Float64 │ String │ -├───────┼──────────┼──────────┼─────────────┤ -│ 3 │ 0.000000 │ 1.000000 │ │ -│ 6 │ 0.289000 │ 0.711000 │ │ -│ 9 │ 0.306710 │ 0.693290 │ │ -└───────┴──────────┴──────────┴─────────────┘ +julia> health_nodes = [["H$i" for i in 1:N]...] +julia> print_state_probabilities(diagram, S_probabilities, health_nodes) + +┌────────┬──────────┬──────────┬─────────────┐ +│ Node │ State 1 │ State 2 │ Fixed state │ +│ String │ Float64 │ Float64 │ String │ +├────────┼──────────┼──────────┼─────────────┤ +│ H1 │ 0.100000 │ 0.900000 │ │ +│ H2 │ 0.270000 │ 0.730000 │ │ +│ H3 │ 0.295300 │ 0.704700 │ │ +│ H4 │ 0.305167 │ 0.694833 │ │ +└────────┴──────────┴──────────┴─────────────┘ + +julia> test_nodes = [["T$i" for i in 1:N-1]...] +julia> print_state_probabilities(diagram, S_probabilities, test_nodes) +┌────────┬──────────┬──────────┬─────────────┐ +│ Node │ State 1 │ State 2 │ Fixed state │ +│ String │ Float64 │ Float64 │ String │ +├────────┼──────────┼──────────┼─────────────┤ +│ T1 │ 0.170000 │ 0.830000 │ │ +│ T2 │ 0.289000 │ 0.711000 │ │ +│ T3 │ 0.306710 │ 0.693290 │ │ +└────────┴──────────┴──────────┴─────────────┘ + +julia> treatment_nodes = [["D$i" for i in 1:N-1]...] +julia> print_state_probabilities(diagram, S_probabilities, treatment_nodes) +┌────────┬──────────┬──────────┬─────────────┐ +│ Node │ State 1 │ State 2 │ Fixed state │ +│ String │ Float64 │ Float64 │ String │ +├────────┼──────────┼──────────┼─────────────┤ +│ D1 │ 0.000000 │ 1.000000 │ │ +│ D2 │ 0.289000 │ 0.711000 │ │ +│ D3 │ 0.306710 │ 0.693290 │ │ +└────────┴──────────┴──────────┴─────────────┘ ``` -### Utility Distribution - -We can also print the utility distribution for the optimal strategy. The selling prices for a healthy and an ill pig are 1000DKK and 300DKK, respectively, while the cost of treatment is 100DKK. We can see that the probability of the pig being ill in the end is the sum of three first probabilities, approximately 30.5%. This matches the probability of state 1 in node 10 in the state probabilities shown above. +### Utility distribution -```julia -udist = UtilityDistribution(S, P, U, Z) -``` +We can also print the utility distribution for the optimal strategy. The selling prices for a healthy and an ill pig are 1000DKK and 300DKK, respectively, while the cost of treatment is 100DKK. We can see that the probability of the pig being ill in the end is the sum of three first probabilities, approximately 30.5%. This matches the probability of state $ill$ in the last node $h_4$ in the state probabilities shown above. ```julia-repl -julia> print_utility_distribution(udist) +julia> print_utility_distribution(U_distribution) ┌─────────────┬─────────────┐ │ Utility │ Probability │ │ Float64 │ Float64 │ @@ -325,7 +295,7 @@ julia> print_utility_distribution(udist) Finally, we print some statistics for the utility distribution. The expected value of the utility is 727DKK, the same as in [^1]. ```julia-repl -julia> print_statistics(udist) +julia> print_statistics(U_distribution) ┌──────────┬────────────┐ │ Name │ Statistics │ │ String │ Float64 │ diff --git a/docs/src/examples/used-car-buyer.md b/docs/src/examples/used-car-buyer.md index f1a0a697..2880c7ec 100644 --- a/docs/src/examples/used-car-buyer.md +++ b/docs/src/examples/used-car-buyer.md @@ -15,130 +15,127 @@ We now add two new features to the problem. A stranger approaches Joe and offers We present the new influence diagram above. The decision node $T$ denotes the decision to accept or decline the stranger's offer, and $R$ is the outcome of the test. We introduce new value nodes $V_1$ and $V_2$ to represent the testing costs and the base profit from purchasing the car. Additionally, the decision node $A$ now can choose to buy with a guarantee. +We start by defining the influence diagram structure. The nodes, as well as their information sets and states, are defined in the first block. Next, the influence diagram parameters consisting of the probabilities and utilities are defined. + + ```julia using JuMP, Gurobi using DecisionProgramming - -const O = 1 # Chance node: lemon or peach -const T = 2 # Decision node: pay stranger for advice -const R = 3 # Chance node: observation of state of the car -const A = 4 # Decision node: purchase alternative -const O_states = ["lemon", "peach"] -const T_states = ["no test", "test"] -const R_states = ["no test", "lemon", "peach"] -const A_states = ["buy without guarantee", "buy with guarantee", "don't buy"] - -S = States([ - (length(O_states), [O]), - (length(T_states), [T]), - (length(R_states), [R]), - (length(A_states), [A]), -]) -C = Vector{ChanceNode}() -D = Vector{DecisionNode}() -V = Vector{ValueNode}() -X = Vector{Probabilities}() -Y = Vector{Consequences}() +diagram = InfluenceDiagram() ``` -We start by defining the influence diagram structure. The decision and chance nodes, as well as their states, are defined in the first block. Next, the influence diagram parameters consisting of the node sets, probabilities, consequences and the state spaces of the nodes are defined. - -### Car's State +### Car's state -The chance node $O$ is defined by its information set $I(O)$ and probability distribution $X_O$. As seen in the influence diagram, the information set is empty and the node is a root node. The probability distribution is thus simply defined over the two states of $O$. +The chance node $O$ is defined by its name, its information set $I(O)$ and its states $lemon$ and $peach$. As seen in the influence diagram, the information set is empty and the node is a root node. ```julia -I_O = Vector{Node}() -X_O = [0.2, 0.8] -push!(C, ChanceNode(O, I_O)) -push!(X, Probabilities(O, X_O)) +add_node!(diagram, ChanceNode("O", [], ["lemon", "peach"])) ``` -### Stranger's Offer Decision +### Stranger's offer decision -A decision node is simply defined by its information state. +A decision node is also defined by its name, its information set and its states. ```julia -I_T = Vector{Node}() -push!(D, DecisionNode(T, I_T)) +add_node!(diagram, DecisionNode("T", [], ["no test", "test"])) ``` -### Test's Outcome +### Test's outcome -The second chance node, $R$, has nodes $O$ and $T$ in its information set, and the probabilities $ℙ(s_j∣𝐬_{I(j)})$ must thus be defined for all combinations of states in $O$, $T$ and $R$. +The second chance node $R$ has nodes $O$ and $T$ in its information set, and three states describing the situations of no test being done, and the test declaring the car to be a lemon or a peach. ```julia -I_R = [O, T] -X_R = zeros(S[O], S[T], S[R]) -X_R[1, 1, :] = [1,0,0] -X_R[1, 2, :] = [0,1,0] -X_R[2, 1, :] = [1,0,0] -X_R[2, 2, :] = [0,0,1] -push!(C, ChanceNode(R, I_R)) -push!(X, Probabilities(R, X_R)) +add_node!(diagram, ChanceNode("R", ["O", "T"], ["no test", "lemon", "peach"])) ``` -### Purchace Decision +### Purchace decision +The purchase decision represented by node $A$ is added as follows. ```julia -I_A = [R] -push!(D, DecisionNode(A, I_A)) +add_node!(diagram, DecisionNode("A", ["R"], ["buy without guarantee", "buy with guarantee", "don't buy"])) ``` +### Testing fee, base profit and repair costs -### Testing Cost +Value nodes are defined by only their names and information sets because they do not have states. Instead, value nodes map their information states to utility values which will be added later on. +```julia +add_node!(diagram, ValueNode("V1", ["T"])) +add_node!(diagram, ValueNode("V2", ["A"])) +add_node!(diagram, ValueNode("V3", ["O", "A"])) +``` -We continue by defining the utilities (consequences) associated with value nodes. The value nodes are defined similarly as the chance nodes, except that instead of probabilities, we define consequences $Y_j(𝐬_{I(j)})$. Value nodes can be named just like the other nodes, e.g. $V1 = 5$, but considering that the index of value nodes is not needed elsewhere (value nodes can't be in information sets), we choose to simply use the index number when creating the node. +### Generate arcs +Now that all of the nodes have been added to our influence diagram we generate the arcs between the nodes. This step automatically orders the nodes, gives them indices and reorganises the information into the appropriate form in the influence diagram structure. +```julia +generate_arcs!(diagram) +``` +### Probabilities +We continue by defining probability distributions for each chance node. + +Node $O$ is a root node and has two states thus, its probability distribution is simply defined over the two states. We can use the `ProbabilityMatrix` structure in creating the probability matrix easily without having to worry about the matrix dimensions. We then set the probability values and add the probabililty matrix to the influence diagram. ```julia -I_V1 = [T] -Y_V1 = [0.0, -25.0] -push!(V, ValueNode(5, I_V1)) -push!(Y, Consequences(5, Y_V1)) +X_O = ProbabilityMatrix(diagram, "O") +set_probability!(X_O, ["peach"], 0.8) +set_probability!(X_O, ["lemon"], 0.2) +add_probabilities!(diagram, "O", X_O) ``` -### Base Profit of Purchase +Node $R$ has two nodes in its information set and three states. The probabilities $P(s_j \mid s_{I(j)})$ must thus be defined for all combinations of states in $O$, $T$ and $R$. We declare the probability distribution over the states of node $R$ for each information state in the following way. More information on defining probability matrices can be found on the [usage page](../usage.md). ```julia -I_V2 = [A] -Y_V2 = [100.0, 40.0, 0.0] -push!(V, ValueNode(6, I_V2)) -push!(Y, Consequences(6, Y_V2)) +X_R = ProbabilityMatrix(diagram, "R") +set_probability!(X_R, ["lemon", "no test", :], [1,0,0]) +set_probability!(X_R, ["lemon", "test", :], [0,1,0]) +set_probability!(X_R, ["peach", "no test", :], [1,0,0]) +set_probability!(X_R, ["peach", "test", :], [0,0,1]) +add_probabilities!(diagram, "R", X_R) ``` -### Repairing Cost -The rows of the consequence matrix Y_V3 correspond to the state of the car, while the columns correspond to the decision made in node $A$. +### Utilities +We continue by defining the utilities associated with the information states of the value nodes. The utilities $Y_j(𝐬_{I(j)})$ are defined and added similarly to the probabilities. + +Value node $V1$ has only node $T$ in its information set and node $T$ only has two states. Therefore, the utility matrix of node $V1$ should hold utility values corresponding to states $test$ and $no \ test$. ```julia -I_V3 = [O, A] -Y_V3 = [-200.0 0.0 0.0; - -40.0 -20.0 0.0] -push!(V, ValueNode(7, I_V3)) -push!(Y, Consequences(7, Y_V3)) +Y_V1 = UtilityMatrix(diagram, "V1") +set_utility!(Y_V1, ["test"], -25) +set_utility!(Y_V1, ["no test"], 0) +add_utilities!(diagram, "V1", Y_V1) ``` -### Validating the Influence Diagram -Validate influence diagram and sort nodes, probabilities and consequences +We then define the utilities associated with the base profit of the purchase in different scenarios. +```julia +Y_V2 = UtilityMatrix(diagram, "V2") +set_utility!(Y_V2, ["buy without guarantee"], 100) +set_utility!(Y_V2, ["buy with guarantee"], 40) +set_utility!(Y_V2, ["don't buy"], 0) +add_utilities!(diagram, "V2", Y_V2) +``` +Finally, we define the utilities corresponding to the repair costs. The rows of the utilities matrix `Y_V3` correspond to the state of the car, while the columns correspond to the decision made in node $A$. Notice that the utility values for the second row are added as a vector, in this case it is important to give the utility values in the correct order. The order of the columns is determined by the order in which the states are given when declaring node $A$. See the [usage page](../usage.md) for more information on the syntax. ```julia -validate_influence_diagram(S, C, D, V) -sort!.((C, D, V, X, Y), by = x -> x.j) +Y_V3 = UtilityMatrix(diagram, "V3") +set_utility!(Y_V3, ["lemon", "buy without guarantee"], -200) +set_utility!(Y_V3, ["lemon", "buy with guarantee"], 0) +set_utility!(Y_V3, ["lemon", "don't buy"], 0) +set_utility!(Y_V3, ["peach", :], [-40, -20, 0]) +add_utilities!(diagram, "V3", Y_V3) ``` -Default path probabilities and utilities are defined as the joint probability of all chance events in the diagram and the sum of utilities in value nodes, respectively. In the [Contingent Portfolio Programming](contingent-portfolio-programming.md) example, we show how to use a user-defined custom path utility function. +### Generate influence diagram +Finally, generate the full influence diagram before defining the decision model. By default this function uses the default path probabilities and utilities, which are defined as the joint probability of all chance events in the diagram and the sum of utilities in value nodes, respectively. In the [Contingent Portfolio Programming](contingent-portfolio-programming.md) example, we show how to use a user-defined custom path utility function. ```julia -P = DefaultPathProbability(C, X) -U = DefaultPathUtility(V, Y) +generate_diagram!(diagram) ``` - -## Decision Model -We then construct the decision model using the DecisionProgramming.jl package, using the expected value as the objective. +## Decision model +We then construct the decision model by declaring a JuMP model and adding decision variables and path compatibility variables to the model. We define the objective function to be the expected value. ```julia model = Model() -z = DecisionVariables(model, S, D) -x_s = PathCompatibilityVariables(model, z, S, P) -EV = expected_value(model, x_s, U, P) +z = DecisionVariables(model, diagram) +x_s = PathCompatibilityVariables(model, diagram, z) +EV = expected_value(model, diagram, x_s) @objective(model, Max, EV) ``` @@ -154,41 +151,36 @@ optimize!(model) ``` -## Analyzing Results -### Decision Strategy -Once the model is solved, we obtain the following decision strategy: +## Analyzing results +Once the model is solved, we extract the results. ```julia Z = DecisionStrategy(z) +S_probabilities = StateProbabilities(diagram, Z) +U_distribution = UtilityDistribution(diagram, Z) ``` -```julia-repl -julia> print_decision_strategy(S, Z) -┌────────┬────┬───┐ -│ Nodes │ () │ 2 │ -├────────┼────┼───┤ -│ States │ () │ 2 │ -└────────┴────┴───┘ -┌────────┬──────┬───┐ -│ Nodes │ (3,) │ 4 │ -├────────┼──────┼───┤ -│ States │ (1,) │ 3 │ -│ States │ (2,) │ 2 │ -│ States │ (3,) │ 1 │ -└────────┴──────┴───┘ -``` - -To start explaining this output, let's take a look at the top table. On the right, we have the decision node 2. We defined earlier that the node $T$ is node number 2. On the left, we have the information set of that decision node, which is empty. The strategy in the first decision node is to choose alternative 2, which we defined to be testing the car. - -In the bottom table, we have node number 4 (node $A$) and its predecessor, node number 3 (node $R$). The first row, where we obtain no test result, is invalid for this strategy since we tested the car. If the car is a lemon, Joe should buy the car with a guarantee (alternative 2), and if it is a peach, buy the car without guarantee (alternative 1). - -### Utility Distribution -```julia -udist = UtilityDistribution(S, P, U, Z) -``` +### Decision strategy +We obtain the following optimal decision strategy: +```julia-repl +julia> print_decision_strategy(diagram, Z, S_probabilities) +┌───────────────┐ +│ Decision in T │ +├───────────────┤ +│ test │ +└───────────────┘ +┌───────────────┬───────────────────────┐ +│ State(s) of R │ Decision in A │ +├───────────────┼───────────────────────┤ +│ lemon │ buy with guarantee │ +│ peach │ buy without guarantee │ +└───────────────┴───────────────────────┘ +``` + +### Utility distribution ```julia-repl -julia> print_utility_distribution(udist) +julia> print_utility_distribution(U_distribution) ┌───────────┬─────────────┐ │ Utility │ Probability │ │ Float64 │ Float64 │ @@ -201,7 +193,7 @@ julia> print_utility_distribution(udist) From the utility distribution, we can see that Joe's profit with this strategy is 15 USD, with a 20% probability (the car is a lemon) and 35 USD with an 80% probability (the car is a peach). ```julia-repl -julia> print_statistics(udist) +julia> print_statistics(U_distribution) ┌──────────┬────────────┐ │ Name │ Statistics │ │ String │ Float64 │ diff --git a/docs/src/figures/2chance_1decision_1value.svg b/docs/src/figures/2chance_1decision_1value.svg new file mode 100644 index 00000000..1b4d05b4 --- /dev/null +++ b/docs/src/figures/2chance_1decision_1value.svg @@ -0,0 +1,3 @@ + + +
C1 \\ \{x, y,...
 D1 \\ \{a,b\}
 C2 \\ \{v,w\}
 VViewer does not support full SVG 1.1 \ No newline at end of file diff --git a/docs/src/figures/chance-node.svg b/docs/src/figures/chance-node.svg deleted file mode 100644 index 019a220b..00000000 --- a/docs/src/figures/chance-node.svg +++ /dev/null @@ -1,3 +0,0 @@ - - -
2 \\ \{1,2,3\...
 1 \\ \{1,2\}
 3 \\ \{1,2\}Viewer does not support full SVG 1.1 \ No newline at end of file diff --git a/docs/src/figures/decision-node.svg b/docs/src/figures/decision-node.svg deleted file mode 100644 index c96f4bab..00000000 --- a/docs/src/figures/decision-node.svg +++ /dev/null @@ -1,3 +0,0 @@ - - -
2 \\ \{1,2,3\...
 1 \\ \{1,2\}
 3 \\ \{1,2\}Viewer does not support full SVG 1.1 \ No newline at end of file diff --git a/docs/src/figures/getting-started.drawio b/docs/src/figures/getting-started.drawio deleted file mode 100644 index e3bc6a09..00000000 --- a/docs/src/figures/getting-started.drawio +++ /dev/null @@ -1 +0,0 @@ -7Vldb5swFP01kbaHSWAIIY8LbVdNm/ZQaZv25oABdwZTY/KxX79rbJLQhjTZmpKoSVBiHxt/3XOObsLACbLFJ4GL9CuPCBsgK1oMnKsBQmPfhk8FLDXgWb4GEkEjDdlr4I7+IQa0DFrRiJStjpJzJmnRBkOe5ySULQwLweftbjFn7VkLnJAnwF2I2VP0B41kqlF/aK3xW0KTtJnZtkxLhpvOZogyxRGfa6ju41wPnEBwLnUpWwSEqbNrzkUPdNPRulqYILnc5wbv+9Rb5veLb78ePpflz+JPtEQfhmZtctlsmESwf1PNeQ5fE8GrPCJqGAtqXMiUJzzH7AvnBYA2gPdEyqWJHq4kByiVGTOtMWUs4IyLegonHqo34KUU/DfZaPHq16qlOXE4q4leqFpd5/6bg+aVCMmOTTc8wiIhckc/ZxUlYDfhGZFiCfcJwrCks/Y6sOFZsuq3DgUUTDQOiIwZd4ZZZWYaII9JdTYFzlsx8x4qRaJJqA/yo1pjMn0HdIALZrc2Su/X3aGU6G8XLjXZYBioSxdGE1hCAHDg1NUr3a9ZBWxKL6QZ5DGPGAONKvrMUyrJXYHrmMzBJdrkwGWhhRvThSLZP7Al5rncwPUL8IIICtEgQs1O88RQqYtcMyIkWeym11M6mBtco3vjc6jxgfmGaxgo3TCMBntx/rhvUdnoHJSNtihbKcvuUmBLfY+C+p/iiobEj9xtAfPR1KkD1qOI3J415LyyBzt7MeDiv9uoszJcwx0b9Uwe2+nDgfcQ9HFzK3dPC7ZHvWZX7iW9Oit5O+jE8iu7l59OvcvbOw95e5cca18l9Z1kaaacaZZ1MDXASm9wRpk6/FvCZkTSEJsGYwK2CjRmNMmhEkJUiThQ8F0E7KDZSyTu/qllX34f9hzHxAvDbdGKRuOpZR3dnsd72jNCvdrz+JJ9nVX25Vknln019vK25L16AnLa8m4GvmRfzyup7+wLbfs38sjZ1yFOK1KeTavyead9wcxqh5i7DOF4mdXIe7XMCqrrx5F128YzXef6Lw== \ No newline at end of file diff --git a/docs/src/figures/usage.drawio b/docs/src/figures/usage.drawio new file mode 100644 index 00000000..8d5d492c --- /dev/null +++ b/docs/src/figures/usage.drawio @@ -0,0 +1 @@ +7Vldb5swFP01kbaHSWASkj4u6ZemTXuo1E17c+AC7gymxvmgv342thNoPppsWUi2pFGA4wu27z3n6FZ0vFE6v+M4T76wEGgHOeG84113EBr0ffmrgFIDPnI0EHMSashdAg/kBQxowyYkhKIRKBijguRNMGBZBoFoYJhzNmuGRYw2Z81xDCvAQ4DpKvqNhCIx2+o5S/weSJzYmV3HjKTYBptHFAkO2UxDVYx30/FGnDGhz9L5CKjKnc2LftDthtHFwjhkYpcb/MexX2ZP868/nj8Vxff8JSzRh55ZmyjthiGU+zeXGcvkYcjZJAtBPcaRV4yLhMUsw/QzY7kEXQk+gRClqR6eCCahRKTUjEaE0hGjjFdTeFFP/Um8EJz9hNqIX30WIzbjMlfD1e3avLIJD2DLHi1tMI9BbInzdJxKQG0Ck8w7YCkIXsoADhQLMm0SBBuexYu4ZSnkianGHpUxz51iOjEzdZBPhcpNjrNGzfzniSLRMNCJ/KjWGI/fSTrIr5zdqZ29X4bLs1gfu/I7clVQb6S++qQ/nOubSn14qcBrHW4XI/em12Of9ZpOlEqpKhbNEiLgIcdVrWbSLJocwUWu9RuRueLab5AmYpmo4foj8Rw4kUUBrmYnWWwYtQfHpsAFzLeywox2jfyN3SFrB7OaeRgoqfmGxQ5Oo+5/IHB0lgJHawSulHW9RohYKnDclN+rqv6husIeDMLuuooN0NirKtaiiroti8g7thejVQpM5U2ziwO/zZ2F5RryuKhl9rheGx68g6IP2mR1dzRht39SLmzXfemzzkvlHjqxRstt5V+pY6vcP1OV+5dma2cptd1tWe6sFmtNWwSyWNFhiyXd7RanhKp03AOdgiABNgNGl65KPaYkzuRFIFMPfE8NbqLEhsIfoqcenFpfNGjDMaMI/CBYV62wfzV2nEM75tWOjonQaTnm1aUvOse+yHdOrC+yLvNPq3zxauLMVL7g16UveltKbfdFljtHs+TH/ayWJywdT4q3rfaAHdYWNW9yhL/XYfX9o3VY8nL55rAaq71+9W5+AQ== \ No newline at end of file diff --git a/docs/src/figures/value-node.svg b/docs/src/figures/value-node.svg deleted file mode 100644 index b5405551..00000000 --- a/docs/src/figures/value-node.svg +++ /dev/null @@ -1,3 +0,0 @@ - - -
2 \\ \{1,2,3\...
 1 \\ \{1,2\}
 3Viewer does not support full SVG 1.1 \ No newline at end of file diff --git a/docs/src/usage.md b/docs/src/usage.md index c63393aa..46d84713 100644 --- a/docs/src/usage.md +++ b/docs/src/usage.md @@ -1,61 +1,265 @@ # Usage -On this page, we demonstrate common patterns for expressing influence diagrams and creating decision models using `DecisionProgramming.jl`. We can import the package with the `using` keyword. +On this page, we demonstrate common patterns for expressing influence diagrams and creating decision models using `DecisionProgramming.jl`. We also discuss the abstraction that is created using the influence diagram structure. We can import the package with the `using` keyword. ```julia using DecisionProgramming ``` -## Chance Nodes -![](figures/chance-node.svg) +## Adding nodes +![](figures/2chance_1decision_1value.svg) -Given the above influence diagram, we can create the `ChanceNode` and `Probabilities` structures for the node `3` as follows: +Given the above influence diagram, we express it as a Decision Programming model as follows. We create `ChanceNode` and `DecisionNode` instances and add them to the influence diagram. Creating a `ChanceNode` or `DecisionNode` requires giving it a unique name, its information set and its states. If the node is a root node, the information set is left empty using square brackets. The order in which nodes are added does not matter. ```julia -S = States([2, 3, 2]) -j = 3 -I_j = Node[1, 2] -X_j = zeros(S[I_j]..., S[j]) -X_j[1, 1, :] = [0.1, 0.9] -X_j[1, 2, :] = [0.0, 1.0] -X_j[1, 3, :] = [0.3, 0.7] -X_j[2, 1, :] = [0.2, 0.8] -X_j[2, 2, :] = [0.4, 0.6] -X_j[2, 3, :] = [1.0, 0.0] -ChanceNode(j, I_j) -Probabilities(j, X_j) +diagram = InfluenceDiagram() +add_node!(diagram, DecisionNode("D1", [], ["a", "b"])) +add_node!(diagram, ChanceNode("C2", ["D1", "C1"], ["v", "w"])) +add_node!(diagram, ChanceNode("C1", [], ["x", "y", "z"])) ``` +Value nodes are added by simply giving it a name and its information set. Value nodes do not have states because their purpose is to map their information state to utility values. -## Decision Nodes -![](figures/decision-node.svg) +```julia +add_node!(diagram, ValueNode("V", ["C2"])) +``` + +Once all the nodes are added, we generate the arcs. This orders the nodes and gives each one a number, so that their predecessors have numbers less than theirs. In effect, the chance and decision nodes are numbered such that $C \cup D = \{1,...,n\}$, where $n = \mid C\mid + \mid D\mid$. The value nodes are numbered $V = \{n+1,..., N\}$, where $N = \mid C\mid + \mid D\mid + \mid V \mid$. For more details on influence diagrams see page [influence diagram](decision-programming/influence-diagram.md). +```julia +generate_arcs!(diagram) +``` +Now the fields `Names`, `I_j`, `States`, `S`, `C`, `D` and `V` in the influence diagram structure have been properly filled. The `Names` field holds the names of all nodes in the order of their numbers. From this we can see that node D1 has been numbered 1, node C1 has been numbered 2 and node C2 has been numbered 3. Field `I_j` holds the information sets of each node. Notice, that the nodes are identified by their numbers. Field `States` holds the names of the states of each node and field `S` holds the number of states each node has. Fields `C`, `D` and `V` contain the chance, decision and value nodes respectively. + +```julia +julia> diagram.Names +4-element Array{String,1}: + "D1" + "C1" + "C2" + "V" + +julia> diagram.I_j +4-element Array{Array{Int16,1},1}: + [] + [] + [1, 2] + [3] + +julia> diagram.States +3-element Array{Array{String,1},1}: + ["a", "b"] + ["x", "y", "z"] + ["v", "w"] + +julia> diagram.S +3-element States: + 2 + 3 + 2 + +julia> diagram.C +2-element Array{Int16,1}: + 2 + 3 + +julia> diagram.D +1-element Array{Int16,1}: + 1 + +julia> diagram.V +1-element Array{Int16,1}: + 4 +``` + +## Probability Matrices +Each chance node needs a probability matrix which describes the probability distribution over its states given an information state. It holds probability values +$$ℙ(X_j=s_j∣X_{I(j)}=𝐬_{I(j)})$$ + +for all $s_j \in S_j$ and $𝐬_{I(j)} \in 𝐒_{I(j)}$. + +Thus, the probability matrix of a chance node needs to have dimensions that correspond to the number of states of the nodes in its information set and number of state of the node itself. + +For example, the node C1 in the influence diagram above has an empty information set and three states $x, y$, and $z$. Therefore its probability matrix needs dimensions (3,1). If the probabilities of events $x, y$, and $z$ occuring are $10\%, 30\%$ and $60\%$, then the probability matrix $X_{C1}$ should be $[0.1 \quad 0.3 \quad 0.6]$. The order of the probability values is determined by the order in which the states are given when the node is added. The states are also stored in this order in the `States` vector. + +In Decision Programming the probability matrix of node C1 can be added in the following way. Note, that probability matrices can only be added after the arcs have been generated. + +```julia +# How C1 was added: add_node!(diagram, ChanceNode("C1", [], ["x", "y", "z"])) +X_C1 = [0.1, 0.3, 0.6] +add_probabilities!(diagram, "C1", X_C1) +``` + +The `add_probabilities!` function adds the probability matrix as a `Probabilities` structure into the influence diagram's `X` field. +```julia +julia> diagram.X +1-element Array{Probabilities,1}: + [0.1, 0.3, 0.6] +``` + + +As another example, we will add the probability matrix of node C2. It has two nodes in its information set: C1 and D1. These nodes have 3 and 2 state respectively. Node C2 itself has 2 states. The question is should the dimensions of the probability matrix be $(|S_{C1}|, |\ S_{D1}|, |\ S_{C2}|) = (3, 2, 2)$ or $(|S_{D1}|, |\ S_{C1}|, \ |S_{C2}|) = (2, 3, 2)$? The answer is that the dimensions should be in ascending order of the nodes' numbers that they correspond to. This is also the order that the information set is in in the field `I_j`. In this case the influence diagram looks like this: +```julia +julia> diagram.Names +4-element Array{String,1}: + "D1" + "C1" + "C2" + "V" + + julia> diagram.I_j +4-element Array{Array{Int16,1},1}: + [] + [] + [1, 2] + [3] + + julia> diagram.S +3-element States: + 2 + 3 + 2 +``` + +Therefore, the probability matrix of node C2 should have dimensions $(|S_{D1}|, |\ S_{C1}|, \ |S_{C2}|) = (2, 3, 2)$. The probability matrix can be added by declaring the matrix and then filling in the probability values as shown below. +```julia +X_C2 = zeros(2, 3, 2) +X_C2[1, 1, 1] = ... +X_C2[1, 1, 2] = ... +X_C2[1, 1, 2] = ... +⋮ +add_probabilities!(diagram, "C2", X_C2) +``` +In order to be able to fill in the probability values, it is crucial to understand what the matrix indices represent. The indices represent a subpath in the influence diagram. The states in the path are referred to with their numbers instead of with their names. The states of a node are numbered according to their positions in the vector of states in field `States`. The order of the states of each node is seen below. From this, we can deduce that for nodes D1, C1, C2 the subpath `(1,1,1)` corresponds to subpath $(a, x, v)$ and subpath `(1, 3, 2)` corresponds to subpath $(a, z, w)$. Therefore, the probability value at `X_C2[1, 3, 2]` should be the probability of the scenario $(a, z, w)$ occuring. +```julia +julia> diagram.States +3-element Array{Array{String,1},1}: + ["a", "b"] + ["x", "y", "z"] + ["v", "w"] +``` +### Helper Syntax +Figuring out the dimensions of a probability matrix and adding the probability values is difficult. Therefore, we have implemented an easier syntax. + +A probability matrix can be initialised with the correct dimensions using the `ProbabilityMatrix` function. It initiliases the probability matrix with zeros. +```julia +julia> X_C2 = ProbabilityMatrix(diagram, "C2") +2×3×2 ProbabilityMatrix{3}: +[:, :, 1] = + 0.0 0.0 0.0 + 0.0 0.0 0.0 -Given the above influence diagram, we can create the `DecisionNode` structure for the node `3` as follows: +[:, :, 2] = + 0.0 0.0 0.0 + 0.0 0.0 0.0 + +julia> size(X_C2) +(2, 3, 2) +``` +A matrix of type `ProbabilityMatrix` can be filled using the names of the states. The states must however be given in the correct order, according to the order of the nodes in the information set vector `I_j`. Notice that if we use the `Colon` (`:`) to indicate several elements of the matrix, the probability values have to be given in the correct order of the states in `States`. +```julia +julia> set_probability!(X_C2, ["a", "z", "w"], 0.25) +0.25 + +julia> set_probability!(X_C2, ["z", "a", "v"], 0.75) +ERROR: DomainError with Node D1 does not have a state called z.: + +julia> set_probability!(X_C2, ["a", "z", "v"], 0.75) +0.75 + +julia> set_probability!(X_C2, ["a", "x", :], [0.3, 0.7]) +2-element Array{Float64,1}: + 0.3 + 0.7 +``` + +A matrix of type `ProbabilityMatrix` can also be filled using the matrix indices if that is more convient. The following achieves the same as what was done above. +```julia +julia> X_C2[1, 3, 2] = 0.25 +0.25 + +julia> X_C2[1, 3, 1] = 0.75 +0.75 + +julia> X_C2[1, 1, :] = [0.3, 0.7] +2-element Array{Float64,1}: + 0.3 + 0.7 +```` + +Now, the probability matrix X_C2 is partially filled. +```julia +julia> X_C2 +2×3×2 ProbabilityMatrix{3}: +[:, :, 1] = + 0.3 0.0 0.75 + 0.0 0.0 0.0 + +[:, :, 2] = + 0.7 0.0 0.25 + 0.0 0.0 0.0 +``` + +The probability matrix can be added to the influence diagram once it has been filled with probability values. The probability matrix of node C2 is added exactly like before, despite X_C2 now being a matrix of type `ProbabilityMatrix`. +```julia +julia> add_probabilities!(diagram, "C2", X_C2) +``` + +## Utility Matrices +Each value node maps its information states to utility values. In Decision Programming the utility values are passed to the influence diagram using utility matrices. Utility matrices are very similar to probability matrices of chance nodes. There are only two important differences. First, the utility matrices hold utility values instead of probabilities, meaning that they do not need to sum to one. Second, since value nodes do not have states, the cardinality of a utility matrix depends only on the number of states of the nodes in the information set. + +As an example, the utility matrix of node V should have dimensions (2,1) because its information set consists of node C2, which has two states. If state $v$ of node C2 yields a utility of -100 and state $w$ yields utility of 400, then the utility matrix of node V can be added in the following way. Note, that utility matrices can only be added after the arcs have been generated. ```julia -S = States([2, 3, 2]) -j = 3 -I_j = Node[1, 2] -DecisionNode(j, I_j) +julia> Y_V = zeros(2) +2-element Array{Float64,1}: + 0.0 + 0.0 + +julia> Y_V[1] = -100 +-100 + +julia> Y_V[2] = 400 +400 + +julia> add_utilities!(diagram, "V", Y_V) ``` -## Value Nodes -![](figures/value-node.svg) +The other option is to add the utility matrix using the `UtilityMatrix` type. This is very similar to the `ProbabilityMatrix` type. The `UtilityMatrix` function initialises the values to `Inf`. Using the `UtilityMatrix` type's functionalities, the utility matrix of node V could also be added like shown below. This achieves the exact same result as we did above with the more abstract syntax. + -Given the above influence diagram, we can create `ValueNode` and `Consequences` structures for node `3` as follows: +```julia +julia> Y_V = UtilityMatrix(diagram, "V") +2-element UtilityMatrix{1}: + Inf + Inf + +julia> set_utility!(Y_V, ["w"], 400) +400 + +julia> set_utility!(Y_V, ["v"], -100) +-100 + +julia> add_utilities!(diagram, "V", Y_V) +``` +The `add_utilities!` function adds the utility matrix as a `Utilities` structure into the influence diagram's `Y` field. ```julia -S = States([2, 3]) -j = 3 -I_j = [1, 2] -Y_j = zeros(S[I_j]...) -Y_j[1, 1] = -1.3 -Y_j[1, 2] = 2.5 -Y_j[1, 3] = 0.1 -Y_j[2, 1] = 0.0 -Y_j[2, 2] = 3.2 -Y_j[2, 3] = -2.7 -ValueNode(j, I_j) -Consequences(j, Y_j) +julia> diagram.Y +1-element Array{Utilities,1}: + [-100.0, 400.0] ``` + +## Generating the influence diagram + +The final part of modeling an influence diagram using the Decision Programming package is generating the full influence diagram. This is done using the `generate_diagram!` function. +```julia +generate_diagram!(diagram) +``` +In this function, first, the probability and utility matrices in fields `X` and `Y` are sorted according to the chance and value nodes' indices. + +Second, the path probability and path utility types are declared and added into fields `P` and `U` respectively. These types define how the path probability $p(𝐬)$ and path utility $\mathcal{U}(𝐬)$ are defined in the model. By default, the function will set them to default path probability and default path utility. See the [influence diagram](decision-programming/influence-diagram.md) for more information on default path probability and utility. + + + + diff --git a/examples/CHD_preventative_care.jl b/examples/CHD_preventative_care.jl index d1c34561..24870234 100644 --- a/examples/CHD_preventative_care.jl +++ b/examples/CHD_preventative_care.jl @@ -6,12 +6,11 @@ using CSV, DataFrames, PrettyTables # Setting subproblem specific parameters -const chosen_risk_level = 13 +const chosen_risk_level = "12%" # Reading tests' technical performance data (dummy data in this case) -data = CSV.read("CHD_preventative_care_data.csv", DataFrame) - +data = CSV.read("examples/CHD_preventative_care_data.csv", DataFrame) # Bayes posterior risk probabilities calculation function @@ -73,13 +72,13 @@ function state_probabilities(risk_p::Array{Float64}, t::Int64, h::Int64, prior:: #if no test is performed, then the probabilities of moving to states (other than the prior risk level) are 0 and to the prior risk element is 1 if t == 3 - state_probabilites = zeros(101, 1) + state_probabilites = zeros(101) state_probabilites[prior] = 1.0 return state_probabilites end # return vector - state_probabilites = zeros(101,1) + state_probabilites = zeros(101) # copying the probabilities of the scores for ease of readability if h == 1 && t == 1 # CHD and TRS @@ -108,157 +107,84 @@ function state_probabilities(risk_p::Array{Float64}, t::Int64, h::Int64, prior:: end -function analysing_results(Z::DecisionStrategy, sprobs::StateProbabilities) - - d = Z.D[1] #taking one of the decision nodes to retrieve the information_set_R - information_set_R = vec(collect(paths(S[d.I_j]))) - results = DataFrame(Information_set = map( x -> string(x) * "%", [0:1:100;])) - # T1 - Z_j = Z.Z_j[1] - probs = map(x -> x > 0 ? 1 : 0, get(sprobs.probs, 1,0)) #these are zeros and ones - dec = [Z_j(s_I) for s_I in information_set_R] - results[!, "T1"] = map(x -> x == 0 ? "" : "$x", probs.*dec) - - # T2 - Z_j = Z.Z_j[2] - probs = map(x -> x > 0 ? 1 : 0, (get(sprobs.probs, 4,0))) #these are zeros and ones - dec = [Z_j(s_I) for s_I in information_set_R] - results[!, "T2"] = map(x -> x == 0 ? "" : "$x", probs.*dec) - - # TD - Z_j = Z.Z_j[3] - probs = map(x -> x > 0 ? 1 : 0, (get(sprobs.probs, 6,0))) #these are zeros and ones - dec = [Z_j(s_I) for s_I in information_set_R] - results[!, "TD"] = map(x -> x == 0 ? "" : "$x", probs.*dec) - - pretty_table(results) -end - - -const R0 = 1 -const H = 2 -const T1 = 3 -const R1 = 4 -const T2 = 5 -const R2 = 6 -const TD = 7 -const TC = 8 -const HB = 9 - +@info("Creating the influence diagram.") +diagram = InfluenceDiagram() const H_states = ["CHD", "no CHD"] const T_states = ["TRS", "GRS", "no test"] const TD_states = ["treatment", "no treatment"] -const R_states = map( x -> string(x) * "%", [0:1:100;]) -const TC_states = ["TRS", "GRS", "TRS & GRS", "no tests"] -const HB_states = ["CHD & treatment", "CHD & no treatment", "no CHD & treatment", "no CHD & no treatment"] +const R_states = [string(x) * "%" for x in [0:1:100;]] -@info("Creating the influence diagram.") -S = States([ - (length(R_states), [R0, R1, R2]), - (length(H_states), [H]), - (length(T_states), [T1, T2]), - (length(TD_states), [TD]) -]) -C = Vector{ChanceNode}() -D = Vector{DecisionNode}() -V = Vector{ValueNode}() -X = Vector{Probabilities}() -Y = Vector{Consequences}() +add_node!(diagram, ChanceNode("R0", [], R_states)) +add_node!(diagram, ChanceNode("R1", ["R0", "H", "T1"], R_states)) +add_node!(diagram, ChanceNode("R2", ["R1", "H", "T2"], R_states)) +add_node!(diagram, ChanceNode("H", ["R0"], H_states)) +add_node!(diagram, DecisionNode("T1", ["R0"], T_states)) +add_node!(diagram, DecisionNode("T2", ["R1"], T_states)) +add_node!(diagram, DecisionNode("TD", ["R2"], TD_states)) -I_R0 = Vector{Node}() -X_R0 = zeros(S[R0]) -X_R0[chosen_risk_level] = 1 -push!(C, ChanceNode(R0, I_R0)) -push!(X, Probabilities(R0, X_R0)) +add_node!(diagram, ValueNode("TC", ["T1", "T2"])) +add_node!(diagram, ValueNode("HB", ["H", "TD"])) -I_H = [R0] -X_H = zeros(S[R0], S[H]) -X_H[:, 1] = data.risk_levels # 1 = "CHD" -X_H[:, 2] = 1 .- X_H[:, 1] # 2 = "no CHD" -push!(C, ChanceNode(H, I_H)) -push!(X, Probabilities(H, X_H)) +generate_arcs!(diagram) +X_R0 = ProbabilityMatrix(diagram, "R0") +set_probability!(X_R0, [chosen_risk_level], 1) +add_probabilities!(diagram, "R0", X_R0) -I_T1 = [R0] -push!(D, DecisionNode(T1, I_T1)) +X_H = ProbabilityMatrix(diagram, "H") +set_probability!(X_H, [:, "CHD"], data.risk_levels) +set_probability!(X_H, [:, "no CHD"], 1 .- data.risk_levels) +add_probabilities!(diagram, "H", X_H) -I_R1 = [R0, H, T1] -X_R1 = zeros(S[I_R1]..., S[R1]) +X_R = ProbabilityMatrix(diagram, "R1") for s_R0 = 1:101, s_H = 1:2, s_T1 = 1:3 - X_R1[s_R0, s_H, s_T1, :] = state_probabilities(update_risk_distribution(s_R0, s_T1), s_T1, s_H, s_R0) + X_R[s_R0, s_H, s_T1, :] = state_probabilities(update_risk_distribution(s_R0, s_T1), s_T1, s_H, s_R0) end -push!(C, ChanceNode(R1, I_R1)) -push!(X, Probabilities(R1, X_R1)) - - -I_T2 = [R1] -push!(D, DecisionNode(T2, I_T2)) +add_probabilities!(diagram, "R1", X_R) +add_probabilities!(diagram, "R2", X_R) -I_R2 = [H, R1, T2] -X_R2 = zeros(S[I_R2]..., S[R2]) -for s_R1 = 1:101, s_H = 1:2, s_T2 = 1:3 - X_R2[s_H, s_R1, s_T2, :] = state_probabilities(update_risk_distribution(s_R1, s_T2), s_T2, s_H, s_R1) -end -push!(C, ChanceNode(R2, I_R2)) -push!(X, Probabilities(R2, X_R2)) - - -I_TD = [R2] -push!(D, DecisionNode(TD, I_TD)) - - -I_TC = [T1, T2] -Y_TC = zeros(S[I_TC]...) cost_TRS = -0.0034645 cost_GRS = -0.004 -cost_forbidden = 0 #the cost of forbidden test combinations is negligible -Y_TC[1 , 1] = cost_forbidden -Y_TC[1 , 2] = cost_TRS + cost_GRS -Y_TC[1, 3] = cost_TRS -Y_TC[2, 1] = cost_GRS + cost_TRS -Y_TC[2, 2] = cost_forbidden -Y_TC[2, 3] = cost_GRS -Y_TC[3, 1] = cost_TRS -Y_TC[3, 2] = cost_GRS -Y_TC[3, 3] = 0 -push!(V, ValueNode(TC, I_TC)) -push!(Y, Consequences(TC, Y_TC)) - - -I_HB = [H, TD] -Y_HB = zeros(S[I_HB]...) -Y_HB[1 , 1] = 6.89713671259061 # sick & treat -Y_HB[1 , 2] = 6.65436854256236 # sick & don't treat -Y_HB[2, 1] = 7.64528451705134 # healthy & treat -Y_HB[2, 2] = 7.70088349200034 # healthy & don't treat -push!(V, ValueNode(HB, I_HB)) -push!(Y, Consequences(HB, Y_HB)) - - -@info("Validate influence diagram.") -validate_influence_diagram(S, C, D, V) -sort!.((C, D, V, X, Y), by = x -> x.j) - -P = DefaultPathProbability(C, X) -U = DefaultPathUtility(V, Y) +forbidden = 0 #the cost of forbidden test combinations is negligible +Y_TC = UtilityMatrix(diagram, "TC") +set_utility!(Y_TC, ["TRS", "TRS"], forbidden) +set_utility!(Y_TC, ["TRS", "GRS"], cost_TRS + cost_GRS) +set_utility!(Y_TC, ["TRS", "no test"], cost_TRS) +set_utility!(Y_TC, ["GRS", "TRS"], cost_TRS + cost_GRS) +set_utility!(Y_TC, ["GRS", "GRS"], forbidden) +set_utility!(Y_TC, ["GRS", "no test"], cost_GRS) +set_utility!(Y_TC, ["no test", "TRS"], cost_TRS) +set_utility!(Y_TC, ["no test", "GRS"], cost_GRS) +set_utility!(Y_TC, ["no test", "no test"], 0) +add_utilities!(diagram, "TC", Y_TC) + +Y_HB = UtilityMatrix(diagram, "HB") +set_utility!(Y_HB, ["CHD", "treatment"], 6.89713671259061) +set_utility!(Y_HB, ["CHD", "no treatment"], 6.65436854256236 ) +set_utility!(Y_HB, ["no CHD", "treatment"], 7.64528451705134) +set_utility!(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, S, D) +z = DecisionVariables(model, diagram) # Defining forbidden paths to include all those where a test is repeated twice -forbidden_tests = ForbiddenPath[([T1,T2], Set([(1,1),(2,2),(3,1), (3,2)]))] +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, z, S, P; fixed = Dict(1 => chosen_risk_level), forbidden_paths = forbidden_tests, probability_cut=false) +x_s = PathCompatibilityVariables(model, diagram, z; fixed = fixed_R0, forbidden_paths = [forbidden_tests], probability_cut=false) -EV = expected_value(model, x_s, U, P, probability_scale_factor = scale_factor) +EV = expected_value(model, diagram, x_s, probability_scale_factor = scale_factor) @objective(model, Max, EV) @info("Starting the optimization process.") @@ -268,25 +194,24 @@ optimizer = optimizer_with_attributes( "MIPGap" => 1e-6, ) set_optimizer(model, optimizer) +GC.enable(false) optimize!(model) +GC.enable(true) @info("Extracting results.") Z = DecisionStrategy(z) +S_probabilities = StateProbabilities(diagram, Z) +U_distribution = UtilityDistribution(diagram, Z) @info("Printing decision strategy using tailor made function:") -sprobs = StateProbabilities(S, P, Z) -analysing_results(Z, sprobs) +print_decision_strategy(diagram, Z, S_probabilities) @info("Printing state probabilities:") # Here we can see that the probability of having a CHD event is exactly that of the chosen risk level -print_state_probabilities(sprobs, [R0, R1, R2]) - -@info("Computing utility distribution.") -udist = UtilityDistribution(S, P, U, Z) +print_state_probabilities(diagram, S_probabilities, ["R0", "R1", "R2"]) @info("Printing utility distribution.") -print_utility_distribution(udist) +print_utility_distribution(U_distribution) -@info("Printing statistics") -print_statistics(udist) +print_statistics(U_distribution) diff --git a/examples/contingent-portfolio-programming.jl b/examples/contingent-portfolio-programming.jl index 8d4baa2b..26f39b19 100644 --- a/examples/contingent-portfolio-programming.jl +++ b/examples/contingent-portfolio-programming.jl @@ -4,43 +4,31 @@ using DecisionProgramming Random.seed!(42) -const dᴾ = 1 # Decision node: range for number of patents -const cᵀ = 2 # Chance node: technical competitiveness -const dᴬ = 3 # Decision node: range for number of applications -const cᴹ = 4 # Chance node: market share -const DP_states = ["0-3 patents", "3-6 patents", "6-9 patents"] -const CT_states = ["low", "medium", "high"] -const DA_states = ["0-5 applications", "5-10 applications", "10-15 applications"] -const CM_states = ["low", "medium", "high"] - -S = States([ - (length(DP_states), [dᴾ]), - (length(CT_states), [cᵀ]), - (length(DA_states), [dᴬ]), - (length(CM_states), [cᴹ]), -]) -C = Vector{ChanceNode}() -D = Vector{DecisionNode}() -V = Vector{ValueNode}() -X = Vector{Probabilities}() -Y = Vector{Consequences}() - -I_DP = Vector{Node}() -push!(D, DecisionNode(dᴾ, I_DP)) - -I_CT = [dᴾ] -X_CT = zeros(S[dᴾ], S[cᵀ]) +function num_states(diagram::InfluenceDiagram, node::Name) + idx = findfirst(isequal(node), diagram.Names) + if isnothing(idx) + throw(DomainError("Name $node not found in the diagram.")) + end + return diagram.S[idx] +end + +@info("Creating the influence diagram.") +diagram = InfluenceDiagram() + +add_node!(diagram, DecisionNode("DP", [], ["0-3 patents", "3-6 patents", "6-9 patents"])) +add_node!(diagram, ChanceNode("CT", ["DP"], ["low", "medium", "high"])) +add_node!(diagram, DecisionNode("DA", ["DP", "CT"], ["0-5 applications", "5-10 applications", "10-15 applications"])) +add_node!(diagram, ChanceNode("CM", ["CT", "DA"], ["low", "medium", "high"])) + +generate_arcs!(diagram) + +X_CT = ProbabilityMatrix(diagram, "CT") X_CT[1, :] = [1/2, 1/3, 1/6] X_CT[2, :] = [1/3, 1/3, 1/3] X_CT[3, :] = [1/6, 1/3, 1/2] -push!(C, ChanceNode(cᵀ, I_CT)) -push!(X, Probabilities(cᵀ, X_CT)) +add_probabilities!(diagram, "CT", X_CT) -I_DA = [dᴾ, cᵀ] -push!(D, DecisionNode(dᴬ, I_DA)) - -I_CM = [cᵀ, dᴬ] -X_CM = zeros(S[cᵀ], S[dᴬ], S[cᴹ]) +X_CM = ProbabilityMatrix(diagram, "CM") X_CM[1, 1, :] = [2/3, 1/4, 1/12] X_CM[1, 2, :] = [1/2, 1/3, 1/6] X_CM[1, 3, :] = [1/3, 1/3, 1/3] @@ -50,23 +38,15 @@ X_CM[2, 3, :] = [1/6, 1/3, 1/2] X_CM[3, 1, :] = [1/3, 1/3, 1/3] X_CM[3, 2, :] = [1/6, 1/3, 1/2] X_CM[3, 3, :] = [1/12, 1/4, 2/3] -push!(C, ChanceNode(cᴹ, I_CM)) -push!(X, Probabilities(cᴹ, X_CM)) +add_probabilities!(diagram, "CM", X_CM) -# Dummy value node -push!(V, ValueNode(5, [cᴹ])) -push!(Y, Consequences(5, zeros(S[cᴹ]))) +generate_diagram!(diagram, default_utility=false) -@info("Validate influence diagram.") -validate_influence_diagram(S, C, D, V) -sort!.((C, D, V, X, Y), by = x -> x.j) -@info("Creating path probability.") -P = DefaultPathProbability(C, X) -@info("Defining DecisionModel") +@info("Creating the decision model.") model = Model() -z = DecisionVariables(model, S, D) +z = DecisionVariables(model, diagram) @info("Creating problem specific constraints and expressions") @@ -85,12 +65,12 @@ O_t = rand(1:3,n_T) # number of patents for each tech project I_a = rand(n_T)*2 # costs of application projects O_a = rand(2:4,n_T) # number of applications for each appl. project -V_A = rand(S[cᴹ], n_A).+0.5 # Value of an application +V_A = rand(num_states(diagram, "CM"), n_A).+0.5 # Value of an application V_A[1, :] .+= -0.5 # Low market share: less value V_A[3, :] .+= 0.5 # High market share: more value -x_T = variables(model, [S[dᴾ]...,n_T]; binary=true) -x_A = variables(model, [S[dᴾ]...,S[cᵀ]...,S[dᴬ]..., n_A]; binary=true) +x_T = variables(model, [num_states(diagram, "DP"), n_T]; binary=true) +x_A = variables(model, [num_states(diagram, "DP"), num_states(diagram, "CT"), num_states(diagram, "DA"), n_A]; binary=true) M = 20 # a large constant ε = 0.5*minimum([O_t O_a]) # a helper variable, allows using ≤ instead of < in constraints (28b) and (29b) @@ -125,10 +105,11 @@ z_dA = z.z[2] @constraint(model, [i=1:3, j=1:3, k=1:3], x_A[i,j,k,2] <= x_T[i,2]) @info("Creating model objective.") -patent_investment_cost = @expression(model, [i=1:S[1]], sum(x_T[i, t] * I_t[t] for t in 1:n_T)) -application_investment_cost = @expression(model, [i=1:S[1], j=1:S[2], k=1:S[3]], sum(x_A[i, j, k, a] * I_a[a] for a in 1:n_A)) -application_value = @expression(model, [i=1:S[1], j=1:S[2], k=1:S[3], l=1:S[4]], sum(x_A[i, j, k, a] * V_A[l, a] for a in 1:n_A)) -@objective(model, Max, sum( sum( P((i,j,k,l)) * (application_value[i,j,k,l] - application_investment_cost[i,j,k]) for j in 1:S[2], k in 1:S[3], l in 1:S[4] ) - patent_investment_cost[i] for i in 1:S[1] )) +patent_investment_cost = @expression(model, [i=1:diagram.S[1]], sum(x_T[i, t] * I_t[t] for t in 1:n_T)) +application_investment_cost = @expression(model, [i=1:diagram.S[1], j=1:diagram.S[2], k=1:diagram.S[3]], sum(x_A[i, j, k, a] * I_a[a] for a in 1:n_A)) +application_value = @expression(model, [i=1:diagram.S[1], j=1:diagram.S[2], k=1:diagram.S[3], l=1:diagram.S[4]], sum(x_A[i, j, k, a] * V_A[l, a] for a in 1:n_A)) +@objective(model, Max, sum( sum( diagram.P((i,j,k,l)) * (application_value[i,j,k,l] - application_investment_cost[i,j,k]) for j in State(1):diagram.S[2], k in State(1):diagram.S[3], l in State(1):diagram.S[4] ) - patent_investment_cost[i] for i in State(1):diagram.S[1] )) + @info("Starting the optimization process.") optimizer = optimizer_with_attributes( @@ -141,29 +122,30 @@ optimize!(model) @info("Extracting results.") Z = DecisionStrategy(z) +S_probabilities = StateProbabilities(diagram, Z) @info("Printing decision strategy:") -print_decision_strategy(S, Z) +print_decision_strategy(diagram, Z, S_probabilities) @info("Extracting path utilities") struct PathUtility <: AbstractPathUtility data::Array{AffExpr} end -Base.getindex(U::PathUtility, i::Int) = getindex(U.data, i) -Base.getindex(U::PathUtility, I::Vararg{Int,N}) where N = getindex(U.data, I...) +Base.getindex(U::PathUtility, i::State) = getindex(U.data, i) +Base.getindex(U::PathUtility, I::Vararg{State,N}) where N = getindex(U.data, I...) (U::PathUtility)(s::Path) = value.(U[s...]) path_utility = [@expression(model, sum(x_A[s[1:3]..., a] * (V_A[s[4], a] - I_a[a]) for a in 1:n_A) - - sum(x_T[s[1], t] * I_t[t] for t in 1:n_T)) for s in paths(S)] -U = PathUtility(path_utility) + sum(x_T[s[1], t] * I_t[t] for t in 1:n_T)) for s in paths(diagram.S)] +diagram.U = PathUtility(path_utility) @info("Computing utility distribution.") -udist = UtilityDistribution(S, P, U, Z) +U_distribution = UtilityDistribution(diagram, Z) @info("Printing utility distribution.") -print_utility_distribution(udist) +print_utility_distribution(U_distribution) @info("Printing statistics") -print_statistics(udist) +print_statistics(U_distribution) diff --git a/examples/figures/simple-id.drawio b/examples/figures/simple-id.drawio index f824a7f8..3303c031 100644 --- a/examples/figures/simple-id.drawio +++ b/examples/figures/simple-id.drawio @@ -1 +1 @@ -7Vhdb5swFP01SNvDJDDYyR6XtFsndQ9THvbsgPnoDKaOKUl//WxsEgghoWoEmTaCgn18ba7v8bm2sNxluv3GcR7/YAGhFrCDreXeWQCg2Uz+K2CnAQ/YGoh4EmjIOQCr5JUYsDYrkoBsWoaCMSqSvA36LMuIL1oY5pyVbbOQ0fZbcxyRDrDyMe2iv5JAxBqdQ/uAP5Akius3O7ZpSXFtbIbYxDhgpYYqG/fecpecMaFL6XZJqIpdHRc90Nee1r1jnGRiSIfd489X9FwG5BU+wALk4ep7/Omz8U3s6gmTQM7fVDOWyceCsyILiBrGljXGRcwilmH6yFguQUeCT0SInWEPF4JJKBYpNa0bwdnvffzkzBdhQumSUcarl7oBJPPA21s2WuZg7SIkW7Sjyrve+deBZgX3yZlJ1+sI84iIM3Zwz5Jc3YSlRPCd7McJxSJ5afuBzTqL9nYHKmTBsPEGZsy4L5gW5k0W8OStQAsu1a0Ls4U0VRVQ1e60WZdXSqVmFJ1lnAiyynEVo1Kqtk0W3uRaSGGyVaQfcxVC9TvFFaou1YNlooHr6+Q6yAlPZLwIV/4kWWRgM3fCBdmep7tLj+ngIiPEOvF4pl42ZGyguKHgGrs6obMOJe/W1duZ6TJwRV2BgbrqIW4cXc1vI+MN4GoCZrwpmQEnMh6iQsUox1mLM/RcqG1z4euwfVE+RusP0FE5sUqGreLHQw9ZivRTZUlnaDLVfshpaVfqMY5W0jsz64BdcMwMCrwby6Buz5YIe7Y8HrN0XWwub3nH2gxD4PuneAjQGsGxeUD2jfFQn82nzqHjnhq9v+HU6F3Ioc3k5143+f3rB0yIppal0+VkkqPN5fR5RVnCgbJ0p5QlHP9o4/0/2vRr15uPtqXK6uEbT9XW+FDm3v8B \ No newline at end of file +5Vhbb5swFP41SNvDJG4m2WNDunZS9zBF2p4NmMtqMHVMgPz62dgkEHJruiZMJZGwv3N8O5/PZwvNctPqgcI8/kEChDVTDyrNmmumaYCJyV8CqRViGAqJaBIobAsskjVSoK7QIgnQsufICMEsyfugT7IM+ayHQUpJ2XcLCe6PmsMIDYCFD/EQ/Z0ELJboFOhb/BElUdyObOjKksLWWXWxjGFASgk1Pta9ZrmUECZLaeUiLKLXxkV29O2AdTMxijJ2ToP66efaeSkDtAaPoDDzcPE9/vJVzY3V7YJRwNevqhnJ+GtGSZEFSHSj8xqhLCYRySB+IiTnoMHBP4ixWrEHC0Y4FLMUK+uSUfK8iR9f+SxMMHYJJrQZ1AoAmgb2xrNjmZqe5TjcMlxuG1dSUB8dWWO7bSCNEDviB6SfCEBnABXMB0RSxGjNHSjCkCWr/gaBap9FG78tFbyg2HgFM6rfFcSFGkkzbf7ngdA14Iq/LExmlWaKSt3U5tJtyCvGPGcEnWWcMLTIYRO0kudtnyy4zGUihUklSN/lKgTit48rp3lEC5KxDi6fvfsgRzTh8UJUzCfJIgUfJHuFKEPVUXqU1XJUIirlsW1VLztprKC4k8Et9s8JnQwoeXNevZ6ZIQOX55V5Zl4Zo8qr6TgU7wyu3p8Ze1TMmHsUz8FMxCiHWY8z56UQx+bMl2G7E3OMvE/AEJrYiGGv+Hnbgpci+RYqeTcUUygbeX0xlfPgy5JTafvY2UlvVNYzTsFrKqhpj0xBrQNH4q8DRx6NSeoVy9NH3m5uhqHp+/t4CBzPAdfmwdFHxkN7N7+1hr7rrdH+L2+N9gkN7YqfOxS/lRS/8hLx++gXTODcOi2NISc3udqcls/L0xKcmZbWqNISXP9qMx9m9zNvhJu2cmyPbpp8+MuOPb3aIcur268+ja3z8cy6/ws= \ No newline at end of file diff --git a/examples/figures/simple-id.svg b/examples/figures/simple-id.svg index 0b5eba7d..2fc5a9ce 100644 --- a/examples/figures/simple-id.svg +++ b/examples/figures/simple-id.svg @@ -1,3 +1,3 @@ -
 2 \\ \{1, 2\}
 1 \\ \{1, 2\}
 5
 3 \\ \{1, 2\}
 4 \\ \{1, 2\}Viewer does not support full SVG 1.1 \ No newline at end of file +
 B \\ \{x, y\}
 A \\ \{a, b\}
 V
 C \\ \{v, w\}
D \\ \{k,l \...Viewer does not support full SVG 1.1 \ No newline at end of file diff --git a/examples/n_monitoring.jl b/examples/n_monitoring.jl index 7ec29b3e..64460e1e 100644 --- a/examples/n_monitoring.jl +++ b/examples/n_monitoring.jl @@ -5,97 +5,69 @@ using DecisionProgramming Random.seed!(13) const N = 4 -const L = [1] -const R_k = [k + 1 for k in 1:N] -const A_k = [(N + 1) + k for k in 1:N] -const F = [2*N + 2] -const T = [2*N + 3] -const L_states = ["high", "low"] -const R_k_states = ["high", "low"] -const A_k_states = ["yes", "no"] -const F_states = ["failure", "success"] const c_k = rand(N) const b = 0.03 fortification(k, a) = [c_k[k], 0][a] @info("Creating the influence diagram.") -S = States([ - (length(L_states), L), - (length(R_k_states), R_k), - (length(A_k_states), A_k), - (length(F_states), F) -]) -C = Vector{ChanceNode}() -D = Vector{DecisionNode}() -V = Vector{ValueNode}() -X = Vector{Probabilities}() -Y = Vector{Consequences}() - -for j in L - I_j = Vector{Node}() - X_j = zeros(S[I_j]..., S[j]) - X_j[1] = rand() - X_j[2] = 1.0 - X_j[1] - push!(C, ChanceNode(j, I_j)) - push!(X, Probabilities(j, X_j)) -end +diagram = InfluenceDiagram() + +add_node!(diagram, ChanceNode("L", [], ["high", "low"])) -for j in R_k - I_j = L - x, y = rand(2) - X_j = zeros(S[I_j]..., S[j]) - X_j[1, 1] = max(x, 1-x) - X_j[1, 2] = 1.0 - X_j[1, 1] - X_j[2, 2] = max(y, 1-y) - X_j[2, 1] = 1.0 - X_j[2, 2] - push!(C, ChanceNode(j, I_j)) - push!(X, Probabilities(j, X_j)) +for i in 1:N + add_node!(diagram, ChanceNode("R$i", ["L"], ["high", "low"])) + add_node!(diagram, DecisionNode("A$i", ["R$i"], ["yes", "no"])) end -for (i, j) in zip(R_k, A_k) - I_j = [i] - push!(D, DecisionNode(j, I_j)) +add_node!(diagram, ChanceNode("F", ["L", ["A$i" for i in 1:N]...], ["failure", "success"])) + +add_node!(diagram, ValueNode("T", ["F", ["A$i" for i in 1:N]...])) + +generate_arcs!(diagram) + +X_L = [rand(), 0] +X_L[2] = 1.0 - X_L[1] +add_probabilities!(diagram, "L", X_L) + +for i in 1:N + x_R, y_R = rand(2) + X_R = ProbabilityMatrix(diagram, "R$i") + set_probability!(X_R, ["high", "high"], max(x_R, 1-x_R)) + set_probability!(X_R, ["high", "low"], 1 - max(x_R, 1-x_R)) + set_probability!(X_R, ["low", "low"], max(y_R, 1-y_R)) + set_probability!(X_R, ["low", "high"], 1-max(y_R, 1-y_R)) + add_probabilities!(diagram, "R$i", X_R) end -for j in F - I_j = L ∪ A_k - x, y = rand(2) - X_j = zeros(S[I_j]..., S[j]) - for s in paths(S[A_k]) - d = exp(b * sum(fortification(k, a) for (k, a) in enumerate(s))) - X_j[1, s..., 1] = max(x, 1-x) / d - X_j[1, s..., 2] = 1.0 - X_j[1, s..., 1] - X_j[2, s..., 1] = min(y, 1-y) / d - X_j[2, s..., 2] = 1.0 - X_j[2, s..., 1] - end - push!(C, ChanceNode(j, I_j)) - push!(X, Probabilities(j, X_j)) + +X_F = ProbabilityMatrix(diagram, "F") +x_F, y_F = rand(2) +for s in paths([State(2) for i in 1:N]) + denominator = exp(b * sum(fortification(k, a) for (k, a) in enumerate(s))) + s1 = [s...] + X_F[1, s1..., 1] = max(x_F, 1-x_F) / denominator + X_F[1, s..., 2] = 1.0 - X_F[1, s..., 1] + X_F[2, s..., 1] = min(y_F, 1-y_F) / denominator + X_F[2, s..., 2] = 1.0 - X_F[2, s..., 1] end +add_probabilities!(diagram, "F", X_F) -for j in T - I_j = A_k ∪ F - Y_j = zeros(S[I_j]...) - for s in paths(S[A_k]) - cost = sum(-fortification(k, a) for (k, a) in enumerate(s)) - Y_j[s..., 1] = cost + 0 - Y_j[s..., 2] = cost + 100 - end - push!(V, ValueNode(j, I_j)) - push!(Y, Consequences(j, Y_j)) + +Y_T = UtilityMatrix(diagram, "T") +for s in paths([State(2) for i in 1:N]) + cost = sum(-fortification(k, a) for (k, a) in enumerate(s)) + Y_T[1, s...] = 0 + cost + Y_T[2, s...] = 100 + cost end +add_utilities!(diagram, "T", Y_T) -validate_influence_diagram(S, C, D, V) -sort!.((C, D, V, X, Y), by = x -> x.j) +generate_diagram!(diagram, positive_path_utility=true) -P = DefaultPathProbability(C, X) -U = DefaultPathUtility(V, Y) -@info("Creating the decision model.") -U⁺ = PositivePathUtility(S, U) model = Model() -z = DecisionVariables(model, S, D) -x_s = PathCompatibilityVariables(model, z, S, P, probability_cut = false) -EV = expected_value(model, x_s, U⁺, P) +z = DecisionVariables(model, diagram) +x_s = PathCompatibilityVariables(model, diagram, z, probability_cut = false) +EV = expected_value(model, diagram, x_s) @objective(model, Max, EV) @info("Starting the optimization process.") @@ -108,22 +80,20 @@ optimize!(model) @info("Extracting results.") Z = DecisionStrategy(z) +U_distribution = UtilityDistribution(diagram, Z) +S_probabilities = StateProbabilities(diagram, Z) @info("Printing decision strategy:") -print_decision_strategy(S, Z) - -@info("Printing state probabilities:") -sprobs = StateProbabilities(S, P, Z) -print_state_probabilities(sprobs, L) -print_state_probabilities(sprobs, R_k) -print_state_probabilities(sprobs, A_k) -print_state_probabilities(sprobs, F) +print_decision_strategy(diagram, Z, S_probabilities) -@info("Computing utility distribution.") -udist = UtilityDistribution(S, P, U, Z) +@info("State probabilities:") +print_state_probabilities(diagram, S_probabilities, ["L"]) +print_state_probabilities(diagram, S_probabilities, [["R$i" for i in 1:N]...]) +print_state_probabilities(diagram, S_probabilities, [["A$i" for i in 1:N]...]) +print_state_probabilities(diagram, S_probabilities, ["F"]) @info("Printing utility distribution.") -print_utility_distribution(udist) +print_utility_distribution(U_distribution) @info("Printing statistics") -print_statistics(udist) +print_statistics(U_distribution) diff --git a/examples/pig_breeding.jl b/examples/pig_breeding.jl index bf26415c..ca6629c1 100644 --- a/examples/pig_breeding.jl +++ b/examples/pig_breeding.jl @@ -3,97 +3,61 @@ using JuMP, Gurobi using DecisionProgramming const N = 4 -const health = [3*k - 2 for k in 1:N] -const test = [3*k - 1 for k in 1:(N-1)] -const treat = [3*k for k in 1:(N-1)] -const cost = [(3*N - 2) + k for k in 1:(N-1)] -const price = [(3*N - 2) + N] -const health_states = ["ill", "healthy"] -const test_states = ["positive", "negative"] -const treat_states = ["treat", "pass"] @info("Creating the influence diagram.") -S = States([ - (length(health_states), health), - (length(test_states), test), - (length(treat_states), treat), -]) -C = Vector{ChanceNode}() -D = Vector{DecisionNode}() -V = Vector{ValueNode}() -X = Vector{Probabilities}() -Y = Vector{Consequences}() - -for j in health[[1]] - I_j = Vector{Node}() - X_j = zeros(S[I_j]..., S[j]) - X_j[1] = 0.1 - X_j[2] = 1.0 - X_j[1] - push!(C, ChanceNode(j, I_j)) - push!(X, Probabilities(j, X_j)) +diagram = InfluenceDiagram() + +add_node!(diagram, ChanceNode("H1", [], ["ill", "healthy"])) +for i in 1:N-1 + # Testing result + add_node!(diagram, ChanceNode("T$i", ["H$i"], ["positive", "negative"])) + # Decision to treat + add_node!(diagram, DecisionNode("D$i", ["T$i"], ["treat", "pass"])) + # Cost of treatment + add_node!(diagram, ValueNode("C$i", ["D$i"])) + # Health of next period + add_node!(diagram, ChanceNode("H$(i+1)", ["H$(i)", "D$(i)"], ["ill", "healthy"])) end - -for (i, k, j) in zip(health[1:end-1], treat, health[2:end]) - I_j = [i, k] - X_j = zeros(S[I_j]..., S[j]) - X_j[2, 2, 1] = 0.2 - X_j[2, 2, 2] = 1.0 - X_j[2, 2, 1] - X_j[2, 1, 1] = 0.1 - X_j[2, 1, 2] = 1.0 - X_j[2, 1, 1] - X_j[1, 2, 1] = 0.9 - X_j[1, 2, 2] = 1.0 - X_j[1, 2, 1] - X_j[1, 1, 1] = 0.5 - X_j[1, 1, 2] = 1.0 - X_j[1, 1, 1] - push!(C, ChanceNode(j, I_j)) - push!(X, Probabilities(j, X_j)) -end - -for (i, j) in zip(health, test) - I_j = [i] - X_j = zeros(S[I_j]..., S[j]) - X_j[1, 1] = 0.8 - X_j[1, 2] = 1.0 - X_j[1, 1] - X_j[2, 2] = 0.9 - X_j[2, 1] = 1.0 - X_j[2, 2] - push!(C, ChanceNode(j, I_j)) - push!(X, Probabilities(j, X_j)) -end - -for (i, j) in zip(test, treat) - I_j = [i] - push!(D, DecisionNode(j, I_j)) +add_node!(diagram, ValueNode("MP", ["H$N"])) + +generate_arcs!(diagram) + +# Add probabilities for node H1 +add_probabilities!(diagram, "H1", [0.1, 0.9]) + +# Declare proability matrix for health nodes H_2, ... H_N-1, which have identical information sets and states +X_H = ProbabilityMatrix(diagram, "H2") +set_probability!(X_H, ["healthy", "pass", :], [0.2, 0.8]) +set_probability!(X_H, ["healthy", "treat", :], [0.1, 0.9]) +set_probability!(X_H, ["ill", "pass", :], [0.9, 0.1]) +set_probability!(X_H, ["ill", "treat", :], [0.5, 0.5]) + +# Declare proability matrix for test result nodes T_1...T_N +X_T = ProbabilityMatrix(diagram, "T1") +set_probability!(X_T, ["ill", "positive"], 0.8) +set_probability!(X_T, ["ill", "negative"], 0.2) +set_probability!(X_T, ["healthy", "negative"], 0.9) +set_probability!(X_T, ["healthy", "positive"], 0.1) + +for i in 1:N-1 + add_probabilities!(diagram, "T$i", X_T) + add_probabilities!(diagram, "H$(i+1)", X_H) end -for (i, j) in zip(treat, cost) - I_j = [i] - Y_j = zeros(S[I_j]...) - Y_j[1] = -100 - Y_j[2] = 0 - push!(V, ValueNode(j, I_j)) - push!(Y, Consequences(j, Y_j)) +for i in 1:N-1 + add_utilities!(diagram, "C$i", [-100.0, 0.0]) end -for (i, j) in zip(health[end], price) - I_j = [i] - Y_j = zeros(S[I_j]...) - Y_j[1] = 300 - Y_j[2] = 1000 - push!(V, ValueNode(j, I_j)) - push!(Y, Consequences(j, Y_j)) -end +add_utilities!(diagram, "MP", [300.0, 1000.0]) -validate_influence_diagram(S, C, D, V) -sort!.((C, D, V, X, Y), by = x -> x.j) +generate_diagram!(diagram, positive_path_utility = true) -P = DefaultPathProbability(C, X) -U = DefaultPathUtility(V, Y) @info("Creating the decision model.") -U⁺ = PositivePathUtility(S, U) model = Model() -z = DecisionVariables(model, S, D) -x_s = PathCompatibilityVariables(model, z, S, P, probability_cut = false) -EV = expected_value(model, x_s, U⁺, P) +z = DecisionVariables(model, diagram) +x_s = PathCompatibilityVariables(model, diagram, z, probability_cut = false) +EV = expected_value(model, diagram, x_s) @objective(model, Max, EV) @info("Starting the optimization process.") @@ -106,30 +70,28 @@ optimize!(model) @info("Extracting results.") Z = DecisionStrategy(z) +S_probabilities = StateProbabilities(diagram, Z) +U_distribution = UtilityDistribution(diagram, Z) + @info("Printing decision strategy:") -print_decision_strategy(S, Z) +print_decision_strategy(diagram, Z, S_probabilities) + +@info("Printing utility distribution.") +print_utility_distribution(U_distribution) + +@info("Printing statistics") +print_statistics(U_distribution) @info("State probabilities:") -sprobs = StateProbabilities(S, P, Z) -print_state_probabilities(sprobs, health) -print_state_probabilities(sprobs, test) -print_state_probabilities(sprobs, treat) +print_state_probabilities(diagram, S_probabilities, [["H$i" for i in 1:N]...]) +print_state_probabilities(diagram, S_probabilities, [["T$i" for i in 1:N-1]...]) +print_state_probabilities(diagram, S_probabilities, [["D$i" for i in 1:N-1]...]) @info("Conditional state probabilities") -node = 1 -for state in 1:2 - sprobs2 = StateProbabilities(S, P, Z, node, state, sprobs) - print_state_probabilities(sprobs2, health) - print_state_probabilities(sprobs2, test) - print_state_probabilities(sprobs2, treat) +for state in ["ill", "healthy"] + S_probabilities2 = StateProbabilities(diagram, Z, "H1", state, S_probabilities) + print_state_probabilities(diagram, S_probabilities2, [["H$i" for i in 1:N]...]) + print_state_probabilities(diagram, S_probabilities2, [["T$i" for i in 1:N-1]...]) + print_state_probabilities(diagram, S_probabilities2, [["D$i" for i in 1:N-1]...]) end - -@info("Computing utility distribution.") -udist = UtilityDistribution(S, P, U, Z) - -@info("Printing utility distribution.") -print_utility_distribution(udist) - -@info("Printing statistics") -print_statistics(udist) diff --git a/examples/used_car_buyer.jl b/examples/used_car_buyer.jl index 5ce895ef..ec188d33 100644 --- a/examples/used_car_buyer.jl +++ b/examples/used_car_buyer.jl @@ -1,76 +1,62 @@ + using Logging using JuMP, Gurobi using DecisionProgramming -const O = 1 # Chance node: lemon or peach -const T = 2 # Decision node: pay stranger for advice -const R = 3 # Chance node: observation of state of the car -const A = 4 # Decision node: purchase alternative -const O_states = ["lemon", "peach"] -const T_states = ["no test", "test"] -const R_states = ["no test", "lemon", "peach"] -const A_states = ["buy without guarantee", "buy with guarantee", "don't buy"] - @info("Creating the influence diagram.") -S = States([ - (length(O_states), [O]), - (length(T_states), [T]), - (length(R_states), [R]), - (length(A_states), [A]), -]) -C = Vector{ChanceNode}() -D = Vector{DecisionNode}() -V = Vector{ValueNode}() -X = Vector{Probabilities}() -Y = Vector{Consequences}() - -I_O = Vector{Node}() -X_O = [0.2, 0.8] -push!(C, ChanceNode(O, I_O)) -push!(X, Probabilities(O, X_O)) - -I_T = Vector{Node}() -push!(D, DecisionNode(T, I_T)) - -I_R = [O, T] -X_R = zeros(S[O], S[T], S[R]) -X_R[1, 1, :] = [1,0,0] -X_R[1, 2, :] = [0,1,0] -X_R[2, 1, :] = [1,0,0] -X_R[2, 2, :] = [0,0,1] -push!(C, ChanceNode(R, I_R)) -push!(X, Probabilities(R, X_R)) - -I_A = [R] -push!(D, DecisionNode(A, I_A)) - -I_V1 = [T] -Y_V1 = [0.0, -25.0] -push!(V, ValueNode(5, I_V1)) -push!(Y, Consequences(5, Y_V1)) - -I_V2 = [A] -Y_V2 = [100.0, 40.0, 0.0] -push!(V, ValueNode(6, I_V2)) -push!(Y, Consequences(6, Y_V2)) - -I_V3 = [O, A] -Y_V3 = [-200.0 0.0 0.0; - -40.0 -20.0 0.0] -push!(V, ValueNode(7, I_V3)) -push!(Y, Consequences(7, Y_V3)) - -validate_influence_diagram(S, C, D, V) -sort!.((C, D, V, X, Y), by = x -> x.j) - -P = DefaultPathProbability(C, X) -U = DefaultPathUtility(V, Y) +diagram = InfluenceDiagram() + +add_node!(diagram, ChanceNode("O", [], ["lemon", "peach"])) +add_node!(diagram, DecisionNode("T", [], ["no test", "test"])) +add_node!(diagram, ChanceNode("R", ["O", "T"], ["no test", "lemon", "peach"])) +add_node!(diagram, DecisionNode("A", ["R"], ["buy without guarantee", "buy with guarantee", "don't buy"])) + +add_node!(diagram, ValueNode("V1", ["T"])) +add_node!(diagram, ValueNode("V2", ["A"])) +add_node!(diagram, ValueNode("V3", ["O", "A"])) + +generate_arcs!(diagram) + +X_O = ProbabilityMatrix(diagram, "O") +set_probability!(X_O, ["peach"], 0.8) +set_probability!(X_O, ["lemon"], 0.2) +add_probabilities!(diagram, "O", X_O) + + +X_R = ProbabilityMatrix(diagram, "R") +set_probability!(X_R, ["lemon", "no test", :], [1,0,0]) +set_probability!(X_R, ["lemon", "test", :], [0,1,0]) +set_probability!(X_R, ["peach", "no test", :], [1,0,0]) +set_probability!(X_R, ["peach", "test", :], [0,0,1]) +add_probabilities!(diagram, "R", X_R) + +Y_V1 = UtilityMatrix(diagram, "V1") +set_utility!(Y_V1, ["test"], -25) +set_utility!(Y_V1, ["no test"], 0) +add_utilities!(diagram, "V1", Y_V1) + + +Y_V2 = UtilityMatrix(diagram, "V2") +set_utility!(Y_V2, ["buy without guarantee"], 100) +set_utility!(Y_V2, ["buy with guarantee"], 40) +set_utility!(Y_V2, ["don't buy"], 0) +add_utilities!(diagram, "V2", Y_V2) + +Y_V3 = UtilityMatrix(diagram, "V3") +set_utility!(Y_V3, ["lemon", "buy without guarantee"], -200) +set_utility!(Y_V3, ["lemon", "buy with guarantee"], 0) +set_utility!(Y_V3, ["lemon", "don't buy"], 0) +set_utility!(Y_V3, ["peach", :], [-40, -20, 0]) +add_utilities!(diagram, "V3", Y_V3) + +generate_diagram!(diagram) + @info("Creating the decision model.") model = Model() -z = DecisionVariables(model, S, D) -x_s = PathCompatibilityVariables(model, z, S, P) -EV = expected_value(model, x_s, U, P) +z = DecisionVariables(model, diagram) +x_s = PathCompatibilityVariables(model, diagram, z) +EV = expected_value(model, diagram, x_s) @objective(model, Max, EV) @info("Starting the optimization process.") @@ -83,15 +69,14 @@ optimize!(model) @info("Extracting results.") Z = DecisionStrategy(z) +S_probabilities = StateProbabilities(diagram, Z) +U_distribution = UtilityDistribution(diagram, Z) @info("Printing decision strategy:") -print_decision_strategy(S, Z) - -@info("Computing utility distribution.") -udist = UtilityDistribution(S, P, U, Z) +print_decision_strategy(diagram, Z, S_probabilities) @info("Printing utility distribution.") -print_utility_distribution(udist) +print_utility_distribution(U_distribution) @info("Printing expected utility.") -print_statistics(udist) +print_statistics(U_distribution) diff --git a/src/DecisionProgramming.jl b/src/DecisionProgramming.jl index af1b4e17..4dcc6ba5 100644 --- a/src/DecisionProgramming.jl +++ b/src/DecisionProgramming.jl @@ -7,6 +7,7 @@ include("analysis.jl") include("printing.jl") export Node, + Name, AbstractNode, ChanceNode, DecisionNode, @@ -15,31 +16,45 @@ export Node, States, Path, paths, + ForbiddenPath, + FixedPath, Probabilities, - Consequences, + Utility, + Utilities, AbstractPathProbability, DefaultPathProbability, AbstractPathUtility, DefaultPathUtility, + validate_influence_diagram, + InfluenceDiagram, + generate_arcs!, + generate_diagram!, + add_node!, + ProbabilityMatrix, + set_probability!, + add_probabilities!, + UtilityMatrix, + set_utility!, + add_utilities!, LocalDecisionStrategy, - DecisionStrategy, - validate_influence_diagram + DecisionStrategy export DecisionVariables, PathCompatibilityVariables, - ForbiddenPath, lazy_probability_cut, - PositivePathUtility, - NegativePathUtility, expected_value, - value_at_risk, conditional_value_at_risk -export random_diagram +export random_diagram!, + random_probabilities!, + random_utilities!, + LocalDecisionStrategy export CompatiblePaths, UtilityDistribution, - StateProbabilities + StateProbabilities, + value_at_risk, + conditional_value_at_risk export print_decision_strategy, print_utility_distribution, diff --git a/src/analysis.jl b/src/analysis.jl index bbb8a4e0..b0540d58 100644 --- a/src/analysis.jl +++ b/src/analysis.jl @@ -1,12 +1,20 @@ +""" + struct CompatiblePaths + S::States + C::Vector{Node} + Z::DecisionStrategy + fixed::FixedPath + end +CompatiblePaths type. +""" struct CompatiblePaths S::States - C::Vector{ChanceNode} + C::Vector{Node} Z::DecisionStrategy - fixed::Dict{Node, State} + fixed::FixedPath function CompatiblePaths(S, C, Z, fixed) - C_j = Set([c.j for c in C]) - if !all(k∈C_j for k in keys(fixed)) + if !all(k∈Set(C) for k in keys(fixed)) throw(DomainError("You can only fix chance states.")) end new(S, C, Z, fixed) @@ -14,9 +22,9 @@ struct CompatiblePaths end """ - CompatiblePaths(S::States, C::Vector{ChanceNode}, Z::DecisionStrategy) + CompatiblePaths(diagram::InfluenceDiagram, Z::DecisionStrategy, fixed::FixedPath=Dict{Node, State}()) -Interface for iterating over paths that are compatible and active given influence diagram and decision strategy. +CompatiblePaths outer construction function. Interface for iterating over paths that are compatible and active given influence diagram and decision strategy. 1) Initialize path `s` of length `n` 2) Fill chance states `s[C]` by generating subpaths `paths(C)` @@ -24,57 +32,58 @@ Interface for iterating over paths that are compatible and active given influenc # Examples ```julia -for s in CompatiblePaths(S, C, Z) +julia> for s in CompatiblePaths(diagram, Z) ... end ``` """ - -function CompatiblePaths(S::States, C::Vector{ChanceNode}, Z::DecisionStrategy) - CompatiblePaths(S, C, Z, Dict{Node, State}()) +function CompatiblePaths(diagram::InfluenceDiagram, Z::DecisionStrategy, fixed::FixedPath=Dict{Node, State}()) + CompatiblePaths(diagram.S, diagram.C, Z, fixed) end -function compatible_path(S::States, C::Vector{ChanceNode}, Z::DecisionStrategy, s_C::Path) - s = Array{Int}(undef, length(S)) +function compatible_path(S::States, C::Vector{Node}, Z::DecisionStrategy, s_C::Path) + s = Array{State}(undef, length(S)) for (c, s_C_j) in zip(C, s_C) - s[c.j] = s_C_j + s[c] = s_C_j end - for (d, Z_j) in zip(Z.D, Z.Z_j) - s[d.j] = Z_j((s[d.I_j]...,)) + for (d, I_d, Z_d) in zip(Z.D, Z.I_d, Z.Z_d) + s[d] = Z_d((s[I_d]...,)) end return (s...,) end -function Base.iterate(a::CompatiblePaths) - C_j = [c.j for c in a.C] - if isempty(a.fixed) - iter = paths(a.S[C_j]) +function Base.iterate(S_Z::CompatiblePaths) + if isempty(S_Z.fixed) + iter = paths(S_Z.S[S_Z.C]) else - ks = sort(collect(keys(a.fixed))) - fixed = Dict{Int, Int}(i => a.fixed[k] for (i, k) in enumerate(ks)) - iter = paths(a.S[C_j], fixed) + ks = sort(collect(keys(S_Z.fixed))) + fixed = Dict{Node, State}(Node(i) => S_Z.fixed[k] for (i, k) in enumerate(S_Z.C) if k in ks) + iter = paths(S_Z.S[S_Z.C], fixed) end next = iterate(iter) if next !== nothing s_C, state = next - return (compatible_path(a.S, a.C, a.Z, s_C), (iter, state)) + return (compatible_path(S_Z.S, S_Z.C, S_Z.Z, s_C), (iter, state)) end end -function Base.iterate(a::CompatiblePaths, gen) +function Base.iterate(S_Z::CompatiblePaths, gen) iter, state = gen next = iterate(iter, state) if next !== nothing s_C, state = next - return (compatible_path(a.S, a.C, a.Z, s_C), (iter, state)) + return (compatible_path(S_Z.S, S_Z.C, S_Z.Z, s_C), (iter, state)) end end Base.eltype(::Type{CompatiblePaths}) = Path -Base.length(a::CompatiblePaths) = prod(a.S[c.j] for c in a.C) +Base.length(S_Z::CompatiblePaths) = prod(S_Z.S[c] for c in S_Z.C) """ - UtilityDistribution + struct UtilityDistribution + u::Vector{Float64} + p::Vector{Float64} + end UtilityDistribution type. @@ -85,23 +94,23 @@ struct UtilityDistribution end """ - UtilityDistribution(S::States, P::AbstractPathProbability, U::AbstractPathUtility, Z::DecisionStrategy) + UtilityDistribution(diagram::InfluenceDiagram, Z::DecisionStrategy) -Constructs the probability mass function for path utilities on paths that are compatible and active. +Construct the probability mass function for path utilities on paths that are compatible with given decision strategy. # Examples ```julia -UtilityDistribution(S, P, U, Z) +julia> UtilityDistribution(diagram, Z) ``` """ -function UtilityDistribution(S::States, P::AbstractPathProbability, U::AbstractPathUtility, Z::DecisionStrategy) +function UtilityDistribution(diagram::InfluenceDiagram, Z::DecisionStrategy) # Extract utilities and probabilities of active paths - S_Z = CompatiblePaths(S, P.C, Z) + S_Z = CompatiblePaths(diagram, Z) utilities = Vector{Float64}(undef, length(S_Z)) probabilities = Vector{Float64}(undef, length(S_Z)) for (i, s) in enumerate(S_Z) - utilities[i] = U(s) - probabilities[i] = P(s) + utilities[i] = diagram.U(s) + probabilities[i] = diagram.P(s) end # Filter zero probabilities @@ -132,84 +141,112 @@ function UtilityDistribution(S::States, P::AbstractPathProbability, U::AbstractP end """ - StateProbabilities + struct StateProbabilities + probs::Dict{Node, Vector{Float64}} + fixed::FixedPath + end StateProbabilities type. """ struct StateProbabilities probs::Dict{Node, Vector{Float64}} - fixed::Dict{Node, State} + fixed::FixedPath end """ - function StateProbabilities(S::States, P::AbstractPathProbability, Z::DecisionStrategy, node::Node, state::State, prev::StateProbabilities) + StateProbabilities(diagram::InfluenceDiagram, Z::DecisionStrategy, node::Node, state::State, prior_probabilities::StateProbabilities) -Associates each node with array of conditional probabilities for each of its states occuring in active paths given fixed states and prior probability. +Associate each node with array of conditional probabilities for each of its states occuring in compatible paths given + fixed states and prior probability. Fix node and state using their indices. # Examples ```julia # Prior probabilities -prev = StateProbabilities(S, P, Z) - -# Select node and fix its state -node = 1 -state = 2 -StateProbabilities(S, P, Z, node, state, prev) +julia> prior_probabilities = StateProbabilities(diagram, Z) +julia> StateProbabilities(diagram, Z, Node(2), State(1), prior_probabilities) ``` """ -function StateProbabilities(S::States, P::AbstractPathProbability, Z::DecisionStrategy, node::Node, state::State, prev::StateProbabilities) - prior = prev.probs[node][state] - fixed = prev.fixed +function StateProbabilities(diagram::InfluenceDiagram, Z::DecisionStrategy, node::Node, state::State, prior_probabilities::StateProbabilities) + prior = prior_probabilities.probs[node][state] + fixed = deepcopy(prior_probabilities.fixed) + push!(fixed, node => state) - probs = Dict(i => zeros(S[i]) for i in 1:length(S)) - for s in CompatiblePaths(S, P.C, Z, fixed), i in 1:length(S) - probs[i][s[i]] += P(s) / prior + probs = Dict(i => zeros(diagram.S[i]) for i in 1:length(diagram.S)) + for s in CompatiblePaths(diagram, Z, fixed), i in 1:length(diagram.S) + probs[i][s[i]] += diagram.P(s) / prior end StateProbabilities(probs, fixed) end """ - function StateProbabilities(S::States, P::AbstractPathProbability, Z::DecisionStrategy) + StateProbabilities(diagram::InfluenceDiagram, Z::DecisionStrategy, node::Name, state::Name, prior_probabilities::StateProbabilities) + +Associate each node with array of conditional probabilities for each of its states occuring in compatible paths given + fixed states and prior probability. Fix node and state using their names. + +# Examples +```julia +# Prior probabilities +julia> prior_probabilities = StateProbabilities(diagram, Z) + +# Select node and fix its state +julia> node = "R" +julia> state = "no test" +julia> StateProbabilities(diagram, Z, node, state, prior_probabilities) +``` +""" +function StateProbabilities(diagram::InfluenceDiagram, Z::DecisionStrategy, node::Name, state::Name, prior_probabilities::StateProbabilities) + node_index = findfirst(j -> j ==node, diagram.Names) + state_index = findfirst(j -> j == state, diagram.States[node_index]) + + return StateProbabilities(diagram, Z, Node(node_index), State(state_index), prior_probabilities) +end + + + + +""" + StateProbabilities(diagram::InfluenceDiagram, Z::DecisionStrategy) -Associates each node with array of probabilities for each of its states occuring in active paths. +Associate each node with array of probabilities for each of its states occuring in compatible paths. # Examples ```julia -StateProbabilities(S, P, Z) +julia> StateProbabilities(diagram, Z) ``` """ -function StateProbabilities(S::States, P::AbstractPathProbability, Z::DecisionStrategy) - probs = Dict(i => zeros(S[i]) for i in 1:length(S)) - for s in CompatiblePaths(S, P.C, Z), i in 1:length(S) - probs[i][s[i]] += P(s) +function StateProbabilities(diagram::InfluenceDiagram, Z::DecisionStrategy) + probs = Dict(i => zeros(diagram.S[i]) for i in 1:length(diagram.S)) + for s in CompatiblePaths(diagram, Z), i in 1:length(diagram.S) + probs[i][s[i]] += diagram.P(s) end StateProbabilities(probs, Dict{Node, State}()) end """ - function value_at_risk(u::Vector{Float64}, p::Vector{Float64}, α::Float64) + value_at_risk(U_distribution::UtilityDistribution, α::Float64) -Value-at-risk. +Calculate value-at-risk. """ -function value_at_risk(u::Vector{Float64}, p::Vector{Float64}, α::Float64) +function value_at_risk(U_distribution::UtilityDistribution, α::Float64) @assert 0 ≤ α ≤ 1 "We should have 0 ≤ α ≤ 1." - i = sortperm(u) - u, p = u[i], p[i] + perm = sortperm(U_distribution.u) + u, p = U_distribution.u[perm], U_distribution.p[perm] index = findfirst(x -> x≥α, cumsum(p)) return if index === nothing; u[end] else u[index] end end """ - function conditional_value_at_risk(u::Vector{Float64}, p::Vector{Float64}, α::Float64) + conditional_value_at_risk(u::Vector{Float64}, p::Vector{Float64}, α::Float64) -Conditional value-at-risk. +Calculate conditional value-at-risk. """ -function conditional_value_at_risk(u::Vector{Float64}, p::Vector{Float64}, α::Float64) - x_α = value_at_risk(u, p, α) +function conditional_value_at_risk(U_distribution::UtilityDistribution, α::Float64) + x_α = value_at_risk(U_distribution, α) if iszero(α) return x_α else - tail = u .≤ x_α - return (sum(u[tail] .* p[tail]) - (sum(p[tail]) - α) * x_α) / α + tail = U_distribution.u .≤ x_α + return (sum(U_distribution.u[tail] .* U_distribution.p[tail]) - (sum(U_distribution.p[tail]) - α) * x_α) / α end end diff --git a/src/decision_model.jl b/src/decision_model.jl index 52d64def..989fdf1e 100644 --- a/src/decision_model.jl +++ b/src/decision_model.jl @@ -1,44 +1,44 @@ using JuMP -function decision_variable(model::Model, S::States, d::DecisionNode, base_name::String="") +function decision_variable(model::Model, S::States, d::Node, I_d::Vector{Node}, base_name::String="") # Create decision variables. - dims = S[[d.I_j; d.j]] - z_j = Array{VariableRef}(undef, dims...) + dims = S[[I_d; d]] + z_d = Array{VariableRef}(undef, dims...) for s in paths(dims) - z_j[s...] = @variable(model, binary=true, base_name=base_name) + z_d[s...] = @variable(model, binary=true, base_name=base_name) end # Constraints to one decision per decision strategy. - for s_I in paths(S[d.I_j]) - @constraint(model, sum(z_j[s_I..., s_j] for s_j in 1:S[d.j]) == 1) + for s_I in paths(S[I_d]) + @constraint(model, sum(z_d[s_I..., s_d] for s_d in 1:S[d]) == 1) end - return z_j + return z_d end struct DecisionVariables - D::Vector{DecisionNode} + D::Vector{Node} + I_d::Vector{Vector{Node}} z::Vector{<:Array{VariableRef}} end """ - DecisionVariables(model::Model, S::States, D::Vector{DecisionNode}; names::Bool=false, name::String="z") + DecisionVariables(model::Model, diagram::InfluenceDiagram; names::Bool=false, name::String="z") Create decision variables and constraints. # Arguments - `model::Model`: JuMP model into which variables are added. -- `S::States`: States structure associated with the influence diagram. -- `D::Vector{DecisionNode}`: Vector containing decicion nodes. +- `diagram::InfluenceDiagram`: Influence diagram structure. - `names::Bool`: Use names or have JuMP variables be anonymous. - `name::String`: Prefix for predefined decision variable naming convention. # Examples ```julia -z = DecisionVariables(model, S, D) +julia> z = DecisionVariables(model, diagram) ``` """ -function DecisionVariables(model::Model, S::States, D::Vector{DecisionNode}; names::Bool=false, name::String="z") - DecisionVariables(D, [decision_variable(model, S, d, (names ? "$(name)_$(d.j)$(s)" : "")) for d in D]) +function DecisionVariables(model::Model, diagram::InfluenceDiagram; names::Bool=false, name::String="z") + DecisionVariables(diagram.D, diagram.I_j[diagram.D], [decision_variable(model, diagram.S, d, I_d, (names ? "$(name)_$(d.j)$(s)" : "")) for (d, I_d) in zip(diagram.D, diagram.I_j[diagram.D])]) end function is_forbidden(s::Path, forbidden_paths::Vector{ForbiddenPath}) @@ -46,7 +46,7 @@ function is_forbidden(s::Path, forbidden_paths::Vector{ForbiddenPath}) end -function path_compatibility_variable(model::Model, z::DecisionVariables, base_name::String="") +function path_compatibility_variable(model::Model, base_name::String="") # Create a path compatiblity variable x = @variable(model, base_name=base_name) @@ -60,6 +60,7 @@ struct PathCompatibilityVariables{N} <: AbstractDict{Path{N}, VariableRef} data::Dict{Path{N}, VariableRef} end +Base.length(x_s::PathCompatibilityVariables) = length(x_s.data) Base.getindex(x_s::PathCompatibilityVariables, key) = getindex(x_s.data, key) Base.get(x_s::PathCompatibilityVariables, key, default) = get(x_s.data, key, default) Base.keys(x_s::PathCompatibilityVariables) = keys(x_s.data) @@ -69,67 +70,63 @@ Base.iterate(x_s::PathCompatibilityVariables) = iterate(x_s.data) Base.iterate(x_s::PathCompatibilityVariables, i) = iterate(x_s.data, i) -function decision_strategy_constraint(model::Model, S::States, d::DecisionNode, D::Vector{DecisionNode}, z::Array{VariableRef}, x_s::PathCompatibilityVariables) +function decision_strategy_constraint(model::Model, S::States, d::Node, I_d::Vector{Node}, D::Vector{Node}, z::Array{VariableRef}, x_s::PathCompatibilityVariables) - # states of nodes in information structure (s_j | s_I(j)) - dims = S[[d.I_j; d.j]] + # states of nodes in information structure (s_d | s_I(d)) + dims = S[[I_d; d]] - # Theoretical upper bound based on number of paths with information structure (s_j | s_I(j)) divided by number of possible decision strategies in other decision nodes - other_decisions = map(d_n -> d_n.j, filter(d_n -> all(d_n.j != i for i in [d.I_j; d.j]), D)) + # Theoretical upper bound based on number of paths with information structure (s_d | s_I(d)) divided by number of possible decision strategies in other decision nodes + other_decisions = filter(j -> all(j != d_set for d_set in [I_d; d]), D) theoretical_ub = prod(S)/prod(dims)/ prod(S[other_decisions]) - # paths that have corresponding path compatibility variable + # paths that have a corresponding path compatibility variable existing_paths = keys(x_s) - for s_j_s_Ij in paths(dims) # iterate through all information states and states of d - # paths with (s_j | s_I(j)) information structure - feasible_paths = filter(s -> s[[d.I_j; d.j]] == s_j_s_Ij, existing_paths) + for s_d_s_Id in paths(dims) # iterate through all information states and states of d + # paths with (s_d | s_I(d)) information structure + feasible_paths = filter(s -> s[[I_d; d]] == s_d_s_Id, existing_paths) - @constraint(model, sum(get(x_s, s, 0) for s in feasible_paths) ≤ z[s_j_s_Ij...] * min(length(feasible_paths), theoretical_ub)) + @constraint(model, sum(get(x_s, s, 0) for s in feasible_paths) ≤ z[s_d_s_Id...] * min(length(feasible_paths), theoretical_ub)) end end """ PathCompatibilityVariables(model::Model, - z::DecisionVariables, - S::States, - P::AbstractPathProbability; + diagram::InfluenceDiagram, + z::DecisionVariables; names::Bool=false, name::String="x", forbidden_paths::Vector{ForbiddenPath}=ForbiddenPath[], - fixed::Dict{Node, State}=Dict{Node, State}(), + fixed::FixedPath=Dict{Node, State}(), probability_cut::Bool=true) Create path compatibility variables and constraints. # Arguments - `model::Model`: JuMP model into which variables are added. +- `diagram::InfluenceDiagram`: Influence diagram structure. - `z::DecisionVariables`: Decision variables from `DecisionVariables` function. -- `S::States`: States structure associated with the influence diagram. -- `P::AbstractPathProbability`: Path probabilities structure for which the function - `P(s)` is defined and returns the path probabilities for path `s`. - `names::Bool`: Use names or have JuMP variables be anonymous. - `name::String`: Prefix for predefined decision variable naming convention. - `forbidden_paths::Vector{ForbiddenPath}`: The forbidden subpath structures. Path compatibility variables will not be generated for paths that include forbidden subpaths. -- `fixed::Dict{Node, State}`: Path compatibility variable will not be generated +- `fixed::FixedPath`: Path compatibility variable will not be generated for paths which do not include these fixed subpaths. - `probability_cut` Includes probability cut constraint in the optimisation model. # Examples ```julia -x_s = PathCompatibilityVariables(model, z, S, P; probability_cut = false) +julia> x_s = PathCompatibilityVariables(model, diagram; probability_cut = false) ``` """ function PathCompatibilityVariables(model::Model, - z::DecisionVariables, - S::States, - P::AbstractPathProbability; + diagram::InfluenceDiagram, + z::DecisionVariables; names::Bool=false, name::String="x", forbidden_paths::Vector{ForbiddenPath}=ForbiddenPath[], - fixed::Dict{Node, State}=Dict{Node, State}(), + fixed::FixedPath=Dict{Node, State}(), probability_cut::Bool=true) if !isempty(forbidden_paths) @@ -137,147 +134,93 @@ function PathCompatibilityVariables(model::Model, end # Create path compatibility variable for each effective path. - N = length(S) + N = length(diagram.S) variables_x_s = Dict{Path{N}, VariableRef}( - s => path_compatibility_variable(model, z, (names ? "$(name)$(s)" : "")) - for s in paths(S, fixed) - if !iszero(P(s)) && !is_forbidden(s, forbidden_paths) + s => path_compatibility_variable(model, (names ? "$(name)$(s)" : "")) + for s in paths(diagram.S, fixed) + if !iszero(diagram.P(s)) && !is_forbidden(s, forbidden_paths) ) x_s = PathCompatibilityVariables{N}(variables_x_s) # Add decision strategy constraints for each decision node for (d, z_d) in zip(z.D, z.z) - decision_strategy_constraint(model, S, d, z.D, z_d, x_s) + decision_strategy_constraint(model, diagram.S, d, diagram.I_j[d], z.D, z_d, x_s) end if probability_cut - @constraint(model, sum(x * P(s) for (s, x) in x_s) == 1.0) + @constraint(model, sum(x * diagram.P(s) for (s, x) in x_s) == 1.0) end x_s end """ - lazy_probability_cut(model::Model, x_s::PathCompatibilityVariables, P::AbstractPathProbability) + lazy_probability_cut(model::Model, diagram::InfluenceDiagram, x_s::PathCompatibilityVariables) -Adds a probability cut to the model as a lazy constraint. +Add a probability cut to the model as a lazy constraint. # Examples ```julia -lazy_probability_cut(model, x_s, P) +julia> lazy_probability_cut(model, diagram, x_s) ``` !!! note Remember to set lazy constraints on in the solver parameters, unless your solver does this automatically. Note that Gurobi does this automatically. """ -function lazy_probability_cut(model::Model, x_s::PathCompatibilityVariables, P::AbstractPathProbability) +function lazy_probability_cut(model::Model, diagram::InfluenceDiagram, x_s::PathCompatibilityVariables) # August 2021: The current implementation of JuMP doesn't allow multiple callback functions of the same type (e.g. lazy) # (see https://github.com/jump-dev/JuMP.jl/issues/2642) # What this means is that if you come up with a new lazy cut, you must replace this # function with a more general function (see discussion and solution in https://github.com/gamma-opt/DecisionProgramming.jl/issues/20) function probability_cut(cb_data) - xsum = sum(callback_value(cb_data, x) * P(s) for (s, x) in x_s) + xsum = sum(callback_value(cb_data, x) * diagram.P(s) for (s, x) in x_s) if !isapprox(xsum, 1.0) - con = @build_constraint(sum(x * P(s) for (s, x) in x_s) == 1.0) + con = @build_constraint(sum(x * diagram.P(s) for (s, x) in x_s) == 1.0) MOI.submit(model, MOI.LazyConstraint(cb_data), con) end end MOI.set(model, MOI.LazyConstraintCallback(), probability_cut) end -# --- Objective Functions --- - -""" - PositivePathUtility(S::States, U::AbstractPathUtility) - -Positive affine transformation of path utility. Always evaluates positive values. - -# Examples -```julia-repl -julia> U⁺ = PositivePathUtility(S, U) -julia> all(U⁺(s) > 0 for s in paths(S)) -true -``` -""" -struct PositivePathUtility <: AbstractPathUtility - U::AbstractPathUtility - min::Float64 - function PositivePathUtility(S::States, U::AbstractPathUtility) - u_min = minimum(U(s) for s in paths(S)) - new(U, u_min) - end -end - -(U::PositivePathUtility)(s::Path) = U.U(s) - U.min + 1 - -""" - NegativePathUtility(S::States, U::AbstractPathUtility) - -Negative affine transformation of path utility. Always evaluates negative values. - -# Examples -```julia-repl -julia> U⁻ = NegativePathUtility(S, U) -julia> all(U⁻(s) < 0 for s in paths(S)) -true -``` -""" -struct NegativePathUtility <: AbstractPathUtility - U::AbstractPathUtility - max::Float64 - function NegativePathUtility(S::States, U::AbstractPathUtility) - u_max = maximum(U(s) for s in paths(S)) - new(U, u_max) - end -end - -(U::NegativePathUtility)(s::Path) = U.U(s) - U.max - 1 - - """ expected_value(model::Model, - x_s::PathCompatibilityVariables, - U::AbstractPathUtility, - P::AbstractPathProbability; + diagram::InfluenceDiagram, + x_s::PathCompatibilityVariables; probability_scale_factor::Float64=1.0) Create an expected value objective. # Arguments - `model::Model`: JuMP model into which variables are added. +- `diagram::InfluenceDiagram`: Influence diagram structure. - `x_s::PathCompatibilityVariables`: Path compatibility variables. -- `S::States`: States structure associated with the influence diagram. -- `P::AbstractPathProbability`: Path probabilities structure for which the function - `P(s)` is defined and returns the path probabilities for path `s`. - `probability_scale_factor::Float64`: Multiplies the path probabilities by this factor. # Examples ```julia -EV = expected_value(model, x_s, U, P) -EV = expected_value(model, x_s, U, P; probability_scale_factor = 10.0) +julia> EV = expected_value(model, diagram, x_s) +julia> EV = expected_value(model, diagram, x_s; probability_scale_factor = 10.0) ``` """ function expected_value(model::Model, - x_s::PathCompatibilityVariables, - U::AbstractPathUtility, - P::AbstractPathProbability; + diagram::InfluenceDiagram, + x_s::PathCompatibilityVariables; probability_scale_factor::Float64=1.0) if probability_scale_factor ≤ 0 throw(DomainError("The probability_scale_factor must be greater than 0.")) end - @expression(model, sum(P(s) * x * U(s) * probability_scale_factor for (s, x) in x_s)) + @expression(model, sum(diagram.P(s) * x * diagram.U(s, diagram.translation) * probability_scale_factor for (s, x) in x_s)) end """ conditional_value_at_risk(model::Model, + diagram, x_s::PathCompatibilityVariables{N}, - U::AbstractPathUtility, - P::AbstractPathProbability, α::Float64; probability_scale_factor::Float64=1.0) where N @@ -285,10 +228,8 @@ Create a conditional value-at-risk (CVaR) objective. # Arguments - `model::Model`: JuMP model into which variables are added. +- `diagram::InfluenceDiagram`: Influence diagram structure. - `x_s::PathCompatibilityVariables`: Path compatibility variables. -- `S::States`: States structure associated with the influence diagram. -- `P::AbstractPathProbability`: Path probabilities structure for which the function - `P(s)` is defined and returns the path probabilities for path `s`. - `α::Float64`: Probability level at which conditional value-at-risk is optimised. - `probability_scale_factor::Float64`: Adjusts conditional value at risk model to be compatible with the expected value expression if the probabilities were scaled there. @@ -297,15 +238,14 @@ Create a conditional value-at-risk (CVaR) objective. # Examples ```julia -α = 0.05 # Parameter such that 0 ≤ α ≤ 1 -CVaR = conditional_value_at_risk(model, x_s, U, P, α) -CVaR = conditional_value_at_risk(model, x_s, U, P, α; probability_scale_factor = 10.0) +julia> α = 0.05 # Parameter such that 0 ≤ α ≤ 1 +julia> CVaR = conditional_value_at_risk(model, x_s, U, P, α) +julia> CVaR = conditional_value_at_risk(model, x_s, U, P, α; probability_scale_factor = 10.0) ``` """ function conditional_value_at_risk(model::Model, + diagram::InfluenceDiagram, x_s::PathCompatibilityVariables{N}, - U::AbstractPathUtility, - P::AbstractPathProbability, α::Float64; probability_scale_factor::Float64=1.0) where N @@ -321,7 +261,7 @@ function conditional_value_at_risk(model::Model, end # Pre-computed parameters - u = collect(Iterators.flatten(U(s) for s in keys(x_s))) + u = collect(Iterators.flatten(diagram.U(s, diagram.translation) for s in keys(x_s))) u_sorted = sort(u) u_min = u_sorted[1] u_max = u_sorted[end] @@ -335,7 +275,7 @@ function conditional_value_at_risk(model::Model, @constraint(model, η ≤ u_max) ρ′_s = Dict{Path{N}, VariableRef}() for (s, x) in x_s - u_s = U(s) + u_s = diagram.U(s, diagram.translation) λ = @variable(model, binary=true) λ′ = @variable(model, binary=true) ρ = @variable(model) @@ -349,14 +289,14 @@ function conditional_value_at_risk(model::Model, @constraint(model, ρ ≤ λ * probability_scale_factor) @constraint(model, ρ′ ≤ λ′* probability_scale_factor) @constraint(model, ρ ≤ ρ′) - @constraint(model, ρ′ ≤ x * P(s) * probability_scale_factor) - @constraint(model, (x * P(s) - (1 - λ))* probability_scale_factor ≤ ρ) + @constraint(model, ρ′ ≤ x * diagram.P(s) * probability_scale_factor) + @constraint(model, (x * diagram.P(s) - (1 - λ))* probability_scale_factor ≤ ρ) ρ′_s[s] = ρ′ end @constraint(model, sum(values(ρ′_s)) == α * probability_scale_factor) # Return CVaR as an expression - CVaR = @expression(model, sum(ρ_bar * U(s) for (s, ρ_bar) in ρ′_s) / α) + CVaR = @expression(model, sum(ρ_bar * diagram.U(s, diagram.translation) for (s, ρ_bar) in ρ′_s) / α) return CVaR end @@ -368,8 +308,8 @@ end Construct decision strategy from variable refs. """ -function LocalDecisionStrategy(j::Node, z::Array{VariableRef}) - LocalDecisionStrategy(j, @. Int(round(value(z)))) +function LocalDecisionStrategy(d::Node, z::Array{VariableRef}) + LocalDecisionStrategy(d, @. Int(round(value(z)))) end """ @@ -383,5 +323,5 @@ Z = DecisionStrategy(z) ``` """ function DecisionStrategy(z::DecisionVariables) - DecisionStrategy(z.D, [LocalDecisionStrategy(d.j, v) for (d, v) in zip(z.D, z.z)]) + DecisionStrategy(z.D, z.I_d, [LocalDecisionStrategy(d, z_var) for (d, z_var) in zip(z.D, z.z)]) end diff --git a/src/influence_diagram.jl b/src/influence_diagram.jl index fc17e339..56a95357 100644 --- a/src/influence_diagram.jl +++ b/src/influence_diagram.jl @@ -4,101 +4,72 @@ using Base.Iterators: product # --- Nodes and States --- """ - Node = Int + Node = Int16 -Primitive type for node. Alias for `Int`. +Primitive type for node index. Alias for `Int16`. """ -const Node = Int +const Node = Int16 """ - abstract type AbstractNode end + Name = String -Node type for directed, acyclic graph. +Primitive type for node names. Alias for `String`. """ -abstract type AbstractNode end +const Name = String -function validate_node(j::Node, I_j::Vector{Node}) - if !allunique(I_j) - throw(DomainError("All information nodes should be unique.")) - end - if !all(i < j for i in I_j) - throw(DomainError("All nodes in the information set must be less than node j.")) - end -end +""" + abstract type AbstractNode end +Node type for directed, acyclic graph. """ - ChanceNode <: AbstractNode +abstract type AbstractNode end -Chance node type. -# Examples -```julia -c = ChanceNode(3, [1, 2]) -``` -""" struct ChanceNode <: AbstractNode - j::Node - I_j::Vector{Node} - function ChanceNode(j::Node, I_j::Vector{Node}) - validate_node(j, I_j) - new(j, I_j) + name::Name + I_j::Vector{Name} + states::Vector{Name} + function ChanceNode(name, I_j, states) + return new(name, I_j, states) end -end - -""" - DecisionNode <: AbstractNode -Decision node type. +end -# Examples -```julia -d = DecisionNode(2, [1]) -``` -""" struct DecisionNode <: AbstractNode - j::Node - I_j::Vector{Node} - function DecisionNode(j::Node, I_j::Vector{Node}) - validate_node(j, I_j) - new(j, I_j) + name::Name + I_j::Vector{Name} + states::Vector{Name} + function DecisionNode(name, I_j, states) + return new(name, I_j, states) end end -""" - ValueNode <: AbstractNode - -Value node type. - -# Examples -```julia -v = ValueNode(4, [1, 3]) -``` -""" struct ValueNode <: AbstractNode - j::Node - I_j::Vector{Node} - function ValueNode(j::Node, I_j::Vector{Node}) - validate_node(j, I_j) - new(j, I_j) + name::Name + I_j::Vector{Name} + function ValueNode(name, I_j) + return new(name, I_j) end end + + """ const State = Int -Primitive type for the number of states. Alias for `Int`. +Primitive type for the number of states. Alias for `Int16`. """ -const State = Int +const State = Int16 """ - States <: AbstractArray{State, 1} + struct States <: AbstractArray{State, 1} States type. Works like `Vector{State}`. # Examples ```julia -S = States([2, 3, 2, 4]) +julia> S = States([2, 3, 2, 4]) ``` """ struct States <: AbstractArray{State, 1} @@ -117,68 +88,47 @@ Base.getindex(S::States, i::Int) = getindex(S.vals, i) Base.length(S::States) = length(S.vals) Base.eltype(S::States) = eltype(S.vals) -""" - function States(states::Vector{Tuple{State, Vector{Node}}}) -Construct states from vector of (state, nodes) tuples. -# Examples -```julia-repl -julia> S = States([(2, [1, 3]), (3, [2, 4, 5])]) -States([2, 3, 2, 3, 3]) -``` +# --- Paths --- + """ -function States(states::Vector{Tuple{State, Vector{Node}}}) - S_j = Vector{State}(undef, sum(length(j) for (_, j) in states)) - for (s, j) in states - S_j[j] .= s - end - States(S_j) -end + const Path{N} = NTuple{N, State} where N +Path type. Alias for `NTuple{N, State} where N`. """ - function validate_influence_diagram(S::States, C::Vector{ChanceNode}, D::Vector{DecisionNode}, V::Vector{ValueNode}) +const Path{N} = NTuple{N, State} where N + -Validate influence diagram. """ -function validate_influence_diagram(S::States, C::Vector{ChanceNode}, D::Vector{DecisionNode}, V::Vector{ValueNode}) - n = length(C) + length(D) - if length(S) != n - throw(DomainError("Each change and decision node should have states.")) - end - if Set(c.j for c in C) ∪ Set(d.j for d in D) != Set(1:n) - throw(DomainError("Union of change and decision nodes should be {1,...,n}.")) - end - if Set(v.j for v in V) != Set((n+1):(n+length(V))) - throw(DomainError("Values nodes should be {n+1,...,n+|V|}.")) - end - I_V = union((v.I_j for v in V)...) - if !(I_V ⊆ Set(1:n)) - throw(DomainError("Each information set I(v) for value node v should be a subset of C∪D.")) - end - # Check for redundant nodes. - leaf_nodes = setdiff(1:n, (c.I_j for c in C)..., (d.I_j for d in D)...) - for i in leaf_nodes - if !(i∈I_V) - @warn("Chance or decision node $i is redundant.") - end - end - for v in V - if isempty(v.I_j) - @warn("Value node $(v.j) is redundant.") - end - end -end + const ForbiddenPath = Tuple{Vector{Node}, Set{Path}} +ForbiddenPath type. + +# Examples +```julia +julia> ForbiddenPath(([1, 2], Set([(1, 2)]))) +julia> ForbiddenPath[ + ([1, 2], Set([(1, 2)])), + ([3, 4, 5], Set([(1, 2, 3), (3, 4, 5)])) +] +``` +""" +const ForbiddenPath = Tuple{Vector{Node}, Set{Path}} -# --- Paths --- """ - const Path{N} = NTuple{N, State} where N + const FixedPath = Dict{Node, State} -Path type. Alias for `NTuple{N, State} where N`. +FixedPath type. + +# Examples +```julia +julia> FixedPath(Dict(1=>1, 2=>3)) +``` """ -const Path{N} = NTuple{N, State} where N +const FixedPath = Dict{Node, State} + """ function paths(states::AbstractVector{State}) @@ -197,18 +147,20 @@ function paths(states::AbstractVector{State}) end """ - function paths(states::AbstractVector{State}, fixed::Dict{Node, State}) + function paths(states::AbstractVector{State}, fixed::FixedPath) Iterate over paths with fixed states in lexicographical order. # Examples ```julia-repl julia> states = States([2, 3]) -julia> vec(collect(paths(states, Dict(1=>2)))) +julia> vec(collect(paths(states, Dict(Node(1) => State(2))))) [(2, 1), (2, 2), (2, 3)] + +julia> vec(collect(paths(states, FixedPath(diagram, Dict("O" => "lemon"))))) ``` """ -function paths(states::AbstractVector{State}, fixed::Dict{Node, State}) +function paths(states::AbstractVector{State}, fixed::FixedPath) iters = collect(UnitRange.(one(eltype(states)), states)) for (i, v) in fixed iters[i] = UnitRange(v, v) @@ -216,21 +168,6 @@ function paths(states::AbstractVector{State}, fixed::Dict{Node, State}) product(iters...) end -""" - const ForbiddenPath = Tuple{Vector{Node}, Set{Path}} - -ForbiddenPath type. - -# Examples -```julia -ForbiddenPath[ - ([1, 2], Set([(1, 2)])), - ([3, 4, 5], Set([(1, 2, 3), (3, 4, 5)])) -] -``` -""" -const ForbiddenPath = Tuple{Vector{Node}, Set{Path}} - # --- Probabilities --- @@ -242,22 +179,22 @@ Construct and validate stage probabilities. # Examples ```julia-repl julia> data = [0.5 0.5 ; 0.2 0.8] -julia> X = Probabilities(2, data) +julia> X = Probabilities(Node(2), data) julia> s = (1, 2) julia> X(s) 0.5 ``` """ struct Probabilities{N} <: AbstractArray{Float64, N} - j::Node + c::Node data::Array{Float64, N} - function Probabilities(j::Node, data::Array{Float64, N}) where N + function Probabilities(c::Node, data::Array{Float64, N}) where N for i in CartesianIndices(size(data)[1:end-1]) if !(sum(data[i, :]) ≈ 1) throw(DomainError("Probabilities should sum to one.")) end end - new{N}(j, data) + new{N}(c, data) end end @@ -278,65 +215,87 @@ Abstract path probability type. # Examples ```julia -struct PathProbability <: AbstractPathProbability - C::Vector{ChanceNode} - # ... +julia> struct PathProbability <: AbstractPathProbability + C::Vector{ChanceNode} + # ... end -(U::PathProbability)(s::Path) = ... +julia> (P::PathProbability)(s::Path) = ... ``` """ abstract type AbstractPathProbability end """ - DefaultPathProbability <: AbstractPathProbability + struct DefaultPathProbability <: AbstractPathProbability Path probability. # Examples ```julia -P = DefaultPathProbability(C, X) -s = (1, 2) -P(s) +julia> P = DefaultPathProbability(diagram.C, diagram.X) +julia> s = (1, 2) +julia> P(s) ``` """ struct DefaultPathProbability <: AbstractPathProbability - C::Vector{ChanceNode} + C::Vector{Node} + I_c::Vector{Vector{Node}} X::Vector{Probabilities} + function DefaultPathProbability(C, I_c, X) + if length(C) == length(I_c) + new(C, I_c, X) + else + throw(DomainError("The number of chance nodes and information sets given to DefaultPathProbability should be equal.")) + end + end + end function (P::DefaultPathProbability)(s::Path) - prod(X(s[[c.I_j; c.j]]) for (c, X) in zip(P.C, P.X)) + prod(X(s[[I_c; c]]) for (c, I_c, X) in zip(P.C, P.I_c, P.X)) end -# --- Consequences --- +# --- Utilities --- + +""" + const Utility = Float32 + +Primitive type for utility. Alias for `Float32`. +""" +const Utility = Float32 """ - Consequences{N} <: AbstractArray{Float64, N} + struct Utilities{N} <: AbstractArray{Utility, N} State utilities. # Examples ```julia-repl julia> vals = [1.0 -2.0; 3.0 4.0] -julia> Y = Consequences(3, vals) +julia> Y = Utilities(3, vals) julia> s = (1, 2) julia> Y(s) -2.0 ``` """ -struct Consequences{N} <: AbstractArray{Float64, N} - j::Node - data::Array{Float64, N} +struct Utilities{N} <: AbstractArray{Utility, N} + v::Node + data::Array{Utility, N} + function Utilities(v::Node, data::Array{Utility, N}) where N + if any(isinf(u) for u in data) + throw(DomainError("A value should be defined for each element of a utility matrix.")) + end + new{N}(v, data) + end end -Base.size(Y::Consequences) = size(Y.data) -Base.IndexStyle(::Type{<:Consequences}) = IndexLinear() -Base.getindex(Y::Consequences, i::Int) = getindex(Y.data, i) -Base.getindex(Y::Consequences, I::Vararg{Int,N}) where N = getindex(Y.data, I...) +Base.size(Y::Utilities) = size(Y.data) +Base.IndexStyle(::Type{<:Utilities}) = IndexLinear() +Base.getindex(Y::Utilities, i::Int) = getindex(Y.data, i) +Base.getindex(Y::Utilities, I::Vararg{Int,N}) where N = getindex(Y.data, I...) -(Y::Consequences)(s::Path) = Y[s...] +(Y::Utilities)(s::Path) = Y[s...] # --- Path Utility --- @@ -348,35 +307,715 @@ Abstract path utility type. # Examples ```julia -struct PathUtility <: AbstractPathUtility - V::Vector{ValueNode} - # ... -end +julia> struct PathUtility <: AbstractPathUtility + V::Vector{ValueNode} + # ... + end -(U::PathUtility)(s::Path) = ... +julia> (U::PathUtility)(s::Path) = ... +julia> (U::PathUtility)(s::Path, translation::Utility) = ... ``` """ abstract type AbstractPathUtility end """ - DefaultPathUtility <: AbstractPathUtility + struct DefaultPathUtility <: AbstractPathUtility Default path utility. # Examples ```julia -U = DefaultPathUtility(V, Y) -s = (1, 2) -U(s) +julia> U = DefaultPathUtility(V, Y) +julia> s = (1, 2) +julia> U(s) + +julia> t = -100.0 +julia> U(s, t) ``` """ struct DefaultPathUtility <: AbstractPathUtility - V::Vector{ValueNode} - Y::Vector{Consequences} + I_v::Vector{Vector{Node}} + Y::Vector{Utilities} end function (U::DefaultPathUtility)(s::Path) - sum(Y(s[v.I_j]) for (v, Y) in zip(U.V, U.Y)) + sum(Y(s[I_v]) for (I_v, Y) in zip(U.I_v, U.Y)) +end + +function (U::DefaultPathUtility)(s::Path, t::Utility) + U(s) + t +end + +# --- Influence diagram --- +""" + mutable struct InfluenceDiagram + Nodes::Vector{AbstractNode} + Names::Vector{Name} + I_j::Vector{Vector{Node}} + States::Vector{Vector{Name}} + S::States + C::Vector{Node} + D::Vector{Node} + V::Vector{Node} + X::Vector{Probabilities} + Y::Vector{Utilities} + P::AbstractPathProbability + U::AbstractPathUtility + translation::Utility + function InfluenceDiagram() + new(Vector{AbstractNode}()) + end + end + +Hold all information related to the influence diagram. + +# Fields +- `Nodes::Vector{AbstractNode}`: Vector of added abstract nodes. +- `Names::Vector{Name}`: Names of nodes in order of their indices. +- `I_j::Vector{Vector{Node}}`: Information sets of nodes in order of their indices. + Nodes of information sets identified by their indices. +- `States::Vector{Vector{Name}}`: States of each node in order of their indices. +- `S::States`: Vector showing the number of states each node has. +- `C::Vector{Node}`: Indices of chance nodes in ascending order. +- `D::Vector{Node}`: Indices of decision nodes in ascending order. +- `V::Vector{Node}`: Indices of value nodes in ascending order. +- `X::Vector{Probabilities}`: Probability matrices of chance nodes in order of chance + nodes in C. +- `Y::Vector{Utilities}`: Utility matrices of value nodes in order of value nodes in V. +- `P::AbstractPathProbability`: Path probabilities. +- `U::AbstractPathUtility`: Path utilities. +- `translation::Utility`: Utility translation for storing the positive or negative + utility translation. + + +# Examples +```julia +julia> diagram = InfluenceDiagram() +``` +""" +mutable struct InfluenceDiagram + Nodes::Vector{AbstractNode} + Names::Vector{Name} + I_j::Vector{Vector{Node}} + States::Vector{Vector{Name}} + S::States + C::Vector{Node} + D::Vector{Node} + V::Vector{Node} + X::Vector{Probabilities} + Y::Vector{Utilities} + P::AbstractPathProbability + U::AbstractPathUtility + translation::Utility + function InfluenceDiagram() + new(Vector{AbstractNode}()) + end +end + + +# --- Adding nodes --- + +function validate_node(diagram::InfluenceDiagram, + name::Name, + I_j::Vector{Name}; + value_node::Bool=false, + states::Vector{Name}=Vector{Name}()) + + if !allunique([map(x -> x.name, diagram.Nodes)..., name]) + throw(DomainError("All node names should be unique.")) + end + + if !allunique(I_j) + throw(DomainError("All nodes in an information set should be unique.")) + end + + if !allunique([name, I_j...]) + throw(DomainError("Node should not be included in its own information set.")) + end + + if !value_node + if length(states) < 2 + throw(DomainError("Each chance and decision node should have more than one state.")) + end + end + + if value_node + if isempty(I_j) + @warn("Value node $name is redundant.") + end + end +end + +""" + function add_node!(diagram::InfluenceDiagram, node::AbstractNode) + +Add node to influence diagram structure. + +# Examples +```julia +julia> add_node!(diagram, ChanceNode("O", [], ["lemon", "peach"])) +``` +""" +function add_node!(diagram::InfluenceDiagram, node::AbstractNode) + if !isa(node, ValueNode) + validate_node(diagram, node.name, node.I_j, states = node.states) + else + validate_node(diagram, node.name, node.I_j, value_node = true) + end + push!(diagram.Nodes, node) +end + + +# --- Adding Probabilities --- +""" + struct ProbabilityMatrix{N} <: AbstractArray{Float64, N} + nodes::Vector{Name} + indices::Vector{Dict{Name, Int}} + matrix::Array{Float64, N} + end + +Construct probability matrix. +""" +struct ProbabilityMatrix{N} <: AbstractArray{Float64, N} + nodes::Vector{Name} + indices::Vector{Dict{Name, Int}} + matrix::Array{Float64, N} +end + +Base.size(PM::ProbabilityMatrix) = size(PM.matrix) +Base.getindex(PM::ProbabilityMatrix, I::Vararg{Int,N}) where N = getindex(PM.matrix, I...) +Base.setindex!(PM::ProbabilityMatrix, p::T, I::Vararg{Int,N}) where {N, T<:Real} = (PM.matrix[I...] = p) +Base.setindex!(PM::ProbabilityMatrix{N}, X::Array{T}, I::Vararg{Any, N}) where {N, T<:Real} = (PM.matrix[I...] .= X) + +""" + function ProbabilityMatrix(diagram::InfluenceDiagram, node::Name) + +Initialise a probability matrix for a given chance node. The matrix is initialised with zeros. + +# Examples +```julia +julia> X_O = ProbabilityMatrix(diagram, "O") +``` +""" +function ProbabilityMatrix(diagram::InfluenceDiagram, node::Name) + if node ∉ diagram.Names + throw(DomainError("Node $node should be added as a node to the influence diagram.")) + end + if node ∉ diagram.Names[diagram.C] + throw(DomainError("Only chance nodes can have probability matrices.")) + end + + # Find the node's indices and it's I_c nodes + c = findfirst(x -> x==node, diagram.Names) + nodes = [diagram.I_j[c]..., c] + names = diagram.Names[nodes] + + indices = Vector{Dict{Name, Int}}() + for j in nodes + states = Dict{Name, Int}(state => i + for (i, state) in enumerate(diagram.States[j]) + ) + push!(indices, states) + end + matrix = fill(0.0, diagram.S[nodes]...) + + return ProbabilityMatrix(names, indices, matrix) +end + +""" + function set_probability!(probability_matrix::ProbabilityMatrix, scenario::Array{String}, probability::Float64) + +Set a single probability value into probability matrix. + +# Examples +```julia +julia> X_O = ProbabilityMatrix(diagram, "O") +julia> set_probability!(X_O, ["peach"], 0.8) +julia> set_probability!(X_O, ["lemon"], 0.2) +``` +""" +function set_probability!(probability_matrix::ProbabilityMatrix, scenario::Array{String}, probability::Real) + index = Vector{Int}() + for (i, s) in enumerate(scenario) + if get(probability_matrix.indices[i], s, 0) == 0 + throw(DomainError("Node $(probability_matrix.nodes[i]) does not have a state called $s.")) + else + push!(index, get(probability_matrix.indices[i], s, 0)) + end + end + + probability_matrix[index...] = probability +end + +""" + function set_probability!(probability_matrix::ProbabilityMatrix, scenario::Array{Any}, probabilities::Array{Float64}) + +Set multiple probability values into probability matrix. + +# Examples +```julia +julia> X_O = ProbabilityMatrix(diagram, "O") +julia> set_probability!(X_O, ["lemon", "peach"], [0.2, 0.8]) +``` +""" +function set_probability!(probability_matrix::ProbabilityMatrix, scenario::Array{Any}, probabilities::Array{T}) where T<:Real + index = Vector{Any}() + for (i, s) in enumerate(scenario) + if isa(s, Colon) + push!(index, s) + elseif get(probability_matrix.indices[i], s, 0) == 0 + throw(DomainError("Node $(probability_matrix.nodes[i]) does not have state $s.")) + else + push!(index, get(probability_matrix.indices[i], s, 0)) + end + end + + probability_matrix[index...] = probabilities +end + +""" + function add_probabilities!(diagram::InfluenceDiagram, node::Name, probabilities::AbstractArray{Float64, N}) where N + +Add probability matrix to influence diagram, specifically to its X vector. + +# Examples +```julia +julia> X_O = ProbabilityMatrix(diagram, "O") +julia> set_probability!(X_O, ["lemon", "peach"], [0.2, 0.8]) +julia> add_probabilities!(diagram, "O", X_O) + +julia> add_probabilities!(diagram, "O", [0.2, 0.8]) + +!!! note +The arcs must be generated before probabilities or utilities can be added to the influence diagram. +``` +""" +function add_probabilities!(diagram::InfluenceDiagram, node::Name, probabilities::AbstractArray{Float64, N}) where N + c = findfirst(x -> x==node, diagram.Names) + + if c ∈ [j.c for j in diagram.X] + throw(DomainError("Probabilities should be added only once for each node.")) + end + + if size(probabilities) == Tuple((diagram.S[j] for j in (diagram.I_j[c]..., c))) + if isa(probabilities, ProbabilityMatrix) + # Check that probabilities sum to one happesn in Probabilities + push!(diagram.X, Probabilities(Node(c), probabilities.matrix)) + else + push!(diagram.X, Probabilities(Node(c), probabilities)) + end + else + throw(DomainError("The dimensions of a probability matrix should match the node's states' and information states' cardinality. Expected $(Tuple((diagram.S[n] for n in (diagram.I_j[c]..., c)))) for node $name, got $(size(probabilities)).")) + end +end + + +# --- Adding Utilities --- + +""" + struct UtilityMatrix{N} <: AbstractArray{Utility, N} + I_v::Vector{Name} + indices::Vector{Dict{Name, Int}} + matrix::Array{Utility, N} + end + +Construct utility matrix. +""" +struct UtilityMatrix{N} <: AbstractArray{Utility, N} + I_v::Vector{Name} + indices::Vector{Dict{Name, Int}} + matrix::Array{Utility, N} +end + +Base.size(UM::UtilityMatrix) = size(UM.matrix) +Base.getindex(UM::UtilityMatrix, I::Vararg{Int,N}) where N = getindex(UM.matrix, I...) +Base.setindex!(UM::UtilityMatrix, y::T, I::Vararg{Int,N}) where {N, T<:Real} = (UM.matrix[I...] = y) +Base.setindex!(UM::UtilityMatrix{N}, Y::Array{T}, I::Vararg{Any, N}) where {N, T<:Real} = (UM.matrix[I...] .= Y) + +""" + function UtilityMatrix(diagram::InfluenceDiagram, node::Name) + +Initialise a utility matrix for a value node. The matrix is initialised with `Inf` values. + +# Examples +```julia +julia> Y_V3 = UtilityMatrix(diagram, "V3") +``` +""" +function UtilityMatrix(diagram::InfluenceDiagram, node::Name) + if node ∉ diagram.Names + throw(DomainError("Node $node should be added as a node to the influence diagram.")) + end + if node ∉ diagram.Names[diagram.V] + throw(DomainError("Only value nodes can have consequence matrices.")) + end + + # Find the node's indexand it's I_v nodes + v = findfirst(x -> x==node, diagram.Names) + I_v = diagram.I_j[v] + names = diagram.Names[I_v] + + indices = Vector{Dict{Name, Int}}() + for j in I_v + states = Dict{Name, Int}(state => i + for (i, state) in enumerate(diagram.States[j]) + ) + push!(indices, states) + end + matrix = Array{Utility}(fill(Inf, diagram.S[I_v]...)) + + return UtilityMatrix(names, indices, matrix) +end + +""" + function set_utility!(utility_matrix::UtilityMatrix, scenario::Array{String}, utility::Real) + +Set a single utility value into utility matrix. + +# Examples +```julia +julia> Y_V3 = UtilityMatrix(diagram, "V3") +julia> set_utility!(Y_V3, ["lemon", "buy without guarantee"], -200) +``` +""" +function set_utility!(utility_matrix::UtilityMatrix, scenario::Array{String}, utility::Real) + index = Vector{Int}() + for (i, s) in enumerate(scenario) + if get(utility_matrix.indices[i], s, 0) == 0 + throw(DomainError("Node $(utility_matrix.I_v[i]) does not have a state called $s.")) + else + push!(index, get(utility_matrix.indices[i], s, 0)) + end + end + + utility_matrix[index...] = utility +end + +""" + function set_utility!(utility_matrix::UtilityMatrix, scenario::Array{Any}, utility::Array{T}) where T<:Real + +Set multiple utility values into utility matrix. + +# Examples +```julia +julia> Y_V3 = UtilityMatrix(diagram, "V3") +julia> set_utility!(Y_V3, ["peach", :], [-40, -20, 0]) +``` +""" +function set_utility!(utility_matrix::UtilityMatrix, scenario::Array{Any}, utility::Array{T}) where T<:Real + index = Vector{Any}() + for (i, s) in enumerate(scenario) + if isa(s, Colon) + push!(index, s) + elseif get(utility_matrix.indices[i], s, 0) == 0 + throw(DomainError("Node $(utility_matrix.I_v[i]) does not have state $s.")) + else + push!(index, get(utility_matrix.indices[i], s, 0)) + end + end + + utility_matrix[index...] = utility +end + + +""" + function add_utilities!(diagram::InfluenceDiagram, node::Name, utilities::AbstractArray{T, N}) where {N,T<:Real} + +Add utility matrix to influence diagram, specifically to its Y vector. + +# Examples +```julia +julia> Y_V3 = UtilityMatrix(diagram, "V3") +julia> set_utility!(Y_V3, ["peach", :], [-40, -20, 0]) +julia> set_utility!(Y_V3, ["lemon", :], [-200, 0, 0]) +julia> add_utilities!(diagram, "V3", Y_V3) + +julia> add_utilities!(diagram, "V1", [0, -25]) +``` +!!! note +The arcs must be generated before probabilities or utilities can be added to the influence diagram. +""" +function add_utilities!(diagram::InfluenceDiagram, node::Name, utilities::AbstractArray{T, N}) where {N,T<:Real} + v = findfirst(x -> x==node, diagram.Names) + + if v ∈ [j.v for j in diagram.Y] + throw(DomainError("Utilities should be added only once for each node.")) + end + if any(u ==Inf for u in utilities) + throw(DomainError("Utility values should be less than infinity.")) + end + + if size(utilities) == Tuple((diagram.S[j] for j in diagram.I_j[v])) + if isa(utilities, UtilityMatrix) + push!(diagram.Y, Utilities(Node(v), utilities.matrix)) + else + # Conversion to Float32 using Utility(), since machine default is Float64 + push!(diagram.Y, Utilities(Node(v), [Utility(u) for u in utilities])) + end + else + throw(DomainError("The dimensions of the utilities matrix should match the node's information states' cardinality. Expected $(Tuple((diagram.S[n] for n in diagram.I_j[v]))) for node $name, got $(size(utilities)).")) + end +end + + +# --- Generating Arcs --- + +function validate_structure(Nodes::Vector{AbstractNode}, C_and_D::Vector{AbstractNode}, n_CD::Int, V::Vector{AbstractNode}, n_V::Int) + # Validating node structure + if n_CD == 0 + throw(DomainError("The influence diagram must have chance or decision nodes.")) + end + if !(union((n.I_j for n in Nodes)...) ⊆ Set(n.name for n in Nodes)) + throw(DomainError("Each node that is part of an information set should be added as a node.")) + end + # Checking the information sets of C and D nodes + if !isempty(union((j.I_j for j in C_and_D)...) ∩ Set(v.name for v in V)) + throw(DomainError("Information sets should not include any value nodes.")) + end + # Checking the information sets of V nodes + if !isempty(V) && !isempty(union((v.I_j for v in V)...) ∩ Set(v.name for v in V)) + throw(DomainError("Information sets should not include any value nodes.")) + end + # Check for redundant chance or decision nodes. + last_CD_nodes = setdiff((j.name for j in C_and_D), (j.I_j for j in C_and_D)...) + for i in last_CD_nodes + if !isempty(V) && i ∉ union((v.I_j for v in V)...) + @warn("Node $i is redundant.") + end + end +end + +""" + function generate_arcs!(diagram::InfluenceDiagram) + +Generate arc structures using nodes added to influence diagram, by ordering nodes, +giving them indices and generating correct values for the vectors Names, I_j, states, +S, C, D, V in the influence digram. Abstraction is created and the names of the nodes +and states are only used in the user interface from here on. + +# Examples +```julia +julia> generate_arcs!(diagram) +``` +""" +function generate_arcs!(diagram::InfluenceDiagram) + + # Chance and decision nodes + C_and_D = filter(x -> !isa(x, ValueNode), diagram.Nodes) + n_CD = length(C_and_D) + # Value nodes + V_nodes = filter(x -> isa(x, ValueNode), diagram.Nodes) + n_V = length(V_nodes) + + validate_structure(diagram.Nodes, C_and_D, n_CD, V_nodes, n_V) + + # Declare vectors for results (final resting place InfluenceDiagram.Names and InfluenceDiagram.I_j) + Names = Vector{Name}(undef, n_CD+n_V) + I_j = Vector{Vector{Node}}(undef, n_CD+n_V) + states = Vector{Vector{Name}}() + S = Vector{State}(undef, n_CD) + C = Vector{Node}() + D = Vector{Node}() + V = Vector{Node}() + + # Declare helper collections + indices = Dict{Name, Node}() + indexed_nodes = Set{Name}() + # Declare index + index = 1 + + + while true + # Index nodes C and D that don't yet have indices but whose I_j have indices + new_nodes = filter(j -> (j.name ∉ indexed_nodes && Set(j.I_j) ⊆ indexed_nodes), C_and_D) + for j in new_nodes + # Update helper collections + push!(indices, j.name => index) + push!(indexed_nodes, j.name) + # Update results + Names[index] = Name(j.name) #TODO datatype conversion happens here, should we use push! ? + I_j[index] = map(x -> Node(indices[x]), j.I_j) + push!(states, j.states) + S[index] = State(length(j.states)) + if isa(j, ChanceNode) + push!(C, Node(index)) + else + push!(D, Node(index)) + end + # Increase index + index += 1 + end + + # If no new nodes were indexed this iteration, terminate while loop + if isempty(new_nodes) + if index < n_CD + throw(DomainError("The influence diagram should be acyclic.")) + else + break + end + end + end + + + # Index value nodes + for v in V_nodes + # Update results + Names[index] = Name(v.name) + I_j[index] = map(x -> Node(indices[x]), v.I_j) + push!(V, Node(index)) + # Increase index + index += 1 + end + + diagram.Names = Names + diagram.I_j = I_j + diagram.States = states + diagram.S = States(S) + diagram.C = C + diagram.D = D + diagram.V = V + # Declaring X and Y + diagram.X = Vector{Probabilities}() + diagram.Y = Vector{Utilities}() +end + + + +# --- Generating Diagram --- +""" +function generate_diagram!(diagram::InfluenceDiagram; + default_probability::Bool=true, + default_utility::Bool=true, + positive_path_utility::Bool=false, + negative_path_utility::Bool=false) + +Generate complete influence diagram with probabilities and utilities as well. + +# Arguments +- `default_probability::Bool=true`: Choice to use default path probabilities. +- `default_utility::Bool=true`: Choice to use default path utilities. +- `positive_path_utility::Bool=false`: Choice to use a positive path utility translation. +- `negative_path_utility::Bool=false`: Choice to use a negative path utility translation. + +# Examples +```julia +julia> generate_diagram!(diagram) +``` + +!!! note +The influence diagram must be generated after probabilities and utilities are added +but before creating the decision model. + +!!! note +If the default probabilities and utilities are not used, define `AbstractPathProbability` +and `AbstractPathUtility` structures and define P(s), U(s) and U(s, t) functions +for them. Add the `AbstractPathProbability` and `AbstractPathUtility` structures +to the influence diagram fields P and U. +""" +function generate_diagram!(diagram::InfluenceDiagram; + default_probability::Bool=true, + default_utility::Bool=true, + positive_path_utility::Bool=false, + negative_path_utility::Bool=false) + + + # Sort probabilities and consequences + sort!(diagram.X, by = x -> x.c) + sort!(diagram.Y, by = x -> x.v) + + + # Declare P and U if defaults are used + if default_probability + diagram.P = DefaultPathProbability(diagram.C, diagram.I_j[diagram.C], diagram.X) + end + if default_utility + diagram.U = DefaultPathUtility(diagram.I_j[diagram.V], diagram.Y) + if positive_path_utility + # Conversion to Float32 using Utility(), since machine default is Float64 + diagram.translation = 1 - minimum(diagram.U(s) for s in paths(diagram.S)) + elseif negative_path_utility + diagram.translation = -1 - maximum(diagram.U(s) for s in paths(diagram.S)) + else + diagram.translation = 0 + end + end + +end + +# --- ForbiddenPath and FixedPath outer construction functions --- +""" + function ForbiddenPath(diagram::InfluenceDiagram, nodes::Vector{Name}, paths::Vector{NTuple{N, Name}}) where N + +ForbiddenPath outer construction function. Create ForbiddenPath variable. + +# Arguments +- `diagram::InfluenceDiagram`: Influence diagram structure +- `nodes::Vector{Name}`: Vector of nodes involved in forbidden paths. Identified by their names. +- `paths`::Vector{NTuple{N, Name}}`: Vector of tuples defining the forbidden combinations of states. States identified by their names. + +# Example +```julia +julia> ForbiddenPath(diagram, ["R1", "R2"], [("high", "low"), ("low", "high")]) +``` +""" +function ForbiddenPath(diagram::InfluenceDiagram, nodes::Vector{Name}, paths::Vector{NTuple{N, Name}}) where N + node_indices = Vector{Node}() + for node in nodes + j = findfirst(i -> i == node, diagram.Names) + if isnothing(j) + throw(DomainError("Node $node does not exist.")) + end + push!(node_indices, j) + end + + path_set = Set{Path}() + for s in paths + s_states = Vector{State}() + for (i, s_i) in enumerate(s) + s_i_index = findfirst(x -> x == s_i, diagram.States[node_indices[i]]) + if isnothing(s_i_index) + throw(DomainError("Node $(nodes[i]) does not have a state called $s_i.")) + end + + push!(s_states, s_i_index) + end + push!(path_set, Path(s_states)) + end + + return ForbiddenPath((node_indices, path_set)) +end + +""" + function FixedPath(diagram::InfluenceDiagram, fixed::Dict{Name, Name}) + +FixedPath outer construction function. Create FixedPath variable. + +# Arguments +- `diagram::InfluenceDiagram`: Influence diagram structure +- `fixed::Dict{Name, Name}`: Dictionary of nodes and their fixed states. Order is node=>state, and both are idefied with their names. + +# Example +```julia +julia> FixedPath(diagram, Dict("R1"=>"high", "R2"=>"high")) +``` +""" +function FixedPath(diagram::InfluenceDiagram, fixed::Dict{Name, Name}) + fixed_paths = Dict{Node, State}() + + for (j, s_j) in fixed + j_index = findfirst(i -> i == j, diagram.Names) + if isnothing(j_index) + throw(DomainError("Node $j does not exist.")) + end + + s_j_index = findfirst(s -> s == s_j, diagram.States[j_index]) + if isnothing(s_j_index) + throw(DomainError("Node $j does not have a state called $s_j.")) + end + push!(fixed_paths, Node(j_index) => State(s_j_index)) + end + + return FixedPath(fixed_paths) end @@ -394,9 +1033,9 @@ Z(s_I) ``` """ struct LocalDecisionStrategy{N} <: AbstractArray{Int, N} - j::Node + d::Node data::Array{Int, N} - function LocalDecisionStrategy(j::Node, data::Array{Int, N}) where N + function LocalDecisionStrategy(d::Node, data::Array{Int, N}) where N if !all(0 ≤ x ≤ 1 for x in data) throw(DomainError("All values x must be 0 ≤ x ≤ 1.")) end @@ -405,7 +1044,7 @@ struct LocalDecisionStrategy{N} <: AbstractArray{Int, N} throw(DomainError("Values should add to one.")) end end - new{N}(j, data) + new{N}(d, data) end end @@ -427,6 +1066,7 @@ end Decision strategy type. """ struct DecisionStrategy - D::Vector{DecisionNode} - Z_j::Vector{LocalDecisionStrategy} + D::Vector{Node} + I_d::Vector{Vector{Node}} + Z_d::Vector{LocalDecisionStrategy} end diff --git a/src/printing.jl b/src/printing.jl index 611ef2ff..b62d05ad 100644 --- a/src/printing.jl +++ b/src/printing.jl @@ -2,38 +2,57 @@ using DataFrames, PrettyTables using StatsBase, StatsBase.Statistics """ - function print_decision_strategy(S::States, Z::DecisionStrategy) + print_decision_strategy(diagram::InfluenceDiagram, Z::DecisionStrategy, state_probabilities::StateProbabilities; show_incompatible_states::Bool = false) Print decision strategy. +# Arguments +- `diagram::InfluenceDiagram`: Influence diagram structure. +- `Z::DecisionStrategy`: Decision strategy structure with optimal decision strategy. +- `state_probabilities::StateProbabilities`: State probabilities structure corresponding to optimal decision strategy. +- `show_incompatible_states::Bool`: Choice to print rows also for incompatible states. + # Examples ```julia -print_decision_strategy(S, Z) +>julia print_decision_strategy(diagram, Z, S_probabilities) ``` """ -function print_decision_strategy(S::States, Z::DecisionStrategy) - for (d, Z_j) in zip(Z.D, Z.Z_j) - a1 = vec(collect(paths(S[d.I_j]))) - a2 = [Z_j(s_I) for s_I in a1] - labels = fill("States", length(a1)) - df = DataFrame(labels = labels, a1 = a1, a2 = a2) - pretty_table(df, ["Nodes", "$((d.I_j...,))", "$(d.j)"]) +function print_decision_strategy(diagram::InfluenceDiagram, Z::DecisionStrategy, state_probabilities::StateProbabilities; show_incompatible_states::Bool = false) + probs = state_probabilities.probs + + for (d, I_d, Z_d) in zip(Z.D, Z.I_d, Z.Z_d) + s_I = vec(collect(paths(diagram.S[I_d]))) + s_d = [Z_d(s) for s in s_I] + + if !isempty(I_d) + informations_states = [join([String(diagram.States[i][s_i]) for (i, s_i) in zip(I_d, s)], ", ") for s in s_I] + decision_probs = [ceil(prod(probs[i][s1] for (i, s1) in zip(I_d, s))) for s in s_I] + decisions = collect(p == 0 ? "--" : diagram.States[d][s] for (s, p) in zip(s_d, decision_probs)) + df = DataFrame(informations_states = informations_states, decisions = decisions) + if !show_incompatible_states + filter!(row -> row.decisions != "--", df) + end + pretty_table(df, ["State(s) of $(join([diagram.Names[i] for i in I_d], ", "))", "Decision in $(diagram.Names[d])"], alignment=:l) + else + df = DataFrame(decisions = diagram.States[d][s_d]) + pretty_table(df, ["Decision in $(diagram.Names[d])"], alignment=:l) + end end end """ - function print_utility_distribution(udist::UtilityDistribution; util_fmt="%f", prob_fmt="%f") + print_utility_distribution(U_distribution::UtilityDistribution; util_fmt="%f", prob_fmt="%f") Print utility distribution. # Examples ```julia -udist = UtilityDistribution(S, P, U, Z) -print_utility_distribution(udist) +>julia U_distribution = UtilityDistribution(diagram, Z) +>julia print_utility_distribution(U_distribution) ``` """ -function print_utility_distribution(udist::UtilityDistribution; util_fmt="%f", prob_fmt="%f") - df = DataFrame(Utility = udist.u, Probability = udist.p) +function print_utility_distribution(U_distribution::UtilityDistribution; util_fmt="%f", prob_fmt="%f") + df = DataFrame(Utility = U_distribution.u, Probability = U_distribution.p) formatters = ( ft_printf(util_fmt, [1]), ft_printf(prob_fmt, [2])) @@ -41,44 +60,46 @@ function print_utility_distribution(udist::UtilityDistribution; util_fmt="%f", p end """ - function print_state_probabilities(sprobs::StateProbabilities, nodes::Vector{Node}; prob_fmt="%f") + print_state_probabilities(sprobs::StateProbabilities, nodes::Vector{Node}; prob_fmt="%f") Print state probabilities with fixed states. # Examples ```julia -sprobs = StateProbabilities(S, P, U, Z) -print_state_probabilities(sprobs, [c.j for c in C]) -print_state_probabilities(sprobs, [d.j for d in D]) +>julia S_probabilities = StateProbabilities(diagram, Z) +>julia print_state_probabilities(S_probabilities, ["R"]) +>julia print_state_probabilities(S_probabilities, ["A"]) ``` """ -function print_state_probabilities(sprobs::StateProbabilities, nodes::Vector{Node}; prob_fmt="%f") - probs = sprobs.probs - fixed = sprobs.fixed +function print_state_probabilities(diagram::InfluenceDiagram, state_probabilities::StateProbabilities, nodes::Vector{Name}; prob_fmt="%f") + node_indices = [findfirst(j -> j==node, diagram.Names) for node in nodes] + + probs = state_probabilities.probs + fixed = state_probabilities.fixed prob(p, state) = if 1≤state≤length(p) p[state] else NaN end fix_state(i) = if i∈keys(fixed) string(fixed[i]) else "" end # Maximum number of states - limit = maximum(length(probs[i]) for i in nodes) + limit = maximum(length(probs[i]) for i in node_indices) states = 1:limit df = DataFrame() df[!, :Node] = nodes for state in states - df[!, Symbol("State $state")] = [prob(probs[i], state) for i in nodes] + df[!, Symbol("State $state")] = [prob(probs[i], state) for i in node_indices] end - df[!, Symbol("Fixed state")] = [fix_state(i) for i in nodes] + df[!, Symbol("Fixed state")] = [fix_state(i) for i in node_indices] pretty_table(df; formatters = ft_printf(prob_fmt, (first(states)+1):(last(states)+1))) end """ -function print_statistics(udist::UtilityDistribution; fmt = "%f") + print_statistics(U_distribution::UtilityDistribution; fmt = "%f") Print statistics about utility distribution. """ -function print_statistics(udist::UtilityDistribution; fmt = "%f") - u = udist.u - w = ProbabilityWeights(udist.p) +function print_statistics(U_distribution::UtilityDistribution; fmt = "%f") + u = U_distribution.u + w = ProbabilityWeights(U_distribution.p) names = ["Mean", "Std", "Skewness", "Kurtosis"] statistics = [mean(u, w), std(u, w, corrected=false), skewness(u, w), kurtosis(u, w)] df = DataFrame(Name = names, Statistics = statistics) @@ -86,14 +107,13 @@ function print_statistics(udist::UtilityDistribution; fmt = "%f") end """ - function print_risk_measures(udist::UtilityDistribution, αs::Vector{Float64}; fmt = "%f") + print_risk_measures(U_distribution::UtilityDistribution, αs::Vector{Float64}; fmt = "%f") Print risk measures. """ -function print_risk_measures(udist::UtilityDistribution, αs::Vector{Float64}; fmt = "%f") - u, p = udist.u, udist.p - VaR = [value_at_risk(u, p, α) for α in αs] - CVaR = [conditional_value_at_risk(u, p, α) for α in αs] +function print_risk_measures(U_distribution::UtilityDistribution, αs::Vector{Float64}; fmt = "%f") + VaR = [value_at_risk(U_distribution, α) for α in αs] + CVaR = [conditional_value_at_risk(U_distribution, α) for α in αs] df = DataFrame(α = αs, VaR = VaR, CVaR = CVaR) pretty_table(df, formatters = ft_printf(fmt)) end diff --git a/src/random.jl b/src/random.jl index 77bc0695..7e5d94c6 100644 --- a/src/random.jl +++ b/src/random.jl @@ -5,9 +5,10 @@ using Random Generates random information sets for chance and decision nodes. """ -function information_set(rng::AbstractRNG, j::Int, n_I::Int) +function information_set(rng::AbstractRNG, j::Node, n_I::Int) m = min(rand(rng, 0:n_I), j-1) - return shuffle(rng, 1:(j-1))[1:m] + I_j = shuffle(rng, 1:(j-1))[1:m] + return sort(I_j) end """ @@ -24,22 +25,36 @@ function information_set(rng::AbstractRNG, leaf_nodes::Vector{Node}, n::Int) else m = rand(rng, 0:l) end - return [leaf_nodes; non_leaf_nodes[1:m]] + I_v = [leaf_nodes; non_leaf_nodes[1:m]] + return sort(I_v) end + """ - function random_diagram(rng::AbstractRNG, n_C::Int, n_D::Int, n_V::Int, m_C::Int, m_D::Int) + function random_diagram!(rng::AbstractRNG, n_C::Int, n_D::Int, n_V::Int, m_C::Int, m_D::Int, states::Vector{Int}) Generate random decision diagram with `n_C` chance nodes, `n_D` decision nodes, and `n_V` value nodes. Parameter `m_C` and `m_D` are the upper bounds for the size of the information set. +# Arguments +- `rng::AbstractRNG`: Random number generator. +- `n_C::Int`: Number of chance nodes. +- `n_D::Int`: Number of decision nodes. +- `n_V::Int`: Number of value nodes. +- `m_C::Int`: Upper bound for size of information set for chance nodes. +- `m_D::Int`: Upper bound for size of information set for decision nodes. +- `states::Vector{State}`: The number of states for each chance and decision node + is randomly chosen from this set of numbers. + + # Examples ```julia rng = MersenneTwister(3) -random_diagram(rng, 5, 2, 3, 2) +diagram = InfluenceDiagram() +random_diagram!(rng, diagram, 5, 2, 3, 2, [2, 4, 5]) ``` """ -function random_diagram(rng::AbstractRNG, n_C::Int, n_D::Int, n_V::Int, m_C::Int, m_D::Int) +function random_diagram!(rng::AbstractRNG, diagram::InfluenceDiagram, n_C::Int, n_D::Int, n_V::Int, m_C::Int, m_D::Int, states::Vector{Int}) n = n_C + n_D n_C ≥ 0 || throw(DomainError("There should be `n_C ≥ 0` chance nodes.")) n_D ≥ 0 || throw(DomainError("There should be `n_D ≥ 0` decision nodes")) @@ -47,44 +62,43 @@ function random_diagram(rng::AbstractRNG, n_C::Int, n_D::Int, n_V::Int, m_C::Int n_V ≥ 1 || throw(DomainError("There should be `n_V ≥ 1` value nodes.")) m_C ≥ 1 || throw(DomainError("Maximum size of information set should be `m_C ≥ 1`.")) m_D ≥ 1 || throw(DomainError("Maximum size of information set should be `m_D ≥ 1`.")) + all(s > 1 for s in states) || throw(DomainError("Minimum number of states possible should be 2.")) # Create node indices U = shuffle(rng, 1:n) - C_j = sort(U[1:n_C]) - D_j = sort(U[(n_C+1):n]) - V_j = collect((n+1):(n+n_V)) + diagram.C = [Node(c) for c in sort(U[1:n_C])] + diagram.D = [Node(d) for d in sort(U[(n_C+1):n])] + diagram.V = [Node(v) for v in collect((n+1):(n+n_V))] + diagram.I_j = Vector{Vector{Node}}(undef, n+n_V) # Create chance and decision nodes - C = [ChanceNode(j, information_set(rng, j, m_C)) for j in C_j] - D = [DecisionNode(j, information_set(rng, j, m_D)) for j in D_j] + for c in diagram.C + diagram.I_j[c] = information_set(rng, c, m_C) + end + for d in diagram.D + diagram.I_j[d] = information_set(rng, d, m_D) + end + # Assign each leaf node to a random value node - leaf_nodes = setdiff(1:n, (c.I_j for c in C)..., (d.I_j for d in D)...) - leaf_nodes_j = Dict(j=>Node[] for j in V_j) - for i in leaf_nodes - k = rand(rng, V_j) - push!(leaf_nodes_j[k], i) + leaf_nodes = setdiff(1:n, (diagram.I_j[c] for c in diagram.C)..., (diagram.I_j[d] for d in diagram.D)...) + leaf_nodes_v = Dict(v=>Node[] for v in diagram.V) + for j in leaf_nodes + v = rand(rng, diagram.V) + push!(leaf_nodes_v[v], j) end # Create values nodes - V = [ValueNode(j, information_set(rng, leaf_nodes_j[j], n)) for j in V_j] - - return C, D, V -end + for v in diagram.V + diagram.I_j[v] = information_set(rng, leaf_nodes_v[v], n) + end -""" - function States(rng::AbstractRNG, states::Vector{State}, n::Int) -Generate `n` random states from `states`. + diagram.S = States(State[rand(rng, states, n)...]) + diagram.X = Vector{Probabilities}(undef, n_C) + diagram.Y = Vector{Utilities}(undef, n_V) -# Examples -```julia -rng = MersenneTwister(3) -S = States(rng, [2, 3], 10) -``` -""" -function States(rng::AbstractRNG, states::Vector{State}, n::Int) - States(rand(rng, states, n)) + return diagram end """ @@ -95,14 +109,19 @@ Generate random probabilities for chance node `c` with `S` states. # Examples ```julia rng = MersenneTwister(3) -c = ChanceNode(2, [1]) -S = States([2, 2]) -Probabilities(rng, c, S) +diagram = InfluenceDiagram() +random_diagram!(rng, diagram, 5, 2, 3, 2, [2, 4, 5]) +c = diagram.C[1] +Probabilities!(rng, diagram, c) ``` """ -function Probabilities(rng::AbstractRNG, c::ChanceNode, S::States; n_inactive::Int=0) - states = S[c.I_j] - state = S[c.j] +function random_probabilities!(rng::AbstractRNG, diagram::InfluenceDiagram, c::Node; n_inactive::Int=0) + if !(c in diagram.C) + throw(DomainError("Probabilities can only be added for chance nodes.")) + end + I_c = diagram.I_j[c] + states = diagram.S[I_c] + state = diagram.S[c] if !(0 ≤ n_inactive ≤ prod([states...; (state - 1)])) throw(DomainError("Number of inactive states must be < prod([S[I_j]...;, S[j]-1])")) end @@ -133,37 +152,48 @@ function Probabilities(rng::AbstractRNG, c::ChanceNode, S::States; n_inactive::I data[s, :] /= sum(data[s, :]) end - Probabilities(c.j, data) + index_c = findfirst(j -> j==c, diagram.C) + diagram.X[index_c] = Probabilities(c, data) end -scale(x::Float64, low::Float64, high::Float64) = x * (high - low) + low +scale(x::Utility, low::Utility, high::Utility) = x * (high - low) + low """ - function Consequences(rng::AbstractRNG, v::ValueNode, S::States; low::Float64=-1.0, high::Float64=1.0) + function Utilities!(rng::AbstractRNG, diagram::InfluenceDiagram, v::Node; low::Float64=-1.0, high::Float64=1.0) -Generate random consequences between `low` and `high` for value node `v` with `S` states. +Generate random utilities between `low` and `high` for value node `v`. # Examples ```julia rng = MersenneTwister(3) -v = ValueNode(3, [1]) -S = States([2, 2]) -Consequences(rng, v, S; low=-1.0, high=1.0) +diagram = InfluenceDiagram() +random_diagram!(rng, diagram, 5, 2, 3, 2, [2, 4, 5]) +v = diagram.V[1] +Utilities!(rng, diagram, v) ``` """ -function Consequences(rng::AbstractRNG, v::ValueNode, S::States; low::Float64=-1.0, high::Float64=1.0) +function random_utilities!(rng::AbstractRNG, diagram::InfluenceDiagram, v::Node; low::Float64=-1.0, high::Float64=1.0) + if !(v in diagram.V) + throw(DomainError("Utilities can only be added for value nodes.")) + end if !(high > low) throw(DomainError("high should be greater than low")) end - data = rand(rng, S[v.I_j]...) - data = scale.(data, low, high) - Consequences(v.j, data) + I_v = diagram.I_j[v] + data = rand(rng, Utility, diagram.S[I_v]...) + data = scale.(data, Utility(low), Utility(high)) + + index_v = findfirst(j -> j==v, diagram.V) + diagram.Y[index_v] = Utilities(v, data) end + + + """ function LocalDecisionStrategy(rng::AbstractRNG, d::DecisionNode, S::States) -Generate random decision strategy for decision node `d` with `S` states. +Generate random decision strategy for decision node `d`. # Examples ```julia @@ -173,13 +203,14 @@ S = States([2, 2]) LocalDecisionStrategy(rng, d, S) ``` """ -function LocalDecisionStrategy(rng::AbstractRNG, d::DecisionNode, S::States) - states = S[d.I_j] - state = S[d.j] +function LocalDecisionStrategy(rng::AbstractRNG, diagram::InfluenceDiagram, d::Node) + I_d = diagram.I_j[d] + states = diagram.S[I_d] + state = diagram.S[d] data = zeros(Int, states..., state) for s in CartesianIndices((states...,)) s_j = rand(rng, 1:state) data[s, s_j] = 1 end - LocalDecisionStrategy(d.j, data) + LocalDecisionStrategy(d, data) end diff --git a/test/decision_model.jl b/test/decision_model.jl index 7a82a782..e6149972 100644 --- a/test/decision_model.jl +++ b/test/decision_model.jl @@ -3,85 +3,84 @@ using DecisionProgramming function influence_diagram(rng::AbstractRNG, n_C::Int, n_D::Int, n_V::Int, m_C::Int, m_D::Int, states::Vector{Int}, n_inactive::Int) - C, D, V = random_diagram(rng, n_C, n_D, n_V, m_C, m_D) - S = States(rng, states, length(C) + length(D)) - X = [Probabilities(rng, c, S; n_inactive=n_inactive) for c in C] - Y = [Consequences(rng, v, S; low=-1.0, high=1.0) for v in V] - - validate_influence_diagram(S, C, D, V) - - s_c = sortperm([c.j for c in C]) - s_d = sortperm([d.j for d in D]) - s_v = sortperm([v.j for v in V]) - C, D, V = C[s_c], D[s_d], V[s_v] - X, Y = X[s_c], Y[s_v] - P = DefaultPathProbability(C, X) - U = DefaultPathUtility(V, Y) - - return D, S, P, U + diagram = InfluenceDiagram() + random_diagram!(rng, diagram, n_C, n_D, n_V, m_C, m_D, states) + for c in diagram.C + random_probabilities!(rng, diagram, c; n_inactive=n_inactive) + end + for v in diagram.V + random_utilities!(rng, diagram, v; low=-1.0, high=1.0) + end + + # Names needed for printing functions only + diagram.Names = ["node$j" for j in 1:n_C+n_D+n_V] + diagram.States = [["s$s" for s in 1:n_s] for n_s in diagram.S] + + diagram.P = DefaultPathProbability(diagram.C, diagram.I_j[diagram.C], diagram.X) + diagram.U = DefaultPathUtility(diagram.I_j[diagram.V], diagram.Y) + + return diagram end -function test_decision_model(D, S, P, U, n_inactive, probability_scale_factor, probability_cut) +function test_decision_model(diagram, n_inactive, probability_scale_factor, probability_cut) model = Model() @info "Testing DecisionVariables" - z = DecisionVariables(model, S, D) + z = DecisionVariables(model, diagram) @info "Testing PathCompatibilityVariables" - x_s = PathCompatibilityVariables(model, z, S, P; probability_cut = probability_cut) - - @info "Testing PositivePathUtility" - U′ = if probability_cut U else PositivePathUtility(S, U) end + x_s = PathCompatibilityVariables(model, diagram, z; probability_cut = probability_cut) @info "Testing probability_cut" - lazy_probability_cut(model, x_s, P) + lazy_probability_cut(model, diagram, x_s) @info "Testing expected_value" if probability_scale_factor > 0 - EV = expected_value(model, x_s, U′, P; probability_scale_factor = probability_scale_factor) + EV = expected_value(model, diagram, x_s; probability_scale_factor = probability_scale_factor) else - @test_throws DomainError expected_value(model, x_s, U′, P; probability_scale_factor = probability_scale_factor) + @test_throws DomainError expected_value(model, diagram, x_s; probability_scale_factor = probability_scale_factor) end @info "Testing conditional_value_at_risk" if probability_scale_factor > 0 - CVaR = conditional_value_at_risk(model, x_s, U′, P, 0.2; probability_scale_factor = probability_scale_factor) + CVaR = conditional_value_at_risk(model, diagram, x_s, 0.2; probability_scale_factor = probability_scale_factor) else - @test_throws DomainError conditional_value_at_risk(model, x_s, U′, P, 0.2; probability_scale_factor = probability_scale_factor) + @test_throws DomainError conditional_value_at_risk(model, diagram, x_s, 0.2; probability_scale_factor = probability_scale_factor) end @test true end -function test_analysis_and_printing(D, S, P, U) +function test_analysis_and_printing(diagram) @info("Creating random decision strategy") - Z_j = [LocalDecisionStrategy(rng, d, S) for d in D] - Z = DecisionStrategy(D, Z_j) + Z_j = [LocalDecisionStrategy(rng, diagram, d) for d in diagram.D] + Z = DecisionStrategy(diagram.D, diagram.I_j[diagram.D], Z_j) @info "Testing CompatiblePaths" - @test all(true for s in CompatiblePaths(S, P.C, Z)) - @test_throws DomainError CompatiblePaths(S, P.C, Z, Dict(D[1].j => 1)) - node, state = (P.C[1].j, 1) - @test all(s[node] == state for s in CompatiblePaths(S, P.C, Z, Dict(node => state))) + @test all(true for s in CompatiblePaths(diagram, Z)) + @test_throws DomainError CompatiblePaths(diagram, Z, Dict(diagram.D[1] => State(1))) + node, state = (diagram.C[1], State(1)) + @test all(s[node] == state for s in CompatiblePaths(diagram, Z, Dict(node => state))) @info "Testing UtilityDistribution" - udist = UtilityDistribution(S, P, U, Z) + U_distribution = UtilityDistribution(diagram, Z) @info "Testing StateProbabilities" - sprobs = StateProbabilities(S, P, Z) + S_probabilities = StateProbabilities(diagram, Z) @info "Testing conditional StateProbabilities" - sprobs2 = StateProbabilities(S, P, Z, node, state, sprobs) - - @info "Testing " - print_decision_strategy(S, Z) - print_utility_distribution(udist) - print_state_probabilities(sprobs, [c.j for c in P.C]) - print_state_probabilities(sprobs, [d.j for d in D]) - print_state_probabilities(sprobs2, [c.j for c in P.C]) - print_state_probabilities(sprobs2, [d.j for d in D]) - print_statistics(udist) - print_risk_measures(udist, [0.0, 0.05, 0.1, 0.2, 1.0]) + S_probabilities2 = StateProbabilities(diagram, Z, node, state, S_probabilities) + + @info "Testing printing functions" + print_decision_strategy(diagram, Z, S_probabilities) + print_decision_strategy(diagram, Z, S_probabilities, show_incompatible_states=true) + print_utility_distribution(U_distribution) + print_state_probabilities(diagram, S_probabilities, [diagram.Names[c] for c in diagram.C]) + print_state_probabilities(diagram, S_probabilities, [diagram.Names[d] for d in diagram.D]) + print_state_probabilities(diagram, S_probabilities2, [diagram.Names[c] for c in diagram.C]) + print_state_probabilities(diagram, S_probabilities2, [diagram.Names[d] for d in diagram.D]) + print_statistics(U_distribution) + print_risk_measures(U_distribution, [0.0, 0.05, 0.1, 0.2, 1.0]) @test true end @@ -89,13 +88,13 @@ end @info "Testing model construction" rng = MersenneTwister(4) for (n_C, n_D, states, n_inactive, probability_scale_factor, probability_cut) in [ - (3, 2, [1, 2, 3], 0, 1.0, true), - (3, 2, [1, 2, 3], 0, -1.0, true), + (3, 2, [2, 3, 4], 0, 1.0, true), + (3, 2, [2, 3], 0, -1.0, true), (3, 2, [3], 1, 100.0, true), - (3, 2, [1, 2, 3], 0, -1.0, false), - (3, 2, [3], 1, 10.0, false) + (3, 2, [2, 3], 0, -1.0, false), + (3, 2, [4], 1, 10.0, false) ] - D, S, P, U = influence_diagram(rng, n_C, n_D, 2, 2, 2, states, n_inactive) - test_decision_model(D, S, P, U, n_inactive, probability_scale_factor, probability_cut) - test_analysis_and_printing(D, S, P, U) + diagram = influence_diagram(rng, n_C, n_D, 2, 2, 2, states, n_inactive) + test_decision_model(diagram, n_inactive, probability_scale_factor, probability_cut) + test_analysis_and_printing(diagram) end diff --git a/test/influence_diagram.jl b/test/influence_diagram.jl index 26bf673d..8333f3ac 100644 --- a/test/influence_diagram.jl +++ b/test/influence_diagram.jl @@ -2,90 +2,186 @@ using Test, Logging, Random using DecisionProgramming @info "Testing ChandeNode" -@test isa(ChanceNode(1, Node[]), ChanceNode) -@test isa(ChanceNode(2, Node[1]), ChanceNode) -@test_throws DomainError ChanceNode(1, Node[1]) -@test_throws DomainError ChanceNode(1, Node[2]) +@test isa(ChanceNode("A", [], ["a", "b"]), ChanceNode) +@test isa(ChanceNode("B", [], ["x", "y", "z"]), ChanceNode) +@test_throws MethodError ChanceNode(1, [], ["x", "y", "z"]) +@test_throws MethodError ChanceNode("B", [1], ["y", "z"]) +@test_throws MethodError ChanceNode("B", ["A"], [1, "y", "z"]) +@test_throws MethodError ChanceNode("B", ["A"]) @info "Testing DecisionNode" -@test isa(DecisionNode(1, Node[]), DecisionNode) -@test isa(DecisionNode(2, Node[1]), DecisionNode) -@test_throws DomainError DecisionNode(1, Node[1]) -@test_throws DomainError DecisionNode(1, Node[2]) +@test isa(DecisionNode("D", [], ["x", "y"]), DecisionNode) +@test isa(DecisionNode("E", ["C"], ["x", "y"]), DecisionNode) +@test_throws MethodError DecisionNode(1, [], ["x", "y", "z"]) +@test_throws MethodError DecisionNode("D", [1], ["y", "z"]) +@test_throws MethodError DecisionNode("D", ["A"], [1, "y", "z"]) +@test_throws MethodError DecisionNode("D", ["A"]) @info "Testing ValueNode" -@test isa(ValueNode(1, Node[]), ValueNode) -@test isa(ValueNode(2, Node[1]), ValueNode) -@test_throws DomainError ValueNode(1, Node[1]) -@test_throws DomainError ValueNode(1, Node[2]) +@test isa(ValueNode("V", []), ValueNode) +@test isa(ValueNode("V", ["E", "D"]), ValueNode) +@test_throws MethodError ValueNode(1, []) +@test_throws MethodError ValueNode("V", [2]) @info "Testing State" -@test isa(States([1, 2, 3]), States) -@test_throws DomainError States([0, 1]) -@test States([(2, [1, 3]), (3, [2, 4, 5])]) == States([2, 3, 2, 3, 3]) - -@info "Testing validate_influence_diagram" -@test_throws DomainError validate_influence_diagram( - States([1]), - [ChanceNode(1, Node[])], - [DecisionNode(2, Node[])], - [ValueNode(3, [1, 2])] -) -@test_throws DomainError validate_influence_diagram( - States([1, 1]), - [ChanceNode(1, Node[])], - [DecisionNode(1, Node[])], - [ValueNode(3, [1, 2])] -) -@test_throws DomainError validate_influence_diagram( - States([1, 1]), - [ChanceNode(1, Node[])], - [DecisionNode(2, Node[])], - [ValueNode(2, [1])] -) -@test_throws DomainError validate_influence_diagram( - States([1, 1]), - [ChanceNode(1, Node[])], - [DecisionNode(2, Node[])], - [ValueNode(3, [2]), ValueNode(4, [3])] -) -# Test redundancy -@test validate_influence_diagram( - States([1, 1]), - [ChanceNode(1, Node[])], - [DecisionNode(2, Node[])], - [ValueNode(3, Node[])] -) === nothing +@test isa(States(State[1, 2, 3]), States) +@test_throws DomainError States(State[0, 1]) @info "Testing paths" -@test vec(collect(paths(States([2, 3])))) == [(1, 1), (2, 1), (1, 2), (2, 2), (1, 3), (2, 3)] -@test vec(collect(paths(States([2, 3]), Dict(1=>2)))) == [(2, 1), (2, 2), (2, 3)] +@test vec(collect(paths(States(State[2, 3])))) == [(1, 1), (2, 1), (1, 2), (2, 2), (1, 3), (2, 3)] +@test vec(collect(paths(States(State[2, 3]), Dict(Node(1)=>State(2))))) == [(2, 1), (2, 2), (2, 3)] +@test vec(collect(paths(States(State[2, 3]), FixedPath(Dict(Node(1)=>State(2)))))) == [(2, 1), (2, 2), (2, 3)] @info "Testing Probabilities" -@test isa(Probabilities(1, [0.4 0.6; 0.3 0.7]), Probabilities) -@test isa(Probabilities(1, [0.0, 0.4, 0.6]), Probabilities) -@test_throws DomainError Probabilities(1, [1.1, 0.1]) +@test isa(Probabilities(Node(1), [0.4 0.6; 0.3 0.7]), Probabilities) +@test isa(Probabilities(Node(1), [0.0, 0.4, 0.6]), Probabilities) +@test_throws DomainError Probabilities(Node(1), [1.1, 0.1]) @info "Testing DefaultPathProbability" P = DefaultPathProbability( - [ChanceNode(1, Node[]), ChanceNode(2, [1])], - [Probabilities(1, [0.4, 0.6]), Probabilities(2, [0.3 0.7; 0.9 0.1])] + [Node(1), Node(2)], + [Node[], [Node(1)]], + [Probabilities(Node(1), [0.4, 0.6]), Probabilities(Node(2), [0.3 0.7; 0.9 0.1])] ) @test isa(P, DefaultPathProbability) -@test P((1, 2)) == 0.4 * 0.7 +@test P((State(1), State(2))) == 0.4 * 0.7 -@info "Testing Consequences" -@test isa(Consequences(1, [-1.1, 0.0, 2.7]), Consequences) -@test isa(Consequences(1, [-1.1 0.0; 2.7 7.0]), Consequences) +@info "Testing Utilities" +@test isa(Utilities(Node(1), Utility[-1.1, 0.0, 2.7]), Utilities) +@test isa(Utilities(Node(1), Utility[-1.1 0.0; 2.7 7.0]), Utilities) @info "Testing DefaultPathUtility" U = DefaultPathUtility( - [ValueNode(3, [2]), ValueNode(4, [1, 2])], - [Consequences(3, [1.0, 1.4]), Consequences(4, [1.0 1.5; 0.6 3.4])] + [Node[2], Node[1, 2]], + [Utilities(Node(3), Utility[1.0, 1.4]), Utilities(Node(4), Utility[1.0 1.5; 0.6 3.4])] ) @test isa(U, DefaultPathUtility) -@test U((2, 1)) == 1.0 + 0.6 +@test U((State(2), State(1))) == Utility(1.0 + 0.6) + +@info "Testing InfluenceDiagram" +diagram = InfluenceDiagram() +@test isa(diagram, InfluenceDiagram) +push!(diagram.Nodes, ChanceNode("A", [], ["a", "b"])) +@test isa(diagram, InfluenceDiagram) + +@info "Testing add_node! and validate_node" +diagram = InfluenceDiagram() +add_node!(diagram, ChanceNode("A", [], ["a", "b"])) +@test isa(diagram.Nodes[1], ChanceNode) +@test_throws DomainError add_node!(diagram, ChanceNode("A", [], ["a", "b"])) +@test_throws DomainError add_node!(diagram, ChanceNode("C", ["B", "B"], ["a", "b"])) +@test_throws DomainError add_node!(diagram, ChanceNode("C", ["C", "B"], ["a", "b"])) +@test_throws DomainError add_node!(diagram, ChanceNode("C", [], ["a"])) +@test_throws DomainError add_node!(diagram, DecisionNode("A", [], ["a", "b"])) +@test_throws DomainError add_node!(diagram, DecisionNode("C", ["B", "B"], ["a", "b"])) +@test_throws DomainError add_node!(diagram, DecisionNode("C", ["C", "B"], ["a", "b"])) +@test_throws DomainError add_node!(diagram, DecisionNode("C", [], ["a"])) +@test_throws DomainError add_node!(diagram, ValueNode("A", [])) +@test_throws DomainError add_node!(diagram, ValueNode("C", ["B", "B"])) +@test_throws DomainError add_node!(diagram, ValueNode("C", ["C", "B"])) +add_node!(diagram, ChanceNode("C", ["A"], ["a", "b"])) +add_node!(diagram, DecisionNode("D", ["A"], ["c", "d"])) +add_node!(diagram, ValueNode("V", ["A"])) +@test length(diagram.Nodes) == 4 +@test isa(diagram.Nodes[3], DecisionNode) +@test isa(diagram.Nodes[4], ValueNode) + +@info "Testing generate_arcs!" +diagram = InfluenceDiagram() +@test_throws DomainError generate_arcs!(diagram) +add_node!(diagram, ChanceNode("A", [], ["a", "b"])) +add_node!(diagram, ChanceNode("C", ["A"], ["a", "b", "c"])) +add_node!(diagram, ValueNode("V", ["A", "C"])) +generate_arcs!(diagram) +@test diagram.Names == ["A", "C", "V"] +@test diagram.I_j == [[], Node[1], Node[1, 2]] +@test diagram.States == [["a", "b"], ["a", "b", "c"]] +@test diagram.S == [State(2), State(3)] +@test diagram.C == Node[1, 2] +@test diagram.D == Node[] +@test diagram.V == Node[3] +@test diagram.X == Probabilities[] +@test diagram.Y == Utilities[] + +#Non-existent node B +diagram = InfluenceDiagram() +add_node!(diagram, ChanceNode("A", [], ["a", "b"])) +add_node!(diagram, ChanceNode("C", ["B"], ["a", "b", "c"])) +add_node!(diagram, ValueNode("V", ["A", "C"])) +@test_throws DomainError generate_arcs!(diagram) + +#Cylic +diagram = InfluenceDiagram() +add_node!(diagram, ChanceNode("A", ["R"], ["a", "b"])) +add_node!(diagram, ChanceNode("R", ["C", "A"], ["a", "b", "c"])) +add_node!(diagram, ChanceNode("C", ["A"], ["a", "b", "c"])) +add_node!(diagram, ValueNode("V", ["A", "C"])) +@test_throws DomainError generate_arcs!(diagram) + +#Value node in I_j +diagram = InfluenceDiagram() +add_node!(diagram, ChanceNode("A", ["R"], ["a", "b"])) +add_node!(diagram, ChanceNode("R", ["C", "A"], ["a", "b", "c"])) +add_node!(diagram, ChanceNode("C", ["A", "V"], ["a", "b", "c"])) +add_node!(diagram, ValueNode("V", ["A"])) +@test_throws DomainError generate_arcs!(diagram) + +@info "Testing ProbabilityMatrix" +diagram = InfluenceDiagram() +add_node!(diagram, ChanceNode("A", [], ["a", "b"])) +add_node!(diagram, ChanceNode("B", ["A"], ["a", "b", "c"])) +add_node!(diagram, DecisionNode("D", ["A"], ["a", "b", "c"])) +generate_arcs!(diagram) +@test ProbabilityMatrix(diagram, "A") == zeros(2) +@test ProbabilityMatrix(diagram, "B") == zeros(2, 3) +@test_throws DomainError ProbabilityMatrix(diagram, "C") +@test_throws DomainError ProbabilityMatrix(diagram, "D") +X_A = ProbabilityMatrix(diagram, "A") +set_probability!(X_A, ["a"], 0.2) +@test X_A == [0.2, 0] +set_probability!(X_A, ["b"], 0.9) +@test X_A == [0.2, 0.9] +@test_throws DomainError add_probabilities!(diagram, "A", X_A) +set_probability!(X_A, ["b"], 0.8) +@test add_probabilities!(diagram, "A", X_A) == [[0.2, 0.8]] +@test_throws DomainError add_probabilities!(diagram, "A", X_A) + +@info "Testing UtilityMatrix" +diagram = InfluenceDiagram() +add_node!(diagram, ChanceNode("A", [], ["a", "b"])) +add_node!(diagram, DecisionNode("D", ["A"], ["a", "b", "c"])) +add_node!(diagram, ValueNode("V", ["A", "D"])) +generate_arcs!(diagram) +@test UtilityMatrix(diagram, "V") == fill(Inf, (2, 3)) +@test_throws DomainError UtilityMatrix(diagram, "C") +@test_throws DomainError UtilityMatrix(diagram, "D") +Y_V = UtilityMatrix(diagram, "V") +@test_throws DomainError add_utilities!(diagram, "V", Y_V) +set_utility!(Y_V, ["a", :], [1, 2, 3]) +set_utility!(Y_V, ["b", "c"], 4) +set_utility!(Y_V, ["b", "a"], 5) +set_utility!(Y_V, ["b", "b"], 6) +@test Y_V == [1 2 3; 5 6 4] +add_utilities!(diagram, "V", Y_V) +@test diagram.Y == [[1 2 3; 5 6 4]] +@test_throws DomainError add_utilities!(diagram, "V", Y_V) + + +@info "Testing generate_diagram!" +diagram = InfluenceDiagram() +add_node!(diagram, ChanceNode("A", [], ["a", "b"])) +add_node!(diagram, DecisionNode("D", ["A"], ["a", "b", "c"])) +add_node!(diagram, ValueNode("V", ["A", "D"])) +generate_arcs!(diagram) +add_utilities!(diagram, "V", [-1 2 3; 5 6 4]) +add_probabilities!(diagram, "A", [0.2, 0.8]) +generate_diagram!(diagram) +@test diagram.translation == Utility(0) -@info "Testing LocalDecisionStrategy" -@test_throws DomainError LocalDecisionStrategy(1, [0, 0, 2]) -@test_throws DomainError LocalDecisionStrategy(1, [0, 1, 1]) +@info "Testing positive and negative path utility translations" +generate_diagram!(diagram, positive_path_utility=true) +@test diagram.translation == Utility(2) +@test all(diagram.U(s, diagram.translation) > 0 for s in paths(diagram.S)) +generate_diagram!(diagram, negative_path_utility=true) +@test diagram.translation == Utility(-7) +@test all(diagram.U(s, diagram.translation) < 0 for s in paths(diagram.S)) diff --git a/test/random.jl b/test/random.jl index 29c31be2..e87e9275 100644 --- a/test/random.jl +++ b/test/random.jl @@ -4,44 +4,58 @@ using DecisionProgramming rng = MersenneTwister(4) @info "Testing random_diagram" -@test_throws DomainError random_diagram(rng, -1, 1, 1, 1, 1) -@test_throws DomainError random_diagram(rng, 1, -1, 1, 1, 1) -@test_throws DomainError random_diagram(rng, 0, 0, 1, 1, 1) -@test_throws DomainError random_diagram(rng, 1, 1, 0, 1, 1) -@test_throws DomainError random_diagram(rng, 1, 1, 1, 0, 1) -@test_throws DomainError random_diagram(rng, 1, 1, 1, 1, 0) +diagram = InfluenceDiagram() +@test_throws DomainError random_diagram!(rng, diagram, -1, 1, 1, 1, 1, [2]) +@test_throws DomainError random_diagram!(rng, diagram, 1, -1, 1, 1, 1, [2]) +@test_throws DomainError random_diagram!(rng, diagram, 0, 0, 1, 1, 1, [2]) +@test_throws DomainError random_diagram!(rng, diagram, 1, 1, 0, 1, 1, [2]) +@test_throws DomainError random_diagram!(rng, diagram, 1, 1, 1, 0, 1, [2]) +@test_throws DomainError random_diagram!(rng, diagram, 1, 1, 1, 1, 0, [2]) +@test_throws DomainError random_diagram!(rng, diagram, 1, 1, 1, 1, 1, [1]) for (n_C, n_D) in [(1, 0), (0, 1)] rng = RandomDevice() - C, D, V = random_diagram(rng, n_C, n_D, 1, 1, 1) - @test isa(C, Vector{ChanceNode}) - @test isa(D, Vector{DecisionNode}) - @test isa(V, Vector{ValueNode}) - @test length(C) == n_C - @test length(D) == n_D - @test length(V) == 1 - @test all(!isempty(v.I_j) for v in V) + diagram = InfluenceDiagram() + random_diagram!(rng, diagram, n_C, n_D, 1, 1, 1, [2]) + @test isa(diagram.C, Vector{Node}) + @test isa(diagram.D, Vector{Node}) + @test isa(diagram.V, Vector{Node}) + @test length(diagram.C) == n_C + @test length(diagram.D) == n_D + @test length(diagram.V) == 1 + @test all(!isempty(I_v) for I_v in diagram.I_j[diagram.V]) + @test isa(diagram.S, States) end -@info "Testing random States" -@test_throws DomainError States(rng, [0], 10) -@test isa(States(rng, [2, 3], 10), States) - @info "Testing random Probabilities" -S = States([2, 3, 2]) -c = ChanceNode(3, [1, 2]) -@test isa(Probabilities(rng, c, S; n_inactive=0), Probabilities) -@test isa(Probabilities(rng, c, S; n_inactive=1), Probabilities) -@test isa(Probabilities(rng, c, S; n_inactive=2*3*(2-1)), Probabilities) -@test_throws DomainError Probabilities(rng, c, S; n_inactive=2*3*(2-1)+1) - -@info "Testing random Consequences" -S = States([2, 3]) -v = ValueNode(3, [1, 2]) -@test isa(Consequences(rng, v, S; low=-1.0, high=1.0), Consequences) -@test_throws DomainError Consequences(rng, v, S; low=1.1, high=1.0) +diagram = InfluenceDiagram() +random_diagram!(rng, diagram, 2, 2, 1, 1, 1, [2]) +@test_throws DomainError random_probabilities!(rng, diagram, diagram.D[1]; n_inactive=0) +random_probabilities!(rng, diagram, diagram.C[1]; n_inactive=0) +@test isa(diagram.X[1], Probabilities) +random_probabilities!(rng, diagram, diagram.C[2]; n_inactive=1) +@test isa(diagram.X[2], Probabilities) + +diagram = InfluenceDiagram() +diagram.C = Node[1,3] +diagram.D = Node[2] +diagram.I_j = [Node[], Node[], Node[1,2]] +diagram.S = States(State[2, 3, 2]) +diagram.X = Vector{Probabilities}(undef, 2) +random_probabilities!(rng, diagram, diagram.C[2]; n_inactive=2*3*(2-1)) +@test isa(diagram.X[2], Probabilities) +@test_throws DomainError random_probabilities!(rng, diagram, diagram.C[2]; n_inactive=2*3*(2-1)+1) + + +@info "Testing random Utilities" +diagram = InfluenceDiagram() +random_diagram!(rng, diagram, 1, 1, 2, 1, 1, [2]) +random_utilities!(rng, diagram, diagram.V[1], low=-1.0, high=1.0) +@test isa(diagram.Y[1], Utilities) +@test_throws DomainError random_utilities!(rng, diagram, diagram.V[1], low=1.1, high=1.0) +@test_throws DomainError random_utilities!(rng, diagram, diagram.C[1], low=-1.0, high=1.0) @info "Testing random LocalDecisionStrategy" -S = States([2, 3, 2]) -d = DecisionNode(3, [1, 2]) -@test isa(LocalDecisionStrategy(rng, d, S), LocalDecisionStrategy) +diagram = InfluenceDiagram() +random_diagram!(rng, diagram, 2, 2, 2, 2, 2, [2]) +@test isa(LocalDecisionStrategy(rng, diagram, diagram.D[1]), LocalDecisionStrategy)