#only run this once include("RateDistortionDecisionMaking.jl") 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 predator-prey example include("PredatorPreyExample.jl") ###### Predator-prey scenario ############# #same utility as used in the main paper w_values, w_strings, a_values, a_strings, p_w, U = setup_predator_prey_example() ###### Mating scenario #################### #alternative utility - mating scenario #w_values, w_strings, a_values, a_strings, p_w, U = setup_predator_prey_example(mating_utility=true) numa = length(a_strings) a_vec = collect(1:numa) numw = length(w_strings) w_vec = collect(1:numw) #pre-compute utility 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="animal size w", ylabel="Action a", legendlabel="U(a,w)") #set up hand-crafted likelihood model and p(x|w) x_values = collect(1:13) x_vec = collect(1:length(x_values)) x_strings = map((x)->string(x), x_values) #slider for selecting λ λ_vals = 0.1:0.1:5 slider_l = slider(λ_vals,label="Perceptual precision λ") #use lift to connect the actual plotting-code to the slider plt_pxgw_sl = lift(λ_vis->begin p_xgw_hc = pogw_handcrafted(x_values, w_values, λ_vis) visualizeBAconditional(p_xgw_hc, x_vec, w_vec, x_strings, w_strings, wlabel="animal size w", alabel="observed size x", legendlabel="p(x|w,λ)") end, slider_l) display(slider_l) display(plt_pxgw_sl) #precision of the hand-crafted perceptual model λ = 1.65 #0.4 #1.65 #the "price" for I(X;W) in the handcrafted case will be β1 (as given below) #to make the two cases comparable. If you want I(X;W) to be the same in both #cases, set β1 and run the sequential case and then tune λ until the mutual #information terms are equal. # #Note that β2 will also be used for the bounded rational decision-maker that #uses the hand-crafted perception. This should allow for easy comparison of #the action-channels in both cases (hand-crafted vs. sequential case) #inverse temperatures for sequential case #β1: perceptual channel -> price for I(X;W) #β2: action channel -> price for I(A;X) #β3=0 sequential case (otherwise the general three-variable case is specified) β1 = 8 #1 β2 = 10 β3 = 0; #compute hand-crafted perception, then feed it into bounded rational decision-maker #with inv. temperature β2 and compute optimal p(a|x) using hand-crafted p(x|w) p_xgw_hc = pogw_handcrafted(x_values, w_values, λ) numx = length(x_vec) numa = length(a_vec) #comptue expected utility under posterior over w given x # #posterior p(w|x) p_wgx = zeros(numw, numx) for k in 1:numx p_wgx[:,k] = p_xgw_hc[k,:]' .* p_w p_wgx[:,k] += eps() #add small epsilon to prevent numerical problems for λ>> p_wgx[:,k] /= sum(p_wgx[:,k]) end #U(a,x)=∑_w p(w|x)U(a,w) U_ax = zeros(numa, numx) for k in 1:numx U_ax[:,k] = U_pre * p_wgx[:,k] end #p(x)=∑_w p(w)p(x|w) p_x_hc = p_xgw_hc * p_w #make sure that p(x) has non-zero entries (otherwise numerical problems arise in I(X;W)) p_x_hc += eps() p_x_hc /= sum(p_x_hc) #solve the decision-making part with one-step Blahut-Arimoto ε = 0.0001 #convergence critetion for BAiterations maxiter = 10000 #maximum number of BA iterations #initialize p(a) uniformly num_acts = length(a_vec) pa_init = rand(num_acts) pa_init /= sum(pa_init) #BA iterations p_agx_hc, p_a_hc, perf_df_hc = BAiterations(pa_init, β2, U_ax, p_x_hc, ε, maxiter, compute_performance=true, performance_as_dataframe=true) ; #suppress output #visualize solution p(x|w), p(a|x), p(a|w) #p(x|w) plt_pxgw_hc = visualizeBAconditional(p_xgw_hc, x_vec, w_vec, x_strings, w_strings, wlabel="animal size w", alabel="observed size x", legendlabel="p(x|w,λ)") #p(a|x) plt_pa_hc, plt_pagx_hc = visualizeBAsolution(p_a_hc, p_agx_hc, a_vec, x_vec, a_strings, x_strings, wlabel="observed size x", alabel="action a", legendlabel_marginal="p(a)", legendlabel_conditional="p(a|x,λ)", suppress_vis=true) #compute p(a|w) #p(a|w) = ∑_x p(x|w)p(a|x) p_agw_hc = p_agx_hc * p_xgw_hc plt_pagw_hc = visualizeBAconditional(p_agw_hc, a_vec, w_vec, a_strings, w_strings, wlabel="animal size w", alabel="action a", legendlabel="p(a|w,λ)") #visualize p(x|w) display(plt_pxgw_hc) #visualize p(a|x) display(plt_pagx_hc) #visualize p(a|w) display(plt_pagw_hc) #compute p(a|x,w) = p(a|x) ∀w p_agxw_hc = zeros(numa, numx, numw) for j in 1:numw p_agxw_hc[:,:,j] = p_agx_hc end #compute performance measures of solution res_hc = analyze_three_var_BAsolution(p_w, p_x_hc, p_a_hc, p_xgw_hc, p_agx_hc, p_agxw_hc, p_agw_hc, U_pre, β1, β2, β3) perf_df_hc = performancemeasures2DataFrame(res_hc...) #visualize p_MI, p_composed, p_perf = plot_three_var_performancemeasures(perf_df_hc, maximum(U_pre), β1, β2, β3) display(hstack(p_MI, p_composed)) display(p_perf) ε = 0.0001 #convergence critetion for BAiterations maxiter = 10000 #maximum number of BA iterations #This function performs Blahut-Arimoto iterations for the three-variable general case px, pa, pxgw, pagx, pagxw, pagw, performance_df = threevarBAiterations(numx, β1, β2, β3, U_pre, p_w, ε, maxiter, compute_performance=true, performance_as_dataframe=true, init_pogw_sparse = true, init_pogw_uniformly = false, init_pagow_uniformly = true); #call the routine that creates plots for all the probability distributions plt_px, plt_pa, plt_pxgw, plt_pagx, plt_pagw, dpdown, plt_pagow_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|w) display(plt_pxgw) #visualize p(a|x) display(plt_pagx) #visualize p(a|w) display(plt_pagw) 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) #compare both solutions #compare performance measures, EU, J, I(A;W) util_all = [perf_df_hc[end,:E_U], performance_df[end,:E_U], perf_df_hc[end,:Objective_value], performance_df[end,:Objective_value]] color_label = ["E[U]_hc", "E[U]", "J_hc", "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("λ=$λ vs. β1=$β1, β2=$β2, β3=$β3"), colscale, Scale.y_continuous(minvalue=0, maxvalue=maximum(U_pre)), BAtheme(key_position = :none)) MI_comp = [perf_df_hc[end,:I_ow], performance_df[end,:I_ow], perf_df_hc[end,:I_ao], performance_df[end,:I_ao]] color_label = ["I(X;W)_hc", "I(X;W)", "I(A;X)_hc", "I(A;X)"] colors = Scale.color_discrete_hue().f(4) colscale = Scale.color_discrete_manual(colors[1],colors[1],colors[2],colors[2]) 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("λ=$λ vs. β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)) #compare final behavior p(a|w) of both cases display(plt_pagw_hc) display(plt_pagw) #compare perceptual model p(x|w) of both cases display(plt_pxgw_hc) display(plt_pxgw) #generate plots for comparison #display(plt_utility) #draw(SVG("Figures/SequentialComparison_Utility.svg", 10cm, 8cm), plt_utility) #uncomment to store figure #plot_final = vstack(hstack(plt_pogw_hc, plt_pagw_hc),hstack(plt_pogw, plt_pagw)) #plot_bars = (vstack(p_perf_comp,p_perf_comp_MI)) #display(plot_final) #display(plot_bars) #draw(SVG("Figures/SequentialComparison_3.svg", 18cm, 12cm), plot_final) #uncomment to store figure# #draw(SVG("Figures/SequentialComparisonBars_3.svg", 8cm, 12cm), plot_bars) #uncomment to store figure