diff --git a/sandbox/IHACRESCMDStateNode.jl b/sandbox/IHACRESCMDStateNode.jl deleted file mode 100644 index 80d3007..0000000 --- a/sandbox/IHACRESCMDStateNode.jl +++ /dev/null @@ -1,279 +0,0 @@ -# IHACRES CMD State-based Node -using Streamfall -using ModelParameters, OnlineStats, Setfield -import Streamfall: run_node!, update_params! - - -Base.@kwdef mutable struct CMDStateNode{P, A<:Real} <: IHACRESNode - Streamfall.@network_node - - # Activate parameters - d::A = 200.0; - d2::A = 2.0 - e::A = 1.0 - f::A = 0.8 - a::A = 0.9 - b::A = 0.1 - storage_coef::A = 2.9 - alpha::A = 0.1 - - # low CMD (wet) - low_d::P = Param(200.0, bounds=(10.0, 550.0)) # flow threshold - low_d2::P = Param(2.0, bounds=(0.0001, 10.0)) # flow threshold2 - low_e::P = Param(1.0, bounds=(0.1, 1.5)) # temperature to PET conversion factor - low_f::P = Param(0.8, bounds=(0.01, 3.0)) # plant stress threshold factor (multiplicative factor of d) - low_a::P = Param(0.9, bounds=(0.1, 10.0)) # quickflow storage coefficient == (1/tau_q) - low_b::P = Param(0.1, bounds=(1e-3, 0.1)) # slowflow storage coefficent == (1/tau_s) - low_storage_coef::P = Param(2.9, bounds=(1e-10, 10.0)) - low_alpha::P = Param(0.1, bounds=(1e-5, 1 - 1/10^9)) - - # mid CMD ("usual") - mid_d::P = Param(200.0, bounds=(10.0, 550.0)) # flow threshold - mid_d2::P = Param(2.0, bounds=(0.0001, 10.0)) # flow threshold2 - mid_e::P = Param(1.0, bounds=(0.1, 1.5)) # temperature to PET conversion factor - mid_f::P = Param(0.8, bounds=(0.01, 3.0)) # plant stress threshold factor (multiplicative factor of d) - mid_a::P = Param(0.9, bounds=(0.1, 10.0)) # quickflow storage coefficient == (1/tau_q) - mid_b::P = Param(0.1, bounds=(1e-3, 0.1)) # slowflow storage coefficent == (1/tau_s) - mid_storage_coef::P = Param(2.9, bounds=(1e-10, 10.0)) - mid_alpha::P = Param(0.1, bounds=(1e-5, 1 - 1/10^9)) - - # high CMD ("dry") - high_d::P = Param(200.0, bounds=(10.0, 550.0)) # flow threshold - high_d2::P = Param(2.0, bounds=(0.0001, 10.0)) # flow threshold2 - high_e::P = Param(1.0, bounds=(0.1, 1.5)) # temperature to PET conversion factor - high_f::P = Param(0.8, bounds=(0.01, 3.0)) # plant stress threshold factor (multiplicative factor of d) - high_a::P = Param(0.9, bounds=(0.1, 10.0)) # quickflow storage coefficient == (1/tau_q) - high_b::P = Param(0.1, bounds=(1e-3, 0.1)) # slowflow storage coefficent == (1/tau_s) - high_storage_coef::P = Param(2.9, bounds=(1e-10, 10.0)) - high_alpha::P = Param(0.1, bounds=(1e-5, 1 - 1/10^9)) - - level_params::Array{P, 1} = [ - Param(-0.01, bounds=(-10.0, -0.01)), # p1 - Param(0.8, bounds=(0.0, 1.5)), # p2 - Param(4.5, bounds=(0.0, 20.0)), # p3 - Param(5.0, bounds=(1.0, 10.0)), # p4 - Param(0.35, bounds=(0.0, 1.0)), # p5 - Param(1.41, bounds=(-2.0, 2.0)), # p6 - Param(-1.45, bounds=(-2.5, 0.0)), # p7 - Param(6.75, bounds=(0.0, 10.0)), # p8 - Param(150.0, bounds=(50.0, 200.0)) # ctf - ] - - storage::Array{A} = [100.0] # CMD - quick_store::Array{A} = [0.0] - slow_store::Array{A} = [0.0] - outflow::Array{A} = [] - effective_rainfall::Array{A} = [] - et::Array{A} = [] - inflow::Array{A} = [] - level::Array{A} = [] - gw_store::Array{A} = [0.0] - - # cache arrays - cache::NamedTuple{(:i_cache, :r_cache), <:Tuple{Vector{Float64}, Vector{Float64}}} = (i_cache=zeros(3), r_cache=zeros(2)) - - # active param records - active_param_set::Array{Int64} = [1] -end - - -"""Generic run_node method that handles any single state variable. - -Develops an "online" running statistic of states from the first N days, -which remain static for the rest of the analysis period. -""" -function run_node!(node::CMDStateNode, climate::Climate; quantiles::Array=[0.1, 0.9], burn_in=1826, releases=nothing, log_transform=false)::Nothing - - # quantiles = [0.0, 0.1, 0.9] - # new quantiles = [0.1, 0.9] - - # Get node parameters - _, x0, __ = param_info(node; with_level=false) - # num_params = length(x0) - - timesteps = length(climate) - active_param_set = zeros(Int, timesteps) - - # state_var = getfield(node, state) - state_var = node.storage - - o = Quantile(quantiles, b=2000) - param_idxs = nothing - thresholds = nothing - for ts in (1:timesteps) - # Update param set based on state - state_val = state_var[ts] - if log_transform - # add constant to avoid log(0) - tmp = log(state_val .+ 1e-2) - else - tmp = state_val - end - - if ts < burn_in - # Update thresholds based on state value so far - fit!(o, tmp) - thresholds = value(o) - end - - active_param_set[ts] = switch_params!(node, tmp, thresholds) - - # param_set_id, node_params, param_idxs = find_state_vars(tmp, thresholds, params, num_params, length(thresholds)) - # if (ts == 1) || (param_set_id != active_param_set[ts-1]) - # update_params!(node, node_params...) - # end - - # record timestep in which this param was active - # active_param_set[ts] = param_set_id - - Streamfall.run_node!(node, climate, ts; extraction=releases) - end - - node.active_param_set = active_param_set - - # return active_param_set, param_idxs - return nothing -end - - -""" -Clunky approach to switching between parameter sets -""" -function switch_params!(node::CMDStateNode, state_val::T, thresholds::Vector{T})::Int64 where {T<:Real} - set_id = first(searchsortedfirst.(Ref(thresholds), state_val)) - if set_id == node.active_param_set[end] - return set_id - end - - if set_id == 1 - node.d = node.low_d.val - node.d2 = node.low_d2.val - node.e = node.low_e.val - node.f = node.low_f.val - node.a = node.low_a.val - node.b = node.low_b.val - node.storage_coef = node.low_storage_coef.val - node.alpha = node.low_alpha.val - elseif set_id == 2 - node.d = node.mid_d.val - node.d2 = node.mid_d2.val - node.e = node.mid_e.val - node.f = node.mid_f.val - node.a = node.mid_a.val - node.b = node.mid_b.val - node.storage_coef = node.mid_storage_coef.val - node.alpha = node.mid_alpha.val - elseif set_id == 3 - node.d = node.high_d.val - node.d2 = node.high_d2.val - node.e = node.high_e.val - node.f = node.high_f.val - node.a = node.high_a.val - node.b = node.high_b.val - node.storage_coef = node.high_storage_coef.val - node.alpha = node.high_alpha.val - else - error("Unknown state id") - end - - return set_id -end - - -""" - -Clunky implementation of updating parameters. -On update, -""" -function update_params!(node::CMDStateNode, params...)::Nothing - node.d = params[1] - node.d2 = params[2] - node.e = params[3] - node.f = params[4] - node.a = params[5] - node.b = params[6] - node.storage_coef = params[7] - node.alpha = params[8] - - node.low_d = Param(params[1], bounds=node.low_d.bounds::Tuple) - node.low_d2 = Param(params[2], bounds=node.low_d2.bounds::Tuple) - node.low_e = Param(params[3], bounds=node.low_e.bounds::Tuple) - node.low_f = Param(params[4], bounds=node.low_f.bounds::Tuple) - node.low_a = Param(params[5], bounds=node.low_a.bounds::Tuple) - node.low_b = Param(params[6], bounds=node.low_b.bounds::Tuple) - node.low_storage_coef = Param(params[7], bounds=node.low_storage_coef.bounds::Tuple) - node.low_alpha = Param(params[8], bounds=node.low_alpha.bounds::Tuple) - - node.mid_d = Param(params[9], bounds=node.mid_d.bounds::Tuple) - node.mid_d2 = Param(params[10], bounds=node.mid_d2.bounds::Tuple) - node.mid_e = Param(params[11], bounds=node.mid_e.bounds::Tuple) - node.mid_f = Param(params[12], bounds=node.mid_f.bounds::Tuple) - node.mid_a = Param(params[13], bounds=node.mid_a.bounds::Tuple) - node.mid_b = Param(params[14], bounds=node.mid_b.bounds::Tuple) - node.mid_storage_coef = Param(params[15], bounds=node.mid_storage_coef.bounds::Tuple) - node.mid_alpha = Param(params[16], bounds=node.mid_alpha.bounds::Tuple) - - node.high_d = Param(params[17], bounds=node.high_d.bounds::Tuple) - node.high_d2 = Param(params[18], bounds=node.high_d2.bounds::Tuple) - node.high_e = Param(params[19], bounds=node.high_e.bounds::Tuple) - node.high_f = Param(params[20], bounds=node.high_f.bounds::Tuple) - node.high_a = Param(params[21], bounds=node.high_a.bounds::Tuple) - node.high_b = Param(params[22], bounds=node.high_b.bounds::Tuple) - node.high_storage_coef = Param(params[23], bounds=node.high_storage_coef.bounds::Tuple) - node.high_alpha = Param(params[24], bounds=node.high_alpha.bounds::Tuple) - - # Below doesn't seem to work yet... - # @set node.low_d.val = params[1] - # @set node.low_d2.val = params[2] - # @set node.low_e.val = params[3] - # @set node.low_f.val = params[4] - # @set node.low_a.val = params[5] - # @set node.low_b.val = params[6] - # @set node.low_storage_coef.val = params[7] - # @set node.low_alpha.val = params[8] - - # @set node.mid_d.val = params[9] - # @set node.mid_d2.val = params[10] - # @set node.mid_e.val = params[11] - # @set node.mid_f.val = params[12] - # @set node.mid_a.val = params[13] - # @set node.mid_b.val = params[14] - # @set node.mid_storage_coef.val = params[15] - # @set node.mid_alpha.val = params[16] - - # @set node.high_d.val = params[17] - # @set node.high_d2.val = params[18] - # @set node.high_e.val = params[19] - # @set node.high_f.val = params[20] - # @set node.high_a.val = params[21] - # @set node.high_b.val = params[22] - # @set node.high_storage_coef.val = params[23] - # @set node.high_alpha.val = params[24] - - node.active_param_set = [1] - - return nothing -end - - -function run_node!(node::CMDStateNode, climate::Climate; inflow=nothing, extraction=nothing, exchange=nothing) - invoke(run_node!, Tuple{IHACRESNode, Climate}, node, climate; inflow, extraction, exchange) -end - - -# function find_state_vars(val, thresholds, params, num_params, n_states) -# # index of parameters (e.g., if 3 partitions, then create 3 arrays of 8 parameters) -# param_idxs = collect(Iterators.partition(1:length(params), num_params)) - -# # searchsortedfirst.(Ref([3, 5, 10]), [2, 2, 2, 12]) - -# # TODO: Check if we need to exclude first bin -# set_id = searchsortedfirst.(Ref(thresholds), val) - - - -# param_set = param_idxs[set_id] -# node_params = params[param_set] - -# return set_id, node_params, param_idxs -# end diff --git a/sandbox/Project.toml b/sandbox/Project.toml deleted file mode 100644 index 81648c0..0000000 --- a/sandbox/Project.toml +++ /dev/null @@ -1 +0,0 @@ -[deps] diff --git a/sandbox/ensemble.jl b/sandbox/ensemble.jl deleted file mode 100644 index 0903b61..0000000 --- a/sandbox/ensemble.jl +++ /dev/null @@ -1,204 +0,0 @@ -using Streamfall -using Statistics, CSV, DataFrames -using Distributed, BlackBoxOptim -using Base.Iterators -using Plots -import Streamfall.Analysis: TemporalCrossSection - - -# Spin up workers if needed -if nprocs() == 1 - addprocs(4, exeflags="--project=..") - @everywhere begin - # Activate environment and precompile - # using Pkg; Pkg.activate(joinpath(@__DIR__, "..")) - # Pkg.instantiate(); Pkg.precompile() - - using Streamfall - using Statistics, CSV, DataFrames - end -end - - -@everywhere begin - -function v_metric(obs, sim; sp=365, metric=Streamfall.NmKGE) - res = Streamfall.naive_split_metric(obs, sim, sp, metric) - m = mean(res) - cv = std(res) / m - - V = (m*0.9) + ((ℯ^-cv)*0.1) # (ℯ .- 1.5) - - return V -end - - -function v_mod_metric(obs, sim; sp=365, metric=Streamfall.NmKGE) - res = Streamfall.naive_split_metric(obs, sim, sp, metric) - m = mean(res) - sd = std(res) - - V = (m*0.9) + ((ℯ^-sd)*0.1) - - return V -end - - -# load observations -HERE = @__DIR__ -DATA_PATH = joinpath(HERE, "../test/data/cotter/") - -date_format = "YYYY-mm-dd" -obs_data = CSV.File(joinpath(DATA_PATH, "climate/CAMELS-AUS_410730.csv"), - comment="#", - dateformat=date_format) |> DataFrame - -hist_streamflow = obs_data[:, ["Date", "410730_Q"]] -hist_streamflow = rename(hist_streamflow, "410730_Q"=>"410730") - -climate_data = obs_data[:, ["Date", "410730_P", "410730_PET"]] -climate = Climate(climate_data, "_P", "_PET") - -# Create nodes/ensemble -ihacres_node = create_node(BilinearNode, "410730", 129.2) -gr4j_node = create_node(GR4JNode, "410730", 129.2) -# ensemble_node = BaseEnsemble([ihacres_node, gr4j_node], weights=[0.5, 0.5]) - - -burn_in = 365 - -split_v = (obs, sim) -> 1.0 - v_metric(obs[burn_in:end], sim[burn_in:end]; sp=14, metric=Streamfall.NnpKGE) -split_mod_v = (obs, sim) -> 1.0 - v_mod_metric(obs[burn_in:end], sim[burn_in:end]; sp=14, metric=Streamfall.NnpKGE) - -end # everywhere block - - -name_nodes = zip(["IHACRES", "GR4J"], [ihacres_node, gr4j_node]) -name_metrics = zip(["split_v", "split_mod_v", "NnpKGE", "NmKGE"], [split_v, split_mod_v, Streamfall.NnpKGE, Streamfall.NmKGE]) -node_metric = [NamedTuple{(:name, :node, :met_name, :metric)}([name, node, met_name, metric]) for ((name, node), (met_name, metric)) in product(name_nodes, name_metrics)] - - -# res, opt = calibrate!(ihacres_node, climate, hist_streamflow; metric=split_v, MaxTime=1200.0, PopulationSize=500) -# res, opt = calibrate!(gr4j_node, climate, hist_streamflow; metric=split_v, MaxTime=1200.0, PopulationSize=500) - - -# p_res = pmap((n, m) -> calibrate!(n, climate, hist_streamflow; metric=m, MaxTime=1200.0)) - -p_res = pmap((x) -> calibrate!(x.node, climate, hist_streamflow; metric=x.metric, MaxTime=600.0, TraceMode=:silent), node_metric) - - -function postprocess(x, ω) - x = copy(x) - idx = (x .+ ω) .>= 0.0 - x[idx] = x[idx] .+ ω[idx] - - # idx2 = (x .+ ω) .< 0.0 - # x[idx2] .= 0.0 - return x -end - - -c_date = climate_data.Date -obs = hist_streamflow[:, "410730"] -figs = [] -outflows = [] -display_timeframe = 13330:(13330+(365*1)) -for (x, (r, o)) in zip(node_metric, p_res) - update_params!(x.node, best_candidate(r)...) - reset!(x.node) - run_node!(x.node, climate) - - t_name = x.name - m_name = x.met_name - push!(outflows, x.node.outflow) - - mae = round(Streamfall.MAE(obs, x.node.outflow), digits=4) - push!(figs, Streamfall.temporal_cross_section(c_date, obs, x.node.outflow; yscale=:log10, title="$t_name ($m_name | $mae)")) - - mod_xsect_res = TemporalCrossSection(c_date, obs, x.node.outflow) - offset_vals = mod_xsect_res.offsets - - pp_outflow = postprocess(x.node.outflow, offset_vals) - mae = round(Streamfall.MAE(obs, pp_outflow), digits=4) - push!(figs, Streamfall.temporal_cross_section(c_date, obs, pp_outflow; yscale=:log10, title="Post-Processed | $mae")) - - t_plot = plot(c_date[display_timeframe], obs[display_timeframe], label="Observed") - plot!(c_date[display_timeframe], x.node.outflow[display_timeframe], label="Modeled") - - push!(figs, t_plot) - - pp_t_plot = plot(c_date[display_timeframe], obs[display_timeframe], label="Observed") - plot!(c_date[display_timeframe], pp_outflow[display_timeframe], label="Post-Processed") - - push!(figs, pp_t_plot) -end - -# ensemble_node = BaseEnsemble([ihacres_node, gr4j_node], weights=[0.5, 0.5]) -# calibrate!(ensemble_node, climate, hist_streamflow; metric=metric, MaxTime=120.0, TraceMode=:silent) - - -plot(figs..., size=((300*4),(250*8)), layout=(8,4), dpi=300) # , link=:y -savefig("metric_calibration_comparison_2022-04-31__90_10_weight_1_year.png") - -# TODO: Try heatmap of mean error over day of year - -# TODO: Reincorporate ensemble - -# TODO: Try out state-based implementation - -# quickplot(obs, outflows[2], climate, "", false; burn_in=13330, limit=13330+3650) - -# # p_res will have pairs of res and opt, where res are the results and opt are the options used -# # bs = best_candidate(res) # extract best parameters -# # update_params!(node, bs...) # update node with parameters - -# # calibrate!(ensemble_node, climate, hist_streamflow; metric=split_v, MaxTime=1500.0, PopulationSize=500) - -# run_node!(ensemble_node, climate) - -# obs = hist_streamflow[:, "410730"] - -# quickplot(obs, ensemble_node, climate; burn_in=365, metric=split_v) -# quickplot(obs, ihacres_node, climate; burn_in=365, metric=split_v) - -# c_date = climate_data.Date - -# ######### -# mod_xsect_res = TemporalCrossSection(c_date, obs, ensemble_node.outflow) -# offset_vals = mod_xsect_res.offsets - -# # fig2 = Streamfall.temporal_cross_section(c_date, obs, ensemble_node.outflow; yscale=:log10, title="NnpKGE") - -# function postprocess(x, ω) -# x = copy(x) -# idx = (x .+ ω) .>= 0.0 -# x[idx] = x[idx] .+ ω[idx] - -# # idx2 = (x .+ ω) .< 0.0 -# # x[idx2] .= 0.0 -# return x -# end - -# ######### - - -# fig1 = Streamfall.temporal_cross_section(c_date, obs, ihacres_node.outflow; yscale=:log10, title="IHACRES") -# fig2 = Streamfall.temporal_cross_section(c_date, obs, gr4j_node.outflow; yscale=:log10, title="GR4J") -# fig3 = Streamfall.temporal_cross_section(c_date, obs, ensemble_node.outflow; yscale=:log10, title="Ensemble (IHACRES-GR4J)") -# fig4 = Streamfall.temporal_cross_section(c_date, obs, postprocess(ensemble_node.outflow, offset_vals); yscale=:log10, title="Post-Processed (IHACRES-GR4J)") - -# plot(fig1, fig2, fig3, fig4, size=(1000,800), layout=(2,2), link=:y, dpi=300) - -# savefig("ensemble_pp_popsize_500_maxtime_900_modded_v_metric.png") - -# quickplot(obs, postprocess(ensemble_node.outflow, offset_vals), c_date) -# savefig("ensemble_qp_pp_popsize_500_modded_v_metric.png") - -# quickplot(obs, ensemble_node.outflow, c_date) -# savefig("ensemble_qp_ensemble_popsize_500_modded_v_metric.png") - -# quickplot(obs, ihacres_node.outflow, c_date) -# savefig("ensemble_qp_ihacres_popsize_500_modded_v_metric.png") - -# quickplot(obs, gr4j_node.outflow, c_date) -# savefig("ensemble_qp_gr4j_popsize_500_modded_v_metric.png") \ No newline at end of file diff --git a/sandbox/ensemble_dev.jl b/sandbox/ensemble_dev.jl deleted file mode 100644 index c294813..0000000 --- a/sandbox/ensemble_dev.jl +++ /dev/null @@ -1,175 +0,0 @@ -using Streamfall -using Statistics, CSV, DataFrames -using Distributed, BlackBoxOptim -using Base.Iterators -using Plots -import Streamfall.Analysis: TemporalCrossSection - - -# Spin up workers if needed -if nprocs() == 1 - addprocs(4, exeflags="--project=..") - @everywhere begin - # Activate environment and precompile - # using Pkg; Pkg.activate(joinpath(@__DIR__, "..")) - # Pkg.instantiate(); Pkg.precompile() - - using Streamfall - using Statistics, CSV, DataFrames - - include("IHACRESCMDStateNode.jl") - end -end - - -@everywhere begin - -function v_metric(obs, sim; sp=365, metric=Streamfall.NmKGE) - res = Streamfall.naive_split_metric(obs, sim, sp, metric) - m = mean(res) - cv = std(res) / m - - V = (m*0.9) + ((ℯ^-cv)*0.1) # (ℯ .- 1.5) - - return V -end - - -function v_mod_metric(obs, sim; sp=365, metric=Streamfall.NmKGE) - res = Streamfall.naive_split_metric(obs, sim, sp, metric) - m = mean(res) - sd = std(res) - - V = (m*0.9) + ((ℯ^-sd)*0.1) - - return V -end - - -# load observations -HERE = @__DIR__ -DATA_PATH = joinpath(HERE, "../test/data/cotter/") - -date_format = "YYYY-mm-dd" -obs_data = CSV.File(joinpath(DATA_PATH, "climate/CAMELS-AUS_410730.csv"), - comment="#", - dateformat=date_format) |> DataFrame - -hist_streamflow = obs_data[:, ["Date", "410730_Q"]] -hist_streamflow = rename(hist_streamflow, "410730_Q"=>"410730") - -climate_data = obs_data[:, ["Date", "410730_P", "410730_PET"]] -climate = Climate(climate_data, "_P", "_PET") - -# Create nodes/ensemble -ihacres_node = create_node(BilinearNode, "410730", 129.2) -gr4j_node = create_node(GR4JNode, "410730", 129.2) - -state_node = create_node(CMDStateNode, "410730", 129.2) - - -burn_in = 1826 - -split_v = (obs, sim) -> 1.0 - v_metric(obs[burn_in:end], sim[burn_in:end]; sp=14, metric=Streamfall.NnpKGE) -split_mod_v = (obs, sim) -> 1.0 - v_mod_metric(obs[burn_in:end], sim[burn_in:end]; sp=14, metric=Streamfall.NnpKGE) - -end # everywhere block - - -name_nodes = zip(["IHACRES", "GR4J", "cmd_state"], [ihacres_node, gr4j_node, state_node]) -name_metrics = zip(["split_v", "split_mod_v", "NnpKGE", "NmKGE"], [split_v, split_mod_v, Streamfall.NnpKGE, Streamfall.NmKGE]) -node_metric = [NamedTuple{(:name, :node, :met_name, :metric)}([name, node, met_name, metric]) for ((name, node), (met_name, metric)) in product(name_nodes, name_metrics)] - -p_res = pmap((x) -> calibrate!(x.node, climate, hist_streamflow; metric=x.metric, MaxTime=1800.0, TraceMode=:silent), node_metric) - - -function postprocess(x, ω) - x = copy(x) - idx = (x .+ ω) .>= 0.0 - x[idx] = x[idx] .+ ω[idx] - - # x .= max.(x .+ ω, 0.0) - return x -end - - -c_date = climate_data.Date -obs = hist_streamflow[:, "410730"] -figs = [] -outflows = [] -display_timeframe = 13330:(13330+(365*1)) -for (x, (r, o)) in zip(node_metric, p_res) - update_params!(x.node, best_candidate(r)...) - reset!(x.node) - run_node!(x.node, climate) - - t_name = x.name - m_name = x.met_name - push!(outflows, x.node.outflow) - - mae = round(Streamfall.MAE(obs, x.node.outflow), digits=4) - push!(figs, Streamfall.temporal_cross_section(c_date, obs, x.node.outflow; yscale=:log10, title="$t_name ($m_name | $mae)")) - - mod_xsect_res = TemporalCrossSection(c_date, obs, x.node.outflow) - offset_vals = mod_xsect_res.offsets - - pp_outflow = postprocess(x.node.outflow, offset_vals) - mae = round(Streamfall.MAE(obs, pp_outflow), digits=4) - push!(figs, Streamfall.temporal_cross_section(c_date, obs, pp_outflow; yscale=:log10, title="Post-Processed | $mae")) - - t_plot = plot(c_date[display_timeframe], obs[display_timeframe], label="Observed") - plot!(c_date[display_timeframe], x.node.outflow[display_timeframe], label="Modeled") - push!(figs, t_plot) - - pp_t_plot = plot(c_date[display_timeframe], obs[display_timeframe], label="Observed") - plot!(c_date[display_timeframe], pp_outflow[display_timeframe], label="Post-Processed") - push!(figs, pp_t_plot) -end - - -# split_mod_v2 = (obs, sim) -> 1.0 - v_mod_metric(obs[burn_in:end], sim[burn_in:end]; sp=14, metric=Streamfall.NnpKGE) - - -## Ensemble modelling section -# ensemble_node = BaseEnsemble([node_metric[3].node, node_metric[4].node], weights=[0.5, 0.5]) -# reset!(ensemble_node) -# e_res, e_opt = calibrate!(ensemble_node, climate, hist_streamflow; metric=split_mod_v, MaxTime=1200.0, TraceMode=:silent) - -# update_params!(ensemble_node, best_candidate(e_res)...) -# reset!(ensemble_node) -# run_node!(ensemble_node, climate) - -# push!(outflows, ensemble_node.outflow) - -# mae = round(Streamfall.MAE(obs, ensemble_node.outflow), digits=4) -# push!(figs, Streamfall.temporal_cross_section(c_date, obs, ensemble_node.outflow; yscale=:log10, title="Ensemble (split_mod_v | $mae)")) - -# mod_xsect_res = TemporalCrossSection(c_date, obs, ensemble_node.outflow) -# offset_vals = mod_xsect_res.offsets - -# pp_outflow = postprocess(ensemble_node.outflow, offset_vals) -# mae = round(Streamfall.MAE(obs, pp_outflow), digits=4) -# push!(figs, Streamfall.temporal_cross_section(c_date, obs, pp_outflow; yscale=:log10, title="Post-Processed | $mae")) - -# t_plot = plot(c_date[display_timeframe], obs[display_timeframe], label="Observed") -# plot!(c_date[display_timeframe], ensemble_node.outflow[display_timeframe], label="Modeled") -# push!(figs, t_plot) - -# pp_t_plot = plot(c_date[display_timeframe], obs[display_timeframe], label="Observed") -# plot!(c_date[display_timeframe], pp_outflow[display_timeframe], label="Post-Processed") -# push!(figs, pp_t_plot) - - -plot(figs..., size=((350*4),(300*12)), layout=(12,4), dpi=300) # , link=:y -savefig("figs/metric_calibration_comparison_2022-05-01b__90_10_weight_1_year.png") - - - - -# TODO: Try heatmap of mean error over day of year - -# TODO: Reincorporate ensemble [done] - -# TODO: Try out state-based implementation [done] - -# TODO: Try out state-based implementation with borg moea, targeting each catchment state diff --git a/sandbox/figs/ensemble_pp_popsize_500.png b/sandbox/figs/ensemble_pp_popsize_500.png deleted file mode 100644 index 78fff32..0000000 Binary files a/sandbox/figs/ensemble_pp_popsize_500.png and /dev/null differ diff --git a/sandbox/figs/ensemble_pp_popsize_500_maxtime_1200_v_metric.png b/sandbox/figs/ensemble_pp_popsize_500_maxtime_1200_v_metric.png deleted file mode 100644 index b0dc8dc..0000000 Binary files a/sandbox/figs/ensemble_pp_popsize_500_maxtime_1200_v_metric.png and /dev/null differ diff --git a/sandbox/figs/ensemble_pp_popsize_500_maxtime_900.png b/sandbox/figs/ensemble_pp_popsize_500_maxtime_900.png deleted file mode 100644 index c5f43eb..0000000 Binary files a/sandbox/figs/ensemble_pp_popsize_500_maxtime_900.png and /dev/null differ diff --git a/sandbox/figs/ensemble_pp_popsize_500_maxtime_900_modded_v_metric.png b/sandbox/figs/ensemble_pp_popsize_500_maxtime_900_modded_v_metric.png deleted file mode 100644 index 05c710d..0000000 Binary files a/sandbox/figs/ensemble_pp_popsize_500_maxtime_900_modded_v_metric.png and /dev/null differ diff --git a/sandbox/figs/ensemble_pp_showcase.png b/sandbox/figs/ensemble_pp_showcase.png deleted file mode 100644 index 0d3d39b..0000000 Binary files a/sandbox/figs/ensemble_pp_showcase.png and /dev/null differ diff --git a/sandbox/figs/ensemble_pp_showcase2.png b/sandbox/figs/ensemble_pp_showcase2.png deleted file mode 100644 index da95227..0000000 Binary files a/sandbox/figs/ensemble_pp_showcase2.png and /dev/null differ diff --git a/sandbox/figs/ensemble_pp_showcase_calibrated.png b/sandbox/figs/ensemble_pp_showcase_calibrated.png deleted file mode 100644 index 58ef3ad..0000000 Binary files a/sandbox/figs/ensemble_pp_showcase_calibrated.png and /dev/null differ diff --git a/sandbox/figs/ensemble_pp_test1.png b/sandbox/figs/ensemble_pp_test1.png deleted file mode 100644 index 32869ca..0000000 Binary files a/sandbox/figs/ensemble_pp_test1.png and /dev/null differ diff --git a/sandbox/figs/ensemble_qp_ensemble_popsize_500_modded_v_metric.png b/sandbox/figs/ensemble_qp_ensemble_popsize_500_modded_v_metric.png deleted file mode 100644 index e359cd3..0000000 Binary files a/sandbox/figs/ensemble_qp_ensemble_popsize_500_modded_v_metric.png and /dev/null differ diff --git a/sandbox/figs/ensemble_qp_gr4j_popsize_500_modded_v_metric.png b/sandbox/figs/ensemble_qp_gr4j_popsize_500_modded_v_metric.png deleted file mode 100644 index 45c1450..0000000 Binary files a/sandbox/figs/ensemble_qp_gr4j_popsize_500_modded_v_metric.png and /dev/null differ diff --git a/sandbox/figs/ensemble_qp_ihacres_popsize_500_modded_v_metric.png b/sandbox/figs/ensemble_qp_ihacres_popsize_500_modded_v_metric.png deleted file mode 100644 index 44b8f27..0000000 Binary files a/sandbox/figs/ensemble_qp_ihacres_popsize_500_modded_v_metric.png and /dev/null differ diff --git a/sandbox/figs/ensemble_qp_pp_popsize_500_modded_v_metric.png b/sandbox/figs/ensemble_qp_pp_popsize_500_modded_v_metric.png deleted file mode 100644 index 0284fb8..0000000 Binary files a/sandbox/figs/ensemble_qp_pp_popsize_500_modded_v_metric.png and /dev/null differ diff --git a/sandbox/figs/metric_calibration_comparison_2022-04-30.png b/sandbox/figs/metric_calibration_comparison_2022-04-30.png deleted file mode 100644 index 6bee1de..0000000 Binary files a/sandbox/figs/metric_calibration_comparison_2022-04-30.png and /dev/null differ diff --git a/sandbox/figs/metric_calibration_comparison_2022-04-30__90_10_weight.png b/sandbox/figs/metric_calibration_comparison_2022-04-30__90_10_weight.png deleted file mode 100644 index 95322ff..0000000 Binary files a/sandbox/figs/metric_calibration_comparison_2022-04-30__90_10_weight.png and /dev/null differ diff --git a/sandbox/figs/metric_calibration_comparison_2022-04-30__90_10_weight_1_year.png b/sandbox/figs/metric_calibration_comparison_2022-04-30__90_10_weight_1_year.png deleted file mode 100644 index 77e1f58..0000000 Binary files a/sandbox/figs/metric_calibration_comparison_2022-04-30__90_10_weight_1_year.png and /dev/null differ diff --git a/sandbox/figs/metric_calibration_comparison_2022-04-30_equal_weight.png b/sandbox/figs/metric_calibration_comparison_2022-04-30_equal_weight.png deleted file mode 100644 index ce2543e..0000000 Binary files a/sandbox/figs/metric_calibration_comparison_2022-04-30_equal_weight.png and /dev/null differ diff --git a/sandbox/figs/metric_calibration_comparison_2022-05-01__90_10_weight_1_year.png b/sandbox/figs/metric_calibration_comparison_2022-05-01__90_10_weight_1_year.png deleted file mode 100644 index 347a9a2..0000000 Binary files a/sandbox/figs/metric_calibration_comparison_2022-05-01__90_10_weight_1_year.png and /dev/null differ diff --git a/sandbox/figs/metric_calibration_comparison_2022-05-01b__90_10_weight_1_year.png b/sandbox/figs/metric_calibration_comparison_2022-05-01b__90_10_weight_1_year.png deleted file mode 100644 index 6d1c4d1..0000000 Binary files a/sandbox/figs/metric_calibration_comparison_2022-05-01b__90_10_weight_1_year.png and /dev/null differ diff --git a/sandbox/figs/t_offset_streamflow_ME.png b/sandbox/figs/t_offset_streamflow_ME.png deleted file mode 100644 index 79da1b7..0000000 Binary files a/sandbox/figs/t_offset_streamflow_ME.png and /dev/null differ diff --git a/sandbox/figs/t_offset_streamflow_ME_log.png b/sandbox/figs/t_offset_streamflow_ME_log.png deleted file mode 100644 index a8ef92b..0000000 Binary files a/sandbox/figs/t_offset_streamflow_ME_log.png and /dev/null differ diff --git a/sandbox/figs/xsection_overview.png b/sandbox/figs/xsection_overview.png deleted file mode 100644 index 5d2ec71..0000000 Binary files a/sandbox/figs/xsection_overview.png and /dev/null differ diff --git a/sandbox/figs/xsection_overview_log.png b/sandbox/figs/xsection_overview_log.png deleted file mode 100644 index 6ddf761..0000000 Binary files a/sandbox/figs/xsection_overview_log.png and /dev/null differ diff --git a/sandbox/perf_improvements.jl b/sandbox/perf_improvements.jl deleted file mode 100644 index 6539a8e..0000000 --- a/sandbox/perf_improvements.jl +++ /dev/null @@ -1,73 +0,0 @@ -using Streamfall -using Statistics, CSV, DataFrames -import .Analysis: TemporalCrossSection -using Plots - -using BenchmarkTools - - -# load observations -HERE = @__DIR__ -DATA_PATH = joinpath(HERE, "../test/data/cotter/") - -date_format = "YYYY-mm-dd" -obs_data = CSV.File(joinpath(DATA_PATH, "climate/CAMELS-AUS_410730.csv"), - comment="#", - dateformat=date_format) |> DataFrame - -hist_streamflow = obs_data[:, ["Date", "410730_Q"]] -hist_streamflow = rename(hist_streamflow, "410730_Q"=>"410730") - -climate_data = obs_data[:, ["Date", "410730_P", "410730_PET"]] -climate = Climate(climate_data, "_P", "_PET") - -obs = hist_streamflow[:, "410730"] -xsect_res = Streamfall.Analysis.TemporalCrossSection(climate_data.Date, obs) -xsect = Streamfall.Analysis.cross_section(xsect_res) - - -ihacres_node = create_node(BilinearNode, "410730", 100.0) -burn_in = 365 - -@benchmark run_node!(ihacres_node, climate) - -# metric = (obs, sim) -> 1.0 - Streamfall.NmKGE(log.(obs[burn_in:end]), log.(sim[burn_in:end])) - -# res, opt = calibrate!(ihacres_node, climate, hist_streamflow; metric=metric, MaxTime=900.0) -# run_node!(ihacres_node, climate) - -# h_dates = climate_data.Date -# burn_dates = h_dates[burn_in:end] -# burn_obs = obs[burn_in:end] -# fig = Streamfall.temporal_cross_section(burn_dates, burn_obs; yscale=:log10) - -# plot(fig, size=(500, 400)) - - -# mod_xsect_res = TemporalCrossSection(burn_dates, burn_obs, ihacres_node.outflow[burn_in:end]) -# offset_vals = mod_xsect_res.offsets - -# fig2 = Streamfall.temporal_cross_section(burn_dates, burn_obs, ihacres_node.outflow[burn_in:end]; yscale=:log10, title="NmKGE") - -# function postprocess(x, ω) -# idx = (x .+ ω) .>= 0.0 -# x[idx] = x[idx] .+ ω[idx] - -# idx2 = (x .+ ω) .< 0 -# x[idx2] .= 0.0 -# return x -# end - -# fig3 = Streamfall.temporal_cross_section(burn_dates, burn_obs, postprocess(ihacres_node.outflow[burn_in:end], offset_vals); yscale=:log10, title="Offset (>= 0.0)") - -# plot(fig2, fig3, size=(1200, 500), link=:y, layout=(1,2), dpi=300) -# savefig("xsection_overview_log.png") - -# off_vals = postprocess(ihacres_node.outflow[burn_in:end], offset_vals) -# orig_score = round(Streamfall.NmKGE(log.(burn_obs), log.(ihacres_node.outflow[burn_in:end])), digits=2) -# off_score = round(Streamfall.NmKGE(log.(burn_obs), log.(off_vals)), digits=2) - -# plot(burn_dates, Streamfall.ME.(burn_obs, ihacres_node.outflow[burn_in:end]), alpha=0.4, label="IHACRES (log NmKGE: $(orig_score))") -# plot!(burn_dates, Streamfall.ME.(burn_obs, off_vals), alpha=0.4, label="Post-Processed (log NmKGE: $(off_score))", dpi=300) -# savefig("t_offset_streamflow_ME_log.png") - diff --git a/sandbox/state_dev.jl b/sandbox/state_dev.jl deleted file mode 100644 index f4a37ec..0000000 --- a/sandbox/state_dev.jl +++ /dev/null @@ -1,67 +0,0 @@ -using Streamfall -using Statistics, CSV, DataFrames -using Distributed, BlackBoxOptim -using Base.Iterators -using Plots -import Streamfall.Analysis: TemporalCrossSection - - -include("IHACRESCMDStateNode.jl") - - -function v_metric(obs, sim; sp=365, metric=Streamfall.NmKGE) - res = Streamfall.naive_split_metric(obs, sim, sp, metric) - m = mean(res) - cv = std(res) / m - - V = (m*0.9) + ((ℯ^-cv)*0.1) # (ℯ .- 1.5) - - return V -end - - -function v_mod_metric(obs, sim; sp=365, metric=Streamfall.NmKGE) - res = Streamfall.naive_split_metric(obs, sim, sp, metric) - m = mean(res) - sd = std(res) - - V = (m*0.9) + ((ℯ^-sd)*0.1) - - return V -end - - -# load observations -HERE = @__DIR__ -DATA_PATH = joinpath(HERE, "../test/data/cotter/") - -date_format = "YYYY-mm-dd" -obs_data = CSV.File(joinpath(DATA_PATH, "climate/CAMELS-AUS_410730.csv"), - comment="#", - dateformat=date_format) |> DataFrame - -hist_streamflow = obs_data[:, ["Date", "410730_Q"]] -hist_streamflow = rename(hist_streamflow, "410730_Q"=>"410730") - -climate_data = obs_data[:, ["Date", "410730_P", "410730_PET"]] -climate = Climate(climate_data, "_P", "_PET") - -state_node = create_node(CMDStateNode, "410730", 129.2) - - -c_date = climate_data.Date -obs = hist_streamflow[:, "410730"] -figs = [] -outflows = [] -display_timeframe = 13330:(13330+(365*1)) - - -burn_in = 1826 - -split_mod_v = (obs, sim) -> 1.0 - v_mod_metric(obs[burn_in:end], sim[burn_in:end]; sp=14, metric=Streamfall.NnpKGE) - - -# calibrate!(state_node, climate, hist_streamflow; metric=split_mod_v, MaxTime=10.0, TraceMode=:silent) - -res, opt = calibrate!(state_node, climate, hist_streamflow; metric=split_mod_v, MaxTime=10.0, TraceMode=:silent); - diff --git a/sandbox/temporal.jl b/sandbox/temporal.jl deleted file mode 100644 index d5973a3..0000000 --- a/sandbox/temporal.jl +++ /dev/null @@ -1,68 +0,0 @@ -using Streamfall -using Statistics, CSV, DataFrames -import .Analysis: TemporalCrossSection -using Plots - - -# load observations -HERE = @__DIR__ -DATA_PATH = joinpath(HERE, "../test/data/cotter/") - -date_format = "YYYY-mm-dd" -obs_data = CSV.File(joinpath(DATA_PATH, "climate/CAMELS-AUS_410730.csv"), - comment="#", - dateformat=date_format) |> DataFrame - -hist_streamflow = obs_data[:, ["Date", "410730_Q"]] -hist_streamflow = rename(hist_streamflow, "410730_Q"=>"410730") - -climate_data = obs_data[:, ["Date", "410730_P", "410730_PET"]] -climate = Climate(climate_data, "_P", "_PET") - -obs = hist_streamflow[:, "410730"] -xsect_res = Streamfall.Analysis.TemporalCrossSection(climate_data.Date, obs) -xsect = Streamfall.Analysis.cross_section(xsect_res) - - -ihacres_node = create_node(BilinearNode, "410730", 100.0) -burn_in = 365 -metric = (obs, sim) -> 1.0 - Streamfall.NmKGE(log.(obs[burn_in:end]), log.(sim[burn_in:end])) - -res, opt = calibrate!(ihacres_node, climate, hist_streamflow; metric=metric, MaxTime=900.0) -run_node!(ihacres_node, climate) - -h_dates = climate_data.Date -burn_dates = h_dates[burn_in:end] -burn_obs = obs[burn_in:end] -fig = Streamfall.temporal_cross_section(burn_dates, burn_obs; yscale=:log10) - -plot(fig, size=(500, 400)) - - -mod_xsect_res = TemporalCrossSection(burn_dates, burn_obs, ihacres_node.outflow[burn_in:end]) -offset_vals = mod_xsect_res.offsets - -fig2 = Streamfall.temporal_cross_section(burn_dates, burn_obs, ihacres_node.outflow[burn_in:end]; yscale=:log10, title="NmKGE") - -function postprocess(x, ω) - idx = (x .+ ω) .>= 0.0 - x[idx] = x[idx] .+ ω[idx] - - idx2 = (x .+ ω) .< 0 - x[idx2] .= 0.0 - return x -end - -fig3 = Streamfall.temporal_cross_section(burn_dates, burn_obs, postprocess(ihacres_node.outflow[burn_in:end], offset_vals); yscale=:log10, title="Offset (>= 0.0)") - -plot(fig2, fig3, size=(1200, 500), link=:y, layout=(1,2), dpi=300) -savefig("xsection_overview_log.png") - -off_vals = postprocess(ihacres_node.outflow[burn_in:end], offset_vals) -orig_score = round(Streamfall.NmKGE(log.(burn_obs), log.(ihacres_node.outflow[burn_in:end])), digits=2) -off_score = round(Streamfall.NmKGE(log.(burn_obs), log.(off_vals)), digits=2) - -plot(burn_dates, Streamfall.ME.(burn_obs, ihacres_node.outflow[burn_in:end]), alpha=0.4, label="IHACRES (log NmKGE: $(orig_score))") -plot!(burn_dates, Streamfall.ME.(burn_obs, off_vals), alpha=0.4, label="Post-Processed (log NmKGE: $(off_score))", dpi=300) -savefig("t_offset_streamflow_ME_log.png") -