open-AIMS / ADRIA.jl

ADRIA: Adaptive Dynamic Reef Intervention Algorithms. A multi-criteria decision support platform for informing reef restoration and adaptation interventions.
MIT License
18 stars 6 forks source link

Use of surrogates for MORDM #255

Open ConnectedSystems opened 1 year ago

ConnectedSystems commented 1 year ago

It is possible to generate a surrogate from a sampled or pre-existing scenario ensemble $\varepsilon$. This could then be used to support Many Objective Robust Decision Making (MORDM) and related analyses.

In fact, much of the necessary infrastructure is already in place.

The example below has not been tested, but gives some idea of what could be done using the Surrogates.jl package.

An alternative is the SurrogateModelOptim which focuses on a Radial Basis Function approach.

using ADRIA
using Surrogates
using DataFrames, Statistics

rs = ADRIA.load_results("some result set")

spec = ADRIA.model_spec(rs)
y = vec(mean(ADRIA.metrics.scenario_total_cover(rs), dims=1))

# Build surrogate for total absolute cover under RCP45
scen_RCP45 = rs.inputs.RCP .== 45
RCP45_scens = Matrix(rs.inputs[scen_RCP45, Not(:RCP)])
y_RCP45 = y[scen_RCP45]

# Build a Wendland surrogate
surrogate = Wendland(RCP45_scens, y_RCP45, spec[:, :lower_bound], spec[:, :upper_bound])
ConnectedSystems commented 1 year ago

The script below using Polynomial Chaos seems to work.

Need to wrap this process with an interface and convenience function.

using DataFrames, Statistics
using Surrogates, SurrogatesPolyChaos, SurrogatesAbstractGPs

using ADRIA

dom = ADRIA.load_domain("path to latest Domain package")
X = ADRIA.sample(dom, 3072)

rs = ADRIA.run_scenarios(X, dom, "45")
spec = ADRIA.model_spec(rs)

y = vec(mean(ADRIA.metrics.scenario_total_cover(rs), dims=1))

# Select only RCP45 scenarios
scen_RCP45 = rs.inputs.RCP .== 45
RCP45_scens = rs.inputs[scen_RCP45, Not(:RCP)]

# Subset down to the factors of interest
# these are a selection of factors which the quantity of interest (total absolute cover) was 
# found to be sensitive to through the analyses done for the Nov 2022 report.
# Considering more factors requires larger samples
foi = [:SRM, :fogging, :a_adapt, :seed_CA, :seed_years, :shade_years, :seed_year_start, :shade_year_start, :n_adapt, :dhw_scenario]
f_subset = Matrix(RCP45_scens[:, foi])
y_RCP45 = y[scen_RCP45]

foi_spec = spec[spec.fieldname .∈ [string.(foi)], :]

# Build surrogate for total absolute cover under RCP45
deg = 4  # try to fit a four degree polynomial
lb = foi_spec[:, :lower_bound]

# Have to build a GaussOrthoPoly for each dimension
mop = SurrogatesPolyChaos.MultiOrthoPoly([SurrogatesPolyChaos.GaussOrthoPoly(deg) for j in 1:length(lb)], deg)
surrogate = PolynomialChaosSurrogate(Tuple.(eachrow(f_subset)), y_RCP45, lb, foi_spec[:, :upper_bound], 
                op=mop)

ADRIA.performance.RMSE(surrogate.(eachrow(f_subset)), y_RCP45)
# 70015.49427951015
# At a rough guesstimate, this is a ~2 - 4% error

Note:

Evaluating the 3072 samples takes about 9mins on my laptop, compared to:

julia> @time surrogate.(eachrow(f_subset))
 10.597233 seconds (166.89 M allocations: 10.681 GiB, 15.10% gc time)
3072-element Vector{Float64}:
 2.1855677667117477e6
 2.2279093276313795e6
 2.124331817734823e6
 2.221577126168572e6
 2.0974871436893954e6
 2.218216664100222e6
 2.2905356431980426e6
 ⋮