Add DEA functions and plotting to ADRIA.analysis

[Rosejoycrocker]

Adds function to perform output-oriented (default) Data Envelopment Analysis (DEA) given inputs X and output metrics Y. DEA is used to measure the performance of entities (scenarios), where inputs are converted to outputs via some process. Each scenario's "efficiency score" is calculated relative to an "efficiency frontier", a region representing scenarios for which outputs
cannot be further increased by changing inputs (scenario settings). Scenarios on the frontier serve as "benchmarks" or "peers", associated with best practice restoration scenarios. Scenarios with efficiencies not equal to 1 can be improved to be more efficient.

Allows usage with generic input functions, such as functions representing cost as a function of intervention scenario.

Flagging @ConnectedSystems as we discussed this previously

Working example is:

using Revise, Infiltrator
using ADRIA
using DataFrames, YAXArrays
using GLMakie, GeoMakie, GraphMakie
using Distributions

dom = ADRIA.load_domain( "domain path", "45")

scens = ADRIA.sample(dom, 2^10)
rs = ADRIA.run_scenarios(dom, scens, "45")

n_scens = size(scens,1)

# Get cost of deploying corals in each scenario, randomised cost example
cost = YAXArray(collect(range(100,stop=1000000, length=n_scens)) .+ rand(Uniform(1000, 2000), n_scens))

# Get mean coral cover and shelter volume for each scenario
s_tac = dropdims(
    mean(ADRIA.metrics.scenario_total_cover(rs); dims=:timesteps); dims=:timesteps
)
s_sv = dropdims(
   mean(ADRIA.metrics.scenario_shelter_volume(rs); dims=:timesteps); dims=:timesteps
    )

# Do output oriented DEA analysis seeking to maximise cover and shelter volume for minimum
# deployment cost.
DEA_scens = ADRIA.economics.data_envelopment_analysis(cost, Array(s_tac), Array(s_sv))
dea_fig = ADRIA.viz.data_envelopment_analysis(rs, DEA_scens)

[codecov]

Codecov Report

Attention: Patch coverage is 84.00000% with 8 lines in your changes missing coverage. Please review.

Project coverage is 51.77%. Comparing base (ad39467) to head (5a339e3).

Files with missing lines	Patch %	Lines
src/analysis/data_envelopment.jl	60.00%	6 Missing ⚠️
src/viz/viz.jl	0.00%	2 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #896      +/-   ##
==========================================
+ Coverage   51.44%   51.77%   +0.33%     
==========================================
  Files          72       74       +2     
  Lines        4782     4832      +50     
==========================================
+ Hits         2460     2502      +42     
- Misses       2322     2330       +8     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[Zapiano]

This PR doesn't add any new test. I think it would be a good idea to add at least some basic tests, at least the "happy case" both for the metrics and for the viz function added.

[Zapiano]

Two comments:

1. Maybe economics is not a good name for this module. We have other analysis "tools" that are also related to economics (like pareto).
2. I'm not sure this needs to be in its own module. Since we already have analysis module, why not just rename this file to data_envelopment.jl or something and include that in analysis module? (this is a real question, do you see a reason why this should be on a separate module? to me it looks like this is too small to deserve its own module, at least at this stage).

[Zapiano]

Suggested change

	function DEAResult(CRS_eff::Vector{Float64}, VRS_eff::Vector{Float64},
	FDH_eff::Vector{Float64}, CRS_peers::DEA.DEAPeers, VRS_peers::DEA.DEAPeers,
	FDH_peers::DEA.DEAPeers, X::Matrix{Float64}, Y::Vector{Float64}
	)::DEAResult
	function DEAResult(
	CRS_eff::Vector{Float64},
	VRS_eff::Vector{Float64},
	FDH_eff::Vector{Float64},
	CRS_peers::DEA.DEAPeers,
	VRS_peers::DEA.DEAPeers,
	FDH_peers::DEA.DEAPeers,
	X::Matrix{Float64},
	Y::Vector{Float64}
	)::DEAResult

[Zapiano]

Make this a single line and update (it doesn't match the actual function signature)

[Zapiano]

Suggested change

	return DEAResult(1 ./ CRS_eff,
	1 ./ VRS_eff,
	1 ./ FDH_eff,
	CRS_peers,
	VRS_peers,
	FDH_peers,
	X,
	Y)
	return DEAResult(
	1 ./ CRS_eff,
	1 ./ VRS_eff,
	1 ./ FDH_eff,
	CRS_peers,
	VRS_peers,
	FDH_peers,
	X,
	Y
	)

[Zapiano]

Please, make this inline and update (doesn't match the function signature)

[Zapiano]

axa? Maybe ax_a? Not only to avoid confusion but also because it seems that we use snake case everywhere in this project.

[Zapiano]

So, I my opinion calling a scatter plot data is not very good. Not only because data can mean almost anything but because this is a figure object (right?). Since you are calling this "data cloud" it could be data_cloud, for example..?

[Zapiano]

So, what's the difference between a data cloud and a scatter plot?

[Rosejoycrocker]

Sorry, what do you mean here?

[Zapiano]

Sorry, I was probably tired and being picky when I wrote this, you can just ignore this comment.

[Zapiano]

Suggested change

	axb = Axis(g[2, 1]; title="Scale efficiency", ylabel=scale_eff_y_lab, axis_opts...)
	ax_b = Axis(g[2, 1]; title="Scale efficiency", ylabel=scale_eff_y_lab, axis_opts...)

[Zapiano]

Suggested change

	axc = Axis(g[3, 1]; title="Technical efficiency", ylabel=tech_eff_y_lab, axis_opts...)
	ax_c = Axis(g[3, 1]; title="Technical efficiency", ylabel=tech_eff_y_lab, axis_opts...)

[Rosejoycrocker]

This PR doesn't add any new test. I think it would be a good idea to add at least some basic tests, at least the "happy case" both for the metrics and for the viz function added.

I've added the plotting and calculation functions to analysis.jl tests

[Rosejoycrocker]

I think I've addressed your comments @Zapiano, let me know if you have others

[Zapiano]

It's almost ready, just a few docstring missing and a refactor suggestion :)

[Zapiano]

Suggested change

	function ADRIA.viz.data_envelopment_analysis(
	rs::ResultSet, DEA_output::DEAResult;
	axis_opts::OPT_TYPE=DEFAULT_OPT_TYPE(), fig_opts::OPT_TYPE=DEFAULT_OPT_TYPE(),
	opts::OPT_TYPE=DEFAULT_OPT_TYPE()
	)
	function ADRIA.viz.data_envelopment_analysis(
	rs::ResultSet,
	DEA_output::DEAResult;
	axis_opts::OPT_TYPE=DEFAULT_OPT_TYPE(),
	fig_opts::OPT_TYPE=DEFAULT_OPT_TYPE(),
	opts::OPT_TYPE=DEFAULT_OPT_TYPE()
	)

[Zapiano]

This docstring doesn't match the function signature (all opts params types are wrong) and it should be one method per line (this is how we do in the resto of this project at least).

[Zapiano]

Suggested change

	function ADRIA.viz.data_envelopment_analysis!(g::Union{GridLayout,GridPosition},
	rs::ResultSet, DEA_output::DEAResult; axis_opts::OPT_TYPE=DEFAULT_OPT_TYPE(),
	opts::OPT_TYPE=DEFAULT_OPT_TYPE()
	)
	function ADRIA.viz.data_envelopment_analysis!(
	g::Union{GridLayout,GridPosition},
	rs::ResultSet,
	DEA_output::DEAResult;
	axis_opts::OPT_TYPE=DEFAULT_OPT_TYPE(),
	opts::OPT_TYPE=DEFAULT_OPT_TYPE()
	)

[Zapiano]

Suggested change

	function ADRIA.viz.data_envelopment_analysis!(g::Union{GridLayout,GridPosition},
	DEA_output::DEAResult; axis_opts::OPT_TYPE=DEFAULT_OPT_TYPE(),
	opts::OPT_TYPE=DEFAULT_OPT_TYPE())
	function ADRIA.viz.data_envelopment_analysis!(
	g::Union{GridLayout,GridPosition},
	DEA_output::DEAResult; axis_opts::OPT_TYPE=DEFAULT_OPT_TYPE(),
	opts::OPT_TYPE=DEFAULT_OPT_TYPE()
	)

[Zapiano]

I think there's a way of making this more interesting, doing the same thing in less lines and extracting part of this code to another function:

If I understand, the reason why you need to save the result of these two following scatter! to variables is to render the Legend. You could also render the Legend in a different way:

marker_colors = [frontier_color, data_color]
marker_els = [MarkerElement(;color=m_color) for m_color in marker_colors]
Legend(g[1,2], marker_els, [frontier_name, data_name]; framevisible=false)

You can even put these three lines inside a separate function and pass opts to it so you don't need to extract these colors and names here.

And if you do that, you don't need to save those scatter! results to variables. Then, you can extract these lines that you repeat three times to an external function:

    # Plot the scale efficiency (ratio of efficiencies assuming CRS vs. assuming VRS)
    ax_b = Axis(g[2, 1]; title="Scale efficiency", ylabel=scale_eff_y_lab, axis_opts...)
    scatter!(ax_b, scale_efficiency; color=data_color)
    scatter!(
        ax_b,
        best_practice_scens,
        scale_efficiency[best_practice_scens];
        color=frontier_color
    )

And, again, you can pass opts to this function directly and extract the relevant variables inside it so you don't need to extract them here).

[Zapiano]

Still about the legend: I've done something similar here:

ADRIA.jl/ext/AvizExt/viz/scenarios.jl

Line 290 in 3acd2a7

function _render_legend(

And a small digression, hope you don't mind:

That makes me think we should have some standardized way of doing this the function I'm suggesting to you and the function in this link are essentially the same except that one generates a Legend with LineElements and the other generates a Legend with MarkerElements - and that's something we could address using multiple dispatch probably.

All that to say: maybe not do what I just said in this PR, but if you just extract this code that generates a legend (the one I suggested in the previous comment) to a function called render_marker_legend, for example, in the future we can make this better and have just a render_legend that can be used for different types of elements.

[Rosejoycrocker]

I'll add render_marker_legend to data_envelopment.jlbut I'm guessing this will eventually be moved to 'viz.jl' when multiple dispatch is incorporated, so it can be used by both yourscenario.jlfunctions and thedata_envelopment.jl` functions?

[Zapiano]

This one doesn't match the function signature.

[Rosejoycrocker]

I think I've addressed your comments @Zapiano, let me know if there is anything else :)