updated rng seeding for better consistency. Added a bash script to ru…

…n all experiments

examples/geothermal_example.jl

@@ -11,8 +11,8 @@ default(framestyle = :box,  color_palette=:seaborn_deep6, fontfamily="Computer M
 
 # Define random seeds
 fix_solve_and_eval_seed = true # Whether the seed is set before each policy gen and evaluation. Seed is the eval index + test set. It is threadsafe. 
-pomdp_gen_seed = 0 # seed used to control the generation of the pomdps
-split_by = :geo # Split the test set for ensuring unique :geo, :econ, or :both
+pomdp_gen_seed = 1234 # seed used to control the generation of the pomdps
+split_by = :both # Split the test set for ensuring unique :geo, :econ, or :both
 
 # Can parse 1 command line argument for the split_by 
 if length(ARGS) > 0 
@@ -125,7 +125,7 @@ pomdps, test_sets = create_pomdps(
     econ_params, 
     obs_actions, 
     Nbins, 
-    rng=MersenneTwister(pomdp_gen_seed), # Set the pomdp random seed 
+    rng_seed=pomdp_gen_seed, # Set the pomdp random seed 
     discount=discount_factor,
     split_by=split_by
 )
@@ -224,16 +224,16 @@ elseif split_by == :econ
     geo_fracs = alls
     econ_fracs = fracs
 elseif split_by == :both
-    geo_fracs = [fracs..., alls...]
-    econ_fracs = [alls..., fracs...]
+    geo_fracs = [fracs..., alls[1:end-1]...]
+    econ_fracs = [alls..., fracs[1:end-1]...]
 end
 
 zip_fracs = [(g,e) for (g, e) in zip(geo_fracs, econ_fracs)]
-pomdps_per_geo, test_sets_per_geo = create_pomdps_with_different_training_fractions(zip_fracs, scenario_csvs, geo_params, econ_params, obs_actions, Nbins; rng=MersenneTwister(0), discount=discount_factor, split_by)
+pomdps_per_geo, test_sets_per_geo = create_pomdps_with_different_training_fractions(zip_fracs, scenario_csvs, geo_params, econ_params, obs_actions, Nbins; rng_seed=pomdp_gen_seed, discount=discount_factor, split_by)
 
 # Solve the policies and evaluate the results #<---- Uncomment the below lines to solve and eval the policies
 results = Dict()
-for (policy, pol_name) in zip(policies, policy_names)
+for (policy, pol_name) in zip(policies[end:end], policy_names[end:end])
     results[pol_name] = Dict(:geo_frac => [], :econ_frac => [], :results =>[])
     for (frac, pomdps, test_sets) in zip(zip_fracs, pomdps_per_geo, test_sets_per_geo)
         geo_frac, econ_frac = frac[1], frac[2]
@@ -246,13 +246,17 @@ end
 JLD2.@save joinpath(savedir, "Nsamples_results.jld2") results
 
 # Alternatively, load from file by uncommenting the following lines
-# results = JLD2.load(joinpath(savedir, "Nstates_results.jld2"))["results"] # <---- Uncomment this line to load the results from file
+# results = JLD2.load(joinpath(savedir, "Nsamples_results.jld2"))["results"] # <---- Uncomment this line to load the results from file
 
 # Show how all metrics across all policies vary with number of subsurface realizations
-train_states_comparison_summary(results)
-savefig(joinpath(savedir, "subsurfaces_realizations_comparison.pdf"))
+# train_states_comparison_summary(results)
+# savefig(joinpath(savedir, "subsurfaces_realizations_comparison.pdf"))
 
 # Print the same for just the reward for just one policy
 reward_vs_ngeolgies(results["SARSOP Policy"], "SARSOP Policy")
-savefig(joinpath(savedir, "SARSOP_subsurface_realizations.pdf")) 
-# savefig(joinpath(savedir, "SARSOP_subsurface_realizations.tex")) 
+savefig(joinpath(savedir, "SARSOP_geological_realizations.pdf")) 
+# savefig(joinpath(savedir, "SARSOP_geological_realizations.tex")) 
+
+reward_vs_necon(results["SARSOP Policy"], "SARSOP Policy")
+savefig(joinpath(savedir, "SARSOP_economic_realizations.pdf")) 
+# savefig(joinpath(savedir, "SARSOP_economic_realizations.tex")) 

runner.sh

@@ -0,0 +1,35 @@
+#!/bin/bash
+
+echo "Determining the operating system..."
+os_name=$(uname)
+
+echo "Operating System detected: $os_name"
+
+# Determine the number of CPU threads based on the operating system
+if [ "$os_name" = "Linux" ]; then
+    num_threads=$(nproc)
+elif [ "$os_name" = "Darwin" ]; then
+    # Darwin is the operating system name for macOS
+    num_threads=$(sysctl -n hw.ncpu)
+else
+    echo "Unsupported Operating System. Defaulting to 1 thread."
+    num_threads=1
+fi
+
+echo "Number of CPU threads available: $num_threads"
+
+export JULIA_NUM_THREADS=$num_threads
+echo "Setting Julia to use $JULIA_NUM_THREADS threads"
+
+echo "Running Julia script with different arguments..."
+
+# Run with first argument "geo"
+julia1.9 examples/geothermal_example.jl both
+
+# Run with second argument "econ"
+julia1.9 examples/geothermal_example.jl geo
+
+# Run with third argument "both"
+julia1.9 examples/geothermal_example.jl econ
+
+echo "Script execution completed."

src/InfoGatheringPOMDPs.jl

@@ -37,7 +37,7 @@ include("metrics.jl")
 
 export policy_results_summary, policy_comparison_summary, train_states_comparison_summary, 
        policy_sankey_diagram, trajectory_regret_metrics, generate_action_table, scenario_returns, 
-       policy_comparison_table, pes_comparison, expected_loss_comparison, reward_vs_ngeolgies
+       policy_comparison_table, pes_comparison, expected_loss_comparison, reward_vs_ngeolgies, reward_vs_necon
 include("plotting.jl")
 
 end

src/plotting.jl

@@ -256,30 +256,43 @@ function policy_comparison_table(policy_results, policy_names)
     end
 end
 
-function train_states_comparison_summary(results)
-    p_reward = plot()
-    p_obs_cost = plot()
-    p_num_obs = plot()
-    p_correct_scenario = plot()
-    p_correct_gonogo = plot()
-    p_legend = plot(legend=false, grid=false, axis=false, ticks=nothing, border=:none, size=(800,400))
-
-    xlab = "No. Subsurface Realizations"
-
-    for (policy_name, pol_results) in results
-        plot!(p_reward, pol_results[:geo_frac]*250, [mean(r[:reward]) for r in pol_results[:results]], xlabel=xlab, ylabel="Mean Discounted Reward", legend = false)
-        plot!(p_obs_cost, pol_results[:geo_frac]*250, [mean(r[:obs_cost]) for r in pol_results[:results]], xlabel=xlab, ylabel="Mean Data Acquisition Cost", legend = false)
-        plot!(p_num_obs, pol_results[:geo_frac]*250, [mean(r[:num_obs]) for r in pol_results[:results]], xlabel=xlab, ylabel="Mean Number of Data Acquisition Actions", legend = false)
-        plot!(p_correct_scenario, pol_results[:geo_frac]*250, [mean(r[:correct_scenario]) for r in pol_results[:results]], xlabel=xlab, ylabel="Mean Correct Scenario", legend = false)
-        plot!(p_correct_gonogo, pol_results[:geo_frac]*250, [mean(r[:correct_gonogo]) for r in pol_results[:results]], xlabel=xlab, ylabel="Mean Correct Go/NoGo", legend = false)
-        plot!(p_legend, [], [], label=policy_name, legend=:topleft)
-    end
-    plot(p_reward, p_obs_cost, p_num_obs, p_correct_scenario, p_correct_gonogo, p_legend, layout=(2,3),size=(1400,800), margin=5mm)
-end
+# function train_states_comparison_summary(results)
+#     p_reward = plot()
+#     p_obs_cost = plot()
+#     p_num_obs = plot()
+#     p_correct_scenario = plot()
+#     p_correct_gonogo = plot()
+#     p_legend = plot(legend=false, grid=false, axis=false, ticks=nothing, border=:none, size=(800,400))
+
+#     xlab = "No. Subsurface Realizations"
+
+#     for (policy_name, pol_results) in results
+#         plot!(p_reward, pol_results[:geo_frac]*250, [mean(r[:reward]) for r in pol_results[:results]], xlabel=xlab, ylabel="Mean Discounted Reward", legend = false)
+#         plot!(p_obs_cost, pol_results[:geo_frac]*250, [mean(r[:obs_cost]) for r in pol_results[:results]], xlabel=xlab, ylabel="Mean Data Acquisition Cost", legend = false)
+#         plot!(p_num_obs, pol_results[:geo_frac]*250, [mean(r[:num_obs]) for r in pol_results[:results]], xlabel=xlab, ylabel="Mean Number of Data Acquisition Actions", legend = false)
+#         plot!(p_correct_scenario, pol_results[:geo_frac]*250, [mean(r[:correct_scenario]) for r in pol_results[:results]], xlabel=xlab, ylabel="Mean Correct Scenario", legend = false)
+#         plot!(p_correct_gonogo, pol_results[:geo_frac]*250, [mean(r[:correct_gonogo]) for r in pol_results[:results]], xlabel=xlab, ylabel="Mean Correct Go/NoGo", legend = false)
+#         plot!(p_legend, [], [], label=policy_name, legend=:topleft)
+#     end
+#     plot(p_reward, p_obs_cost, p_num_obs, p_correct_scenario, p_correct_gonogo, p_legend, layout=(2,3),size=(1400,800), margin=5mm)
+# end
 
 function reward_vs_ngeolgies(pol_results, policy_name; p=plot())
-    xlab = "No. Subsurface Realizations"
-    plot!(p, pol_results[:geo_frac]*250, [mean(r[:reward]) for r in pol_results[:results]], xlabel=xlab, ylabel="Mean Discounted Reward", legend = false, title=policy_name)
+    xlab = "No. Geological Realizations"
+    max_econ = maximum(pol_results[:econ_frac])
+    xs = 250*pol_results[:geo_frac][pol_results[:econ_frac] .== max_econ]
+    perm = sortperm(xs)
+    ys = [mean(r[:reward]) for (r, e) in zip(pol_results[:results], pol_results[:econ_frac]) if e == max_econ]
+    plot!(p, xs[perm], ys[perm], xlabel=xlab, ylabel="Mean Discounted Reward", legend = false, title=policy_name)
+end
+
+function reward_vs_necon(pol_results, policy_name; p=plot())
+    xlab = "No. Economic Realizations"
+    max_geo = maximum(pol_results[:geo_frac])
+    xs = 50*pol_results[:econ_frac][pol_results[:geo_frac] .== max_geo]
+    perm = sortperm(xs)
+    ys = [mean(r[:reward]) for (r, g) in zip(pol_results[:results], pol_results[:geo_frac]) if g == max_geo]
+    plot!(p, xs[perm], ys[perm], xlabel=xlab, ylabel="Mean Discounted Reward", legend = false, title=policy_name)
 end
 
 function generate_action_table(pomdp, var_desc)

src/pomdp.jl

@@ -115,12 +115,12 @@ function parseCSV(
 end
 
 # Generate k folds of the data
-function kfolds(states, nfolds; split_by=:geo, geo_train_frac = 1.0-1.0/nfolds, econ_train_frac = 1.0-1.0/nfolds, rng=Random.GLOBAL_RNG)
+function kfolds(states, nfolds; split_by=:geo, geo_train_frac = 1.0-1.0/nfolds, econ_train_frac = 1.0-1.0/nfolds)
     geo_indices = unique([s[:GeoID] for s in states])
     econ_indices = unique([s[:EconID] for s in states])
 
-    shuffled_geo_indices = shuffle(rng, geo_indices)
-    shuffled_econ_indices = shuffle(rng, econ_indices)
+    shuffled_geo_indices = shuffle(geo_indices)
+    shuffled_econ_indices = shuffle(econ_indices)
 
     Ngeo = length(shuffled_geo_indices)
     Necon = length(shuffled_econ_indices)
@@ -286,11 +286,12 @@ function POMDPs.gen(m::InfoGatheringPOMDP, s, a, rng)
     return (;sp, o, r)
 end
 
-function create_pomdps(scenario_csvs, geo_params, econ_params, obs_actions, Nbins; split_by=:geo, Nsamples_per_bin=10, nfolds=5, geo_train_frac = 1.0 - 1.0/nfolds, econ_train_frac = 1.0 - 1.0/nfolds, discount=0.9, rng=Random.GLOBAL_RNG)
+function create_pomdps(scenario_csvs, geo_params, econ_params, obs_actions, Nbins; split_by=:geo, Nsamples_per_bin=10, nfolds=5, geo_train_frac = 1.0 - 1.0/nfolds, econ_train_frac = 1.0 - 1.0/nfolds, discount=0.9, rng_seed=0)
     # Parse all of the states from the csv files
     statevec = parseCSV(scenario_csvs, geo_params, econ_params)
 
-    train_sets, test_sets = kfolds(statevec, nfolds; geo_train_frac, econ_train_frac, rng, split_by)
+    Random.seed!(rng_seed)
+    train_sets, test_sets = kfolds(statevec, nfolds; geo_train_frac, econ_train_frac, split_by)
     pomdps = Array{Any}(undef, nfolds)
     p = Progress(nfolds, 1, "Creating POMDPs...")
     Threads.@threads for i in 1:nfolds
@@ -316,11 +317,11 @@ function create_pomdps(scenario_csvs, geo_params, econ_params, obs_actions, Nbin
     return pomdps, test_sets
 end
 
-function create_pomdps_with_different_training_fractions(train_fracs, scenario_csvs, geo_params, econ_params, obs_actions, Nbins; split_by=:geo, Nsamples_per_bin=10, nfolds=5, discount=0.9, rng=Random.GLOBAL_RNG)
+function create_pomdps_with_different_training_fractions(train_fracs, scenario_csvs, geo_params, econ_params, obs_actions, Nbins; split_by=:geo, Nsamples_per_bin=10, nfolds=5, discount=0.9, rng_seed=0)
     pomdps_per_ngeo = []
     test_sets_per_frac = []
     for (geo_train_frac, econ_train_frac) in train_fracs
-        pomdps, test_sets = create_pomdps(scenario_csvs, geo_params, econ_params, obs_actions, Nbins; split_by, Nsamples_per_bin, nfolds, geo_train_frac, econ_train_frac, discount, rng=deepcopy(rng))
+        pomdps, test_sets = create_pomdps(scenario_csvs, geo_params, econ_params, obs_actions, Nbins; split_by, Nsamples_per_bin, nfolds, geo_train_frac, econ_train_frac, discount, rng_seed)
         push!(pomdps_per_ngeo, pomdps)
         push!(test_sets_per_frac, test_sets)
     end