text:
# External packages
using Revise, DelimitedFiles, SymPy

# Add Exasim to Julia search path
cdir = pwd(); ii = findlast("Exasim", cdir);
include(cdir[1:ii[end]] * "/Installation/setpath.jl");

# Exasim packages
using Preprocessing, Mesh, Gencode, Postprocessing

# create pde structure and mesh structure
pde, mesh = Preprocessing.initializeexasim();

# Define PDE model: governing equations, initial solutions, and boundary conditions
pde.model = "ModelC";            # ModelC, ModelD, ModelW
include("pdemodel.jl");          # include the PDE model file

# Set discretization parameters, physical parameters, and solver parameters
pde.porder = 4;          # polynomial degree
pde.torder = 3;          # time-stepping order of accuracy
pde.nstage = 3;          # time-stepping number of stages
pde.dt = 0.05*ones(200);   # time step sizes
pde.soltime = collect(1:length(pde.dt)); # steps at which solution are collected
pde.visdt = 0.05; # visualization timestep size

gam = 1.4;                      # specific heat ratio
M_ref = sqrt(1/gam);            # Mach number
pde.physicsparam = [gam M_ref];
pde.tau = [1+1/M_ref];          # DG stabilization parameter

# Choose computing platform and set number of processors
#pde.platform = "gpu";           # choose this option if NVIDIA GPUs are available
pde.mpiprocs = 1;                # number of MPI processors

# create a linear mesh for a square domain
mesh.p, mesh.t = Mesh.SquareMesh(10,10,1); # a mesh of 8 by 8 quadrilaterals
mesh.p = 10.0*mesh.p .- 5.0;
# expressions for disjoint boundaries
mesh.boundaryexpr = [p -> (p[2,:] .< -5+1e-3), p -> (p[1,:] .> 5-1e-3), p -> (p[2,:] .> 5-1e-3), p -> (p[1,:] .< -5+1e-3)];
mesh.boundarycondition = [1 1 1 1]; # Set boundary condition for each disjoint boundary
mesh.periodicexpr = [2 p->p[2,:] 4 p->p[2,:]; 1 p->p[1,:] 3 p->p[1,:]];

# call exasim to generate and run C++ code to solve the PDE model
sol, pde, mesh,~,~,~,~  = Postprocessing.exasim(pde,mesh);

# visualize the numerical solution of the PDE model using Paraview
pde.visscalars = ["density", 1, "energy", 4];  # list of scalar fields for visualization
pde.visvectors = ["momentum", [2, 3]]; # list of vector fields for visualization
Postprocessing.vis(sol,pde,mesh); # visualize the numerical solution
print("Done!");


text:
function make_initial_load(cur_package::Module, all_files::AbstractArray; verbose::Bool=false, is_test::Bool=false)

  cur_cabinet = FileCabinet(cur_package)

  verbose && println("\ninitial load:\n")

  for cur_file in all_files

    verbose && print(cur_file * " ")

    cur_shards = parse_file(cur_file)

    cur_info = FileInfo(cur_file, cur_cabinet)

    bad_index = 0

    for (cur_index, cur_shard) in enumerate(cur_shards)

      did_load = attempt_shard_load!(
        cur_package, cur_info, cur_shard, is_test
      )

      if !did_load
        bad_index = cur_index
        break
      end

      push!(cur_info.loaded_shards, cur_shard)

    end

    verbose && println( iszero(bad_index) ? "✓" : "X" )

    if iszero(bad_index)
      push!(cur_cabinet.load_order, cur_file)
      continue
    end

    cur_info.unloaded_shards = cur_shards[bad_index:end]

    for (cur_index, cur_shard) in enumerate(cur_info.loaded_shards)

      cur_info.loaded_shards[cur_index] =
        clean_shard(cur_package, cur_shard)

    end

    push!(cur_cabinet.file_infos, cur_info)

  end

  cur_cabinet

end

ext: jl


text:
module RcSetup
using Parameters
import DataStructures: OrderedDict
export Setup, Port, all_ports, OrderedDict

module Port
export PressPorts, all_ports
@enum PressPorts in1 in2 in3 in4 out1 out2 out3 out4 ice1 ice2 ice3 ice4
const all_ports = [in1, in2, in3, in4, ice1, ice2, ice3, ice4, out1, out2, out3, out4]
end
"Distance between fixed installed ports"
const def_portdists = [0.88, 0.44, 0.44, 1.68, 0.34, 0.34, 0.68]

# TODO: figure something out for pressure ports in ice:  The change depending on time

const P2C = OrderedDict{Port.PressPorts,Vector{Int}}
@with_kw struct Setup
    drawing_pdf::String # path to pdf
    port2channels::P2C # mapping physical port to channels
                                                            # keep the ports in order
    @assert length(used_channels(port2channels))==length(unique(used_channels(port2channels)))
    portdists::Vector{Float64}=def_portdists # needs to be in order of keys(port2channels)
end
"Return the used measurement channels"
used_channels(di::P2C) = sort(vcat(values(di)...))
used_channels(setup::Setup) = used_channels(setup.port2channels)
"Return the used measurement ports"
used_ports(di::P2C) = collect(keys(di))
used_ports(setup::Setup) = used_ports(setup.port2channels)

end


text:
function value(obj::RealFieldBivariateFunction, arg0::RealFieldElement, arg1::RealFieldElement)
    return jcall(obj, "value", RealFieldElement, (RealFieldElement, RealFieldElement), arg0, arg1)
end


text:
"""
	kml2gmt(cmd0::String="", arg1=nothing, kwargs...)

kml2gmt - Extract GMT table data from Google Earth KML files

Full option list at [`kml2gmt`]($(GMTdoc)kml2gmt.html)

Parameters
----------

- **F** | **select** :: [Type => Str]        ``Arg = s|l|p``

    Specify a particular feature type to output. Choose from points (s), line (l), or polygon (p).
    By default we output all geometries.
    ($(GMTdoc)kml2gmt.html#f)
- **Z** | **altitudes** :: [Type => Bool]

    Output the altitude coordinates as GMT z coordinates [Default will output just longitude and latitude].
    ($(GMTdoc)kml2gmt.html#z)
- $(GMT.opt_V)
- $(GMT.opt_bo)
- $(GMT.opt_do)
- $(GMT.opt_swap_xy)
"""
function kml2gmt(cmd0::String="", arg1=nothing; kwargs...)

	length(kwargs) == 0 && occursin(" -", cmd0) && return monolitic("kml2gmt", cmd0, arg1)

	d = KW(kwargs)
	help_show_options(d)		# Check if user wants ONLY the HELP mode
	cmd, = parse_common_opts(d, "", [:V_params :bo :do :yx])
	cmd  = parse_these_opts(cmd, d, [[:F :select], [:Z :altitudes], [:E :extended]])
	common_grd(d, cmd0, cmd, "kml2gmt ", arg1)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
kml2gmt(arg1, cmd0::String=""; kw...) = kml2gmt(cmd0, arg1; kw...)


text:
immutable VCJL <: Version end

Base.show(io::IO, v::VCJL) = print(io, "CJL")

"""
    vCJL()

Identifier for the DICE-CJL version of the model.

# Examples
```jldoctest
julia> vCJL()
CJL
```
"""
function vCJL()
    VCJL()
end

export vCJL

@extend immutable OptionsVCJL <: Options
    δσ₂::Float64 #Quadratic term in decarbonization
    σ₀::Float64 #CO₂-equivalent emissions-GNP ratio 2005
    backrat::Float64 #atio initial to final backstop cost
    partfract1::Float64 #Fraction of emissions under control regime 2005
    partfract2::Float64 #Fraction of emissions under control regime 2015
    partfract21::Float64 #Fraction of emissions under control regime 2205
    dpartfract::Float64 #Decline rate of participation
end

function options(version::VCJL;
    N::Int = 600, #Number of years to calculate (from 2005 onwards)
    tstep::Int = 1, #Years per Period (UNUSED IN CJL)
    α::Float64 = 2.0, #Elasticity of marginal utility of consumption
    ρ::Float64 = 0.015, #Initial rate of social time preference per year
    γₑ::Float64 = 0.3, #Capital elasticity in production function
    pop₀::Int = 6514, #2005 world population (millions)
    popadj::Float64 = 0.035, #->gpop₀ Growth rate of population per year
    popasym::Int = 8600, #Asymptotic population (millions)
    δk::Float64 = 0.1, #Depreciation rate on capital (per year)
    q₀::Float64 = 61.1, #Initial world gross output (trill 2005 USD)
    k₀::Float64 = 137.0, #Initial capital value (trill 2005 USD)
    a₀::Float64 = 0.02722, #Initial level of total factor productivity
    ga₀::Float64 = 0.092, #Initial growth rate for TFP per year
    δₐ::Float64 = 0.001, #Decline rate of TFP per year
    gσ₁::Float64 = -0.00730, #Initial growth of sigma (continuous per year)
    δσ::Float64 = 0.003, #Decline rate of decarbonization per period
    eland₀::Float64 = 1.1, #Carbon emissions from land 2005 (GtCO2 per year)
    deland::Float64 = 0.01, #Decline rate of land emissions (per period)
    mat₀::Float64 = 808.9, #Initial Concentration in atmosphere 2005 (GtC)
    mu₀::Float64 = 1255.0, #Initial Concentration in upper strata 2005 (GtC)
    ml₀::Float64 = 18365.0, #Initial Concentration in lower strata 2005 (GtC)
    mateq::Float64 = 587.473, #Equilibrium concentration atmosphere  (GtC)
    mueq::Float64 = 1143.894, #Equilibrium concentration in upper strata (GtC)
    mleq::Float64 = 18340.0, #Equilibrium concentration in lower strata (GtC)
    ϕ₁₂::Float64 = 0.0189288, #Carbon cycle transition matrix coefficient
    ϕ₂₃::Float64 = 0.005, #Carbon cycle transition matrix coefficient
    t2xco2::Float64 = 3.0, #Equilibrium temp impact (oC per doubling CO2)
    fₑₓ0::Float64 = -0.06, #2000 forcings of non-CO2 GHG (Wm-2)
    fₑₓ1::Float64 = 0.3, #2000 forcings of non-CO2 GHG (Wm-2)
    tocean₀::Float64 = 0.0068, #2000 lower stratum temp change (C from 1900)
    tatm₀::Float64 = 0.7307, #2000 atmospheric temp change (C from 1900)
    ξ₁::Float64 = 0.0220, #Climate equation coefficient for upper level
    ξ₃::Float64 = 0.3, #Transfer coefficient upper to lower stratum
    ξ₄::Float64 = 0.005, #Transfer coefficient for lower level
    η::Float64 = 3.8, #Forcings of equilibrium CO2 doubling (Wm-2)
    ψ₁::Float64 = 0.0, #Damage intercept
    ψ₂::Float64 = 0.0028388, #Damage quadratic term
    ψ₃::Float64 = 2.0, #Damage exponent
    θ₂::Float64 = 2.8, #Exponent of control cost function
    pback::Float64 = 1.17, #Cost of backstop 2005$ per tCO2 2005
    gback::Float64 = 0.005, #Initial cost decline backstop cost per year
    limμ::Float64 = 1.0, #Upper limit on control rate after 2150
    fosslim::Float64 = 6000.0, #Maximum cumulative extraction fossil fuels (GtC)
    scale1::Float64 = 194., #Multiplicative scaling coefficient
    scale2::Float64 = 381800., #Additive scaling coefficient
    δσ₂::Float64 = 0.0, #Quadratic term in decarbonization
    σ₀::Float64 = 0.13418, #CO₂-equivalent emissions-GNP ratio 2005
    backrat::Float64 = 2.0, #atio initial to final backstop cost
    partfract1::Float64 = 1.0, #Fraction of emissions under control regime 2005
    partfract2::Float64 = 1.0, #Fraction of emissions under control regime 2015
    partfract21::Float64 = 1.0, #Fraction of emissions under control regime 2205
    dpartfract::Float64 = 0.0) #Decline rate of participation

    OptionsVCJL(N,tstep,α,ρ,γₑ,pop₀,popadj,popasym,δk,q₀,k₀,a₀,ga₀,δₐ,gσ₁,δσ,eland₀,deland,mat₀,mu₀,ml₀,mateq,mueq,mleq,ϕ₁₂,ϕ₂₃,t2xco2,fₑₓ0,fₑₓ1,tocean₀,tatm₀,ξ₁,ξ₃,ξ₄,η,ψ₁,ψ₂,ψ₃,θ₂,pback,gback,limμ,fosslim,scale1,scale2,δσ₂,σ₀,backrat,partfract1,partfract2,partfract21,dpartfract)
end

function Base.show(io::IO, ::MIME"text/plain", opt::OptionsVCJL)
    println(io, "Options for DICE-CJL");
    println(io, "N: $(opt.N)");
    println(io, "Preferences");
    println(io, "α: $(opt.α), ρ: $(opt.ρ)");
    println(io, "Population and Technology");
    println(io, "γₑ: $(opt.γₑ), pop₀: $(opt.pop₀), popadj: $(opt.popadj), popasym: $(opt.popasym), δk: $(opt.δk)");
    println(io, "q₀: $(opt.q₀), k₀: $(opt.k₀), a₀: $(opt.a₀), ga₀: $(opt.ga₀), δₐ: $(opt.δₐ)");
    println(io, "Emissions Parameters");
    println(io, "σ₀: $(opt.σ₀), gσ₁: $(opt.gσ₁), δσ: $(opt.δσ), δσ₂: $(opt.δσ₂), eland₀: $(opt.eland₀)");
    println(io, "Carbon Cycle");
    println(io, "mat₀: $(opt.mat₀), mu₀: $(opt.mu₀), ml₀: $(opt.ml₀)");
    println(io, "mateq: $(opt.mateq), mueq: $(opt.mueq), mleq: $(opt.mleq)");
    println(io, "Flow Parameters");
    println(io, "ϕ₁₂: $(opt.ϕ₁₂), ϕ₂₃: $(opt.ϕ₂₃)");
    println(io, "Climate Model Parameters");
    println(io, "t2xco2: $(opt.t2xco2), fₑₓ0: $(opt.fₑₓ0), fₑₓ1: $(opt.fₑₓ1)");
    println(io, "tocean₀: $(opt.tocean₀), tatm₀: $(opt.tatm₀), ξ₁: $(opt.ξ₁)");
    println(io, "ξ₃: $(opt.ξ₃), ξ₄: $(opt.ξ₄), η: $(opt.η)");
    println(io, "Climate Damage Parameters");
    println(io, "ψ₁: $(opt.ψ₁), ψ₂: $(opt.ψ₂), ψ₃: $(opt.ψ₃)");
    println(io, "Abatement Cost");
    println(io, "θ₂: $(opt.θ₂), pback: $(opt.pback), gback: $(opt.gback), backrat: $(opt.backrat), limμ: $(opt.limμ)");
    println(io, "Fossil Fuel Availability");
    println(io, "fosslim: $(opt.fosslim)");
    println(io, "Scaling Parameters");
    println(io, "scale1: $(opt.scale1), scale2: $(opt.scale2)");
    println(io, "Emissions control");
    print(io, "partfract1: $(opt.partfract1), partfract2: $(opt.partfract2), partfract21: $(opt.partfract21), dpartfract: $(opt.dpartfract)");
end

@extend immutable ParametersVCJL <: Parameters
    rr::Array{Float64,1} # Average utility social discount rate
    partfract::Array{Float64,1} # Fraction of emissions in control regime
end

function generate_parameters(c::OptionsVCJL)
    ϕ₁₁::Float64 = 1 - c.ϕ₁₂; # Carbon cycle transition matrix coefficient
    ϕ₂₁::Float64 = c.mateq*c.ϕ₁₂/c.mueq; # Carbon cycle transition matrix coefficient
    ϕ₂₂::Float64 = 1 - ϕ₂₁ - c.ϕ₂₃; # Carbon cycle transition matrix coefficient
    ϕ₃₂::Float64 = c.mueq*c.ϕ₂₃/c.mleq; # Carbon cycle transition matrix coefficient
    ϕ₃₃::Float64 = 1 - ϕ₃₂; # Carbon cycle transition matrix coefficient
    σ₀::Float64 = c.σ₀; # CO2-equivalent emissions-GNP ratio 2005
    λ::Float64 = c.η/c.t2xco2; # Climate model parameter

    # Growth factor population
    #NOTE: This is only a helper function for L, so we don't include it in the Narrative anywhere.
    gfacpop = Array{Float64}(c.N);
    # Level of population and labor
    L = Array{Float64}(c.N);
    # Growth rate of productivity from 0 to N
    gₐ = Array{Float64}(c.N);
    # Change in sigma (cumulative improvement of energy efficiency)
    gσ = Array{Float64}(c.N);
    # Emissions from deforestation
    Etree = Array{Float64}(c.N);
    # Average utility social discount rate
    rr = Array{Float64}(c.N);

    for i in 1:c.N
        gfacpop[i] = (exp.(c.popadj*(i-1))-1)/exp.(c.popadj*(i-1));
        L[i] = c.pop₀*(1-gfacpop[i])+gfacpop[i]*c.popasym;
        gₐ[i] = c.ga₀*exp.(-c.δₐ*(i-1));
        gσ[i] = c.gσ₁*exp.(-c.δσ*(i-1))-c.δσ₂*((i-1)^2);
        Etree[i] = c.eland₀*(1-c.deland)^(i-1);
        rr[i] = 1/((1+c.ρ)^(i-1));
    end

    # Initial conditions and offset required
    # Level of total factor productivity
    A = Array{Float64}(c.N);
    A[1] = c.a₀;
    # CO2-equivalent-emissions output ratio
    σ = Array{Float64}(c.N);
    σ[1] = σ₀;

    for i in 1:c.N-1
        A[i+1] = A[i]/(1-gₐ[i]);
        σ[i+1] = σ[i]/(1-gσ[i+1]);
    end

    # Adjusted cost for backstop
    θ₁ = Array{Float64}(c.N);
    # Exogenous forcing for other greenhouse gases
    fₑₓ = Array{Float64}(c.N);
    # Fraction of emissions in control regime
    partfract = Array{Float64}(c.N);

    for i in 1:c.N
        θ₁[i] = (c.pback*σ[i]/c.θ₂)*((c.backrat-1+exp.(-c.gback*(i-1)))/c.backrat);
        fₑₓ[i] = if i < 101
                        c.fₑₓ0+0.01*(c.fₑₓ1-c.fₑₓ0)*(i-1)
                    else
                        0.36
                    end;
        partfract[i] = if i <= 259
                            c.partfract21+(c.partfract2-c.partfract21)*exp.(-c.dpartfract*(i-2))
                       else
                            c.partfract21
                       end;
    end
    partfract[1] = 0.25372; #NOTE: This is overriden in the setup, so partfract1 is never really used...
    ParametersVCJL(ϕ₁₁,ϕ₂₁,ϕ₂₂,ϕ₃₂,ϕ₃₃,σ₀,λ,gₐ,Etree,L,A,gσ,σ,θ₁,fₑₓ,rr,partfract)
end

function Base.show(io::IO, ::MIME"text/plain", opt::ParametersVCJL)
    println(io, "Calculated Parameters for DICE-CJL");
    println(io, "Carbon cycle transition matrix coefficients");
    println(io, "ϕ₁₁: $(opt.ϕ₁₁), ϕ₂₁: $(opt.ϕ₂₁), ϕ₂₂: $(opt.ϕ₂₂), ϕ₃₂: $(opt.ϕ₃₂), ϕ₃₃: $(opt.ϕ₃₃)");
    println(io, "2015 Carbon intensity: $(opt.σ₀)");
    println(io, "Climate model parameter: $(opt.λ)");
    println(io, "Growth rate of productivity: $(opt.gₐ)");
    println(io, "Emissions from deforestation: $(opt.Etree)");
    println(io, "Avg utility social discout rate: $(opt.rr)");
    println(io, "Population and labour: $(opt.L)");
    println(io, "Total factor productivity: $(opt.A)");
    println(io, "Cumulative improvement of energy efficiency: $(opt.gσ)");
    println(io, "CO₂-equivalent-emissions output ratio: $(opt.σ)");
    println(io, "Adjusted cost for backstop: $(opt.θ₁)");
    println(io, "Exogenious forcing: $(opt.fₑₓ)");
    print(io, "Fraction of emissions in control regime: $(opt.partfract)");
end

@extend struct VariablesVCJL <: Variables
    MₐₜAV::Array{JuMP.Variable,1} # Average concentrations
    PCY::Array{JuMP.Variable,1} # Per capita income thousands US dollars
end

function model_vars(version::VCJL, model::JuMP.Model, N::Int64, cca_ubound::Float64, μ_ubound::Float64)
    @variable(model, 0.0 <= μ[i=1:N] <= μ_ubound); # Emission control rate GHGs
    @variable(model, FORC[1:N]); # Increase in radiative forcing (watts per m2 from 1900)
    @variable(model, 0.0 <= Tₐₜ[1:N] <= 20.0); # Increase temperature of atmosphere (degrees C from 1900)
    @variable(model, -1.0 <= Tₗₒ[1:N] <= 20.0); # Increase temperatureof lower oceans (degrees C from 1900)
    @variable(model, Mₐₜ[1:N] >= 10.0); # Carbon concentration increase in atmosphere (GtC from 1750)
    @variable(model, Mᵤₚ[1:N] >= 100.0); # Carbon concentration increase in shallow oceans (GtC from 1750)
    @variable(model, Mₗₒ[1:N] >= 1000.0); # Carbon concentration increase in lower oceans (GtC from 1750)
    @variable(model, E[1:N] >= 0.0); # Total CO2 emissions (GtCO2 per year)
    @variable(model, C[1:N] >= 20.0); # Consumption (trillions 2005 US dollars per year)
    @variable(model, K[1:N] >= 100.0); # Capital stock (trillions 2005 US dollars)
    @variable(model, CPC[1:N] >= 0.01); #  Per capita consumption (thousands 2005 USD per year)
    @variable(model, I[1:N] >= 0.0); # Investment (trillions 2005 USD per year)
    @variable(model, S[1:N]); # Gross savings rate as fraction of gross world product
    @variable(model, RI[1:N]); # Real interest rate (per annum)
    @variable(model, Y[1:N] >= 0.0); # Gross world product net of abatement and damages (trillions 2005 USD per year)
    @variable(model, YGROSS[1:N] >= 0.0); # Gross world product GROSS of abatement and damages (trillions 2005 USD per year)
    @variable(model, YNET[1:N]); # Output net of damages equation (trillions 2005 USD per year)
    @variable(model, DAMAGES[1:N]); # Damages (trillions 2005 USD per year)
    @variable(model, MCABATE[1:N]); #->ABATECOST Marginal cost of abatement (2005$ per ton CO2)
    @variable(model, 0.0 <= CCA[1:N] <= cca_ubound); # Cumulative industrial carbon emissions (GTC)
    @variable(model, PERIODU[1:N]); # One period utility function
    @variable(model, UTILITY); # Welfare function
    @variable(model, MₐₜAV[1:N] >= 0.0); # Average concentrations
    @variable(model, PCY[1:N]); # Per capita income thousands US dollars
    VariablesVCJL(μ,FORC,Tₐₜ,Tₗₒ,Mₐₜ,Mᵤₚ,Mₗₒ,E,C,K,CPC,I,S,RI,Y,YGROSS,YNET,DAMAGES,MCABATE,CCA,PERIODU,UTILITY,MₐₜAV,PCY)
end

@extend struct EquationsVCJL <: Equations
    kk::Array{JuMP.ConstraintRef,1} # Capital balance equation
end

function model_eqs(model::JuMP.Model, config::OptionsVCJL, params::ParametersVCJL, vars::VariablesVCJL)
    N = config.N;
    # Equations #
    # Emissions Equation
    eeq = @NLconstraint(model, [i=1:N], vars.E[i] == params.σ[i] * (1-vars.μ[i])*params.A[i]*params.L[i]^(1-config.γₑ)*vars.K[i]^config.γₑ + params.Etree[i]);
    # Radiative forcing equation
    @NLconstraint(model, [i=1:N], vars.FORC[i] == config.η * (log((vars.MₐₜAV[i]+.000001)/596.4)/log(2)) + params.fₑₓ[i]);
    # Average concentrations equation
    @constraint(model, [i=1:N], vars.MₐₜAV[i] == vars.Mₐₜ[i]);
    # Output gross equation
    @NLconstraint(model, [i=1:N], vars.YGROSS[i] == params.A[i]*params.L[i]^(1-config.γₑ)*vars.K[i]^config.γₑ);
    # Damage equation
    @NLconstraint(model, [i=1:N], vars.DAMAGES[i] == vars.YGROSS[i]-vars.YGROSS[i]/(1+config.ψ₁*vars.Tₐₜ[i]+config.ψ₂*vars.Tₐₜ[i]^config.ψ₃));
    # Output net of damages equation
    @NLconstraint(model, [i=1:N], vars.YNET[i] == vars.YGROSS[i]/(1+config.ψ₁*vars.Tₐₜ[i]+config.ψ₂*vars.Tₐₜ[i]^config.ψ₃));
    # ->ABATECOST Abatement cost
    @NLconstraint(model, [i=1:N], vars.MCABATE[i] == params.partfract[i]^(1-config.θ₂)*vars.YGROSS[i]*(params.θ₁[i]*(vars.μ[i]^config.θ₂)));
    # Output net equation
    @NLconstraint(model, [i=1:N], vars.Y[i] == vars.YGROSS[i]*(1-(params.partfract[i]^(1-config.θ₂))*params.θ₁[i]*(vars.μ[i]^config.θ₂))/(1+config.ψ₁*vars.Tₐₜ[i]+config.ψ₂*vars.Tₐₜ[i]^config.ψ₃));
    # Savings rate equation
    @NLconstraint(model, [i=1:N], vars.S[i] == vars.I[i]/(0.001 + vars.Y[i]));
    # Interest rate equation
    @NLconstraint(model, [i=1:N], vars.RI[i] == config.γₑ*vars.Y[i]/vars.K[i]-config.δk);
    # Consumption equation
    @constraint(model, [i=1:N], vars.C[i] == vars.Y[i] - vars.I[i]);
    # Per capita consumption definition
    @constraint(model, [i=1:N], vars.CPC[i] == vars.C[i] * 1000.0 / params.L[i]);
    # Per capita income definition
    @constraint(model, [i=1:N], vars.PCY[i] == vars.Y[i] * 1000.0 / params.L[i]);
    # Instantaneous utility function equation
    @NLconstraint(model, [i=1:N], vars.PERIODU[i] == ((vars.C[i]/params.L[i])^(1-config.α)-1)/(1-config.α));

    # Equations (offset) #
    # Cumulative carbon emissions
    @constraint(model, [i=1:N-1], vars.CCA[i+1] == vars.CCA[i] + vars.E[i]);
    # Capital balance equation
    kk = @constraint(model, [i=1:N-1], vars.K[i+1] <= (1-config.δk) * vars.K[i] + vars.I[i]);
    # Atmospheric concentration equation
    @constraint(model, [i=1:N-1], vars.Mₐₜ[i+1] == vars.Mₐₜ[i]*params.ϕ₁₁ + vars.Mᵤₚ[i]*params.ϕ₂₁ + vars.E[i]);
    # Lower ocean concentration
    @constraint(model, [i=1:N-1], vars.Mₗₒ[i+1] == vars.Mₗₒ[i]*params.ϕ₃₃ + vars.Mᵤₚ[i]*config.ϕ₂₃);
    # Shallow ocean concentration
    @constraint(model, [i=1:N-1], vars.Mᵤₚ[i+1] == vars.Mₐₜ[i]*config.ϕ₁₂ + vars.Mᵤₚ[i]*params.ϕ₂₂ + vars.Mₗₒ[i]*params.ϕ₃₂);
    # Temperature-climate equation for atmosphere
    @constraint(model, [i=1:N-1], vars.Tₐₜ[i+1] == vars.Tₐₜ[i] + config.ξ₁*(vars.FORC[i]-params.λ*vars.Tₐₜ[i]-config.ξ₃*(vars.Tₐₜ[i]-vars.Tₗₒ[i])));
    # Temperature-climate equation for lower oceans
    @constraint(model, [i=1:N-1], vars.Tₗₒ[i+1] == vars.Tₗₒ[i] + config.ξ₄*(vars.Tₐₜ[i]-vars.Tₗₒ[i]));

    # Terminal condition for capital
    @constraint(model, 0.02*vars.K[end] <= vars.I[end]);
    # Fix savings assumption for standardization if needed
    #@constraint(model, vars.S[i=1:N] .== 0.22); #NOTE: This kills all savings??
    # Initial conditions
    @constraint(model, vars.CCA[1] == 0.0);
    @constraint(model, vars.K[1] == config.k₀);
    @constraint(model, vars.Mₐₜ[1] == config.mat₀);
    @constraint(model, vars.Mᵤₚ[1] == config.mu₀);
    @constraint(model, vars.Mₗₒ[1] == config.ml₀);
    @constraint(model, vars.Tₐₜ[1] == config.tatm₀);
    @constraint(model, vars.Tₗₒ[1] == config.tocean₀);
    #First period predetermined by Kyoto Protocol
    @constraint(model, vars.μ[1] == 0.005);

    @constraint(model, vars.UTILITY == sum(params.rr[i]*params.L[i]*vars.PERIODU[i]/config.scale1 for i=1:N) + config.scale2);

    # Objective function
    @objective(model, Max, vars.UTILITY);

    EquationsVCJL(eeq,kk)
end

# Just put the default into OptimalPrice for the moment.
function assign_scenario(s::OptimalPriceScenario, model::JuMP.Model, config::OptionsVCJL, params::ParametersVCJL, vars::VariablesVCJL)
end

function assign_scenario(s::Scenario, model::JuMP.Model, config::OptionsVCJL, params::ParametersVCJL, vars::VariablesVCJL)
    error("$(s) is not a valid scenario for DICE-CJL");
end

@extend struct ResultsVCJL <: Results
    MCABATE::Array{Float64,1} # Marginal cost of abatement (2005$ per ton CO2)
    MₐₜAV::Array{Float64,1} # Average concentrations
    PCY::Array{Float64,1} # Per capita income thousands US dollars
    mcemis::Array{Float64,1} #Unsure: TODO
end

function model_results(model::JuMP.Model, config::OptionsVCJL, params::ParametersVCJL, vars::VariablesVCJL, eqs::EquationsVCJL)
    years = 2005+(1:config.N); # TODO: check this
    Mₐₜ = getvalue(vars.Mₐₜ);
    MₐₜAV = getvalue(vars.MₐₜAV);
    Mₐₜppm = getvalue(vars.Mₐₜ)/2.13;
    Mᵤₚ = getvalue(vars.Mᵤₚ);
    Mₗₒ = getvalue(vars.Mₗₒ);
    CCA = getvalue(vars.CCA);
    CCAratio = getvalue(vars.CCA)/config.fosslim;
    Tₐₜ = getvalue(vars.Tₐₜ);
    FORC = getvalue(vars.FORC);
    Tₗₒ = getvalue(vars.Tₗₒ);
    YGROSS = getvalue(vars.YGROSS);
    DAMAGES = getvalue(vars.DAMAGES);
    YNET = getvalue(vars.YNET);
    MCABATE = getvalue(vars.MCABATE);
    Y = getvalue(vars.Y);
    E = getvalue(vars.E);
    I = getvalue(vars.I);
    K = getvalue(vars.K);
    PCY = getvalue(vars.PCY)
    MPK = config.γₑ.*YGROSS./K;
    C = getvalue(vars.C);
    CPC = getvalue(vars.CPC);
    PERIODU = getvalue(vars.PERIODU);
    UTILITY = getvalue(vars.UTILITY);
    S = getvalue(vars.S);
    μ = getvalue(vars.μ);
    RI = getvalue(vars.RI);
    scc = -1000.*getdual(eqs.eeq)./(getdual(eqs.kk)+.00000000001);
    mcemis = config.θ₂.*params.θ₁.*μ.^(config.θ₂-1)./params.σ.*1000.;
    ResultsVCJL(years,Mₐₜ,Mₐₜppm,Mᵤₚ,Mₗₒ,CCA,CCAratio,Tₐₜ,FORC,Tₗₒ,YGROSS,DAMAGES,YNET,
               Y,E,I,K,MPK,C,CPC,PERIODU,UTILITY,S,μ,RI,scc,MCABATE,MₐₜAV,PCY,mcemis)
end

function solve(scenario::Scenario, version::VCJL;
    config::OptionsVCJL = options(version),
    solver = IpoptSolver(print_level=5, max_iter=99900,print_frequency_iter=50,sb="yes"))

    model = JuMP.Model(solver = solver);

    params = generate_parameters(config);

    variables = model_vars(version, model, config.N, config.fosslim, config.limμ);

    assign_scenario(scenario, model, config, params, variables);

    equations = model_eqs(model, config, params, variables);

    JuMP.solve(model);
    #JuMP.solve(model);
    #JuMP.solve(model);
    #JuMP.solve(model);
    #JuMP.solve(model);
    #JuMP.solve(model);

    #results = model_results(model, config, params, variables, equations);

    #DICENarrative(config,params,model,scenario,version,variables,equations,results)
    println("Done");
end


text:
using Pidoh
using Random

@testset "Fixed Cooling" begin
    constant = 0.1234
    fc = FixedCooling(constant)
    @test temperature(fc, 1) == constant
    @test temperature(fc, 1000) == constant
end

max_stars_repo_name: DTUComputeTONIA/SDRLS


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %% [markdown]
# https://discourse.julialang.org/t/how-to-reduce-large-memory-allocation-by-lu-with-sparse-matrices/64791/7

# %%
VERSION

# %%
using BenchmarkTools
using Random; Random.seed!(4649373)
using LinearAlgebra
using SparseArrays
using SparseArrays: getcolptr
using SuiteSparse: decrement
using SuiteSparse.UMFPACK: UmfpackLU, UMFVTypes, UMFITypes, umfpack_numeric!

"""https://github.com/JuliaLang/SuiteSparse.jl/blob/master/src/umfpack.jl#L220-L278"""
function lu_org!(F::UmfpackLU, S::SparseMatrixCSC{<:UMFVTypes,<:UMFITypes}; check::Bool=true)
    zerobased = getcolptr(S)[1] == 0
    F.m = size(S, 1)
    F.n = size(S, 2)
    F.colptr = zerobased ? copy(getcolptr(S)) : decrement(getcolptr(S))
    F.rowval = zerobased ? copy(rowvals(S)) : decrement(rowvals(S))
    F.nzval = copy(nonzeros(S))

    umfpack_numeric!(F, reuse_numeric = false)
    check && (issuccess(F) || throw(LinearAlgebra.SingularException(0)))
    return F
end

"""A modification of https://github.com/JuliaLang/SuiteSparse.jl/blob/master/src/umfpack.jl#L220-L278"""
function lu_dot!(F::UmfpackLU, S::SparseMatrixCSC{<:UMFVTypes,<:UMFITypes}; check::Bool=true)
    zerobased = getcolptr(S)[1] == 0
    F.m = size(S, 1)
    F.n = size(S, 2)
    if zerobased
        F.colptr .= getcolptr(S)
        F.rowval .= rowvals(S)
    else
        F.colptr .= getcolptr(S) .- oneunit(eltype(S))
        F.rowval .= rowvals(S) .- oneunit(eltype(S))
    end
    F.nzval .= nonzeros(S)

    umfpack_numeric!(F, reuse_numeric = false)
    check && (issuccess(F) || throw(LinearAlgebra.SingularException(0)))
    return F
end

n = 2^10
X = I + sprandn(n, n, 1e-3)

F = lu(X)
G = lu(X)
@show lu(X).L == lu!(F, X).L
@show lu(X).U == lu!(F, X).U
@show lu!(F, X).L == lu_dot!(G, X).L
@show lu!(F, X).U == lu_dot!(G, X).U
println()
print("lu(X):        ")
@btime lu($X)
print("lu!(F, X):    ")
@btime lu!($F, $X)
print("lu_dot!(F, X):")
@btime lu_dot!($F, $X);

# %%
print("umfpack_numeric!(F, reuse_numeric = false):")
let F = F, S = X
    zerobased = getcolptr(S)[1] == 0
    F.m = size(S, 1)
    F.n = size(S, 2)
    if zerobased
        F.colptr .= getcolptr(S)
        F.rowval .= rowvals(S)
    else
        F.colptr .= getcolptr(S) .- oneunit(eltype(S))
        F.rowval .= rowvals(S) .- oneunit(eltype(S))
    end
    F.nzval .= nonzeros(S)

    @btime umfpack_numeric!(F, reuse_numeric = false)
end;

# %%
@which umfpack_numeric!(F, reuse_numeric = false)

# %% [markdown]
# https://github.com/JuliaLang/SuiteSparse.jl/blob/master/src/umfpack.jl#L384-L400

# %%
mutable struct Factor F::UmfpackLU end
factor = Factor(F)
factor.F === F

# %%

max_stars_repo_path: 0012/lu!(F, X) example.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0012/lu!(F, X) example.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0012/lu!(F, X) example.jl

max_forks_repo_name: genkuroki/public


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: license_manager
using AWS.Compat
using AWS.UUIDs

"""
    CreateLicenseConfiguration()

Creates a license configuration. A license configuration is an abstraction of a customer license agreement that can be consumed and enforced by License Manager. Components include specifications for the license type (licensing by instance, socket, CPU, or vCPU), allowed tenancy (shared tenancy, Dedicated Instance, Dedicated Host, or all of these), license affinity to host (how long a license must be associated with a host), and the number of licenses purchased and used.

# Required Parameters
- `LicenseCountingType`: Dimension used to track the license inventory.
- `Name`: Name of the license configuration.

# Optional Parameters
- `Description`: Description of the license configuration.
- `LicenseCount`: Number of licenses managed by the license configuration.
- `LicenseCountHardLimit`: Indicates whether hard or soft license enforcement is used. Exceeding a hard limit blocks the launch of new instances.
- `LicenseRules`: License rules. The syntax is #name=value (for example, #allowedTenancy=EC2-DedicatedHost). The available rules vary by dimension, as follows.    Cores dimension: allowedTenancy | licenseAffinityToHost | maximumCores | minimumCores     Instances dimension: allowedTenancy | maximumCores | minimumCores | maximumSockets | minimumSockets | maximumVcpus | minimumVcpus     Sockets dimension: allowedTenancy | licenseAffinityToHost | maximumSockets | minimumSockets     vCPUs dimension: allowedTenancy | honorVcpuOptimization | maximumVcpus | minimumVcpus    The unit for licenseAffinityToHost is days and the range is 1 to 180. The possible values for allowedTenancy are EC2-Default, EC2-DedicatedHost, and EC2-DedicatedInstance. The possible values for honorVcpuOptimization are True and False.
- `ProductInformationList`: Product information.
- `Tags`: Tags to add to the license configuration.
"""
create_license_configuration(LicenseCountingType, Name; aws_config::AWSConfig=global_aws_config()) = license_manager("CreateLicenseConfiguration", Dict{String, Any}("LicenseCountingType"=>LicenseCountingType, "Name"=>Name); aws_config=aws_config)
create_license_configuration(LicenseCountingType, Name, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = license_manager("CreateLicenseConfiguration", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("LicenseCountingType"=>LicenseCountingType, "Name"=>Name), args)); aws_config=aws_config)

"""
    DeleteLicenseConfiguration()

Deletes the specified license configuration. You cannot delete a license configuration that is in use.

# Required Parameters
- `LicenseConfigurationArn`: ID of the license configuration.

"""
delete_license_configuration(LicenseConfigurationArn; aws_config::AWSConfig=global_aws_config()) = license_manager("DeleteLicenseConfiguration", Dict{String, Any}("LicenseConfigurationArn"=>LicenseConfigurationArn); aws_config=aws_config)
delete_license_configuration(LicenseConfigurationArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = license_manager("DeleteLicenseConfiguration", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("LicenseConfigurationArn"=>LicenseConfigurationArn), args)); aws_config=aws_config)

"""
    GetLicenseConfiguration()

Gets detailed information about the specified license configuration.

# Required Parameters
- `LicenseConfigurationArn`: Amazon Resource Name (ARN) of the license configuration.

"""
get_license_configuration(LicenseConfigurationArn; aws_config::AWSConfig=global_aws_config()) = license_manager("GetLicenseConfiguration", Dict{String, Any}("LicenseConfigurationArn"=>LicenseConfigurationArn); aws_config=aws_config)
get_license_configuration(LicenseConfigurationArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = license_manager("GetLicenseConfiguration", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("LicenseConfigurationArn"=>LicenseConfigurationArn), args)); aws_config=aws_config)

"""
    GetServiceSettings()

Gets the License Manager settings for the current Region.

"""
get_service_settings(; aws_config::AWSConfig=global_aws_config()) = license_manager("GetServiceSettings"; aws_config=aws_config)
get_service_settings(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = license_manager("GetServiceSettings", args; aws_config=aws_config)

"""
    ListAssociationsForLicenseConfiguration()

Lists the resource associations for the specified license configuration. Resource associations need not consume licenses from a license configuration. For example, an AMI or a stopped instance might not consume a license (depending on the license rules).

# Required Parameters
- `LicenseConfigurationArn`: Amazon Resource Name (ARN) of a license configuration.

# Optional Parameters
- `MaxResults`: Maximum number of results to return in a single call.
- `NextToken`: Token for the next set of results.
"""
list_associations_for_license_configuration(LicenseConfigurationArn; aws_config::AWSConfig=global_aws_config()) = license_manager("ListAssociationsForLicenseConfiguration", Dict{String, Any}("LicenseConfigurationArn"=>LicenseConfigurationArn); aws_config=aws_config)
list_associations_for_license_configuration(LicenseConfigurationArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = license_manager("ListAssociationsForLicenseConfiguration", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("LicenseConfigurationArn"=>LicenseConfigurationArn), args)); aws_config=aws_config)

"""
    ListFailuresForLicenseConfigurationOperations()

Lists the license configuration operations that failed.

# Required Parameters
- `LicenseConfigurationArn`: Amazon Resource Name of the license configuration.

# Optional Parameters
- `MaxResults`: Maximum number of results to return in a single call.
- `NextToken`: Token for the next set of results.
"""
list_failures_for_license_configuration_operations(LicenseConfigurationArn; aws_config::AWSConfig=global_aws_config()) = license_manager("ListFailuresForLicenseConfigurationOperations", Dict{String, Any}("LicenseConfigurationArn"=>LicenseConfigurationArn); aws_config=aws_config)
list_failures_for_license_configuration_operations(LicenseConfigurationArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = license_manager("ListFailuresForLicenseConfigurationOperations", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("LicenseConfigurationArn"=>LicenseConfigurationArn), args)); aws_config=aws_config)

"""
    ListLicenseConfigurations()

Lists the license configurations for your account.

# Optional Parameters
- `Filters`: Filters to scope the results. The following filters and logical operators are supported:    licenseCountingType - The dimension on which licenses are counted. Possible values are vCPU | Instance | Core | Socket. Logical operators are EQUALS | NOT_EQUALS.    enforceLicenseCount - A Boolean value that indicates whether hard license enforcement is used. Logical operators are EQUALS | NOT_EQUALS.    usagelimitExceeded - A Boolean value that indicates whether the available licenses have been exceeded. Logical operators are EQUALS | NOT_EQUALS.
- `LicenseConfigurationArns`: Amazon Resource Names (ARN) of the license configurations.
- `MaxResults`: Maximum number of results to return in a single call.
- `NextToken`: Token for the next set of results.
"""
list_license_configurations(; aws_config::AWSConfig=global_aws_config()) = license_manager("ListLicenseConfigurations"; aws_config=aws_config)
list_license_configurations(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = license_manager("ListLicenseConfigurations", args; aws_config=aws_config)

"""
    ListLicenseSpecificationsForResource()

Describes the license configurations for the specified resource.

# Required Parameters
- `ResourceArn`: Amazon Resource Name (ARN) of a resource that has an associated license configuration.

# Optional Parameters
- `MaxResults`: Maximum number of results to return in a single call.
- `NextToken`: Token for the next set of results.
"""
list_license_specifications_for_resource(ResourceArn; aws_config::AWSConfig=global_aws_config()) = license_manager("ListLicenseSpecificationsForResource", Dict{String, Any}("ResourceArn"=>ResourceArn); aws_config=aws_config)
list_license_specifications_for_resource(ResourceArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = license_manager("ListLicenseSpecificationsForResource", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ResourceArn"=>ResourceArn), args)); aws_config=aws_config)

"""
    ListResourceInventory()

Lists resources managed using Systems Manager inventory.

# Optional Parameters
- `Filters`: Filters to scope the results. The following filters and logical operators are supported:    account_id - The ID of the AWS account that owns the resource. Logical operators are EQUALS | NOT_EQUALS.    application_name - The name of the application. Logical operators are EQUALS | BEGINS_WITH.    license_included - The type of license included. Logical operators are EQUALS | NOT_EQUALS. Possible values are sql-server-enterprise | sql-server-standard | sql-server-web | windows-server-datacenter.    platform - The platform of the resource. Logical operators are EQUALS | BEGINS_WITH.    resource_id - The ID of the resource. Logical operators are EQUALS | NOT_EQUALS.
- `MaxResults`: Maximum number of results to return in a single call.
- `NextToken`: Token for the next set of results.
"""
list_resource_inventory(; aws_config::AWSConfig=global_aws_config()) = license_manager("ListResourceInventory"; aws_config=aws_config)
list_resource_inventory(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = license_manager("ListResourceInventory", args; aws_config=aws_config)

"""
    ListTagsForResource()

Lists the tags for the specified license configuration.

# Required Parameters
- `ResourceArn`: Amazon Resource Name (ARN) of the license configuration.

"""
list_tags_for_resource(ResourceArn; aws_config::AWSConfig=global_aws_config()) = license_manager("ListTagsForResource", Dict{String, Any}("ResourceArn"=>ResourceArn); aws_config=aws_config)
list_tags_for_resource(ResourceArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = license_manager("ListTagsForResource", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ResourceArn"=>ResourceArn), args)); aws_config=aws_config)

"""
    ListUsageForLicenseConfiguration()

Lists all license usage records for a license configuration, displaying license consumption details by resource at a selected point in time. Use this action to audit the current license consumption for any license inventory and configuration.

# Required Parameters
- `LicenseConfigurationArn`: Amazon Resource Name (ARN) of the license configuration.

# Optional Parameters
- `Filters`: Filters to scope the results. The following filters and logical operators are supported:    resourceArn - The ARN of the license configuration resource. Logical operators are EQUALS | NOT_EQUALS.    resourceType - The resource type (EC2_INSTANCE | EC2_HOST | EC2_AMI | SYSTEMS_MANAGER_MANAGED_INSTANCE). Logical operators are EQUALS | NOT_EQUALS.    resourceAccount - The ID of the account that owns the resource. Logical operators are EQUALS | NOT_EQUALS.
- `MaxResults`: Maximum number of results to return in a single call.
- `NextToken`: Token for the next set of results.
"""
list_usage_for_license_configuration(LicenseConfigurationArn; aws_config::AWSConfig=global_aws_config()) = license_manager("ListUsageForLicenseConfiguration", Dict{String, Any}("LicenseConfigurationArn"=>LicenseConfigurationArn); aws_config=aws_config)
list_usage_for_license_configuration(LicenseConfigurationArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = license_manager("ListUsageForLicenseConfiguration", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("LicenseConfigurationArn"=>LicenseConfigurationArn), args)); aws_config=aws_config)

"""
    TagResource()

Adds the specified tags to the specified license configuration.

# Required Parameters
- `ResourceArn`: Amazon Resource Name (ARN) of the license configuration.
- `Tags`: One or more tags.

"""
tag_resource(ResourceArn, Tags; aws_config::AWSConfig=global_aws_config()) = license_manager("TagResource", Dict{String, Any}("ResourceArn"=>ResourceArn, "Tags"=>Tags); aws_config=aws_config)
tag_resource(ResourceArn, Tags, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = license_manager("TagResource", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ResourceArn"=>ResourceArn, "Tags"=>Tags), args)); aws_config=aws_config)

"""
    UntagResource()

Removes the specified tags from the specified license configuration.

# Required Parameters
- `ResourceArn`: Amazon Resource Name (ARN) of the license configuration.
- `TagKeys`: Keys identifying the tags to remove.

"""
untag_resource(ResourceArn, TagKeys; aws_config::AWSConfig=global_aws_config()) = license_manager("UntagResource", Dict{String, Any}("ResourceArn"=>ResourceArn, "TagKeys"=>TagKeys); aws_config=aws_config)
untag_resource(ResourceArn, TagKeys, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = license_manager("UntagResource", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ResourceArn"=>ResourceArn, "TagKeys"=>TagKeys), args)); aws_config=aws_config)

"""
    UpdateLicenseConfiguration()

Modifies the attributes of an existing license configuration.

# Required Parameters
- `LicenseConfigurationArn`: Amazon Resource Name (ARN) of the license configuration.

# Optional Parameters
- `Description`: New description of the license configuration.
- `LicenseConfigurationStatus`: New status of the license configuration.
- `LicenseCount`: New number of licenses managed by the license configuration.
- `LicenseCountHardLimit`: New hard limit of the number of available licenses.
- `LicenseRules`: New license rule. The only rule that you can add after you create a license configuration is licenseAffinityToHost.
- `Name`: New name of the license configuration.
- `ProductInformationList`: New product information.
"""
update_license_configuration(LicenseConfigurationArn; aws_config::AWSConfig=global_aws_config()) = license_manager("UpdateLicenseConfiguration", Dict{String, Any}("LicenseConfigurationArn"=>LicenseConfigurationArn); aws_config=aws_config)
update_license_configuration(LicenseConfigurationArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = license_manager("UpdateLicenseConfiguration", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("LicenseConfigurationArn"=>LicenseConfigurationArn), args)); aws_config=aws_config)

"""
    UpdateLicenseSpecificationsForResource()

Adds or removes the specified license configurations for the specified AWS resource. You can update the license specifications of AMIs, instances, and hosts. You cannot update the license specifications for launch templates and AWS CloudFormation templates, as they send license configurations to the operation that creates the resource.

# Required Parameters
- `ResourceArn`: Amazon Resource Name (ARN) of the AWS resource.

# Optional Parameters
- `AddLicenseSpecifications`: ARNs of the license configurations to add.
- `RemoveLicenseSpecifications`: ARNs of the license configurations to remove.
"""
update_license_specifications_for_resource(ResourceArn; aws_config::AWSConfig=global_aws_config()) = license_manager("UpdateLicenseSpecificationsForResource", Dict{String, Any}("ResourceArn"=>ResourceArn); aws_config=aws_config)
update_license_specifications_for_resource(ResourceArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = license_manager("UpdateLicenseSpecificationsForResource", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ResourceArn"=>ResourceArn), args)); aws_config=aws_config)

"""
    UpdateServiceSettings()

Updates License Manager settings for the current Region.

# Optional Parameters
- `EnableCrossAccountsDiscovery`: Activates cross-account discovery.
- `OrganizationConfiguration`: Enables integration with AWS Organizations for cross-account discovery.
- `S3BucketArn`: Amazon Resource Name (ARN) of the Amazon S3 bucket where the License Manager information is stored.
- `SnsTopicArn`: Amazon Resource Name (ARN) of the Amazon SNS topic used for License Manager alerts.
"""
update_service_settings(; aws_config::AWSConfig=global_aws_config()) = license_manager("UpdateServiceSettings"; aws_config=aws_config)
update_service_settings(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = license_manager("UpdateServiceSettings", args; aws_config=aws_config)


text:
# tests for the dynamics type and related functions.

using Base.Test
# using PyPlot

const OUTPUTS_DIRECTORY = "data/plots/"

# Dynamics type tests
function test_dyn_constructor()
    m = MersenneTwister(1)
    d = Dynamics(3, 1., 4., 5., 5., false, m)
    @test d.dt == 3
    @test d.dt_step == 1.
    @test d.pid_omega == 4
    @test d.pid_b_std == 5
end

# pid controller tests
function test_build_pid_policy()
    m = MersenneTwister(1)
    d = Dynamics(3, 1., 4., 5., 5., false, m)
    action = deg2rad(10)
    uav_state_vector = collect(1.:8)
    pid = build_pid_policy(d, action)
    actual_action = pid(uav_state_vector)
    @test actual_action == -173.20747319680908
end

# step tests
function test_multi_uav_step()
    m = MersenneTwister(1)
    omega = .3
    d = Dynamics(1, 1., omega, 0., 0., false, m)
    uav_state_own = UAVState(ones(Float64, 8)...)
    uav_state_int = UAVState(ones(Float64, 8)...)
    state = State([uav_state_own; uav_state_int])
    actions = [deg2rad(10); deg2rad(10)]
    next_state = step(d, state, actions)
end

function test_single_uav_step()
    m = MersenneTwister(1)
    omega = .3
    d = Dynamics(0, 0., omega, 0., 0., false, m)
    action = deg2rad(10)
    pid = build_pid_policy(d, action)
    uav_state = UAVState(ones(Float64, 8)...)
    next_uav_state = step(uav_state, pid, 1, 1.)
end

tic()
# dynamics tests
test_dyn_constructor()

# # pid controller tests
test_build_pid_policy()

# step tests
test_multi_uav_step()
test_single_uav_step()
toc()


text:
module PlatformEngineTests
import ..Pkg # ensure we are using the correct Pkg

using Test, Pkg.PlatformEngines, Pkg.BinaryPlatforms, SHA

# Explicitly probe platform engines in verbose mode to get coverage and make
# CI debugging easier
probe_platform_engines!(;verbose=true)

@testset "Tarball listing parsing" begin
    fake_7z_output = """
	7-Zip [64] 9.20  Copyright (c) 1999-2010 Igor Pavlov  2010-11-18
	p7zip Version 9.20 (locale=en_US.UTF-8,Utf16=on,HugeFiles=on,4 CPUs)
	Listing archive:
	--
	Path =
	Type = tar
	   Date      Time    Attr         Size   Compressed  Name
	------------------- ----- ------------ ------------  ------------------------
	2017-04-10 14:45:00 D....            0            0  bin
	2017-04-10 14:44:59 .....          211          512  bin/socrates
	------------------- ----- ------------ ------------  ------------------------
									   211          512  1 files, 1 folders
	"""
	@test parse_7z_list(fake_7z_output) == ["bin/socrates"]

	fake_tar_output = """
	bin/
	bin/socrates
	"""
	@test parse_tar_list(fake_tar_output) == normpath.(["bin/socrates"])
end

@testset "Packaging" begin
    # Gotta set this guy up beforehand
    tarball_path = nothing
    tarball_hash = nothing

    mktempdir() do prefix
        # Create random files
        mkpath(joinpath(prefix, "bin"))
        mkpath(joinpath(prefix, "lib"))
        mkpath(joinpath(prefix, "etc"))
        bar_path = joinpath(prefix, "bin", "bar.sh")
        open(bar_path, "w") do f
            write(f, "#!/bin/sh\n")
            write(f, "echo yolo\n")
        end
        baz_path = joinpath(prefix, "lib", "baz.so")
        open(baz_path, "w") do f
            write(f, "this is not an actual .so\n")
        end

        qux_path = joinpath(prefix, "etc", "qux.conf")
        open(qux_path, "w") do f
            write(f, "use_julia=true\n")
        end

        # Next, package it up as a .tar.gz file
        mktempdir() do output_dir
            tarball_path =  joinpath(output_dir, "foo.tar.gz")
            package(prefix, tarball_path)
            @test isfile(tarball_path)

            # Test that we can inspect the contents of the tarball
            contents = list_tarball_files(tarball_path)
            @test joinpath("bin", "bar.sh") in contents
            @test joinpath("lib", "baz.so") in contents
            @test joinpath("etc", "qux.conf") in contents
        end
    end

end

@testset "Verification" begin
    mktempdir() do prefix
        foo_path = joinpath(prefix, "foo")
        open(foo_path, "w") do file
            write(file, "test")
        end
        foo_hash = bytes2hex(sha256("test"))

        # Check that verifying with the right hash works
        @test_logs (:info, r"No hash cache found") match_mode=:any begin
            ret, status = verify(foo_path, foo_hash; verbose=true, report_cache_status=true)
            @test ret == true
            @test status == :hash_cache_missing
        end

        # Check that it created a .sha256 file
        @test isfile("$(foo_path).sha256")

        # Check that it verifies the second time around properly
        @test_logs (:info, r"Hash cache is consistent") match_mode=:any begin
            ret, status = verify(foo_path, foo_hash; verbose=true, report_cache_status=true)
            @test ret == true
            @test status == :hash_cache_consistent
        end

        # Sleep for imprecise filesystems
        sleep(2)

        # Get coverage of messing with different parts of the verification chain
        touch(foo_path)
        @test_logs (:info, r"File has been modified") match_mode=:any begin
            ret, status = verify(foo_path, foo_hash; verbose=true, report_cache_status=true)
            @test ret == true
            @test status == :file_modified
        end

        # Ensure that we print an error when verification fails
        rm("$(foo_path).sha256"; force=true)
        @test_logs (:error, r"Hash Mismatch!") match_mode=:any begin
            @test !verify(foo_path, "0"^64; verbose=true)
        end

        # Ensure that incorrect lengths cause an exception
        @test_throws ErrorException verify(foo_path, "0"^65; verbose=true)

        # Ensure that messing with the hash file works properly
        touch(foo_path)
        @test verify(foo_path, foo_hash; verbose=true)
        open("$(foo_path).sha256", "w") do file
            write(file, "this is not the right hash")
        end
        @test_logs (:info, r"hash cache invalidated") match_mode=:any begin
            ret, status = verify(foo_path, foo_hash; verbose=true, report_cache_status=true)
            @test ret == true
            @test status == :hash_cache_mismatch
        end

        # Ensure that messing with the actual file works properly
        open("$(foo_path)", "w") do file
            write(file, "this is not the right content")
        end

        # Delete hash cache file to force re-verification
        rm("$(foo_path).sha256"; force=true)
        @test_logs (:error, r"Hash Mismatch!") match_mode=:any begin
            ret, status = verify(foo_path, foo_hash; verbose=true, report_cache_status=true)
            @test ret == false
            @test status == :hash_mismatch
        end
    end
end

const socrates_urls = [
    "https://github.com/staticfloat/small_bin/raw/f1a92f5eafbd30a0c6a8efb6947485b0f6d1bec3/socrates.tar.gz" =>
    "e65d2f13f2085f2c279830e863292312a72930fee5ba3c792b14c33ce5c5cc58",
    "https://github.com/staticfloat/small_bin/raw/f1a92f5eafbd30a0c6a8efb6947485b0f6d1bec3/socrates.tar.bz2" =>
    "13fc17b97be41763b02cbb80e9d048302cec3bd3d446c2ed6e8210bddcd3ac76",
    "https://github.com/staticfloat/small_bin/raw/f1a92f5eafbd30a0c6a8efb6947485b0f6d1bec3/socrates.tar.xz" =>
    "61bcf109fcb749ee7b6a570a6057602c08c836b6f81091eab7aa5f5870ec6475",
]
const socrates_hash = "adcbcf15674eafe8905093183d9ab997cbfba9056fc7dde8bfa5a22dfcfb4967"

@testset "Downloading" begin
    for (url, hash) in socrates_urls
        mktempdir() do prefix
            tarball_path = joinpath(prefix, "download_target.tar$(splitext(url)[2])")

            target_dir = joinpath(prefix, "target")
            download_verify_unpack(url, hash, target_dir; tarball_path=tarball_path, verbose=true)

            # Test downloading a second time, to get the "already exists" path
            download_verify_unpack(url, hash, target_dir; tarball_path=tarball_path, verbose=true)

            # And a third time, after corrupting it, to get the "redownloading" path
            open(tarball_path, "w") do io
                println(io, "corruptify")
            end
            download_verify_unpack(url, hash, target_dir; tarball_path=tarball_path, verbose=true, force=true)

            # Test that it has the contents we expect
            socrates_path = joinpath(target_dir, "bin", "socrates")
            @test isfile(socrates_path)
            unpacked_hash = open(socrates_path) do f
                bytes2hex(sha256(f))
            end
            @test unpacked_hash == socrates_hash
        end
    end
end

const collapse_url = "https://github.com/staticfloat/small_bin/raw/master/collapse_the_symlink/collapse_the_symlink.tar.gz"
const collapse_hash = "956c1201405f64d3465cc28cb0dec9d63c11a08cad28c381e13bb22e1fc469d3"
@testset "Copyderef unpacking" begin
    withenv("BINARYPROVIDER_COPYDEREF" => "true") do
        mktempdir() do prefix
            target_dir = joinpath(prefix, "target")
            download_verify_unpack(collapse_url, collapse_hash, target_dir; verbose=true)

            # Test that we get the files we expect
            @test isfile(joinpath(target_dir, "collapse_the_symlink", "foo"))
            @test isfile(joinpath(target_dir, "collapse_the_symlink", "foo.1"))
            @test isfile(joinpath(target_dir, "collapse_the_symlink", "foo.1.1"))

            # Test that these are definitely not links
            @test !islink(joinpath(target_dir, "collapse_the_symlink", "foo"))
            @test !islink(joinpath(target_dir, "collapse_the_symlink", "foo.1.1"))

            # Test that broken symlinks get transparently dropped
            @test !ispath(joinpath(target_dir, "collapse_the_symlink", "broken"))
        end
    end
end

@testset "Download GitHub API #88" begin
    mktempdir() do tmp
        PlatformEngines.download("https://api.github.com/repos/JuliaPackaging/BinaryProvider.jl/tarball/c2a4fc38f29eb81d66e3322e585d0199722e5d71", joinpath(tmp, "BinaryProvider"); verbose=true)
        @test isfile(joinpath(tmp, "BinaryProvider"))
    end
end

end # module


text:
module DemoLoopVectorization

export comparemethods
using LoopVectorization
using BenchmarkTools

#Entry point when called as an executable
function julia_main()::Cint
    try
        comparemethods()
        return (0)
    catch e
        @error "$e"
        return (1)
    end
end

"""
    threemeans(x::Array{Float64,2})
Returns a tuple of column means of x, column means of positive part of x,
column means of negative part of x.
"""
function threemeans(x::Array{Float64,2})
    nr, nc = size(x)
    epes = Vector{Float64}(undef, nc)
    enes = Vector{Float64}(undef, nc)
    ees = Vector{Float64}(undef, nc)
    @avx for j = 1:nc
        epe = 0.0
        ene = 0.0
        ee = 0.0
        for i = 1:nr
            ee += x[i, j]
            epe += max(0.0, x[i, j])
            ene += min(0.0, x[i, j])
        end
        epes[j] = epe / nr
        enes[j] = ene / nr
        ees[j] = ee / nr
    end
    epes, enes, ees
end

function threemeans_noavx(x::Array{Float64,2})
    nr, nc = size(x)
    epes = Vector{Float64}(undef, nc)
    enes = Vector{Float64}(undef, nc)
    ees = Vector{Float64}(undef, nc)
    for j = 1:nc
        epe = ene = ee = 0.0
        for i = 1:nr
            ee += x[i, j]
            if x[i, j] > 0
                epe += x[i, j]
            else
                ene += x[i, j]
            end
        end
        epes[j] = epe / nr
        enes[j] = ene / nr
        ees[j] = ee / nr
    end
    epes, enes, ees
end

function threemeans_naive(x::Array{Float64,2})
    epe = vec(mean(max.(x, 0.0), dims = 1))
    ene = vec(mean(min.(x, 0.0), dims = 1))
    ee = vec(mean(x, dims = 1))
    epe, ene, ee
end

#= 13 Dec 2021
julia> comparemethods()
threemeans  475.800 μs (3 allocations: 23.81 KiB)
threemeans_noavx  1.537 ms (3 allocations: 23.81 KiB)
threemeans_naive  6.343 ms (31 allocations: 15.28 MiB)
res1[1] ≈ res2[1] && (res1[2] ≈ res2[2] && res1[3] ≈ res2[3]) = true
res1[1] ≈ res3[1] && (res1[2] ≈ res3[2] && res1[3] ≈ res3[3]) = true
=#
function comparemethods()
    x = rand(1000, 1000)
    print("threemeans")
    res1 = @btime(threemeans($x))
    print("threemeans_noavx")
    res2 = @btime(threemeans_noavx($x))
    print("threemeans_naive")
    res3 = @btime(threemeans_naive($x))

    @show res1[1] ≈ res2[1] && res1[2] ≈ res2[2] && res1[3] ≈ res2[3]
    @show res1[1] ≈ res3[1] && res1[2] ≈ res3[2] && res1[3] ≈ res3[3]
    nothing
end

end #module


text:
# Chapter 8.1 Maximum-likelihood Estimation
## Visualizing the likelihood function
using Plots
N = 50
S = range(1., N, step=0.1)
θ = range(0.1, 0.9,length=100)
L(S, θ) = S*log(θ) + (N-S)*log(1. - θ)
# the surface plot
plotlyjs() # this backend allows rotating with mouse
p1 = surface(S, θ, (S, θ) -> L(S, θ), color=:jet, xlabel="S", ylabel="θ",
    title="ℒ(θ|S)")
# the bird's eye view
gr()
p2 = heatmap(S, θ, (S, θ) -> L(S, θ), color=:jet, xlabel="S", ylabel="θ",
    title="bird's eye view")
SS = 25 # set the value for the slice here
vline!([SS], label=false, color=:black)
# slice through the log likelihood
p3 = plot(θ, θ -> L(SS, θ),label=false, xlabel="θ", title="ℒ(θ|S=$SS)")
plot(p2, p3)

## ML estimation for single-photon imaging
using Images, Distributions
using ImageView:imshow
download("https://probability4datascience.com/data/cameraman.tif", "./cameraman.tif")
λ  = Float64.(load("cameraman.tif"))
T = 100
x = [rand.(Poisson.(λ)) for _=1:T]
y = [x[i] .>= 1. for i=1:T]
λhat = -log.(1. .- mean(y))
imshow(x[1]) # a single sample image
imshow(λhat) # the ML recovered image

## Chapter 8.2 Properties of the ML estimation
# Visualizing the invariance principle
using Plots
N = 50
S = 20
θ = range(0.1,0.9,length=1000)
η  = -log.(1. .- θ) # the reparameterization transformation
L₁(θ) = S*log(θ) + (N-S)*log(1. - θ) # log likelihood in terms of θ
L₂(η) = S*log(1. - exp(-η)) - (N-S)*η # log likelihood in terms of η
# plot the log likelihood with θ parameterization
p4 = plot(θ, θ -> L₁(θ), linewidth=5, color=:blue, label=false)
xlabel!("θ")
title!("ℒ₁(θ|S=$S)")
vline!([0.4], label=false, color=:red)
# plot the transformation
p2 = plot(θ, η, linewidth=5, color=:black, label=false)
xlabel!("θ")
ylabel!("η")
title!("η(θ)")
vline!([0.4], label=false, color=:red)
hline!([0.5], label=false, color=:green)
# plot the log likelihood with η parameterization
p1 = plot(L₂.(η), η,  linewidth=5, color=:blue, label=false, xflip=true)
ylabel!("η")
hline!([0.5], label=false, color=:green)
title!("ℒ₂(η|S=$S)")
# blank place holder plot
p3 = plot(grid=false, xaxis=false, yaxis=false, xticks=false, yticks=false)
# plot the whole set in the final image
plot(p1, p2, p3, p4)

## Chapter 8.3 Maximum-a-Posteriori Estimation
# Influence of the priors
# likelihood is N(θ,1)
using Plots, Distributions
N = 5       # sample size
μ₀  = 0.0   # prior mean
σ₀  = 1.    # prior std. dev.
θ = 5.      # true mean of sample data
x = rand(Normal(θ,1), N)    # the sample data
xbar  = mean(x) # mean of sample data
t = range(-3.,7.,length=1000)
θhat = (σ₀^2. * xbar + μ₀/N)/(σ₀^2. + 1. /N) # MAP
σhat = sqrt(1. / (1. / σ₀^2. + N))  # s.d. posterior
p0 = pdf(Normal(xbar,1.0), t)       # likelihood
p1 = pdf(Normal(θhat, σhat), t)  # posterior
prior = pdf(Normal(μ₀,σ₀),t)     # prior, divided by 10 to
plot(t, [p0, p1, prior], label=["Likelihood" "Posterior" "prior"], linewidth=3, legend=:topleft)
plot!([x], [0.1*ones(N)], line=:stem, label="data", marker=:circle, linewidth=3, title="N = $N")

## Conjugate priors
using Plots, Distributions
σ₀ = 0.25
μ₀ = 0.0
μ  = 1.
σ  = 0.25
Nset = [0 1 2 5 8 12 20]
x0 = σ .* randn(100)
posterior = zeros(1000,7)
t = range(-1.,1.5,length=1000)
p0 = pdf(Normal(0.,1.), t)
for i=1:7
    N = Nset[i]
    x = x0[1:N]
    θ  = μ * (N*σ₀^2.)/(N*σ₀^2. + σ ^2.) + μ₀*(σ^2.)/(N*σ₀^2. + σ^2.)
    σN= sqrt(1. /(1. /σ₀^2. + N/σ^2.))
    posterior[:,i] = pdf(Normal(θ, σN), t)
end
plot(t, posterior, linewidth=2.,
    label=["N = 0" "N = 1" "N = 2" "N = 5" "N = 8" "N = 12" "N = 20"],
    legend=:topleft)


text:
@reexport module Synapse

using DataStructures: Queue, enqueue!, dequeue!, empty!
using DataStructures: CircularBuffer, fill!, push!, empty!
using Adapt
using CUDA
using LoopVectorization

import ..SpikingNN: excite!, spike!, evaluate!, reset!, isactive

export AbstractSynapse, QueuedSynapse, DelayedSynapse,
       excite!, spike!, reset!, isactive

"""
    AbstractSynapse

Inherit from this type to create a concrete synapse.
"""
abstract type AbstractSynapse end

"""
    excite!(synapse::AbstractSynapse, spikes::Vector{<:Integer})
    excite!(synapses::AbstractArray{<:AbstractSynapse}, spikes::Vector{<:Integer})

Excite a `synapse` with a vector of spikes by calling `excite!(synapse, spike) for spike in spikes`.
"""
excite!(synapse::AbstractSynapse, spikes::Vector{<:Integer}) = map(x -> excite!(synapse, x), spikes)
excite!(synapses::AbstractArray{<:AbstractSynapse}, spikes::Vector{<:Integer}) = map(x -> excite!(synapses, x), spikes)

"""
    spike!(synapse::AbstractSynapse, spike::Integer; dt::Real = 1.0)
    spike!(synapse::AbstractArray{<:AbstractSynapse}, spikes::AbstractArray{<:Integer}; dt::Real = 1.0)

Notify a synapse that the post-synaptic neuron has released a spike.
The default implmentation is to do nothing. Override this behavior by dispatching on your synapse type.
"""
spike!(synapse::AbstractSynapse, spike::Integer; dt::Real = 1.0) = nothing
spike!(synapses::AbstractArray{<:AbstractSynapse}, spikes; dt::Real = 1.0) = nothing

"""
    QueuedSynapse{ST<:AbstractSynapse, IT<:Integer}

A `QueuedSynapse` excites its internal synapse when the timestep matches the head of the queue.
Wrapping a synapse in this type allows you to pre-load several spike excitation times, and the
  internal synapse will be excited as those time stamps are evaluated.
This can be useful for cases where it is more efficient to load all the input spikes before simulation.

*Note: currently only supported on CPU.*
"""
struct QueuedSynapse{ST<:AbstractSynapse, IT<:Integer} <: AbstractSynapse
    core::ST
    queue::Queue{IT}
end
QueuedSynapse{IT}(synapse) where {IT<:Integer} = QueuedSynapse{typeof(synapse), IT}(synapse, Queue{IT}())
QueuedSynapse(synapse) = QueuedSynapse{typeof(synapse), Int}(synapse, Queue{Int}())

_ispending(queue, t) = !isempty(queue) && first(queue) <= t
function _shiftspike!(queue, lastspike, t)
    while _ispending(queue, t)
        lastspike = dequeue!(queue)
    end

    return lastspike
end
function _shiftspike!(queues::AbstractArray, lastspikes, t)
    pending = map(x -> _ispending(x, t), queues)
    while any(pending)
        @. lastspikes[pending] = dequeue!(queues[pending])
        pending = map(x -> _ispending(x, t), queues)
    end

    return lastspikes
end

"""
    excite!(synapse::QueuedSynapse, spike::Integer)
    excite!(synapses::AbstractArray{<:QueuedSynapse}, spike::Integer)

Excite `synapse` with a `spike` (`spike` == time step of spike) by pushing
  `spike` onto `synapse.queue`.
"""
excite!(synapse::QueuedSynapse, spike::Integer) = enqueue!(synapse.queue, spike)
excite!(synapses::T, spike::Integer) where T<:AbstractArray{<:QueuedSynapse} =
    map(x -> enqueue!(x, spike), synapses.queue)

isactive(synapse::QueuedSynapse, t::Integer; dt::Real = 1.0) = _ispending(synapse.queue, t) || isactive(synapse.core, t; dt = dt)
isactive(synapses::T, t::Integer; dt::Real = 1.0) where T<:AbstractArray{<:QueuedSynapse} =
    any(map(x -> _ispending(x, t), synapses.queue)) || isactive(synapses.core, t; dt = dt)

"""
    evaluate!(synapse::QueuedSynapse, t::Integer; dt::Real = 1.0)
    (synapse::QueuedSynapse)(t::Integer; dt::Real = 1.0)
    evaluate!(synapses::AbstractArray{<:QueuedSynapse}, t::Integer; dt::Real = 1.0)

Evaluate `synapse` at time `t` by first exciting `synapse.core` with a spike if
  there is one to process, then evaluating `synapse.core`.
"""
function evaluate!(synapse::QueuedSynapse, t::Integer; dt::Real = 1.0)
    excite!(synapse.core, _shiftspike!(synapse.queue, 0, t))

    return synapse.core(t; dt = dt)
end
(synapse::QueuedSynapse)(t::Integer; dt::Real = 1.0) = evaluate!(synapse, t; dt = dt)
function evaluate!(synapses::T, t::Integer; dt::Real = 1.0) where T<:AbstractArray{<:QueuedSynapse}
    lastspikes = _shiftspike!(synapses.queue, zeros(Int, size(synapses)), t)
    @inbounds for i in eachindex(synapses)
        excite!(view(synapses.core, i), lastspikes[i])
    end

    return evaluate!(synapses.core, t; dt = dt)
end

"""
    reset!(synapse::QueuedSynapse)
    reset!(synapses::AbstractArray{<:QueuedSynapse})

Clear `synapse.queue` and reset `synapse.core`.
"""
function reset!(synapse::QueuedSynapse)
    empty!(synapse.queue)
    reset!(synapse.core)
end
function reset!(synapses::T) where T<:AbstractArray{<:QueuedSynapse}
    empty!.(synapses.queue)
    reset!(synapses.core)
end

"""
    DelayedSynapse

A `DelayedSynapse` adds a fixed delay to spikes when exciting its internal synapse.
"""
struct DelayedSynapse{T<:Real, ST<:AbstractSynapse} <: AbstractSynapse
    core::ST
    delay::T
end

"""
    excite!(synapse::DelayedSynapse, spike::Integer)
    excite!(synapses::AbstractArray{<:DelayedSynapse}, spike::Integer)

Excite `synapse.core` with a `spike` + `synapse.delay` (`spike` == time step of spike).
"""
excite!(synapse::DelayedSynapse, spike::Integer) = excite!(synapse.core, spike + synapse.delay)
function excite!(synapses::T, spike::Integer) where T<:AbstractArray{<:DelayedSynapse}
    delayedspikes = adapt(Array{eltype(synapses.delay), ndims(synapses)}, spike .+ synapses.delay)
    if spike > 0
        @inbounds for i in eachindex(synapses)
            excite!(view(synapses.core, i), delayedspikes[i])
        end
    end
end

isactive(synapse::DelayedSynapse, t::Integer; dt::Real = 1.0) = isactive(synapse.core, t; dt = dt)
isactive(synapses::T, t::Integer; dt::Real = 1.0) where T<:AbstractArray{<:DelayedSynapse} =
    isactive(synapses.core, t; dt = dt)

"""
    evaluate!(synapse::DelayedSynapse, t::Integer; dt::Real = 1.0)
    (synapse::DelayedSynapse)(t::Integer; dt::Real = 1.0)
    evaluate!(synapses::AbstractArray{<:DelayedSynapse}, t::Integer; dt::Real = 1.0)

Evaluate `synapse.core` at time `t`.
"""
evaluate!(synapse::DelayedSynapse, t::Integer; dt::Real = 1.0) = evaluate!(synapse.core, t; dt = dt)
(synapse::DelayedSynapse)(t::Integer; dt::Real = 1.0) = evaluate!(synapse, t; dt = dt)
evaluate!(synapses::T, t::Integer; dt::Real = 1.0) where T<:AbstractArray{<:DelayedSynapse} =
    evaluate!(synapses.core, t; dt = dt)
evaluate!(current, synapses::T, t::Integer; dt::Real = 1.0) where T<:AbstractArray{<:DelayedSynapse} =
    evaluate!(current, synapses.core, t; dt = dt)

"""
    reset!(synapse::DelayedSynapse)
    reset!(synapses::AbstractArray{<:DelayedSynapse})

Reset `synapse.core`.
"""
reset!(synapse::DelayedSynapse) = reset!(synapse.core)
reset!(synapses::T) where T<:AbstractArray{<:DelayedSynapse} = reset!(synapses.core)

"""
    Delta{IT<:Integer, VT<:Real}

A synapse representing a Dirac-delta at `lastspike` with amplitude `q`.
"""
mutable struct Delta{IT<:Integer, VT<:Real} <: AbstractSynapse
    lastspike::VT
    q::VT
end

"""
    Delta{IT, VT}(;q::Real = 1)
    Delta(;q::Real = 1)

Create a Dirac-delta synapse with amplitude `q`.
"""
Delta{IT, VT}(;q::Real = 1) where {IT<:Integer, VT<:Real} = Delta{IT, VT}(-Inf, q)
Delta(;q::Real = 1) = Delta{Int, Float32}(q = q)

"""
    excite!(synapse::Delta, spike::Integer)
    excite!(synapses::AbstractArray{<:Delta}, spike::Integer)

Excite `synapse` with a `spike` (`spike` == time step of spike).
"""
excite!(synapse::Delta, spike::Integer) = (spike > 0) && (synapse.lastspike = spike)
excite!(synapses::T, spike::Integer) where T<:AbstractArray{<:Delta} = (spike > 0) && (synapses.lastspike .= spike)

isactive(synapse::Delta, t::Integer; dt::Real = 1.0) = (t * dt == synapse.lastspike)
isactive(synapses::T, t::Integer; dt::Real = 1.0) where T<:AbstractArray{<:Delta} = any(t * dt .== synapses.lastspike)

"""
    delta(t::Real, lastspike, q)
    delta(t::Real, lastspike::AbstractArray{<:Real}, q::AbstractArray{<:Real})
    delta(t::Real, lastspike::CuVecOrMat{<:Real}, q::CuVecOrMat{<:Real})

Evaluate a Dirac-delta synapse.
Use `CuVector` instead of `Vector` for GPU support.

# Fields
- `t`: current time
- `lastspike`: last pre-synaptic spike time
- `q`: amplitude
"""
delta(t::Real, q) = (t ≈ 0) ? q : zero(q)
delta(t::AbstractArray{<:Real}, q::AbstractArray{<:Real}) = delta.(t, q)

"""
    evaluate!(synapse::Delta, t::Integer; dt::Real = 1.0)
    (synapse::Delta)(t::Integer; dt::Real = 1.0)
    evaluate!(synapses::AbstractArray{<:Delta}, t::Integer; dt::Real = 1.0)

Return `synapse.q` if `t == synapse.lastspike` otherwise return zero.
"""
evaluate!(synapse::Delta, t::Integer; dt::Real = 1.0) = delta((t - synapse.lastspike) * dt, synapse.q)
(synapse::Delta)(t::Integer; dt::Real = 1.0) = evaluate!(synapse, t; dt = dt)
evaluate!(synapses::T, t::Integer; dt::Real = 1.0) where T<:AbstractArray{<:Delta} =
    delta(t * dt, synapses.lastspike * dt, synapses.q)

"""
    reset!(synapse::Delta)
    reset!(synapses::AbstractArray{<:Delta})

Reset `synapse`.
"""
reset!(synapse::Delta) = (synapse.lastspike = -Inf)
reset!(synapses::T) where T<:AbstractArray{<:Delta} = (synapses.lastspike .= -Inf)

"""
    Alpha{IT<:Integer, VT<:Real}

Synapse that returns `(t - lastspike) * (q / τ) * exp(-(t - lastspike - τ) / τ) Θ(t - lastspike)`
(where `Θ` is the Heaviside function).
"""
mutable struct Alpha{IT<:Integer, VT<:Real} <: AbstractSynapse
    lastspike::VT
    q::VT
    τ::VT
end

"""
    Alpha{IT, VT}(;q::Real = 1, τ::Real = 1)
    Alpha(;q::Real = 1, τ::Real = 1)

Create an alpha synapse with amplitude `q` and time constant `τ`.
"""
Alpha{IT, VT}(;q::Real = 1, τ::Real = 1) where {IT<:Integer, VT<:Real} = Alpha{IT, VT}(-Inf, q, τ)
Alpha(;q::Real = 1, τ::Real = 1) = Alpha{Int, Float32}(q = q, τ = τ)

"""
    excite!(synapse::Alpha, spike::Integer)
    excite!(synapses::AbstractArray{<:Alpha}, spike::Integer)

Excite `synapse` with a `spike` (`spike` == time step of spike).
"""
excite!(synapse::Alpha, spike::Integer) = (spike > 0) && (synapse.lastspike = spike)
excite!(synapses::T, spike::Integer) where T<:AbstractArray{<:Alpha} = (spike > 0) && (synapses.lastspike .= spike)

isactive(synapse::Alpha, t::Real; dt::Real = 1.0) = dt * (t - synapse.lastspike) <= 10 * synapse.τ
isactive(synapses::T, t::Integer; dt::Real = 1.0) where T<:AbstractArray{<:Alpha} =
    any(dt .* (t .- synapses.lastspike) .<= 10 .* synapses.τ)

"""
    alpha(t::Real, lastspike, q, tau)
    alpha(t::Real, lastspike::AbstractArray{<:Real}, q::AbstractArray{<:Real}, tau::AbstractArray{<:Real})
    alpha(t::Real, lastspike::CuVecOrMat{<:Real}, q::CuVecOrMat{<:Real}, tau::CuVecOrMat{<:Real})

Evaluate an alpha synapse. Modeled as `(t - lastspike) * (q / τ) * exp(-(t - lastspike - τ) / τ) Θ(t - lastspike)`
  (where `Θ` is the Heaviside function).
Use `CuVector` instead of `Vector` for GPU support.

# Fields
- `t`: current time
- `lastspike`: last pre-synaptic spike time
- `q`: amplitude
- `tau`: time constant
"""
function alpha(t::Real, lastspike, q, tau)
    Δ = t - lastspike

    return (Δ >= 0 && Δ < Inf) * Δ * (q / tau) * exp(-(Δ - tau) / tau)
end
function alpha(t::Real, lastspike::AbstractArray{<:Real}, q::AbstractArray{<:Real}, tau::AbstractArray{<:Real})
    Δ = t .- lastspike
    I = @avx @. Δ * (q / tau) * exp(-(Δ - tau) / tau)

    return map((δ, i) -> (δ >= 0) && (δ < Inf) ? δ * i : zero(i), Δ, I)
end
function alpha(t::Real, lastspike::CuVecOrMat{<:Real}, q::CuVecOrMat{<:Real}, tau::CuVecOrMat{<:Real})
    Δ = t .- lastspike

    return @. (Δ >= 0) * (Δ < Inf) * Δ * (q / tau) * exp(-(Δ - tau) / tau)
end

"""
    evaluate!(synapse::Alpha, t::Integer; dt::Real = 1.0)
    (synapse::Alpha)(t::Integer; dt::Real = 1.0)
    evaluate!(synapses::AbstractArray{<:Alpha}, t::Integer; dt::Real = 1.0)

Evaluate an alpha synapse. See [`Synapse.Alpha`](@ref).
"""
evaluate!(synapse::Alpha, t::Integer; dt::Real = 1.0) = alpha(t * dt, synapse.lastspike * dt, synapse.q, synapse.τ)
(synapse::Alpha)(t::Integer; dt::Real = 1.0) = evaluate!(synapse, t; dt = dt)
evaluate!(synapses::T, t::Integer; dt::Real = 1.0) where T<:AbstractArray{<:Alpha} =
    alpha(t * dt, synapses.lastspike * dt, synapses.q, synapses.τ)

"""
    reset!(synapse::Alpha)
    reset!(synapses::AbstractArray{<:Alpha})

Reset `synapse`.
"""
reset!(synapse::Alpha) = (synapse.lastspike = -Inf)
reset!(synapses::T) where T<:AbstractArray{<:Alpha}= (synapses.lastspike .= -Inf)

"""
    EPSP{T<:Real}

Synapse that returns `(ϵ₀ / τm - τs) * (exp(-Δ / τm) - exp(-Δ / τs)) Θ(Δ)`
(where `Θ` is the Heaviside function and `Δ = t - lastspike`).

Specifically, this is the EPSP time course for the SRM0 model introduced by Gerstner.
Details: [Spiking Neuron Models: Single Neurons, Populations, Plasticity]
         (https://icwww.epfl.ch/~gerstner/SPNM/node27.html#SECTION02323400000000000000)
"""
mutable struct EPSP{IT<:Integer, VT<:Real} <: AbstractSynapse
    spikes::CircularBuffer{VT}
    ϵ₀::VT
    τm::VT
    τs::VT
end

"""
    EPSP{IT, VT}(;ϵ₀::Real = 1, τm::Real = 1, τs::Real = 1, N = 100)
    EPSP(;ϵ₀::Real = 1, τm::Real = 1, τs::Real = 1, N = 100)

Create an EPSP synapse with amplitude `ϵ₀`, rise time `τs`, and fall time `τm`.
Specify `N` to adjust how many pre-synaptic spikes are remembered between post-synaptic spikes.
"""
EPSP{IT, VT}(;ϵ₀::Real = 1, τm::Real = 1, τs::Real = 1, N = 100) where {IT<:Integer, VT<:Real} =
    EPSP{IT, VT}(fill!(CircularBuffer{VT}(N), -Inf), ϵ₀, τm, τs)
EPSP(;ϵ₀::Real = 1, τm::Real = 1, τs::Real = 1, N = 100) = EPSP{Int, Float32}(ϵ₀ = ϵ₀, τm = τm, τs = τs, N = N)

"""
    excite!(synapse::EPSP, spike::Integer)
    excite!(synapses::AbstractArray{<:EPSP}, spike::Integer)

Excite `synapse` with a `spike` (`spike` == time step of spike).
"""
excite!(synapse::EPSP, spike::Integer) = (spike > 0) && push!(synapse.spikes, spike)
excite!(synapses::T, spike::Integer) where T<:AbstractArray{<:EPSP} = (spike > 0) && push!.(synapses.spikes, spike)

"""
    spike!(synapse::EPSP, spike::Integer; dt::Real = 1.0)
    spike!(synapses::AbstractArray{<:EPSP}, spikes; dt::Real = 1.0)

Reset `synapse` when the post-synaptic neuron spikes.
"""
spike!(synapse::EPSP, spike::Integer; dt::Real = 1.0) = reset!(synapse)
spike!(synapses::T, spikes; dt::Real = 1.0) where T<:AbstractArray{<:EPSP} = reset!(synapses)

isactive(synapse::EPSP, t::Integer; dt::Real) = dt * (t - first(synapse.spikes)) <= synapse.τs + 8 * synapse.τm
isactive(synapses::T, t::Integer; dt::Real = 1.0) where T<:AbstractArray{<:EPSP} =
    any(dt .* (t .- first.(synapses.spikes)) .<= synapses.τs .+ 8 .* synapses.τm)

"""
    epsp(t::Real, lastspike, q, taum, taus)
    epsp(t::Real, lastspike::AbstractArray{<:Real}, q::AbstractArray{<:Real}, taum::AbstractArray{<:Real}, taus::AbstractArray{<:Real})
    epsp(t::Real, lastspike::CuVecOrMat{<:Real}, q::CuVecOrMat{<:Real}, taum::CuVecOrMat{<:Real}, taus::CuVecOrMat{<:Real})

Evaluate an EPSP synapse. Modeled as `(ϵ₀ / τm - τs) * (exp(-Δ / τm) - exp(-Δ / τs)) Θ(Δ)`
  (where `Θ` is the Heaviside function and `Δ = t - lastspike`).
Specifically, this is the EPSP time course for the SRM0 model introduced by Gerstner.
Details: [Spiking Neuron Models: Single Neurons, Populations, Plasticity]
         (https://icwww.epfl.ch/~gerstner/SPNM/node27.html#SECTION02323400000000000000)

Use `CuVector` instead of `Vector` for GPU support.

# Fields
- `t`: current time
- `lastspike`: last pre-synaptic spike time
- `q`: amplitude
- `taum`: rise time constant
- `taus`: fall time constant
"""
function epsp(t::Real, lastspike, q, taum, taus)
    Δ = t - lastspike

    return (Δ >= 0 && Δ < Inf && taus != taum) * q / (1 - taus / taum) * (exp(-Δ / taum) - exp(-Δ / taus))
end
function epsp(t::Real, lastspike::AbstractArray{<:Real}, q::AbstractArray{<:Real}, taum::AbstractArray{<:Real}, taus::AbstractArray{<:Real})
    Δ = t .- lastspike
    I = @. q / (1 - taus / taum) * (exp(-Δ / taum) - exp(-Δ / taus))

    return map((ts, tm, δ, i) -> (δ >= 0) && (δ < Inf) && (ts != tm) ? i : zero(i), taus, taum, Δ, I)
end
function epsp(t::Real, lastspike::CuVecOrMat{<:Real}, q::CuVecOrMat{<:Real}, taum::CuVecOrMat{<:Real}, taus::CuVecOrMat{<:Real})
    Δ = t .- lastspike

    return @. (Δ >= 0) * (Δ < Inf) * (taus != taum) * q / (1 - taus / taum) * (exp(-Δ / taum) - exp(-Δ / taus))
end

"""
    evaluate!(synapse::EPSP, t::Integer; dt::Real = 1.0)
    (synapse::EPSP)(t::Integer; dt::Real = 1.0)
    evaluate!(synapses::AbstractArray{<:EPSP}, t::Integer; dt::Real = 1.0)

Evaluate an EPSP synapse. See [`Synapse.EPSP`](@ref).
"""
evaluate!(synapse::EPSP, t::Integer; dt::Real = 1.0) =
    mapreduce(tf -> epsp(t * dt, tf * dt, synapse.ϵ₀, synapse.τm, synapse.τs), +, synapse.spikes)
(synapse::EPSP)(t::Integer; dt::Real = 1.0) = evaluate!(synapse, t; dt = dt)
function evaluate!(synapses::T, t::Integer; dt::Real = 1.0) where T<:AbstractArray{<:EPSP}
    N = length(synapses.spikes[1])
    return mapreduce(i -> epsp(t * dt, adapt(typeof(synapses.ϵ₀), getindex.(synapses.spikes, i) * dt), synapses.ϵ₀, synapses.τm, synapses.τs), +, 1:N)
end

"""
    reset!(synapse::EPSP)
    reset!(synapses::AbstractArray{<:EPSP})

Reset `synapse`.
"""
reset!(synapse::EPSP) = fill!(empty!(synapse.spikes), -Inf)
reset!(synapses::T) where T<:AbstractArray{<:EPSP}= fill!.(empty!.(synapses.spikes), -Inf)

end


text:
"""
What happens if you simulate a consumer-resource from the same
initial condition with different models of functional reponse?

Are the resulting trajectories from differnet functional repsones
similar or do they vary widely?
"""

params = ()

funcresps = [LotkaVolterra(α), HollingTypeII(α), HollingTypeII()

traetments =

funcresps = [LotkaVolterra(), HollingTypeII(), HollingTypeII()


text:
# # Tutorial 01 Rijke Tube

# ## Introduction
#
# A Rijke tube is the most simple thermo-acoustic configuration. It comprises a
# tube with an unsteady heat source somewhere inside the tube. This example will
# walk you through the basic steps of setting up a Rijke tube in a
# Helmholtz-solver stability analysis.

# ### The Helmholtz equation
# The Helmholtz equation is the Fourier transform of the acoustic wave equation.
# Given that there is a heat source, it reads:
#
# ∇⋅(c² ∇ p̂) + ω² p̂ = -iω(γ-1)q̂
#
# Here c is speed-of-sound-field, p̂ the (complex) pressure fluctuation amplitude
# ω the (complex) frequency of the problem, i the imaginary unit, γ the ratio of
# specifc heats, and q̂ the amplitude of the heat release fluctuation.
#
# (Note that the Fourier transform here follows a f'(t) --> f̂(ω)exp(+iωt) convention.)
#
# Together with some boundary conditions the Helmholtz equation models
# thermo-acoustic resonance in a cavity. The minimum required inputs to specify
# a problem are
#
# 1. the shape of the domain
# 2. the speed-of-sound field
# 3. the boundary conditions
#
# In case of an active flame (q̂≠0). We will also need
#
# 4. some additional gas porperties and
# 5. a flame response model linking the heat release fluctuations q̂ to the pressure fluctuations p̂
#
# Once you completed this tutorial you will know the basic steps of how to
# conduct a thermo-acoustic stability analysis

## #jl

# ## Header
# First you will need to load the Helmholtz solver. The following line brings all
# necessary tools into scope:
using WavesAndEigenvalues.Helmholtz

## #jl
# ## Mesh
# Now we can load the mesh file. It is the essential piece of information
# defining the domain shape . This tutorial's mesh has been specified in mm
# which is why we scale it by `0.001` to load the coordinates as m:
mesh=Mesh("Rijke_mm.msh",scale=0.001)

# You can have a look at some basic mesh data by printing it
print(mesh)
# In an interactive session, is suffices to just type the mesh's variable name
# without the enclosing `print` command.
#
#
# This info tells us that the mesh features `1006` points as vertices
# `1562` (addressable) triangles on its surface and `3380`tetrahedra forming
# the tube. Note, that there are currently no lines stored. This is because
# line information is only needed when using second-order finite elements. To
# save memory this information is only assembled and stored when explicitly
# requested. We will come back to this aspect in a later tutorial.
#
# You can also see that the mesh features several named domains such as
# `"Interior"`,`"Flame"`, `"Inlet"` and `"Outlet"`. These domains are used to
# specify certain conditions for our model.
#
# A model descriptor is always given as a dictionairy featureing domain names
# as keys and tuples as values. The first tuple entry is a Julia symbol
# specifying the operator to be build on the associated domain, while the second
# entry is again a tuple holding information that is specific to the chosen
# operator.

## #src
# ## Model set-up
# For instance, we certainly want to build the wave operator on the entire mesh.
# So first, we initialize an empty dictionairy
dscrp=Dict()

# And then specify that the wave operator should be build everywhere. For the
# given mesh the domain-specifier to address the entire resononant cavity is
# `"Interior"` and in general the symbol to create the wave operator is
# `:interior`. It requires no options so the options tuple remains empty.
dscrp["Interior"]=(:interior, ())

## src
# ## Boundary Conditions
# The tube should have a pressure node at its outlet, in order to model an
# open-end boundary condition. Boundary conditions are specified in terms of
# admittance values. A pressure note corresponds to an infinite admittance.
# To represent this on a computer we just give it a crazily high value like
# `1E15`. We will also need to specify a variable name for our admmittance value.
# This will allow to quickly change the value after discretization of the
# model by addressing it by this very name. This feature is one of the core
# concepts of the solver. As will be demonstrated later.
#
# The complete description of our boundary condition reads
dscrp["Outlet"]=(:admittance, (:Y,1E15))

# You may wonder why we do not specify conditions for the other boundaries of
# the model. This is because the discretization is based on Bubnov-Galerkin
# finite elements. All unspecified, boundaries will therefore *naturally* be
# discretized as solid walls.

## #src
# ## Flame
# The Rijke tube's main feauture is a domain of heat release. In the current
# mesh there is a designated domain `"Flame"` addressing a thin volume at the
# center of the tube. We can use this key to specify a flame with simple
# n-tau-dynamics. Therefore, we first need to define some basic gas properties.

γ=1.4 #ratio of specific heats
ρ=1.225 #density at the reference location upstream to the flame in kg/m^3
Tu=300.0    #K unburnt gas temperature
Tb=1200.0    #K burnt gas temperature
P0=101325.0 # ambient pressure in Pa
A=pi*0.025^2 # cross sectional area of the tube
Q02U0=P0*(Tb/Tu-1)*A*γ/(γ-1) #the ratio of mean heat release to mean velocity Q02U0

# We will also need to provide the position where the reference velocity has
# been taken and the direction of that velocity
x_ref=[0.0; 0.0; -0.00101]
n_ref=[0.0; 0.0; 1.00] #must be a unit vector
# And of course, we need some values for n and tau
n=0.0 #interaction index
τ=0.001 #time delay
# With these values the specification of the flame reads
dscrp["Flame"]=(:flame,(γ,ρ,Q02U0,x_ref,n_ref,:n,:τ,n,τ))
# Note that the order of the values in the options tuple is important. Also note
# that we assign the symbols `:n`and `:τ` for later analysis.

## #jl
# ## Speed of Sound
# The description of the Rijke tube model is nearly complete. We just need to
# specify the speed of sound field.  For this example, the field is fairly
# simple and can be specified analytically, using the `generate_field` function
# and a custom three-parameter function.
R=287.05 # J/(kg*K) specific gas constant (air)
function speedofsound(x,y,z)
    if z<0.
        return sqrt(γ*R*Tu)#m/s
    else
        return sqrt(γ*R*Tb)#m/s
    end
end
c=generate_field(mesh,speedofsound)
# Note that in more elaborate models, you may read the field `c` from a file
# containing simulation or experimental data, rather than specifying it
# analytically.

## #src
# ## Model Discretization
# Now we have all ingredients together and we can discretize the model.
L=discretize(mesh,dscrp,c)
# The return value `L` here is a family of linear operators. You can display
# an algebraic representation of the family plus a list of associated parameters
# by just printing it
print(L)
# (In an interactive session the enclosing print function is not necessary.)
#
# You might notice that the list contains all parameters that have been
# specified during the model description (`n`,`τ`,`Y`). However, there are also
# two parameters that were added by default: `ω` and `λ`. `ω` is the complex
# frequency of the model and `λ` an auxiliary value that is important for the
# eigenfrequency computation.
## #src
# ## Solving the Model
# ### Global Solver
# We can solve the model for some eigenvalues using one of the eigenvalue
# solvers. The package provides you with two types of eigenvalue solvers. A global
# contour-integration-based solver. That finds you all eigenvalues inside of a
# specified contour Γ in the complex plane.
Γ=[150.0+5.0im, 150.0-5.0im, 1000.0-5.0im, 1000.0+5.0im].*2*pi #corner points for the contour (in this case a rectangle)
Ω, P = beyn(L,Γ,l=10,N=64, output=true)

# The found eigenvalues are stored in the array `Ω`. The corresponding
# eigenvectors are the columns of `P`.
# The huge advantage of the global eigenvalue solver is that it finds you
# multiple eigenvalues. Nonetheless, its accuracy may be low and sometimes it
# provides you with outright wrong solutions.
# However, for the current case the method works as we can varify that in our
# search window there are two eigenmodes oscilating at 272 and 695 Hz
# respectively:
for ω in Ω
    println("Frequency: $(ω/2/pi)")
end
## #src
# ### Local Solver
# To get high accuracy eigenvalues , there is also local iterative eigenvalue
# solver. Based on an initial guess, it only finds you one eigenvalue at a time
# but with machine-precise accuracy.
sol,nn,flag=householder(L,250*2*pi,output=true);

# The return values of this solver are a solution object `sol`, the number of
# iterations `nn` performed by the solver and an error flag `flag` providing
# information on the quality of the solution. If `flag>0` the solution has
# converged. If `flag<0`something went wrong. And if `flag==0`... well, probably
# the solution is fine but you should check it.
#
# In the current case the flag is -1 indicating that the maximum number of
# iterations has been reached. This is because we haven't specified a stopping
# criterion and the iteration just ran for the maximum number of iterations.
# Nevertheless, the 272 Hz mode is correct. Indeed, as you can see from the
# printed output it already converged to machine-precision after 5 iterations.

## #src
# ## Changing a specified parameter.
#
# Remember that we have specified the parameters `Y`,`n`, and `τ`?. We can change
# these easily without rediscretizing our model. For instance currently `n==0.0`
# so the solution is purely accoustic. The effect of the flame just enters the
# problem via the speed of sound field (this is known as passive flame). By
# setting the interaction index to `1.0` we activate the flame.
L.params[:n]=1
# Now we can just solve the model again to get the active solutions
sol_actv,nn,flag=householder(L,245*2*pi-82im*2*pi,output=true,order=3);
# The eigenfrequency is now complex:
sol_actv.params[:ω]/2/pi
# with a growth rate `-imag(sol_actv.params[:ω]/2/pi)≈  59.22`

# Instead of recomputing the eigenvalue by one of the solvers. We can also
# approximate the eigenvalue as a function of one of the model parameters
# utilizing high-order perturbation theory. For instance this gives you
# the 30th order diagonal Padé estimate expanded from the passive solution.

perturb_fast!(sol,L,:n,16) #compute the coefficients
freq(n)=sol(:n,n,8,8)/2/pi #create some function for convenience
freq(1) #evaluate the function at any value for `n` you are interested in.

# The method is slow when only computing one eigenvalue. But its computational
# costs are almost constant when computing multiple eigenvalues. Therefore,
# for larger parametric studies this is clearly the method of choice.

## #src
# ## VTK Output for Paraview
#
# To visualizte your solutions you can store them as `"*.vtu"` files to open
# them with paraview. Just, create a dictionairy, where your modes are stored as
# fields.
data=Dict()
data["speed_of_sound"]=c
data["abs"]=abs.(sol_actv.v)/maximum(abs.(sol_actv.v)) #normalize so that max=1
data["phase"]=angle.(sol_actv.v)

vtk_write("tutorial_01", mesh, data) # Write the solution to paraview

# The `vtk_write` function automatically sorts your data by its interpolation
# scheme. Data that is constant on a tetrahedron will be written to a file
# `*_const.vtu`, data that is linearly interpolated on a tetrahedron will go
# to `*_lin.vtu`, and data that is quadratically interpolated to `*_quad.vtu`.
# The current example uses constant speed of sound on a tetrahedron and linear
# finite elements for the discretization of p. Therefore, two files are
# generated, namely `tutorial_01_const.vtu` containing the speed-of-sound field
# and `tutorial_01_lin.vtu` containing the mode shape. Open them with paraview
# and have a look!.

## src
# ## Summary
#
# You learnt how to set-up a simple Rijke-tube study. This  already introduced
# all basic steps that are needed to work with the Helmholtz solver. However,
# there are a lot of details you can fine tune and even features that weren't
# mentioned in this tutorial. The next tutorials will introduce these aspects in
# more detail.


text:
using DiffEqSensitivity, SafeTestsets
using Test

const GROUP = get(ENV, "GROUP", "All")
const is_APPVEYOR = Sys.iswindows() && haskey(ENV,"APPVEYOR")
const is_TRAVIS = haskey(ENV,"TRAVIS")

@time begin
if GROUP == "All" || GROUP == "Core" || GROUP == "Downstream"
    @safetestset "Local Sensitivity" begin include("local.jl") end
    @safetestset "Adjoint Sensitivity" begin include("adjoint.jl") end
    @safetestset "Morris Method" begin include("morris_method.jl") end
    @safetestset "Sobol Method" begin include("sobol_method.jl") end
    @safetestset "DGSM Method" begin include("DGSM.jl") end
end

if GROUP == "DiffEqFlux"
    using Pkg
    if is_TRAVIS
      using Pkg
      Pkg.add("DiffEqFlux")
    end
    Pkg.test("DiffEqFlux")
end
end


text:
# publish hists over ZMQ
function make_coadded_hist(masschannels)
	coadded_hist = Histogram(0:1:20000)
	n = 0
	for (ch,mc) in masschannels
		if haskey(mc,:energy_hist)
			coadded_hist+=mc[:energy_hist]
			n+=1
		end
	end
	coadded_hist,n
end

function publish_hists_zmq(masschannels, zmqexitchannel)
	context=Context()
	pub=Socket(context, PUB)
	ZMQ.bind(pub, "tcp://*:5555")
	try
		while !isready(zmqexitchannel)
			coadded_hist,n = make_coadded_hist(masschannels)
			io = IOBuffer()
			write(io,"0")
			serialize(io, (coadded_hist,n))
			ZMQ.send(pub, Message(io))
			sleep(1)
		end
	catch ex
		isa(ex, InterruptException) && return
		println(ex)
	finally
		ZMQ.close(pub)
		ZMQ.close(context)
	end
end

zmqexitchannel = Channel()
@schedule publish_hists_zmq(masschannels, zmqexitchannel)


text:
# visualization service calls (based on previous similar Collections and Director visualizations)

"""
    $(SIGNATURES)

Get all the bearing range constraints and measurement values from in memory object (has Graff twin).

Usage: ReprojectBearingRange plugin.
"""
function getBearingRangesOnPose(fgl::FactorGraph, psym::Symbol)
    tags = Dict()

    error("getBearingRangesOnPose not implemented yet")

    return tags
end

#


text:
"""
    Attribute

An `Attribute` represents a lazily evaluated field in a log `Record`.

# Fields
* `f::Function`: A function to evaluate in order to get a value if one is not set.
* `x::Union{Some{T}, Nothing}`: A value that may or may not exist yet.
"""
mutable struct Attribute{T}
    f::Function
    x::Union{Some{T}, Nothing}
end

"""
    Attribute(f::Function) -> Attribute{Any}
    Attribute{T}(f::Function) -> Attribute{T}

Creates an `Attribute` with the given function.
"""
Attribute{T}(f::Function) where {T} = Attribute{T}(f, nothing)
Attribute(f::Function) = Attribute{Any}(f, nothing)

"""
    Attribute(x)

Simply wraps the value `x` in `Some` and sticks that in an `Attribute` with an
empty `Function`.
"""
Attribute(x::T) where {T} = Attribute{T}(() -> x)
Attribute{T}(x::T) where {T} = Attribute{T}(() -> x)

"""
    get(attr::Attribute{T}) -> T

Run set `attr.x` to the output of `attr.f` if `attr.x` is not already set.
We then return the value stored in `attr.x`
"""
function Base.get(attr::Attribute{T}) where T
    if attr.x === nothing
        attr.x = Some{T}(attr.f())
    end

    return something(attr.x)::T
end

if VERSION < v"1.2.0-DEV.272"
    hasfield(T::DataType, name::Symbol) = Base.fieldindex(T, name, false) > 0
end

"""
    Record

A dictionary-like container with `Symbol` keys used to store information about a log events
including the msg, date, level, stacktrace, etc. `Formatter`s use `Records` to format log
message strings.

You can access the properties of a `Record` by using `getproperty` (ie: record.msg).

Subtypes of `Record` should implement `getproperty(::MyRecord, ::Symbol)` and key-value pair
iteration.
"""
abstract type Record <: AbstractDict{Symbol, Any} end

"""
    getlevel(::Record) -> AbstractString

Returns the record level.
"""
getlevel(rec::Record) = rec.level

"""
    AttributeRecord <: Record

A `Record` which stores its properties as `Attribute`s for lazy evaluation.

Calling `getproperty` or iterating will evaluate and cache the properties accessed.
"""
abstract type AttributeRecord <: Record end

@deprecate getindex(rec::T, attr::Symbol) where {T <: Record} getproperty(rec, attr)
Base.haskey(rec::T, attr::Symbol) where {T <: Record} = hasfield(T, attr)
Base.keys(rec::T) where {T <: Record} = (fieldname(T, i) for i in 1:fieldcount(T))

function Base.iterate(rec::T, state=0) where T <: Record
    state >= fieldcount(T) && return nothing
    state += 1
    return (fieldname(T, state) => getfield(rec, state), state)
end

function Base.iterate(rec::T, state=0) where T <: AttributeRecord
    state >= fieldcount(T) && return nothing
    state += 1
    return (fieldname(T, state) => get(getfield(rec, state)), state)
end

Base.getproperty(rec::AttributeRecord, attr::Symbol) = get(getfield(rec, attr))

"""
    Dict(rec::Record)

Extracts the `Record` and its properties into a `Dict`

!!! warn

    On `AttributeRecord`s this may be an expensive operation, so you probably don't want to
    do this for every log record unless you're planning on using every `Attribute`.
"""
Base.Dict(::Record)

"""
    DefaultRecord <: AttributeRecord

Stores the most common logging event information.
NOTE: if you'd like more logging attributes you can:

1. add them to DefaultRecord and open a pull request if the new attributes are applicable to most applications.
2. make a custom `Record` type.

# Fields
* `date::Attribute{ZonedDateTime}`: timestamp of log event
* `level::Attribute{AbstractString}`: log level
* `levelnum::Attribute{Int}`: integer value for log level
* `msg::Attribute{AbstractString}`: the log message itself
* `name::Attribute{AbstractString}`: the name of the source logger
* `pid::Attribute{Int}`: the pid of where the log event occured
* `lookup::Attribute{StackFrame}`: the top StackFrame
* `stacktrace::Attribute{StackTrace}`: a stacktrace
"""
struct DefaultRecord <: AttributeRecord
    date::Attribute
    level::Attribute
    levelnum::Attribute
    msg::Attribute
    name::Attribute
    pid::Attribute
    lookup::Attribute
    stacktrace::Attribute
end

"""
    DefaultRecord(name::AbstractString, level::AbstractString, levelnum::Int, msg)

Takes a few initial log record arguments and creates a `DefaultRecord`.

# Arguments
* `name::AbstractString`: the name of the source logger.
* `level::AbstractString`: the log level.
* `levelnum`::Int: the log level number.
* `msg`: the message being logged.
"""
function DefaultRecord(name::AbstractString, level::AbstractString, levelnum::Int, msg)
    trace = Attribute{StackTrace}(get_trace)

    DefaultRecord(
        Attribute{ZonedDateTime}(() -> Dates.now(tz"UTC")),
        Attribute(level),
        Attribute(levelnum),
        Attribute{AbstractString}(msg),
        Attribute(name),
        Attribute(getpid()),
        Attribute{Union{StackFrame, Nothing}}(get_lookup(trace)),
        trace,
    )
end

"""
    get_trace()

Returns the `StackTrace` with `StackFrame`s from the `Memento` module filtered out.
"""
function get_trace()
    trace = StackTraces.stacktrace()
    return filter!(trace) do frame
        !Base.StackTraces.from(frame, Memento)
    end
end

"""
    get_lookup(trace::Attribute{StackTrace})

Returns the top `StackFrame` for `trace` if it isn't empty.
"""
function get_lookup(trace::Attribute{StackTrace})
    inner() = isempty(get(trace)) ? nothing : first(get(trace))
    return inner
end