Fix rule extraction to handle constant factors

docs/src/usage/analysis.md

@@ -405,18 +405,18 @@ scenario, it is possible to use a Rule Induction algorithm (SIRUS) and plot each
 rule as a scatter graph:
 
 ```julia
-# Select only desired features
-fields_iv = ADRIA.component_params(rs, [Intervention, CriteriaWeights]).fieldname
-scenarios_iv = scens[:, fields_iv]
+# Select features of interest
+foi = ADRIA.component_params(rs, [Intervention, SeedCriteriaWeights]).fieldname
 
 # Use SIRUS algorithm to extract rules
 max_rules = 10
-rules_iv = ADRIA.analysis.cluster_rules(target_clusters, scenarios_iv, max_rules)
+rules_iv = ADRIA.analysis.cluster_rules(dom, target_clusters, scens, foi, max_rules)
+
 
 # Plot scatters for each rule highlighting the area selected them
 rules_scatter_fig = ADRIA.viz.rules_scatter(
     rs,
-    scenarios_iv,
+    scens,
     target_clusters,
     rules_iv;
     fig_opts=fig_opts,

ext/AvizExt/outcome_metadata.jl

@@ -0,0 +1,35 @@
+function outcome_title(outcomes::YAXArray)::String
+    return get(outcomes.properties, :metric_name, "")
+end
+
+function set_plot_opts!(
+    outcomes::AbstractArray,
+    opts::OPT_TYPE,
+    opts_key::Symbol;
+    metadata_key::Symbol=:metric_feature,
+    label_case=titlecase
+)
+    if !haskey(opts, opts_key) && (outcomes isa YAXArray)
+        opts[opts_key] = outcome_label(
+            outcomes; metadata_key=metadata_key, label_case=label_case
+        )
+    end
+end
+
+function outcome_label(
+    outcomes::YAXArray; metadata_key::Symbol=:metric_feature, label_case=titlecase
+)::String
+    outcome_metadata = outcomes.properties
+
+    return if all(haskey.([outcome_metadata], [metadata_key, :metric_unit]))
+        _metric_feature = label_case(outcome_metadata[metadata_key])
+        _metric_label = if !isempty(outcome_metadata[:metric_unit])
+            "[$(outcome_metadata[:metric_unit])]"
+        else
+            ""
+        end
+        "$(_metric_feature) $(_metric_label)"
+    else
+        ""
+    end
+end

ext/AvizExt/viz/clustering.jl

@@ -79,6 +79,10 @@ function ADRIA.viz.clustered_scenarios(
     fig_opts::OPT_TYPE=DEFAULT_OPT_TYPE(),
     axis_opts::OPT_TYPE=DEFAULT_OPT_TYPE()
 )::Figure
+    if !haskey(axis_opts, :title)
+        axis_opts[:title] = "$(outcome_title(outcomes)) Clusters"
+    end
+
     return ADRIA.viz.scenarios(
         outcomes, clusters; opts=opts, fig_opts=fig_opts, axis_opts=axis_opts
     )
@@ -101,7 +105,7 @@ Visualize clustered time series for each location and map.
 
 # Arguments
 - `rs` : ResultSet
-- `data` : Vector of summary statistics data for each location
+- `loc_outcomes` : Vector of summary statistics data for each location
 - `clusters` : Vector of numbers corresponding to clusters
 - `opts` : Options specific to this plotting method
     - `highlight` : Vector of colors indicating cluster membership for each location.
@@ -113,22 +117,29 @@ Figure
 """
 function ADRIA.viz.map(
     rs::Union{Domain,ResultSet},
-    data::AbstractArray{<:Real},
+    loc_outcomes::AbstractVector{<:Real},
     clusters::Union{BitVector,AbstractVector{Int64}};
     opts::OPT_TYPE=DEFAULT_OPT_TYPE(),
     fig_opts::OPT_TYPE=DEFAULT_OPT_TYPE(),
-    axis_opts::OPT_TYPE=DEFAULT_OPT_TYPE()
+    axis_opts::OPT_TYPE=set_axis_defaults(DEFAULT_OPT_TYPE())
 )::Figure
+    # Setup fig_opts size before the other default fig_opts
+    fig_opts[:size] = get(fig_opts, :size, (800, 700))
+    set_figure_defaults(fig_opts)
+
     f = Figure(; fig_opts...)
     g = f[1, 1] = GridLayout()
-    ADRIA.viz.map!(g, rs, data, clusters; opts=opts, axis_opts=axis_opts)
+
+    set_plot_opts!(loc_outcomes, opts, :colorbar_label)
+
+    ADRIA.viz.map!(g, rs, loc_outcomes, clusters; opts=opts, axis_opts=axis_opts)
 
     return f
 end
 function ADRIA.viz.map!(
     g::Union{GridLayout,GridPosition},
     rs::Union{Domain,ResultSet},
-    data::AbstractVector{<:Real},
+    loc_outcomes::AbstractVector{<:Real},
     clusters::Union{BitVector,Vector{Int64}};
     opts::OPT_TYPE=DEFAULT_OPT_TYPE(),
     axis_opts::OPT_TYPE=DEFAULT_OPT_TYPE()
@@ -137,7 +148,7 @@ function ADRIA.viz.map!(
     loc_groups::Dict{Symbol,BitVector} = ADRIA.analysis.scenario_clusters(clusters)
     group_colors::Dict{Symbol,Union{Symbol,RGBA{Float32}}} = colors(loc_groups)
 
-    legend_params::Tuple = _cluster_legend_params(data, loc_groups, group_colors)
+    legend_params::Tuple = _cluster_legend_params(loc_outcomes, loc_groups, group_colors)
 
     _colors::Vector{Union{Symbol,RGBA{Float32}}} = Vector{Union{Symbol,RGBA{Float32}}}(
         undef, length(clusters)
@@ -150,7 +161,7 @@ function ADRIA.viz.map!(
     opts[:highlight] = get(opts, :highlight, _colors)
     opts[:legend_params] = get(opts, :legend_params, legend_params)
 
-    ADRIA.viz.map!(g, rs, data; opts=opts, axis_opts=axis_opts)
+    ADRIA.viz.map!(g, rs, loc_outcomes; opts=opts, axis_opts=axis_opts)
 
     return g
 end
@@ -187,7 +198,7 @@ function _cluster_legend_params(
 
     legend_labels =
         labels(group_keys) .* ": " .* ADRIA.to_scientific.(label_means, digits=2)
-    legend_title = "Cluster mean"
+    legend_title = "Cluster mean $(outcome_label(data; label_case=lowercase))"
 
     return (legend_entries, legend_labels, legend_title)
 end

ext/AvizExt/viz/rule_extraction.jl

@@ -34,12 +34,18 @@ function ADRIA.viz.rules_scatter(
     fig_opts::OPT_TYPE=DEFAULT_OPT_TYPE(),
     axis_opts::OPT_TYPE=DEFAULT_OPT_TYPE()
 )::Figure
-    f = Figure(; fig_opts...)
-    g = f[1, 1] = GridLayout()
-
     # For now we only plot conditions with two clauses
     rules = filter(r -> length(r.condition) == 2, rules)
 
+    # Setup figure size
+    n_rows = ceil(10 / 4)
+    base_width = 1200
+    base_height = n_rows * base_width / 5
+    fig_opts[:size] = get(fig_opts, :size, (base_width, base_height))
+
+    f = Figure(; fig_opts...)
+    g = f[1, 1] = GridLayout()
+
     ADRIA.viz.rules_scatter!(
         g, rs, scenarios, clusters, rules; opts=opts, axis_opts=axis_opts
     )
@@ -141,6 +147,11 @@ function ADRIA.viz.rules_scatter!(
         margin=(5, 5, 5, 5)
     )
 
+    Label(
+        g[1, :, Top()], "SIRUS extracted rules"; padding=(0, 0, 50, 0), font=:bold,
+        valign=:bottom
+    )
+
     return g
 end
 

ext/AvizExt/viz/scenarios.jl

@@ -129,6 +129,11 @@ function ADRIA.viz.scenarios!(
     # Ensure last year is always shown in x-axis
     xtick_vals = get(axis_opts, :xticks, _time_labels(timesteps(outcomes)))
     xtick_rot = get(axis_opts, :xticklabelrotation, 2 / π)
+
+    if !haskey(axis_opts, :title)
+        axis_opts[:title] = outcome_title(outcomes)
+    end
+
     ax = Axis(g[1, 1]; xticks=xtick_vals, xticklabelrotation=xtick_rot, axis_opts...)
 
     _scenarios = copy(scenarios[1:end .∈ [outcomes.scenarios], :])
@@ -180,8 +185,7 @@ function ADRIA.viz.scenarios!(
     end
 
     ax.xlabel = "Year"
-    # ax.ylabel = metric_label(metric)
-
+    ax.ylabel = outcome_label(outcomes)
     return g
 end
 
@@ -198,7 +202,7 @@ function _confints(
     agg_dim = symdiff(axes_names(outcomes), [:timesteps])[1]
     for (idx, group) in enumerate(group_names)
         confints[:, idx, :] = series_confint(
-            outcomes[:, scen_groups[group]]; agg_dim=agg_dim
+            outcomes.data[:, scen_groups[group]]; agg_dim=agg_dim
         )
     end
 

ext/AvizExt/viz/spatial.jl

@@ -182,6 +182,8 @@ function ADRIA.viz.map(
     fig_opts::OPT_TYPE=set_figure_defaults(DEFAULT_OPT_TYPE()),
     axis_opts::OPT_TYPE=set_axis_defaults(DEFAULT_OPT_TYPE())
 )
+    set_plot_opts!(y, opts, :colorbar_label; metadata_key=:metric_name)
+
     f = Figure(; fig_opts...)
     g = f[1, 1] = GridLayout()
 

ext/AvizExt/viz/taxa_dynamics.jl

@@ -20,12 +20,12 @@ ADRIA.viz.taxonomy(scenarios, relative_taxa_cover)
 # Arguments
 - `rs` : ADRIA result set
 - `scenarios` : Scenario specification
-- `relative_taxa_cover` : YAXArray of dimensions [timesteps ⋅ taxa ⋅ scenarios]
+- `relative_taxa_cover` : YAXArray of dimensions [timesteps ⋅ species ⋅ scenarios]
 - `opts` : Aviz options
     - `by_RCP` : Split plots by RCP otherwise split by scenario type. Defaults to false.
     - `by_functional_groups` : If true, split plots by scenario types, otherwise split by taxonomy. Defaults to true.
     - `show_confints` : Show confidence intervals around series. Defaults to true.
-    - `colors` : Colormap for each taxonomy or scenario type. Defaults to Set1_5 for taxa and ADRIA defaults for scenario type.
+    - `colors` : Colormap for each taxonomy or scenario type. Defaults to Set1_5 for species and ADRIA defaults for scenario type.
 - `axis_opts` : Additional options to pass to adjust Axis attributes.
   See: https://docs.makie.org/v0.19/api/index.html#Axis
 - `series_opts` : Additional options to pass to adjust Series attributes
@@ -103,7 +103,7 @@ function ADRIA.viz.taxonomy!(
     by_functional_groups::Bool = get(opts, :by_functional_groups, true)
     if by_functional_groups
         # Create colors
-        n_functional_groups::Int64 = length(relative_taxa_cover.taxa)
+        n_functional_groups::Int64 = length(relative_taxa_cover.species)
         default_color = Symbol("Set1_" * string(n_functional_groups))
         color = get(opts, :colors, default_color)
         _colors = categorical_colors(color, n_functional_groups)
@@ -139,6 +139,10 @@ function ADRIA.viz.taxonomy!(
         )
     end
 
+    Label(
+        g[1, :, Top()], "Taxa dynamics"; padding=(0, 0, 30, 0), font=:bold, valign=:bottom
+    )
+
     return g
 end
 
@@ -193,12 +197,13 @@ function taxonomy_by_intervention!(
     series_opts::OPT_TYPE=DEFAULT_OPT_TYPE()
 )::Nothing where {T<:Float32}
     n_timesteps::Int64 = length(relative_taxa_cover.timesteps)
-    n_functional_groups::Int64 = length(relative_taxa_cover.taxa)
+    functional_groups::Vector{Int64} = ADRIA.axis_labels(relative_taxa_cover, :species)
+    n_functional_groups::Int64 = length(functional_groups)
 
     # Plot and calculate confidence intervals
     confints = zeros(n_timesteps, n_functional_groups, 3)
-    for (idx, taxa) in enumerate(relative_taxa_cover.taxa)
-        confints[:, idx, :] = series_confint(relative_taxa_cover[taxa=At(taxa)])
+    for (idx, group) in enumerate(functional_groups)
+        confints[:, idx, :] = series_confint(relative_taxa_cover[species=At(group)])
         show_confints ?
         band!(
             ax, 1:n_timesteps, confints[:, idx, 1], confints[:, idx, 3];
@@ -247,7 +252,7 @@ function intervention_by_taxonomy!(
 
         intervention_by_taxonomy!(
             ax,
-            relative_taxa_cover[taxa=idx],
+            relative_taxa_cover[species=idx],
             colors,
             scen_groups;
             show_confints=show_confints,

ext/AvizExt/viz/viz.jl

@@ -85,6 +85,7 @@ function _calc_gridsize(n_factors::Int64; max_cols::Int64=4)::Tuple{Int64,Int64}
     return n_rows, n_cols
 end
 
+include("../outcome_metadata.jl")
 include("scenarios.jl")
 include("sensitivity.jl")
 include("clustering.jl")

src/analysis/rule_extraction.jl

@@ -103,16 +103,18 @@ function print_rules(rules::Vector{Rule{Vector{Vector},Vector{Float64}}})::Nothi
 end
 
 """
-    cluster_rules(clusters::Vector{T}, X::DataFrame, max_rules::T; seed::Int64=123, kwargs...) where {T<:Integer,F<:Real}
-    cluster_rules(clusters::Union{BitVector,Vector{Bool}}, X::DataFrame, max_rules::T; kwargs...) where {T<:Int64}
+    cluster_rules(domain::Domain, clusters::Vector{T}, scenarios::DataFrame, factors::Vector{Symbol}, max_rules::T; seed::Int64=123, kwargs...) where {T<:Integer,F<:Real}
+    cluster_rules(domain::Domain, clusters::Union{BitVector,Vector{Bool}}, scenarios::DataFrame, factors::Vector{Symbol}, max_rules::T; kwargs...) where {T<:Int64}
 
 Use SIRUS package to extract rules from time series clusters based on some summary metric
 (default is median). More information about the keyword arguments accepeted can be found in
 MLJ's doc (https://juliaai.github.io/MLJ.jl/dev/models/StableRulesClassifier_SIRUS/).
 
 # Arguments
+- `domain` : Domain
 - `clusters` : Vector of cluster indexes for each scenario outcome
-- `X` : Features to be used as input by SIRUS
+- `scenarios` : Scenarios DataFrame
+- `factors` : Vector of factors of interest
 - `max_rules` : Maximum number of rules, to be used as input by SIRUS
 - `seed` : Seed to be used by RGN
 
@@ -129,8 +131,24 @@ A StableRules object (implemented by SIRUS).
    Electron. J. Statist. 15 (1) 427 - 505.
    https://doi.org//10.1214/20-EJS1792
 """
-function cluster_rules(clusters::Vector{T}, X::DataFrame, max_rules::T;
-    seed::Int64=123, kwargs...) where {T<:Int64}
+function cluster_rules(
+    domain::ADRIA.Domain,
+    clusters::Vector{T},
+    scenarios::DataFrame,
+    factors::Vector{Symbol},
+    max_rules::T;
+    seed::Int64=123,
+    kwargs...
+) where {T<:Int64}
+    ms = ADRIA.model_spec(domain)
+    variable_factors_filter::BitVector = .!ms[ms.fieldname .∈ [factors], :is_constant]
+    variable_factors::Vector{Symbol} = factors[variable_factors_filter]
+
+    if isempty(variable_factors)
+        throw(ArgumentError("Factors of interest have to be non constant."))
+    end
+
+    X = scenarios[:, variable_factors]
     # Set seed and Random Number Generator
     rng = StableRNG(seed)
 
@@ -150,9 +168,17 @@ function cluster_rules(clusters::Vector{T}, X::DataFrame, max_rules::T;
 
     return rules(mach.fitresult)
 end
-function cluster_rules(clusters::Union{BitVector,Vector{Bool}}, X::DataFrame, max_rules::T;
-    kwargs...) where {T<:Int64}
-    return cluster_rules(convert.(Int64, clusters), X, max_rules; kwargs...)
+function cluster_rules(
+    domain::ADRIA.Domain,
+    clusters::Union{BitVector,Vector{Bool}},
+    scenarios::DataFrame,
+    factors::Vector{Symbol},
+    max_rules::T;
+    kwargs...
+) where {T<:Int64}
+    return cluster_rules(
+        domain, convert.(Int64, clusters), scenarios, factors, max_rules; kwargs...
+    )
 end
 
 """