Add CHD example comparison

gamma-opt · Aug 16, 2024 · 27d28f7 · 27d28f7
1 parent 459e7e1
commit 27d28f7
Show file tree

Hide file tree

Showing 5 changed files with 362 additions and 0 deletions.
diff --git a/experiments/CHD.jl b/experiments/CHD.jl
@@ -0,0 +1,41 @@
+using DelimitedFiles, Plots, Statistics
+t1 = zeros(101)
+t2 = zeros(101)
+val1 = zeros(101)
+val2 = zeros(101)
+for i in 0:100
+    try
+        result = readdlm((@__DIR__)*"/results/CHD_"*string(i)*".csv", ',')
+        t1[i+1] = result[2]
+        val1[i+1] = result[1]
+    catch e
+        println("Missing file "*(@__DIR__)*"/results/CHD_"*string(i)*".csv")
+    end
+    try
+        result = readdlm((@__DIR__)*"/../../DecProg_0_1_0/experiments/results/CHD_"*string(i+1)*".csv", ',')
+        t2[i+1] = result[2]
+        val2[i+1] = result[1]
+        # println(val1[i+1]*1000/result[1])
+    catch e
+        t2[i+1] = 7200
+        val2[i+1] = Inf
+    end
+end
+
+#println(sort(t1))
+#println(sort(t2))
+
+# log-axis, all data
+# scatter(t1, t2, xlabel="v1.2.0", ylabel="v0.1.0", markershape=:x, label="Solution time (seconds)", xlim=[0.001,1800], ylim=[0.001,7200], xaxis=:log, yaxis=:log)
+# plot!([0.001,7200], [0.001,7200], ls=:dash, lc=:black, label=false)
+
+# log-axis, >1s
+# scatter(t1, t2, xlabel="v1.2.0", ylabel="v0.1.0", markershape=:x, label="Solution time (seconds)", xlim=[1,1800], ylim=[1,7200], xaxis=:log, yaxis=:log)
+# plot!([1,7200], [1,7200], ls=:dash, lc=:black, label=false)
+
+# "regular" axis
+scatter(t1[t2.<7200], t2[t2.<7200], xlabel="v1.2.0", ylabel="v0.1.0", markershape=:x, label="Solution time (seconds)", xlim=[0,900], ylim=[0,7200], legend=:right)
+scatter!(t1[t2.>=7200], t2[t2.>=7200], markershape=:x, label=false)
+plot!([0,7200], [0,7200], ls=:dash, lc=:black, label=false)
+
+Plots.pdf("CHD")
diff --git a/experiments/CHD.pdf b/experiments/CHD.pdf
diff --git a/experiments/CHD_preventative_care_data.csv b/experiments/CHD_preventative_care_data.csv
@@ -0,0 +1,102 @@
+risk_levels,prior_risk_distribution,TRS_if_healthy,TRS_if_sick,GRS_if_healthy,GRS_if_sick
+0,0.11293990700000,0.39740624000000,0.04893543700000,0.02252611500000,0.00283870400000
+0.01,0.12465484400000,0.14375880600000,0.03547850500000,0.07311748200000,0.03579809800000
+0.02,0.09890158000000,0.10569106600000,0.05107539100000,0.12684417900000,0.10940458400000
+0.03,0.07045429200000,0.06429586300000,0.05294659100000,0.18990717200000,0.11918963900000
+0.04,0.04390572300000,0.05056367600000,0.05803457100000,0.11897423000000,0.13241767500000
+0.05,0.06083044000000,0.03052270900000,0.04297088900000,0.11373319200000,0.13705796900000
+0.06,0.03031388200000,0.02095879400000,0.04511271500000,0.08843425800000,0.09294875300000
+0.07,0.03966468200000,0.03275560100000,0.03669923600000,0.04883144600000,0.08072231700000
+0.08,0.02749771500000,0.03116940200000,0.02995658000000,0.03414356500000,0.06469767600000
+0.09,0.03349420200000,0.02003049500000,0.03945800100000,0.02638107900000,0.03036662900000
+0.1,0.01894239500000,0.00655804300000,0.03828682200000,0.01679492900000,0.03121447300000
+0.11,0.01655726700000,0.00880561100000,0.02568766200000,0.01527908700000,0.02907518800000
+0.12,0.02527272300000,0.00632351800000,0.03223956000000,0.02132732600000,0.01200694000000
+0.13,0.02489687600000,0.00062280900000,0.02773266400000,0.01853267700000,0.01152809400000
+0.14,0.01810201800000,0.00678644900000,0.03350329900000,0.00744960200000,0.01785592900000
+0.15,0.00000000000000,0.00833391300000,0.02556680800000,0.00742368700000,0.01787462400000
+0.16,0.01774679900000,0.00060849200000,0.01627451000000,0.01190778400000,0.01298614400000
+0.17,0.02056393700000,0.00288562700000,0.02217757000000,0.00953420300000,0.01666081400000
+0.18,0.01324023100000,0.00270084800000,0.02124150300000,0.00576241100000,0.01480715600000
+0.19,0.01531385100000,0.00589751000000,0.02479449000000,0.01144990600000,0.01115178600000
+0.2,0.01733562400000,0.00335190500000,0.02195192200000,0.00876195700000,0.01413693600000
+0.21,0.02116504900000,0.00402482500000,0.01956580500000,0.00895109400000,0.00413693600000
+0.22,0.01000000000000,0.00167715800000,0.01787645700000,0.00831543600000,0.00112293600000
+0.23,0.01601149200000,0.00307826100000,0.02319705100000,0.00539759000000,0.00000000000000
+0.24,0.00000000000000,0.00264528400000,0.01020576700000,0.00014600000000,0.00000000000000
+0.25,0.01846514000000,0.00356238400000,0.01657776100000,0.00003414600000,0.00000000000000
+0.26,0.00000000000000,0.00206032000000,0.01563507400000,0.00001677000000,0.00000000000000
+0.27,0.00694367200000,0.00181872700000,0.01312182500000,0.00002267700000,0.00000000000000
+0.28,0.00297060000000,0.00063827400000,0.01557198900000,0.00000000000000,0.00000000000000
+0.29,0.00469651100000,0.00374629400000,0.01719464600000,0.00000000000000,0.00000000000000
+0.3,0.00544617600000,0.00124890500000,0.01135377200000,0.00000000000000,0.00000000000000
+0.31,0.00123335500000,0.00149568700000,0.01212877300000,0.00000000000000,0.00000000000000
+0.32,0.00378659900000,0.00041196500000,0.01390704100000,0.00000000000000,0.00000000000000
+0.33,0.00830991800000,0.00142893800000,0.01750051600000,0.00000000000000,0.00000000000000
+0.34,0.00000000000000,0.00043471200000,0.00608215400000,0.00000000000000,0.00000000000000
+0.35,0.00939161600000,0.00056847000000,0.00733665200000,0.00000000000000,0.00000000000000
+0.36,0.00000000000000,0.00016719700000,0.00514977400000,0.00000000000000,0.00000000000000
+0.37,0.00000000000000,0.00020063600000,0.00276073400000,0.00000000000000,0.00000000000000
+0.38,0.00000000000000,0.00016719700000,0.00109835500000,0.00000000000000,0.00000000000000
+0.39,0.00161001700000,0.00253184700000,0.00279437100000,0.00000000000000,0.00000000000000
+0.4,0.00000000000000,0.00087930000000,0.00135650700000,0.00000000000000,0.00000000000000
+0.41,0.00000000000000,0.00122641700000,0.00150422700000,0.00000000000000,0.00000000000000
+0.42,0.00624891500000,0.00088858100000,0.00147217300000,0.00000000000000,0.00000000000000
+0.43,0.00000000000000,0.00013375800000,0.00251989500000,0.00000000000000,0.00000000000000
+0.44,0.00000000000000,0.00086195100000,0.00391024100000,0.00000000000000,0.00000000000000
+0.45,0.00597112500000,0.00235978300000,0.00372377600000,0.00000000000000,0.00000000000000
+0.46,0.00000000000000,0.00124009800000,0.00280890000000,0.00000000000000,0.00000000000000
+0.47,0.00176904300000,0.00013375800000,0.00000000000000,0.00000000000000,0.00000000000000
+0.48,0.00000000000000,0.00026751500000,0.00147103500000,0.00000000000000,0.00000000000000
+0.49,0.00323970300000,0.00004749700000,0.00148041900000,0.00000000000000,0.00000000000000
+0.5,0.00000000000000,0.00010031800000,0.00104513800000,0.00000000000000,0.00000000000000
+0.51,0.00000000000000,0.00190945300000,0.00116669500000,0.00000000000000,0.00000000000000
+0.52,0.00000000000000,0.00198581600000,0.00399394800000,0.00000000000000,0.00000000000000
+0.53,0.00000000000000,0.00006690000000,0.00129514400000,0.00000000000000,0.00000000000000
+0.54,0.00000000000000,0.00005890000000,0.00009660000000,0.00000000000000,0.00000000000000
+0.55,0.00430848100000,0.00029373000000,0.00000000000000,0.00000000000000,0.00000000000000
+0.56,0.00402460400000,0.00003340000000,0.00025752700000,0.00000000000000,0.00000000000000
+0.57,0.00000000000000,0.00217202900000,0.00269709800000,0.00000000000000,0.00000000000000
+0.58,0.00000000000000,0.00000000000000,0.00009690000000,0.00000000000000,0.00000000000000
+0.59,0.00000000000000,0.00205186700000,0.00107007000000,0.00000000000000,0.00000000000000
+0.6,0.00119385400000,0.00000000000000,0.00015545500000,0.00000000000000,0.00000000000000
+0.61,0.00000000000000,0.00027794900000,0.00016539900000,0.00000000000000,0.00000000000000
+0.62,0.00000000000000,0.00005370000000,0.00018467900000,0.00000000000000,0.00000000000000
+0.63,0.00346794400000,0.00003350000000,0.00003030000000,0.00000000000000,0.00000000000000
+0.64,0.00000000000000,0.00003340000000,0.00000000000000,0.00000000000000,0.00000000000000
+0.65,0.00658831300000,0.00005210000000,0.00378400900000,0.00000000000000,0.00000000000000
+0.66,0.00000000000000,0.00000000000000,0.00000000000000,0.00000000000000,0.00000000000000
+0.67,0.00000000000000,0.00004440000000,0.00000000000000,0.00000000000000,0.00000000000000
+0.68,0.00000000000000,0.00047234300000,0.00000000000000,0.00000000000000,0.00000000000000
+0.69,0.00340850900000,0.00000000000000,0.00000000000000,0.00000000000000,0.00000000000000
+0.7,0.00197460100000,0.00000000000000,0.00000000000000,0.00000000000000,0.00000000000000
+0.71,0.00000000000000,0.00003750000000,0.00000000000000,0.00000000000000,0.00000000000000
+0.72,0.00000000000000,0.00001110000000,0.00000000000000,0.00000000000000,0.00000000000000
+0.73,0.00000000000000,0.00001170000000,0.00000000000000,0.00000000000000,0.00000000000000
+0.74,0.00000000000000,0.00003340000000,0.00000000000000,0.00000000000000,0.00000000000000
+0.75,0.00386885800000,0.00000139000000,0.00000000000000,0.00000000000000,0.00000000000000
+0.76,0.00000000000000,0.00003340000000,0.00000000000000,0.00000000000000,0.00000000000000
+0.77,0.00000000000000,0.00000000000000,0.00000529000000,0.00000000000000,0.00000000000000
+0.78,0.00000000000000,0.00000000000000,0.00010941900000,0.00000000000000,0.00000000000000
+0.79,0.00000000000000,0.00000000000000,0.00000000000000,0.00000000000000,0.00000000000000
+0.8,0.00334460600000,0.00000000000000,0.00000000000000,0.00000000000000,0.00000000000000
+0.81,0.00000000000000,0.00003430000000,0.00000000000000,0.00000000000000,0.00000000000000
+0.82,0.00000000000000,0.00004220000000,0.00132149500000,0.00000000000000,0.00000000000000
+0.83,0.00000000000000,0.00005170000000,0.00128714500000,0.00000000000000,0.00000000000000
+0.84,0.00000000000000,0.00003340000000,0.00028930000000,0.00000000000000,0.00000000000000
+0.85,0.00000000000000,0.00003340000000,0.00029201000000,0.00000000000000,0.00000000000000
+0.86,0.00425751100000,0.00003340000000,0.00063885500000,0.00000000000000,0.00000000000000
+0.87,0.00000000000000,0.00002340000000,0.00000000000000,0.00000000000000,0.00000000000000
+0.88,0.00000000000000,0.00000000000000,0.00058710800000,0.00000000000000,0.00000000000000
+0.89,0.00000000000000,0.00000000000000,0.00000000000000,0.00000000000000,0.00000000000000
+0.9,0.00000000000000,0.00000658000000,0.00000000000000,0.00000000000000,0.00000000000000
+0.91,0.00000000000000,0.00000020900000,0.00000000000000,0.00000000000000,0.00000000000000
+0.92,0.00000000000000,0.00000000000000,0.00000000000000,0.00000000000000,0.00000000000000
+0.93,0.00238340700000,0.00000000000000,0.00000000000000,0.00000000000000,0.00000000000000
+0.94,0.00000000000000,0.00000000000000,0.00000000000000,0.00000000000000,0.00000000000000
+0.95,0.00290569600000,0.00000000000000,0.00000000000000,0.00000000000000,0.00000000000000
+0.96,0.00000000000000,0.00000000000000,0.00000000000000,0.00000000000000,0.00000000000000
+0.97,0.00000000000000,0.00000000000000,0.00000000000000,0.00000000000000,0.00000000000000
+0.98,0.00000000000000,0.00000000000000,0.00000000000000,0.00000000000000,0.00000000000000
+0.99,0.00038569700000,0.00000000000000,0.00000000000000,0.00000000000000,0.00000000000000
+1,0.00000000000000,0.00000000000000,0.00000000000000,0.00000000000000,0.00000000000000
diff --git a/experiments/array_script_CHD.sh b/experiments/array_script_CHD.sh
@@ -0,0 +1,10 @@
+#!/bin/bash
+#SBATCH --time=0:40:00
+#SBATCH --mem=16G
+#SBATCH --cpus-per-task=8
+#SBATCH --output=results/CHD_%a.out
+#SBATCH --error=/dev/null
+#SBATCH --array=0-100
+
+module load julia
+srun julia slurmjob_CHD.jl $SLURM_ARRAY_TASK_ID
diff --git a/experiments/slurmjob_CHD.jl b/experiments/slurmjob_CHD.jl
@@ -0,0 +1,209 @@
+using Pkg
+Pkg.activate((@__DIR__)*"/..")
+using Logging
+using JuMP, Gurobi
+using DecisionProgramming
+using CSV, DataFrames, PrettyTables
+using DelimitedFiles
+
+
+
+# Setting subproblem specific parameters
+idx = parse(Int64, ARGS[1])
+chosen_risk_level = string(idx)*"%" 
+
+
+
+# Reading tests' technical performance data (dummy data in this case)
+data = CSV.read("CHD_preventative_care_data.csv", DataFrame)
+
+
+# Bayes posterior risk probabilities calculation function
+# prior = prior risk level for which the posterior risk distribution is calculated for,
+# t = test done
+# returns a 100x1 vector with the probabilities of getting CHD given the prior risk level and test result
+# for no test done (i.e. t = 3) returns a zero vector
+function update_risk_distribution(prior::Int64, t::Int64)
+    if t == 1 # the test is TRS
+        # P(TRS = result | sick) = P(TRS_if_sick = result) * P(sick) = P(TRS_if_sick = result) * P(prior_risk)
+        numerators = data.TRS_if_sick .* data.risk_levels[prior]
+
+        # P(TRS = result) = P(TRS_if_sick = result) * P(sick) + P(TRS_if_healthy = result) * P(healthy)
+        denominators = data.TRS_if_sick .* data.risk_levels[prior]  + data.TRS_if_healthy .* (1 - data.risk_levels[prior])
+
+        posterior_risks = numerators./denominators
+
+        # if the denominator is zero, post_risk is NaN, changing those to 0
+        for i = 1:101
+            if isnan(posterior_risks[i])
+                posterior_risks[i] = 0
+            end
+        end
+
+        return posterior_risks
+
+
+    elseif t == 2 #the test is GRS
+        numerators = (data.GRS_if_sick .* data.risk_levels[prior])
+        denominators = data.GRS_if_sick .* data.risk_levels[prior]  + data.GRS_if_healthy .* (1 .- data.risk_levels[prior])
+
+        posterior_risks =  numerators./denominators
+
+        # if the denominator is zero, post_risk is NaN, changing those to 0
+        for i = 1:101
+            if isnan(posterior_risks[i])
+                posterior_risks[i] = 0
+            end
+        end
+
+        return posterior_risks
+
+
+    else # no test performed
+        risks_unchanged = zeros(100,1)
+
+
+        return risks_unchanged
+
+    end
+end
+
+# State probabilites calculation function
+# risk_p = the resulting array from update_risk_distribution
+# t = test done
+# h = CHD or no CHD
+# returns the probability distribution in 101x1 vector for the states of the R node given the prior risk level (must be same as to function update_risk_distribution), test t and health h
+function state_probabilities(risk_p::Array{Float64}, t::Int64, h::Int64, prior::Int64)
+
+    #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)
+        state_probabilites[prior] = 1.0
+        return state_probabilites
+    end
+
+    # return vector
+    state_probabilites = zeros(101)
+
+    # copying the probabilities of the scores for ease of readability
+    if h == 1 && t == 1    # CHD and TRS
+        p_scores = data.TRS_if_sick
+    elseif t ==1    # no CHD and TRS
+        p_scores = data.TRS_if_healthy
+    elseif h == 1 && t == 2 # CHD and GRS
+        p_scores = data.GRS_if_sick
+    else # no CHD and GRS
+        p_scores = data.GRS_if_healthy
+    end
+
+    for i = 1:101 #iterating through all risk levels 0%, 1%, ..., 99%, 100% in data.risk_levels
+        for j = 1:101 #iterates through all risk estimates in risk_p
+            #finding all risk estimates risk_p[j] within risk level i
+            # risk_level[i] <= risk_p < risk_level[i]
+            if i < 101 && data.risk_levels[i] <= risk_p[j] && risk_p[j] < data.risk_levels[i+1]
+                state_probabilites[i] += p_scores[j]
+            elseif i == 101 && data.risk_levels[i] <= risk_p[j] #special case: the highest risk level[101] = 100%
+                state_probabilites[i] += p_scores[j]
+            end
+        end
+    end
+
+    return state_probabilites
+end
+
+
+@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 = [string(x) * "%" for x in [0:1:100;]]
+
+
+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))
+
+add_node!(diagram, ValueNode("TC", ["T1", "T2"]))
+add_node!(diagram, ValueNode("HB", ["H", "TD"]))
+
+
+generate_arcs!(diagram)
+
+X_R0 = ProbabilityMatrix(diagram, "R0")
+X_R0[chosen_risk_level] = 1
+add_probabilities!(diagram, "R0", X_R0)
+
+
+X_H = ProbabilityMatrix(diagram, "H")
+X_H[:, "CHD"] = data.risk_levels
+X_H[:, "no CHD"] = 1 .- data.risk_levels
+add_probabilities!(diagram, "H", X_H)
+
+X_R = ProbabilityMatrix(diagram, "R1")
+for s_R0 = 1:101, s_H = 1:2, s_T1 = 1:3
+    X_R[s_R0, s_H,  s_T1, :] =  state_probabilities(update_risk_distribution(s_R0, s_T1), s_T1, s_H, s_R0)
+end
+add_probabilities!(diagram, "R1", X_R)
+add_probabilities!(diagram, "R2", X_R)
+
+
+cost_TRS = -0.0034645
+cost_GRS = -0.004
+forbidden = 0     #the cost of forbidden test combinations is negligible
+Y_TC = UtilityMatrix(diagram, "TC")
+Y_TC["TRS", "TRS"] = forbidden
+Y_TC["TRS", "GRS"] = cost_TRS + cost_GRS
+Y_TC["TRS", "no test"] = cost_TRS
+Y_TC["GRS", "TRS"] = cost_TRS + cost_GRS
+Y_TC["GRS", "GRS"] = forbidden
+Y_TC["GRS", "no test"] = cost_GRS
+Y_TC["no test", "TRS"] = cost_TRS
+Y_TC["no test", "GRS"] = cost_GRS
+Y_TC["no test", "no test"] = 0
+add_utilities!(diagram, "TC", Y_TC)
+
+Y_HB = UtilityMatrix(diagram, "HB")
+Y_HB["CHD", "treatment"] = 6.89713671259061
+Y_HB["CHD", "no treatment"] = 6.65436854256236
+Y_HB["no CHD", "treatment"] = 7.64528451705134
+Y_HB["no CHD", "no treatment"] = 7.70088349200034
+add_utilities!(diagram, "HB", Y_HB)
+
+generate_diagram!(diagram)
+
+
+@info("Creating the decision model.")
+model = Model()
+z = DecisionVariables(model, diagram)
+
+# Defining forbidden paths to include all those where a test is repeated twice
+forbidden_tests = ForbiddenPath(diagram, ["T1","T2"], [("TRS", "TRS"),("GRS", "GRS"),("no test", "TRS"), ("no test", "GRS")])
+fixed_R0 = FixedPath(diagram, Dict("R0" => chosen_risk_level))
+scale_factor = 10000.0
+x_s = PathCompatibilityVariables(model, diagram, z; fixed = fixed_R0, forbidden_paths = [forbidden_tests], probability_cut=false)
+
+EV = expected_value(model, diagram, x_s)
+@objective(model, Max, EV)
+
+@info("Starting the optimization process.")
+optimizer = optimizer_with_attributes(
+    () -> Gurobi.Optimizer(Gurobi.Env()),
+    "TimeLimit" => 1800,
+    "IntFeasTol"=> 1e-9,
+    "MIPGap" => 1e-6
+)
+set_optimizer(model, optimizer)
+
+optimize!(model)
+
+if termination_status(model) == OPTIMAL
+    writedlm((@__DIR__)*"/results/CHD_"*string(idx)*".csv", [objective_value(model), solve_time(model)], ',')
+end
+