This notebook is part of the supplementary material for:
Genewein T., Leibfried F., Grau-Moya J., Braun D.A. (2015) Bounded rationality, abstraction and hierarchical decision-making: an information-theoretic optimality principle, Frontiers in Robotics and AI.
More information on how to run the notebook on the accompanying github repsitory where you can also find updated versions of the code and notebooks.
This notebook corresponds to Sections 4 and 5 and reproduces Figure 11 of the paper.
This notebook allows to compare the solutions of any three-variable case (general, parallel, serial) against the solution of the one-step rate-distortion case (Section 2 of the paper). By default, it illustrates how information processing in the parallel hierarchical case can be split up on different information processing pathways (see Section 4.3 of the paper).
Below we show (in the following order)
#only run this once
include("RateDistortionDecisionMaking.jl")
WARNING: replacing module RateDistortionDecisionMaking
RateDistortionDecisionMaking
using RateDistortionDecisionMaking, Distances, DataFrames, Colors, Gadfly, Distributions, Interact, Reactive
#make the default plot size a bit larger
set_default_plot_size(15cm, 12cm)
#set up taxonomy example
#include("TaxonomyExample.jl")
#w_vec, w_strings, a_vec, a_strings, p_w, U = setuptaxonomy()
#a_values = a_vec
#w_values = w_vec
#set up medical system example
include("MedicalSystemExample.jl")
w_values, w_strings, a_values, a_strings, p_w, U = setup_medical_example(uniform_w=true)
num_acts = length(a_values)
a_vec = collect(1:num_acts)
num_worldstates = length(w_values)
w_vec = collect(1:num_worldstates)
#pre-compute utilities, find maxima
U_pre, Umax = setuputilityarrays(a_values,w_values,U)
#visualize utility
plt_utility = visualizeMatrix(U_pre, w_values, a_values, w_strings, a_strings, xlabel="Disease type w",
ylabel="Treatment a", legendlabel="U(a,w)")
Make sure to re-evaluate all the subsequent cells after changing parameters here. The one-step case corresponds to Section 2 of the paper. The general case is found in Section 5 and by choosing certain settings of the inverse temperatures (see Table 2 in the paper), the parallel, serial or even the one-step case can be recovered.
#one-step Blahut-Arimoto
β = 10 #inverse temperature
#general case
#β1=β3 and β2=∞ is identical to the one-variable case
#β1 = β #inverse temperature for p(x)→p(x|w)
#β2 = Inf #inverse temperature for p(a)→p(a|x)
#β3 = β #inverse temperature for p(a|x)→p(a|x,w)
#parallel case: β2=∞
β1 = 15 #inverse temperature for p(x)→p(x|w)
β2 = Inf #17.5 #inverse temperature for p(a)→p(a|x)
β3 = 14.99 #3.999 #inverse temperature for p(a|x)→p(a|x,w)
#sequential case: β3=0
#β1 = 1 #inverse temperature for p(o)→p(x|w)
#β2 = 8 #inverse temperature for p(a)→p(a|x)
#β3 = 0 #inverse temperature for p(a|x)→p(a|x,w)
#cardinality of percept
num_x = 3#length(a_vec) #cardinality of X is equal to action-space
ε = 0.0001 #convergence critetion for BAiterations
maxiter = 2000 #maximum number of BA iterations
2000
#solve the example with one-step Blahut-Arimoto
#initialize p(a) uniformly
num_acts = length(a_vec)
pa_init = ones(num_acts)/num_acts
#Blahut-Arimotot iterations
p_agw_one, p_a_one, perf_df_one = BAiterations(pa_init, β, U_pre, p_w, ε, maxiter,
compute_performance=true, performance_as_dataframe=true)
#visualize solution
plt_marg, plt_cond = visualizeBAsolution(p_a_one, p_agw_one, a_vec, w_vec, a_strings, w_strings,
wlabel="Disease type w", alabel="Treatment a",
legendlabel_marginal="p(a)", legendlabel_conditional="p*(a|w)",
suppress_vis=true);
display(plt_utility)
display(plt_cond)
x_vec = collect(1:num_x)
x_strings = map((x)->"x"*string(x), x_vec)
num_worldstates = length(w_vec)
#This function performs Blahut-Arimoto iterations for the three-variable general case
px, pa, pxgw, pagx, pagxw, pagw, performance_df = threevarBAiterations(num_x, β1, β2, β3, U_pre, p_w, ε, maxiter,
compute_performance=true, performance_per_iteration=true,
performance_as_dataframe=true,
init_pogw_sparse = true, init_pogw_uniformly = false,
init_pagow_uniformly = true);
! Note that in the legends below $O$ is used instead of $X$
#plot convergence (i.e. evolution across BA iterations)
plot_convergence_MI, plot_convergence_Util = plot_three_var_BA_convergence(performance_df);
Note that $X$ is labeled here with the general term "Observation" but in case of the medical example it indicates the first diagnosis (or model)
#call the routine that creates plots for all the probability distributions
plt_px, plt_pa, plt_pxgw, plt_pagx, plt_pagw, dpdown, plt_pagxw_vis = visualize_three_var_BAsolution(px, pa,
pxgw, pagx, pagxw, pagw,
x_vec, a_vec, w_vec,
x_strings, a_strings, w_strings,
olabel_string="x", alabel_string="a", wlabel_string="w")
#visualize p(x) and p(a)
display(hstack(plt_px,plt_pa))
#visualize p(x|w)
display(plt_pxgw)
#visualize p(a|x)
display(plt_pagx)
#visualize p(a|w)
display(plt_pagw)
#visualize p(a|x,w)
display(dpdown) #show the dropdown box to interactively select wk in p(a|x,w=wk)
display(plt_pagxw_vis)
If you see the drowpdown-box but not the plot and "Javascript error adding output!"
, try selecting another $w$ with the dropdwon once, then the plot should appear.
If this does not work (i.e. kernel is busy with this cell but nothing ever happens), try ckecking out the most current version of Interact with: Pkg.checkout("Interact")
in a Julia console. You can undo this and go back to the latest released version with Pkg.free("Interact")
.
! Note that in the legends below $O$ is used instead of $X$
p_MI, p_composed, p_perf=plot_three_var_performancemeasures(performance_df, maximum(U_pre), β1, β2, β3)
display(hstack(p_MI, p_composed))
display(p_perf)
Whenever the same quantity is shown twice in any of the plots, the left quantity corresponds to the one-step case (denoted by subindex _one
and the right quantity corresponds to the general case.
#compare solution against one-step case
#compare performance measures, EU, J, I(A;W)
util_all = [perf_df_one[end,:E_U], performance_df[end,:E_U],
perf_df_one[end,:RD_obj], performance_df[end,:Objective_value]]
color_label = ["E[U]_one", "E[U]", "J_one", "J"]
colors = BAdiscretecolorscale(4).f(4)
colscale = Scale.color_discrete_manual(colors[1],colors[1],colors[end],colors[end])
p_perf_comp = plot(x=color_label, y=util_all, color=color_label, Geom.bar(),
Guide.ylabel("[utils]"), Guide.xlabel(""), Guide.colorkey(""),
Guide.title("One-step solution β=$β vs. \n multi-step β1=$β1, β2=$β2, β3 = $β3"),
colscale, Scale.y_continuous(minvalue=0, maxvalue=maximum(U_pre)),
BAtheme(key_position = :none))
#compare mutual informations I(A;W)
MI_comp = [perf_df_one[end,:I_aw], performance_df[end,:I_aw]]
color_label = ["I(A;W)_one", "I(A;W)"]
colors = Scale.color_discrete_hue().f(4)
colscale = Scale.color_discrete_manual(colors[end],colors[end])
p_perf_comp_MI = plot(x=color_label, y=MI_comp, color=color_label, Geom.bar(),
Guide.ylabel("[bits]"), Guide.xlabel(""), Guide.colorkey(""),
Guide.title("I(A;W) for β=$β vs. \n multi-step β1=$β1, β2=$β2, β3 = $β3"),
colscale, Scale.y_continuous(minvalue=0), BAtheme(key_position = :none))
display(hstack(p_perf_comp,p_perf_comp_MI))
Top: one-step case Bottom: general case
#compare solutions p(a|w)
plt_cond_one = visualizeBAconditional(p_agw_one, a_vec, w_vec, a_strings, w_strings,
wlabel="w", alabel="a",legendlabel="p*(a|w)_one")
plt_cond_multi = visualizeBAconditional(pagw, a_vec, w_vec, a_strings, w_strings,
wlabel="w", alabel="a",legendlabel="p*(a|w)")
display(plt_cond_one)
display(plt_cond_multi)
display("Difference: norm(p*(a|w)_one - p*(a|w))=$(norm(p_agw_one-pagw))")
plt_cond_diff = visualizeMatrix(p_agw_one - pagw, w_vec, a_vec, w_strings, a_strings,
xlabel="Disease type w", ylabel="Treatment a",legendlabel="p*(a|w)_one - p*(a|w)")
"Difference: norm(p*(a|w)_one - p*(a|w))=6.118579364182079e-5"
! Note that in the legends below $O$ is used instead of $X$
#for the paper
plt_final = vstack((hstack(plt_cond_multi, p_perf_comp_MI )), hstack(p_MI, p_perf_comp))
display(plt_final)
w = 18cm
h = 15cm
#draw(SVG("Figures/Hierarchical_OneStep_Comparison.svg", w, h), plt_final) #uncomment to store figure#
Compose.Measure{Compose.MeasureNil,Compose.MeasureNil}(150.0,Compose.MeasureNil(),Compose.MeasureNil(),0.0,0.0)