Add counterfactual minus guided/unguided diff at location level

ext/AvizExt/viz/spatial.jl

@@ -48,6 +48,7 @@ end
         highlight::Union{Vector,Tuple,Nothing},
         show_colorbar::Bool=true,
         colorbar_label::String="",
+        color_map::$COLORMAP_TYPE_DOCSTRING,
         legend_params::Union{Tuple,Nothing}=nothing,
         axis_opts::Dict{Symbol, <:Any}=set_axis_defaults(Dict{Symbol,Any}())
     )
@@ -71,13 +72,13 @@ function create_map!(
     f::Union{GridLayout,GridPosition},
     geodata::Vector{<:GeoMakie.MultiPolygon},
     data::Observable,
-    highlight::Union{Vector,Tuple,Nothing},
+    highlight::Union{Vector,Tuple,Nothing};
     show_colorbar::Bool=true,
     colorbar_label::String="",
-    color_map::Union{Symbol,Vector{Symbol},RGBA{Float32},Vector{RGBA{Float32}}}=:grayC,
+    color_map::COLORMAP_TYPE(T)=:grayC,
     legend_params::Union{Tuple,Nothing}=nothing,
     axis_opts::OPT_TYPE=set_axis_defaults(DEFAULT_OPT_TYPE())
-)
+)::Union{GridLayout,GridPosition} where {T<:Real}
     spatial = GeoAxis(
         f[1, 1];
         axis_opts...
@@ -86,10 +87,11 @@ function create_map!(
     # spatial.xticklabelsize = 14
     # spatial.yticklabelsize = 14
 
+    min_val = @lift(minimum($data))
     max_val = @lift(maximum($data))
 
     # Plot geodata polygons using data as internal color
-    color_range = (0.0, max_val[])
+    color_range = min_val[] < 0 ? (min_val[], max_val[]) : (0, max_val[])
 
     poly!(
         spatial,
@@ -242,12 +244,30 @@ function ADRIA.viz.map!(
         g,
         geodata,
         data,
-        highlight,
-        show_colorbar,
-        c_label,
-        color_map,
-        legend_params,
-        axis_opts
+        highlight;
+        show_colorbar=show_colorbar,
+        colorbar_label=c_label,
+        color_map=color_map,
+        legend_params=legend_params,
+        axis_opts=axis_opts
+    )
+end
+
+function ADRIA.viz.diff_map(
+    rs::ResultSet,
+    diff_outcome::YAXArray{T};
+    opts::Dict{Symbol,<:Any}=Dict{Symbol,Any}(),
+    fig_opts::Dict{Symbol,<:Any}=Dict{Symbol,Any}(),
+    axis_opts::Dict{Symbol,<:Any}=Dict{Symbol,Any}()
+) where {T}
+    # TODO hande the cases where thes only positive or only negative values
+    min_val, max_res = extrema(diff_outcome)
+    mid_val = -min_val / (max_res - min_val)
+
+    div_cmap::Vector{RGB{Float64}} = diverging_palette(10, 200; mid=mid_val)
+    opts[:color_map] = div_cmap
+    return ADRIA.viz.map(
+        rs, diff_outcome; axis_opts=axis_opts, opts=opts, fig_opts=fig_opts
     )
 end
 

ext/AvizExt/viz/viz.jl

@@ -13,6 +13,13 @@ using .AvizExt
 const OPT_TYPE = Dict{Symbol,<:Any}
 const DEFAULT_OPT_TYPE = Dict{Symbol,Any}
 
+const COLORMAP_TYPE_DOCSTRING = replace("""
+Union{
+Symbol,RGB{T},RGBA{T},Vector{Symbol},Vector{RGBA{T}},Vector{RGB{T}}
+}
+""", "\n" => "")
+const COLORMAP_TYPE = eval(Meta.parse("f(T)=$COLORMAP_TYPE_DOCSTRING"))
+
 """
     _time_labels(labels)
 

src/metrics/spatial.jl

@@ -1,5 +1,8 @@
 """Functions and methods to produce location-level summaries."""
 
+using Bootstrap
+using Random
+
 """
     per_loc(metric, data::YAXArray{D,T,N,A})::YAXArray where {D,T,N,A}
 
@@ -122,3 +125,105 @@ function summarize(
 )::YAXArray where {D,T,N,A}
     return summarize(data[timesteps=timesteps], alongs_axis, metric)
 end
+
+"""
+    ensemble_loc_difference(outcome::YAXArray{T,3}, scens::DataFrame; agg_metric::Union{Function,AbstractFloat}=median, diff_target=:guided, conf::Float64=0.95, rand_seed=1234)::YAXArray where {T}
+
+Mean bootstrapped difference (counterfactual - target) between some outcome aggregated for
+each location.
+
+# Arguments
+- `outcome` : Metric outcome with dimensions (:timesteps, :locations, :scenarios).
+- `scens` : Scenarios DataFrame.
+- `agg_metric` : Metric used to aggregate scenarios when comparing between counterfactual and
+target. If it is an `AbstractFloat` between 0 and 1, it uses the `bs_metric`-th quantile.
+Defaults to `median`.
+- `diff_target` : Target group of scenarios to compare with. Valid options are `:guided` and
+`:unguided`. Defaults to `:guided`
+- `conf` : Percentile used for the confidence interval. Defaults to 0.95.
+- `rand_seed` : Seed used for random number generator. Defaults to 1234.
+
+# Example
+```
+# Load domain
+dom = ADRIA.load_domain(path_to_domain)
+
+# Create scenarios
+num_scens = 2^6
+scens = ADRIA.sample(dom, num_scens)
+
+# Run model
+rs = ADRIA.run_scenarios(dom, scens, "45")
+
+# Calculate difference to the counterfactual for given metric
+_relative_cover = metrics.relative_cover(rs)
+
+# Compute difference between guided and counterfactual using the 0.6-th quantile
+gd_res = metrics.ensemble_loc_difference(r_cover, scens; agg_metric=0.6)
+
+# Compute difference between unguided and counterfactual using the median
+ug_res = metrics.cf_difference_loc(r_cover, scens; diff_target=:unguided)
+
+# Plot maps of difference to the counterfactual
+ADRIA.viz.diff_map(rs, gd_res[2, :])
+ADRIA.viz.diff_map(rs, ug_res[2, :])
+```
+
+# Returns
+Vector with bootstrapped difference (counterfactual - guided) for each location.
+"""
+function ensemble_loc_difference(
+    outcome::YAXArray{T,3},
+    scens::DataFrame;
+    agg_metric::Union{Function,AbstractFloat}=median,
+    diff_target=:guided,
+    conf::Float64=0.95,
+    rand_seed=1234
+)::YAXArray where {T}
+    is_quantile_metric = isa(agg_metric, AbstractFloat)
+    if is_quantile_metric && !(0 <= agg_metric <= 1)
+        error("When metric is a number, it must be within the interval [0,1]")
+    end
+
+    Random.seed!(rand_seed)
+
+    # Mean over all timesteps
+    outcomes_agg = dropdims(mean(outcome; dims=:timesteps); dims=:timesteps)
+
+    # Counterfactual, target outcomes
+    cf_outcomes = outcomes_agg[scenarios=scens.guided .== -1]
+    target_outcomes = if diff_target == :guided
+        outcomes_agg[scenarios=scens.guided .> 0]
+    elseif diff_target == :unguided
+        outcomes_agg[scenarios=scens.guided .== 0]
+    else
+        error("Invalid diff_target value. Valid values are :guided and :unguided.")
+    end
+
+    # Find the smallest set of outcomes because they might not have the same size
+    n_cf_outcomes = size(cf_outcomes, :scenarios)
+    n_target_outcomes = size(target_outcomes, :scenarios)
+    min_n_outcomes = min(n_cf_outcomes, n_target_outcomes)
+
+    # Build (counterfactual-guided) and (counterfactual-unguided) result DataCubes
+    _locations = axis_labels(outcome, :locations)
+    results = ZeroDataCube(;
+        T=T, value=[:lower_bound, :value, :upper_bound], locations=_locations
+    )
+
+    n_locs = length(_locations)
+    for loc in 1:n_locs
+        cf_shuf_set::Vector{Int64} = shuffle(1:n_cf_outcomes)[1:min_n_outcomes]
+        target_shuf_set::Vector{Int64} = shuffle(1:n_target_outcomes)[1:min_n_outcomes]
+
+        @views target_diff = collect(
+            target_outcomes[loc, target_shuf_set] .- cf_outcomes[loc, cf_shuf_set]
+        )
+
+        bootstrap_func(x) = is_quantile_metric ? quantile(x, agg_metric) : agg_metric(x)
+        cf_target_bootstrap = bootstrap(bootstrap_func, target_diff, BalancedSampling(100))
+        results[[2, 1, 3], loc] .= confint(cf_target_bootstrap, PercentileConfInt(conf))[1]
+    end
+
+    return results
+end

src/viz/viz.jl

@@ -42,6 +42,7 @@ function ranks_to_frequencies!() end
 # Spatial
 function map() end
 function map!() end
+function diff_map() end
 function connectivity() end
 function connectivity!() end
 

test/analysis.jl

@@ -7,9 +7,10 @@ if !@isdefined(TEST_RS)
     const TEST_DOM, TEST_N_SAMPLES, TEST_SCENS, TEST_RS = test_rs()
 end
 
+Makie.inline!(false)
+
 """Test larger scenario run with figure creation"""
 function test_rs_w_fig(rs::ADRIA.ResultSet, scens::ADRIA.DataFrame)
-    Makie.inline!(false)
 
     # Visualize results (in terms of absolute coral cover)
     s_tac = ADRIA.metrics.scenario_total_cover(rs)

test/metrics/spatial.jl

@@ -0,0 +1,27 @@
+using ADRIA: metrics, metrics.cf_difference_loc
+using ADRIA: YAXArray
+
+if !@isdefined(TEST_RS)
+    const TEST_DOM, TEST_N_SAMPLES, TEST_SCENS, TEST_RS = test_rs()
+end
+
+@testset "cf_difference_loc" begin
+    test_metrics = [
+        metrics.relative_cover,
+        metrics.juvenile_indicator,
+        metrics.absolute_shelter_volume,
+        metrics.relative_juveniles,
+        metrics.coral_evenness,
+        metrics.relative_shelter_volume
+    ]
+
+    for metric in test_metrics
+        gd_diff::YAXArray{<:Real,2}, ug_diff::YAXArray{<:Real,2} = cf_difference_loc(
+            metric(TEST_RS), TEST_SCENS
+        )
+
+        for diff in [gd_diff, ug_diff]
+            @test all(diff[1, :] .<= diff[2, :] .<= diff[3, :])
+        end
+    end
+end

test/runtests.jl

@@ -22,9 +22,11 @@ include("mcda.jl")
 include("metrics/test_metrics_helper.jl")
 include("metrics/scenario.jl")
 include("metrics/metrics.jl")
+include("metrics/spatial.jl")
 include("utils/scale.jl")
 include("utils/text_display.jl")
 include("viz/taxa_dynamics.jl")
+include("viz/spatial.jl")
 
 # TODO Fix spatial_clustering and site_selection tests
 # include("site_selection.jl")

test/viz/spatial.jl

@@ -0,0 +1,22 @@
+using WGLMakie, GeoMakie, GraphMakie
+using ADRIA: viz
+using ADRIA: YAXArray
+using ADRIA: metrics, metrics.cf_difference_loc
+
+if !@isdefined(TEST_RS)
+    const TEST_DOM, TEST_N_SAMPLES, TEST_SCENS, TEST_RS = test_rs()
+end
+
+Makie.inline!(false)
+@testset "spatial" begin
+    fig_opts = Dict(:size => (1600, 800))
+
+    @testset "diff_map" begin
+        gd_diff::YAXArray{<:Real,2}, ug_diff::YAXArray{<:Real,2} = cf_difference_loc(
+            metrics.relative_cover(TEST_RS), TEST_SCENS
+        )
+
+        viz.diff_map(TEST_RS, gd_diff[2, :]; fig_opts=fig_opts)
+        viz.diff_map(TEST_RS, ug_diff[2, :]; fig_opts=fig_opts)
+    end
+end