text:
module LambertW

export
    lambertw

end # module

@doc """
```
lambertw(x[, k, prec])
```

The `k`'th branch of Lambert's W function evaluated at `x`.

- `k` must be `0` (default) or `-1`.
- `prec` defaults to `eps()`.
"""->
function lambertw(x::Real, k::Int=0, prec=eps())
	@assert k == 0 || k == -1 "The branch k must be either 0 or -1"
	@assert prec >= eps()

	# Check argument
	const minx = -1/e
	if x < minx
		return NaN
	end

	if k == -1 && x >= 0
		return NaN
	end

	# First approximation
	W = lambertw(x, Val{k})

	# Compute residual using logarithms to avoid numerical overflow
	# When x == 0, r = NaN, but here the approximation is exact and
	# the while loop below is unnecessary
	r = abs( W - log(abs(x)) + log(abs(W)) )

	# Apply Fritsch’s method to increase precision
	n = 1
	while r > prec && n <= 5
		z = log(x/W) - W
		q = 2*(1 + W)*(1 + W + 2/3*z)
		epsilon = z*(q - z) / ((1 + W)*(q - 2*z))
		W *= 1 + epsilon

		r = abs( W - log(abs(x)) + log(abs(W)) )
		n += 1
	end

	return W
end

function lambertw{T<:Real}(x::Array{T}, k::Int=0, prec=eps())
	W = Array(Float64, size(x))
	for idx in eachindex(x)
		W[idx] = lambertw( x[idx], k, prec )
	end
	return W
end

# ------------------------------------------------------------
# Initial approximations

# For branch 0
function lambertw(x::Real, ::Type{Val{0}})
	if x <= 1
		const sqrt2 = sqrt(2)
		sqeta = sqrt( 2 + 2*e*x )

		N2 = 3*sqrt2 + 6 - (((2237+1457*sqrt2)*e - 4108*sqrt2 - 5764)*sqeta)/((215+199*sqrt2)*e - 430*sqrt2-796);
		N1 = (1-1/sqrt2)*(N2+sqrt2);

		W = -1 + sqeta/(1 + N1*sqeta/(N2 + sqeta));
	else
		W = log( 6*x/(5*log( 12/5*(x/log(1+12*x/5)) )) )
	end

	return W
end

# For branch -1
function lambertw(x::Real, ::Type{Val{-1}})
	const M1 = 0.3361
	const M2 = -0.0042
	const M3 = -0.0201

	sigma = -1 - log(-x)
	W = -1 - sigma - 2/M1*(1 - 1/(1 + (M1*sqrt(sigma/2)) / (1 + M2*sigma*exp(M3*sqrt(sigma)))))

	return W
end


text:
# Autogenerated wrapper script for GStreamer_jll for i686-w64-mingw32
export gst_inspect, gst_launch, gst_stats, gst_tester, gst_typefind, libgstbase, libgstcheck, libgstcontroller, libgstnet, libgstreamer

using Glib_jll
using LibUnwind_jll
using Elfutils_jll
using GMP_jll
using GSL_jll
using libcap_jll
JLLWrappers.@generate_wrapper_header("GStreamer")
JLLWrappers.@declare_executable_product(gst_inspect)
JLLWrappers.@declare_executable_product(gst_launch)
JLLWrappers.@declare_executable_product(gst_stats)
JLLWrappers.@declare_executable_product(gst_tester)
JLLWrappers.@declare_executable_product(gst_typefind)
JLLWrappers.@declare_library_product(libgstbase, "libgstbase-1.0-0.dll")
JLLWrappers.@declare_library_product(libgstcheck, "libgstcheck-1.0-0.dll")
JLLWrappers.@declare_library_product(libgstcontroller, "libgstcontroller-1.0-0.dll")
JLLWrappers.@declare_library_product(libgstnet, "libgstnet-1.0-0.dll")
JLLWrappers.@declare_library_product(libgstreamer, "libgstreamer-1.0-0.dll")
function __init__()
    JLLWrappers.@generate_init_header(Glib_jll, LibUnwind_jll, Elfutils_jll, GMP_jll, GSL_jll, libcap_jll)
    JLLWrappers.@init_executable_product(
        gst_inspect,
        "bin\\gst-inspect-1.0.exe",
    )

    JLLWrappers.@init_executable_product(
        gst_launch,
        "bin\\gst-launch-1.0.exe",
    )

    JLLWrappers.@init_executable_product(
        gst_stats,
        "bin\\gst-stats-1.0.exe",
    )

    JLLWrappers.@init_executable_product(
        gst_tester,
        "bin\\gst-tester-1.0.exe",
    )

    JLLWrappers.@init_executable_product(
        gst_typefind,
        "bin\\gst-typefind-1.0.exe",
    )

    JLLWrappers.@init_library_product(
        libgstbase,
        "bin\\libgstbase-1.0-0.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgstcheck,
        "bin\\libgstcheck-1.0-0.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgstcontroller,
        "bin\\libgstcontroller-1.0-0.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgstnet,
        "bin\\libgstnet-1.0-0.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgstreamer,
        "bin\\libgstreamer-1.0-0.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()

ext: jl


text:
using JuLIP
using JuLIP.Potentials
using JuLIP.Solve
using JuLIP.Constraints
import JuLIP.Preconditioners: Exp

println("===================================================")
println("          TEST SOLVE ")
println("===================================================")

println("-----------------------------------------------------------------")
println("Testing `minimise!` with equilibration with LJ calculator to lattice")
println("-----------------------------------------------------------------")
calc = lennardjones(r0=JuLIP.ASE.rnn("Al"))
at = Atoms("Al", cubic=true, repeatcell=(3,3,3), pbc=(true,true,true))
X0 = positions(at) |> mat
at = rattle!(at, 0.02)
set_calculator!(at, calc)
set_constraint!(at, FixedCell(at))
minimise!(at)
X1 = positions(at) |> mat
X0 .-= X0[:, 1]
X1 .-= X1[:, 1]
F = X1 / X0
println("check that the optimiser really converged to a lattice")
@show vecnorm(F'*F - eye(3), Inf)
@show vecnorm(F*X0 - X1, Inf)

println("-------------------------------------------------")
println("same test but large and with Exp preconditioner")
println("-------------------------------------------------")

at = Atoms("Al", repeatcell=(20,20,2), pbc=(true,true,true), cubic=true)
at = rattle!(at, 0.02)
set_calculator!(at, calc)
set_constraint!(at, FixedCell(at))
minimise!(at, precond = Exp(at))


text:
handle("completions") do data
  @destruct [path || nothing,
             mod || "Main",
             editorContent || "",
             lineNumber || 1,
             startLine || 0,
             column || 1,
             line, force] = data

  withpath(path) do
    m = getmodule(mod)

    cs, pre = basecompletionadapter(line, m, force, lineNumber - startLine, column, editorContent)

    Dict(:completions => cs,
         :prefix      => string(pre))
  end
end

using REPL.REPLCompletions

function basecompletionadapter(line, mod, force, lineNumber, column, text)
  comps, replace, shouldcomplete = try
    completions(line, lastindex(line), mod)
  catch err
    # might error when e.g. type inference fails
    [], 1:0, false
  end

  # Suppress completions if there are too many of them unless activated manually
  # @TODO: Checking whether `line` is a valid text to be completed in atom-julia-client
  #        in advance and drop this check
  if !force && length(comps) > MAX_COMPLETIONS
    comps = []
    replace = 1:0
  end

  pre = line[replace]
  d = []
  for c in comps
    if REPLCompletions.afterusing(line, Int(first(replace))) # need `Int` for correct dispatch on x86
      c isa REPLCompletions.PackageCompletion || continue
    end
    try
      push!(d, completion(mod, c))
    catch err
      continue
    end
  end

  # completions from the local code block:
  for c in localcompletions(text, lineNumber, column)
    if (force || !isempty(pre)) && startswith(c[:text], pre)
      pushfirst!(d, c)
    end
  end

  d, pre
end

const MAX_COMPLETIONS = 500

function completion(mod, c)
  return Dict(:type               => completiontype(c),
              :icon               => completionicon(c),
              :rightLabel         => completionmodule(mod, c),
              :leftLabel          => completionreturntype(c),
              :text               => completiontext(c),
              :description        => completionsummary(mod, c),
              :descriptionMoreURL => completionurl(c))
end

completiontext(c) = completion_text(c)
completiontext(c::REPLCompletions.MethodCompletion) = begin
  ct = completion_text(c)
  m = match(r"^(.*) in .*$", ct)
  m isa Nothing ? ct : m[1]
end
completiontext(c::REPLCompletions.DictCompletion) = rstrip(completion_text(c), [']', '"'])
completiontext(c::REPLCompletions.PathCompletion) = rstrip(completion_text(c), '"')

using JuliaInterpreter: sparam_syms

completionreturntype(c) = ""
completionreturntype(c::REPLCompletions.MethodCompletion) = begin
  m = c.method
  atypes = m.sig
  sparams = Core.svec(sparam_syms(m)...)
  wa = Core.Compiler.Params(typemax(UInt))  # world age
  inf = try
    Core.Compiler.typeinf_type(m, atypes, sparams, wa)
  catch err
    nothing
  end
  inf in (nothing, Any, Union{}) && return ""
  shortstr(inf)
end
completionreturntype(c::REPLCompletions.PropertyCompletion) =
  shortstr(typeof(getproperty(c.value, c.property)))
completionreturntype(c::REPLCompletions.FieldCompletion) =
  shortstr(fieldtype(c.typ, c.field))
completionreturntype(c::REPLCompletions.DictCompletion) =
  shortstr(valtype(c.dict))
completionreturntype(::REPLCompletions.PathCompletion) = "Path"

using Base.Docs

completionsummary(mod, c) = ""
completionsummary(mod, c::REPLCompletions.ModuleCompletion) = begin
  m, word = c.parent, c.mod
  cangetdocs(m, word) || return ""
  docs = getdocs(m, word, mod)
  description(docs)
end
completionsummary(mod, c::REPLCompletions.MethodCompletion) = begin
  ct = Symbol(c.func)
  cangetdocs(mod, ct) || return ""
  docs = try
    Docs.doc(Docs.Binding(mod, ct), Base.tuple_type_tail(c.method.sig))
  catch err
    ""
  end
  description(docs)
end
completionsummary(mod, c::REPLCompletions.KeywordCompletion) =
  description(getdocs(mod, c.keyword))

using Markdown

description(docs) = ""
description(docs::Markdown.MD) = begin
  md = CodeTools.flatten(docs).content
  for part in md
    if part isa Markdown.Paragraph
      desc = Markdown.plain(part)
      occursin("No documentation found.", desc) && return ""
      return strlimit(desc, 200)
    end
  end
  return ""
end

completionurl(c) = ""
completionurl(c::REPLCompletions.ModuleCompletion) = begin
  mod, name = c.parent, c.mod
  val = getfield′(mod, name)
  if val isa Module # module info
    urimoduleinfo(parentmodule(val) == val || val ∈ (Base, Core) ? name : "$mod.$name")
  else
    uridocs(mod, name)
  end
end
completionurl(c::REPLCompletions.MethodCompletion) = uridocs(c.method.module, c.method.name)
completionurl(c::REPLCompletions.PackageCompletion) = urimoduleinfo(c.package)
completionurl(c::REPLCompletions.KeywordCompletion) = uridocs("Main", c.keyword)

completionmodule(mod, c) = shortstr(mod)
completionmodule(mod, c::REPLCompletions.ModuleCompletion) = shortstr(c.parent)
completionmodule(mod, c::REPLCompletions.MethodCompletion) = shortstr(c.method.module)
completionmodule(mod, c::REPLCompletions.FieldCompletion) = shortstr(c.typ) # predicted type
completionmodule(mod, ::REPLCompletions.KeywordCompletion) = ""
completionmodule(mod, ::REPLCompletions.PathCompletion) = ""

completiontype(c) = "variable"
completiontype(c::REPLCompletions.ModuleCompletion) = begin
  ct = completion_text(c)
  ismacro(ct) && return "snippet"
  ct == "Vararg" && return ""
  mod, name = c.parent, Symbol(ct)
  val = getfield′(mod, name)
  wstype(mod, name, val)
end
completiontype(::REPLCompletions.MethodCompletion) = "method"
completiontype(::REPLCompletions.PackageCompletion) = "import"
completiontype(::REPLCompletions.PropertyCompletion) = "property"
completiontype(::REPLCompletions.FieldCompletion) = "property"
completiontype(::REPLCompletions.DictCompletion) = "property"
completiontype(::REPLCompletions.KeywordCompletion) = "keyword"
completiontype(::REPLCompletions.PathCompletion) = "path"

completionicon(c) = ""
completionicon(c::REPLCompletions.ModuleCompletion) = begin
  ismacro(c.mod) && return "icon-mention"
  mod, name = c.parent, Symbol(c.mod)
  val = getfield′(mod, name)
  wsicon(mod, name, val)
end
completionicon(::REPLCompletions.DictCompletion) = "icon-key"
completionicon(::REPLCompletions.PathCompletion) = "icon-file"

function localcompletions(text, line, col)
  ls = locals(text, line, col)
  reverse!(ls)
  map(localcompletion, ls)
end

function localcompletion(l)
  return Dict(
    :type        => l[:type] == "variable" ? "attribute" : l[:type],
    :icon        => l[:icon] == "v" ? "icon-chevron-right" : l[:icon],
    :rightLabel  => l[:root],
    :leftLabel   => "",
    :text        => l[:name],
    :description => ""
  )
end

handle("cacheCompletions") do mod
  # m = getthing(mod)
  # m = isa(m, Module) ? m : Main
  # CodeTools.completions(m)
end


text:
mutable struct DEStats
  nf::Int
  nw::Int
  nsolve::Int
  naccept::Int
  nreject::Int
  maxeig::Float64
end

DEStats() = DEStats(ntuple(_->0, 5)..., 0.0)

function Base.show(io::IO, s::DEStats)
  println(io, summary(s))
  println(io, "Number of function evaluations: $(s.nf)")
  println(io, "Number of W matrix evaluations: $(s.nw)")
  println(io, "Number of linear solves:        $(s.nsolve)")
  println(io, "Number of accepted steps:       $(s.naccept)")
  print(io,   "Number of rejected steps:       $(s.nreject)")
  iszero(s.maxeig) || print(io, "\nMaximum eigenvalue recorded:    $(s.maxeig)")
end


text:
module LinearFitXYerrors

using Statistics
using Distributions
using Printf
using Plots; gr()    # TODO: define Plot recipe instead

include("linearfitxy.jl")

export linearfitxy

end


text:
const OccipitalPole = BilateralStructure("occipital pole")

const PreoccipitalNotch = BilateralStructure("preoccipital notch")

const SuperiorOccipitalGyrus = BilateralStructure("superior occipital gyrus")

const MiddleOccipitalGyrus = BilateralStructure("middle occipital gyrus")

const InferiorOccipitalGyrus = BilateralStructure("inferior occipital gyrus")

const LingualGyrus = BilateralStructure("lingual gyrus")

# TODO also in temporal lobe
#const FUSIFORM GYRUS = BilateralStructure(:fusiform gyrus) # (Brodmann area 37)
#parent(::BilateralStructure{:fusiform gyrus}) = OccipitalLobe

const LateralOccipitalCortex = BilateralStructure("lateral occipital cortex")

const CalcarineCortex = BilateralStructure("calcarine cortex") # visual cortex

const Cuneus = BilateralStructure("cuneus")

const OccipitalLobe = BilateralStructure("occipital lobe")


text:
module DiffEqFlux

using Adapt, Base.Iterators, ConsoleProgressMonitor, DataInterpolations,
    DiffEqBase, DiffEqSensitivity, DiffResults, Distributions, DistributionsAD,
    ForwardDiff, Optimization, OptimizationPolyalgorithms, LinearAlgebra,
    Logging, LoggingExtras, Printf, ProgressLogging, Random, RecursiveArrayTools,
    Reexport, SciMLBase, StaticArrays, TerminalLoggers, Zygote, ZygoteRules

@reexport using DiffEqSensitivity
@reexport using Zygote

# deprecate

import OptimizationFlux
import NNlib
import Lux
using Requires
using Cassette
@reexport using Flux
@reexport using OptimizationOptimJL

gpu_or_cpu(x) = Array

# ForwardDiff integration

ZygoteRules.@adjoint function ForwardDiff.Dual{T}(x, ẋ::Tuple) where T
  @assert length(ẋ) == 1
  ForwardDiff.Dual{T}(x, ẋ), ḋ -> (ḋ.partials[1], (ḋ.value,))
end

ZygoteRules.@adjoint ZygoteRules.literal_getproperty(d::ForwardDiff.Dual{T}, ::Val{:partials}) where T =
  d.partials, ṗ -> (ForwardDiff.Dual{T}(ṗ[1], 0),)

ZygoteRules.@adjoint ZygoteRules.literal_getproperty(d::ForwardDiff.Dual{T}, ::Val{:value}) where T =
  d.value, ẋ -> (ForwardDiff.Dual{T}(0, ẋ),)

ZygoteRules.@adjoint ZygoteRules.literal_getproperty(A::Tridiagonal, ::Val{:dl}) = A.dl,y -> Tridiagonal(dl,zeros(length(d)),zeros(length(du)),)
ZygoteRules.@adjoint ZygoteRules.literal_getproperty(A::Tridiagonal, ::Val{:d}) = A.d,y -> Tridiagonal(zeros(length(dl)),d,zeros(length(du)),)
ZygoteRules.@adjoint ZygoteRules.literal_getproperty(A::Tridiagonal, ::Val{:du}) = A.dl,y -> Tridiagonal(zeros(length(dl)),zeros(length(d),du),)
ZygoteRules.@adjoint Tridiagonal(dl, d, du) = Tridiagonal(dl, d, du), p̄ -> (diag(p̄[2:end,1:end-1]),diag(p̄),diag(p̄[1:end-1,2:end]))

include("ffjord.jl")
include("train.jl")
include("fast_layers.jl")
include("neural_de.jl")
include("require.jl")
include("spline_layer.jl")
include("tensor_product_basis.jl")
include("tensor_product_layer.jl")
include("collocation.jl")
include("hnn.jl")
include("multiple_shooting.jl")

Flux.device(::FastLayer) = @warn "device(f::FastLayer) is a no-op: to move FastChain computations to a GPU, apply gpu(x) to the weight vector"
Flux.gpu(::FastLayer) = @warn "device(f::FastLayer) is a no-op: to move FastChain computations to a GPU, apply gpu(x) to the weight vector"
Flux.cpu(::FastLayer) = @warn "device(f::FastLayer) is a no-op: to move FastChain computations to a CPU, apply cpu(x) to the weight vector"

export DeterministicCNF, FFJORD, NeuralODE, NeuralDSDE, NeuralSDE, NeuralCDDE, NeuralDAE, NeuralODEMM, TensorLayer, AugmentedNDELayer, SplineLayer, NeuralHamiltonianDE
export HamiltonianNN
export ChebyshevBasis, SinBasis, CosBasis, FourierBasis, LegendreBasis, PolynomialBasis
export FastDense, StaticDense, FastChain, initial_params
export FFJORDDistribution

export EpanechnikovKernel, UniformKernel, TriangularKernel, QuarticKernel
export TriweightKernel, TricubeKernel, GaussianKernel, CosineKernel
export LogisticKernel, SigmoidKernel, SilvermanKernel
export collocate_data

export multiple_shoot

end


text:
module RegService

using Sockets
using Distributed
using Pkg
using Logging

import RegistryTools: register
using ..Messaging

const CONFIG = Dict{String,Any}()

include("management.jl")

"""
    service(zsock::ReplySocket)
Start a server that calls `RegistryTools.register` on incoming
registration requests from `zsock`. Registration requests
must be serialized `RegistryTools.RegisterParams` objects. Close
the socket to halt the server.

Parameters:
- `zsock::ReplySocket`

Returns:
nothing
"""
function service(zsock::ReplySocket)
    while true
        ret = recvsend(zsock) do regp
            if isempty(regp.gitconfig)
                haskey(CONFIG, "user") && (regp.gitconfig["user.name"] = CONFIG["user"])
                haskey(CONFIG, "email") && (regp.gitconfig["user.email"] = CONFIG["email"])
            end
            register(regp)
        end
        ret || return
    end
    nothing
end

function main(config::AbstractString=isempty(ARGS) ? "config.toml" : first(ARGS))
    merge!(CONFIG, Pkg.TOML.parsefile(config)["regservice"])
    global_logger(SimpleLogger(stdout, get_log_level(CONFIG["log_level"])))
    zsock = ReplySocket(get(CONFIG, "port", 5555))

    @info("Starting registration service...")
    t = @async status_monitor(CONFIG["stop_file"], zsock)
    service(zsock)
    wait(t)

    # The !stopped! part is grep'ed by restart.sh, do not change
    @warn("Registration service !stopped!")
end

end    # module


text:
#=
MIT License

Optimal Testing and Containment Strategies for Universities in Mexico amid COVID-19

Copyright © 2021 Test and Contain. Luis Benavides-Vázquez, Héctor Alonso Guzmán-Gutiérrez, Jakob Jonnerby, Philip Lazos, Edwin Lock, Francisco J. Marmolejo-Cossío, Ninad Rajgopal,and José Roberto Tello-Ayala. https://www.testandcontain.com/

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
=#
using JSON

using ProgressMeter

include("iterators.jl")

struct Population
    n::Vector{Int}
    p::Vector{Float64}
    q::Vector{Float64}
    v::Vector{Float64}
    d::Array{Float64,2}
    π::Array{Float64,2}
end

mutable struct Config
    budget::Int  # number of tests available
    stepsize::Int  # size at which to bundle tests
    bucketsizes::Union{Nothing, Vector{Float64}}  # bucket size for all objectives (1 health + k quarantine)
    groupsizes::Vector{Int}
    guaranteed::Vector{Int}
end

Solutions = Vector{Vector{Float64}}  # for semantic convenience

"""Newer json reader, to ease automatization"""

function json_automatization(filename)
    data = JSON.parsefile(filename)

    info = data[1]["es"]
    populations = Array{Population}(undef, length(info))
    configs = Array{Config}(undef, length(info))
    file_names = Array{String}(undef, length(info))

    for campus in keys(info)
        n = [collect(values(i))[2] for i in values(info[campus]["categories"])]
        p = convert(Array{Float64, 1}, info[campus]["p"])
        v = convert(Array{Float64, 1}, info[campus]["v"])
        d = convert(Array{Float64, 2}, to_mat(info[campus]["d"]))
        π = convert(Array{Float64, 2}, to_mat(info[campus]["pi"]))
        q = 1 .- p
        populations[parse(Int,campus[7])] = Population(n, p, q, v, d, π)

        budget = info[campus]["budget"]
        stepsize = "stepsize" in keys(info[campus]) ? info[campus]["stepsize"] : 1
        groupsizes = "groupsizes" in keys(info[campus]) ? info[campus]["groupsizes"] : [1, 3, 5, 10]
        k = length(p)
        bucketsize = "bucketsizes" in keys(info[campus]) ? info[campus]["bucketsizes"] : nothing
        guaranteed = "guaranteed" in keys(info[campus]) ? info[campus]["guaranteed"] : zeros(k)
        configs[parse(Int,campus[7])] = Config(budget, stepsize, bucketsize, groupsizes, guaranteed)

        file_names[parse(Int,campus[7])] =  info[campus]["file"]
    end

    return populations, configs, file_names
end

"""Loads population and configuration data from json file."""
function load_from_json(filename)
    data = JSON.parsefile(filename)
    # load population
    n = convert(Array{Int, 1}, data["n"])
    p = convert(Array{Float64, 1}, data["p"])
    v = convert(Array{Float64, 1}, data["v"])
    d = convert(Array{Float64, 2}, to_mat(data["d"]))
    π = convert(Array{Float64, 2}, to_mat(data["pi"]))
    q = 1 .- p
    pop = Population(n, p, q, v, d, π)
    # load config
    budget = data["budget"]
    stepsize = data["stepsize"]
    bucketsize = "bucketsizes" in keys(data) ? data["bucketsizes"] : nothing
    groupsizes = data["groupsizes"]
    k = length(p)
    guaranteed = "guaranteed" in keys(data) ? data["guaranteed"] : zeros(k)
    config = Config(budget, stepsize, bucketsize, groupsizes, guaranteed)
    return pop, config
end
function to_mat(arrs) # for lists-of-lists parsed by JSON
    return Matrix(hcat(Vector{Float64}.(arrs)...)')
end

"""
Evaluate the health and quarantine objectives for given group sizes and number
of tests for each category.

NOTE: health objective is stored negated for simplicity.
"""
function objectives!(pop, prealloc, g, t, do_nothing)
    k = length(g)
    q = pop.q
    n = pop.n
    # Initialise objectives array
    prealloc.obj[1] = h_objective(pop, prealloc, g, t, do_nothing)
    prealloc.obj[2:end] = [q_objective(g[i], t[i], q[i], n[i]) for i in 1:k]
    return nothing
end

"""
Compute quarantine objective for a single category.

Note: Here q, t, q are numbers, not vectors.
"""
@inline q_objective(g, t, q, n) = (n==0 ? 0 : t*g*(q-q^g))

"""
Compute healthcare objective.

Note: the value is negated w.r.t. the IJCAI submission.
"""
@inline function h_objective(pop, prealloc, g, t, do_nothing)
    criticals = critical(pop, prealloc, g, t)
    return criticals - do_nothing
end

"""
Compute number of critical cases that occur when t tests of size g
are applied to each category.
"""
@inline function critical(pop, prealloc, g, t)
    n, p, q, v, π, d = pop.n, pop.p, pop.q, pop.v, pop.π, pop.d

    mask = (n.==0)
    n_to_divide = n
    n_to_divide = replace(n, 0=>1)

    z, temp = prealloc.z, prealloc.temp
    k = length(g)
    @. z = (n-t*g)*q/n_to_divide + t*g/n_to_divide*q^g
    @. temp = p*(n-t*g)/n_to_divide

    z[mask] .= 1
    temp[mask] .= 1

    @inbounds cases = sum(
        n[i]*v[i]*z[i]*(1-prod((1-π[i,j]*temp[j])^d[i,j] for j in 1:k))
            for i in 1:k)
    return cases
end

"""
Takes a single value or array and rounds each value to the nearest
multiple of bucketsize `bsize`. If no bucket size is provided, it returns
the original input.
"""
@inline function bucket(a, bsize=nothing)
    bsize === nothing && return a
    return round(a / bsize) * bsize
end

"""
Compute the Pareto frontier, output as an array.

Creates list of Pareto-optimal solutions, each consisting of a
Vector{Float64} with 3k+1 entries.
First 2k entries are for group sizes and number of tests.
Last k+1 entries are for health and quarantine objectives.

NOTE: The health objective is stored negated in order to simplify the code.
"""
function frontier(pop, config, nobuckets)
    k = length(pop.p)
    budget, step, guaranteed = config.budget, config.stepsize, config.guaranteed
    if nobuckets == true
        gsize, bsizes = config.groupsizes, nothing
    else
        gsize, bsizes = config.groupsizes, config.bucketsizes
    end
    solutions = Solutions()  # Create solutions list
    prealloc = preallocate(k)  # Preallocate memory
    # Compute number of critical cases if we apply no tests
    do_nothing = critical(pop, prealloc, ones(Int, k), zeros(Int, k))
    # Iterate over all possible solutions (groups, tests)
    @showprogress for groups in GroupIterator(gsize, k)
        for tests in TestIterator(budget, step, guaranteed, pop.n .÷ groups, k)
            objectives!(pop, prealloc, groups, tests, do_nothing)
            s = Float64[groups; tests; prealloc.obj]
            update!(solutions, s, bsizes)
        end
    end
    return solutions
end

"""
Given a new solution s, updates the list of solutions as follows.
If s is  by any solution in `solutions`, don't do anything.
Otherwise add s to `solutions` and remove all solutions s' in `solutions`
dominated by s.
"""
@inline function update!(solutions::Solutions, s::Vector{Float64}, bsizes=nothing)
    to_delete = Int[]  # collect solutions to be deleted
    for (i, other) in enumerate(solutions)
        # Check if sol dominates, or is dominated by, the current solution
        dominates, is_dominated = compare(s, other, bsizes)
        is_dominated && return nothing  # `s` is dominated by `other`
        dominates && push!(to_delete, i)  # store index for deletion
    end
    deleteat!(solutions, to_delete)  # delete all dominated solutions
    push!(solutions, s)  # add s to solutions
    return nothing
end

"""
Compares solutions `s` and `other` up to bucketing parameters (default=nothing).
Returns tuple (dominates, is_dominated) that expresses whether `s` dominates or
is_dominated by `other`.

NOTE: By default, function implements 'weak' dominance. This deviates from the
IJCAI paper, which stipulates that s dominates s' iff s[i] is strictly greater
than s'[i] for some i, and s[j] >= s[j] for all j.
Set `strong=true` to force `compare()` to follow the latter definition instead.
"""
@inline function compare(s, other, bsizes=nothing; strong=false)
    @assert length(s) === length(other)
    k = (length(s)-1) ÷ 3  # number of categories
    has_less, has_greater = false, false
    for i in 1:k+1
        bsize = bsizes===nothing ? nothing : bsizes[i]
        @inbounds comp = bucket(s[2k+i], bsize) - bucket(other[2k+i], bsize)
        if comp < 0
            has_less = true
        elseif comp > 0
            has_greater = true
        end
    end
    if !strong  # weak dominance parameter is true
        return !has_less, !has_greater
    else  # here we want to return truth values for strong dominance
        return !has_less && has_greater, !has_greater && has_less
    end
end

"""
Merges two solution lists into one. Specifically, it updates `solutions' to
contain the merged solutions and returns nothing.
"""
function merge!(solutions::Solutions, other::Solutions, bsizes=nothing)
    for s in other
        update!(solutions, s, bsizes)
    end
    return solutions
end

"""Preallocates memory for vectors that we compute frequently."""
function preallocate(k)
    s = zeros(Float64, 3k+1)  # new solution
    obj = zeros(k+1)  # objective of new solution
    z = zeros(Float64, k)
    temp = zeros(Float64, k)
    dominated = Int[]
    return (s=s, obj=obj, z=z, temp=temp, dominated=dominated)
end

max_stars_repo_name: hguzmang/eaamo

max_issues_repo_name: hguzmang/eaamo

max_forks_repo_name: hguzmang/eaamo


text:
using Manifolds
using Test
using LinearAlgebra

@testset "Testing @defVariable" begin
##

struct NotAManifold end

@test_throws AssertionError  @defVariable(MyVar, NotAManifold(), zeros(3,3))

##

@defVariable(TestVarType1, Euclidean(3), zeros(3))
@defVariable(TestVarType2, SpecialEuclidean(3), ProductRepr(zeros(3), diagm(ones(3))))

##

@test getManifold( TestVarType1) == Euclidean(3)
@test getManifold( TestVarType2) == SpecialEuclidean(3)

@test getDimension(TestVarType1) === 3
@test getDimension(TestVarType2) === 6

@test getPointType(TestVarType1) == Vector{Float64}
@test getPointType(TestVarType2) == ProductRepr{Tuple{Vector{Float64}, Matrix{Float64}}}

@test getPointIdentity(TestVarType1) == zeros(3)
@test getPointIdentity(TestVarType2).parts[1] == ProductRepr(zeros(3), diagm(ones(3))).parts[1]
@test getPointIdentity(TestVarType2).parts[2] == ProductRepr(zeros(3), diagm(ones(3))).parts[2]

##

@test getManifold( TestVarType1()) == Euclidean(3)
@test getManifold( TestVarType2()) == SpecialEuclidean(3)

@test getDimension(TestVarType1()) === 3
@test getDimension(TestVarType2()) === 6

@test getPointType(TestVarType1()) == Vector{Float64}
@test getPointType(TestVarType2()) == ProductRepr{Tuple{Vector{Float64}, Matrix{Float64}}}

@test getPointIdentity(TestVarType1()) == zeros(3)
@test getPointIdentity(TestVarType2()).parts[1] == ProductRepr(zeros(3), diagm(ones(3))).parts[1]
@test getPointIdentity(TestVarType2()).parts[2] == ProductRepr(zeros(3), diagm(ones(3))).parts[2]

end


text:
using TravelingSalesmanExact, GLPK, Test
set_default_optimizer!(with_optimizer(GLPK.Optimizer))

function test_valid_tour(tour, L)
    @test length(tour) == L
    @test isperm(tour)
    @test tour isa Vector{Int}
end

function tour_cost(tour, cost::AbstractMatrix)
    n = length(tour)
    inc(a) = a == n ? one(a) : a + 1
    s = zero(eltype(cost))
    for (index, city) in enumerate(tour)
        nextcity = tour[inc(index)]
        s += cost[city, nextcity]
    end
    return s
end

function tour_cost(
    tour,
    cities::AbstractVector;
    distance = TravelingSalesmanExact.euclidean_distance,
)
    n = length(cities)
    cost = [distance(cities[i], cities[j]) for i = 1:n, j = 1:n]
    return tour_cost(tour, cost)
end

function test_tour(input, opt = TravelingSalesmanExact.get_default_optimizer(); kwargs...)
    t, c = get_optimal_tour(input, opt; kwargs...)
    test_valid_tour(t, size(input, 1))
    if :distance in keys(kwargs)
        @test c ≈ tour_cost(t, input; distance = kwargs[:distance])
    else
        @test c ≈ tour_cost(t, input)
    end
    return t, c
end

@testset "Small random asymmetric" begin
    # chosen via rand(5,5)
    cost = 10 *
           [
        0.474886 0.350983 0.262651 0.138455 0.904042;
        0.683586 0.922968 0.278874 0.408406 0.0224372;
        0.513651 0.778167 0.140392 0.981211 0.891122;
        0.20529 0.976361 0.784706 0.98504 0.385203;
        0.489131 0.783738 0.538762 0.998821 0.0324331
    ]
    t1, c1 = test_tour(cost; verbose = true)
    t2, c2 = test_tour(cost; verbose = false)
    @test c1 ≈ c2

    # incorrect `symmetric` should give the wrong answers
    t3, c3 = get_optimal_tour(cost; verbose = false, symmetric = true)
    test_valid_tour(t3, 5)
    @test !(c2 ≈ c3)
    @test !(c3 ≈ tour_cost(t3, cost))
end

@testset "Small random symmetric" begin
    cost = 10 *
           [
        0.474886 0.350983 0.262651 0.138455 0.904042;
        0.683586 0.922968 0.278874 0.408406 0.0224372;
        0.513651 0.778167 0.140392 0.981211 0.891122;
        0.20529 0.976361 0.784706 0.98504 0.385203;
        0.489131 0.783738 0.538762 0.998821 0.0324331
    ]
    cost_sym = cost + transpose(cost)
    t4, c4 = test_tour(cost_sym; verbose = false)
    t5, c5 = test_tour(cost_sym; symmetric = true, verbose = false)
    t6, c6 = test_tour(cost_sym; symmetric = false, verbose = false)
    @test c4 ≈ c5
    @test c5 ≈ c6
end

@testset "Small random cities" begin
    # cities = [ 100*rand(2) for _ in 1:5]
    cities = Array{Float64,1}[
        [48.8885, 41.0517],
        [35.6635, 12.1844],
        [95.6122, 15.9847],
        [67.5772, 9.54407],
        [16.6325, 51.9001],
    ]
    t7, c7 = test_tour(cities; verbose = true)
    t8, c8 = @inferred test_tour(cities; verbose = false, symmetric = false)
    cost = [TravelingSalesmanExact.euclidean_distance(c1, c2) for c1 in cities, c2 in cities]
    t9, c9 = test_tour(cost; verbose = false, symmetric = false)
    t10, c10 = test_tour(cost; verbose = false, symmetric = true)
    @test c7 ≈ c8
    @test c8 ≈ c9
    @test c9 ≈ c10
    @test c8 isa Float64

end

@testset "Exceptions" begin
    cost = rand(5, 4)
    @test_throws ArgumentError test_tour(cost)

    cost = rand(5, 5)
    TravelingSalesmanExact.reset_default_optimizer!()
    @test_throws ArgumentError test_tour(cost)
    set_default_optimizer!(with_optimizer(GLPK.Optimizer))
end

@testset "att48.tsp" begin
    cities = simple_parse_tsp(joinpath(@__DIR__, "att48.tsp"))
    sym_tour, sym_cost = test_tour(
        cities;
        distance = TravelingSalesmanExact.ATT,
        verbose = true,
    )
    @test sym_cost ≈ 10628

    asym_tour, asym_cost = test_tour(
        cities;
        distance = TravelingSalesmanExact.ATT,
        symmetric = false,
    )
    @test asym_cost ≈ 10628
end


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

using AstroLib
using Test, Dates

include("utils-tests.jl")
include("misc-tests.jl")

# Dummy calls to "show" for new data types, just to increase code coverage.
show(devnull, AstroLib.planets["mercury"])
show(devnull, AstroLib.observatories["ca"])
show(devnull, AstroLib.observatories["vbo"])


text:
export graph

@reexport using LightGraphs
@reexport using MetaGraphs
function graph(m::Model)
    vars = variables(m)
    g = MetaDiGraph(length(vars))
    for (n,v) in enumerate(vars)
        set_prop!(g, n, :name, v)
    end
    set_indexing_prop!(g, :name)

    for (v, v_parents) in dependencies(m)
        for p in v_parents
            add_edge!(g, (g[p,:name],g[v,:name]))
        end
    end
    g
end

export graphEdges
function graphEdges(m::Model)
    g = graph(m)
    [(g[e.src,:name] => g[e.dst,:name]) for e in edges(g)]
end

# TODO: order statements using topological_sort_by_dfs(graph(m))

# EXAMPLE

# julia> g = graph(lda)
# {10, 14} directed Int64 metagraph with Float64 weights defined by :weight (default weight 1.0)

# julia> [(g[e.src,:name] => g[e.dst,:name]) for e in edges(g)]
# 14-element Array{Pair{Symbol,Symbol},1}:
#  :w => :N
#  :w => :M
#  :w => :z
#  :w => :β
#  :M => :N
#  :z => :N
#  :z => :M
#  :z => :θ
#  :β => :K
#  :β => :V
#  :β => :η
#  :θ => :α
#  :θ => :M
#  :θ => :K


text:
using DashExamples

app = DashExamples.make_app()

port = haskey(ENV, "PORT") ? parse(Int64, ENV["PORT"]) : 8050

DashExamples.Dash.run_server(app, "0.0.0.0", port)


text:
using TableTraits

import Gadfly
import DataFrames

function Gadfly.plot(data_source, elements::Gadfly.ElementOrFunctionOrLayers...; mapping...)
    return Gadfly.plot(DataFrames.DataFrame(data_source), elements...; mapping...)
end


text:
using
    MathOptInterface,
    BenchmarkTools,
    LinearAlgebra,
    Statistics,
    QuantEcon,
    Interpolations,
    NLsolve,
    Optim,
    Random,
    IterTools,
    JuMP,
    Ipopt,
    GaussianProcesses,
    Distances,
    Plots,
    Random;

numSectors = 2
numPoints1D = 35
gridMax = 50
gridMin = 10
grid = Vector{Float64}[]
for p in product(LinRange(gridMin,gridMax,numPoints1D),LinRange(gridMin,gridMax,numPoints1D))
    push!(grid, collect(p))
end
A = 1
l = 0.1
σ = 0.8
Id = 1.0*Matrix(I, numPoints1D^2, numPoints1D^2)
kappa(a, b) = A * exp(-.5 * norm(a - b)^2/l)
Kxx = Float64[]
for n in 1:length(grid)
    for m in 1:length(grid)
        x = grid[n]
        y = grid[m]
        push!(Kxx, kappa(x, y))
    end
end
Kxx = reshape(Kxx,numPoints1D^2,numPoints1D^2)
invKxx = inv(Kxx*Kxx + σ^2 * Id)
#the following needs wTarg to have been generated somewhere (at this stage I have via quantecon-constraints.jl in mind)
f = wTarg

xstar = grid[6]
Kxs = Float64[]
for n in 1:length(grid)
    x = grid[n]
    push!(Kxs, kappa(x, xstar))
end

mustar = Kxs'*inv(Kxx*Kxx)*f

function mu(y)
Kxs = Float64[]
for n in 1:length(grid)
    x = grid[n]
    push!(Kxs, kappa(x, y))
end
return Kxs'*invKxx*wTarg
end
numPlot = 30
plotGrid = Vector{Float64}[]

for p in product(LinRange(gridMin, gridMax,numPlot), LinRange(gridMin, gridMax,numPlot))
    push!(plotGrid, collect(p))
end
attemptGrid = zeros(numPlot,numPlot)
for i in 1:length(plotGrid)
    attemptGrid[Int(trunc((i-1)/numPlot))+1, i-Int(trunc((i-1)/numPlot))*numPlot] = mu(plotGrid[i])
end
trueGrid = zeros(numPlot,numPlot)
for i in 1:length(plotGrid)
    trueGrid[Int(trunc((i-1)/numPlot))+1, i-Int(trunc((i-1)/numPlot))*numPlot] = fstar(plotGrid[i])[1]
end

plot(trueGrid, st=:surface, label = "True Function")
plot!(attemptGrid, st=:surface, label = "GP Attempt")

max_stars_repo_name: uq-aibe/aim

max_issues_repo_name: uq-aibe/aim

max_forks_repo_name: uq-aibe/aim


text:
using Dierckx

"""
Bellman equation for the continuation of the Lucas-Stokey Problem
"""
mutable struct BellmanEquation_Recursive{TP <: Model, TI <: Integer, TR <: Real}
    model::TP
    S::TI
    xbar::Array{TR}
    time_0::Bool
    z0::Array{Array}
    cFB::Vector{TR}
    nFB::Vector{TR}
    xFB::Vector{TR}
    zFB::Vector{Vector{TR}}
end

"""
Compute the planner's allocation by solving Bellman
equation.
"""
struct RecursiveAllocation{TP <: Model,
                           TI <: Integer,
                           TVg <: AbstractVector,
                           TT <: Tuple}
    model::TP
    mc::MarkovChain
    S::TI
    T::BellmanEquation_Recursive
    μgrid::TVg
    xgrid::TVg
    Vf::Array
    policies::TT
end

"""
Initializes the type from the calibration Model
"""
function RecursiveAllocation(model::Model, μgrid::AbstractArray)
    G = model.G
    S = size(model.Π, 1)             # number of states
    mc = MarkovChain(model.Π)
    # now find the first best allocation
    Vf, policies, T, xgrid = solve_time1_bellman(model, μgrid)
    T.time_0 = true                  # Bellman equation now solves time 0 problem
    return RecursiveAllocation(model, mc, S, T, μgrid, xgrid, Vf, policies)
end

"""
Solve the time  1 Bellman equation for calibration Model and initial grid μgrid
"""
function solve_time1_bellman{TR <: Real}(model::Model{TR}, μgrid::AbstractArray)
    Π = model.Π
    S = size(model.Π, 1)

    # First get initial fit from lucas stockey solution.
    # Need to change things to be ex_ante
    PP = SequentialAllocation(model)

    function incomplete_allocation(PP::SequentialAllocation,
                                   μ_::AbstractFloat,
                                   s_::Integer)
        c, n, x, V = time1_value(PP, μ_)
        return c, n, dot(Π[s_, :], x), dot(Π[s_, :], V)
    end

    cf = Array{Function}(S, S)
    nf = Array{Function}(S, S)
    xprimef = Array{Function}(S, S)
    Vf = Vector{Function}(S)
    xgrid = Array{TR}(S, length(μgrid))

    for s_ in 1:S
        c = Array{TR}(length(μgrid), S)
        n = Array{TR}(length(μgrid), S)
        x = Array{TR}(length(μgrid))
        V = Array{TR}(length(μgrid))
        for (i_μ, μ) in enumerate(μgrid)
            c[i_μ, :], n[i_μ, :], x[i_μ], V[i_μ] =
                incomplete_allocation(PP, μ, s_)
        end
        xprimes = repmat(x, 1, S)
        xgrid[s_, :] = x
        for sprime = 1:S
            splc = Spline1D(x[end:-1:1], c[:, sprime][end:-1:1], k=3)
            spln = Spline1D(x[end:-1:1], n[:, sprime][end:-1:1], k=3)
            splx = Spline1D(x[end:-1:1], xprimes[:, sprime][end:-1:1], k=3)
            cf[s_, sprime] = y -> splc(y)
            nf[s_, sprime] = y -> spln(y)
            xprimef[s_, sprime] = y -> splx(y)
            # cf[s_, sprime] = LinInterp(x[end:-1:1], c[:, sprime][end:-1:1])
            # nf[s_, sprime] = LinInterp(x[end:-1:1], n[:, sprime][end:-1:1])
            # xprimef[s_, sprime] = LinInterp(x[end:-1:1], xprimes[:, sprime][end:-1:1])
        end
        splV = Spline1D(x[end:-1:1], V[end:-1:1], k=3)
        Vf[s_] = y -> splV(y)
        # Vf[s_] = LinInterp(x[end:-1:1], V[end:-1:1])
    end

    policies = [cf, nf, xprimef]

    # Create xgrid
    xbar = [maximum(minimum(xgrid)), minimum(maximum(xgrid))]
    xgrid = linspace(xbar[1], xbar[2], length(μgrid))

    # Now iterate on Bellman equation
    T = BellmanEquation_Recursive(model, xgrid, policies)
    diff = 1.0
    while diff > 1e-4
        PF = (i_x, x, s) -> get_policies_time1(T, i_x, x, s, Vf, xbar)
        Vfnew, policies = fit_policy_function(T, PF, xgrid)

        diff = 0.0
        for s=1:S
            diff = max(diff, maximum(abs, (Vf[s].(xgrid) - Vfnew[s].(xgrid)) ./ Vf[s].(xgrid)))
        end

        println("diff = $diff")
        Vf = copy(Vfnew)
    end

    return Vf, policies, T, xgrid
end

"""
Fits the policy functions
"""
function fit_policy_function{TF<:AbstractFloat}(T::BellmanEquation_Recursive,
                                                PF::Function,
						xgrid::AbstractVector{TF})
    S = T.S
    # preallocation
    PFvec = Array{TF}(4S + 1, length(xgrid))
    cf = Array{Function}(S, S)
    nf = Array{Function}(S, S)
    xprimef = Array{Function}(S, S)
    TTf = Array{Function}(S, S)
    Vf = Vector{Function}(S)
    # fit policy fuctions
    for s_ in 1:S
        for (i_x, x) in enumerate(xgrid)
            PFvec[:, i_x] = PF(i_x, x, s_)
        end
        splV = Spline1D(xgrid, PFvec[1,:], k=3)
        Vf[s_] = y -> splV(y)
        # Vf[s_] = LinInterp(xgrid, PFvec[1, :])
        for sprime=1:S
            splc = Spline1D(xgrid, PFvec[1 + sprime, :], k=3)
            spln = Spline1D(xgrid, PFvec[1 + S + sprime, :], k=3)
            splxprime = Spline1D(xgrid, PFvec[1 + 2S + sprime, :], k=3)
            splTT = Spline1D(xgrid, PFvec[1 + 3S + sprime, :], k=3)
            cf[s_, sprime] = y -> splc(y)
            nf[s_, sprime] = y -> spln(y)
            xprimef[s_, sprime] = y -> splxprime(y)
            TTf[s_, sprime] = y -> splTT(y)
        end
    end
    policies = (cf, nf, xprimef, TTf)
    return Vf, policies
end

"""
Computes Tau given c,n
"""
function Tau(pab::RecursiveAllocation,
             c::AbstractArray,
	     n::AbstractArray)
    model = pab.model
    Uc, Un = model.Uc(c, n), model.Un(c, n)
    return 1 + Un ./ (model.Θ .* Uc)
end

Tau(pab::RecursiveAllocation, c::Real, n::Real) = Tau(pab, [c], [n])

"""
Finds the optimal allocation given initial government debt B_ and state s_0
"""
function time0_allocation(pab::RecursiveAllocation, B_::Real, s0::Integer)
    T, Vf = pab.T, pab.Vf
    xbar = T.xbar
    z0 = get_policies_time0(T, B_, s0, Vf, xbar)

    c0, n0, xprime0, T0 = z0[2], z0[3], z0[4], z0[5]
    return c0, n0, xprime0, T0
end

"""
Simulates planners policies for `T` periods
"""
function simulate{TF <: AbstractFloat}(pab::RecursiveAllocation,
                                       B_::TF, s_0::Integer, T::Integer,
                                       sHist::Vector=simulate(pab.mc, T, init=s_0))
    model, mc, Vf, S = pab.model, pab.mc, pab.Vf, pab.S
    Π, Uc = model.Π, model.Uc
    cf, nf, xprimef, TTf = pab.policies

    cHist = Array{TF}(T)
    nHist = Array{TF}(T)
    Bhist = Array{TF}(T)
    xHist = Array{TF}(T)
    TauHist = Array{TF}(T)
    THist = Array{TF}(T)
    μHist = Array{TF}(T)

    #time0
    cHist[1], nHist[1], xHist[1], THist[1]  = time0_allocation(pab, B_, s_0)
    TauHist[1] = Tau(pab, cHist[1], nHist[1])[s_0]
    Bhist[1] = B_
    μHist[1] = Vf[s_0](xHist[1])

    #time 1 onward
    for t in 2:T
        s_, x, s = sHist[t-1], xHist[t-1], sHist[t]
        c = Array{TF}(S)
        n = Array{TF}(S)
        xprime = Array{TF}(S)
        TT = Array{TF}(S)
        for sprime=1:S
            c[sprime], n[sprime], xprime[sprime], TT[sprime] =
                cf[s_, sprime](x), nf[s_, sprime](x),
                xprimef[s_, sprime](x), TTf[s_, sprime](x)
        end

        Tau_val = Tau(pab, c, n)[s]
        u_c = Uc(c, n)
        Eu_c = dot(Π[s_, :], u_c)

        μHist[t] = Vf[s](xprime[s])

        cHist[t], nHist[t], Bhist[t], TauHist[t] = c[s], n[s], x/Eu_c, Tau_val
        xHist[t], THist[t] = xprime[s], TT[s]
    end
    return cHist, nHist, Bhist, xHist, TauHist, THist, μHist, sHist
end

"""
Initializes the class from the calibration Model
"""
function BellmanEquation_Recursive{TF <: AbstractFloat}(model::Model{TF},
                                                        xgrid::AbstractVector{TF},
                                                        policies0::Array)

    S = size(model.Π, 1)                                # number of states
    xbar = [minimum(xgrid), maximum(xgrid)]
    time_0 = false
    z0 = Array{Array}(length(xgrid), S)
    cf, nf, xprimef = policies0[1], policies0[2], policies0[3]
    for s in 1:S
        for (i_x, x) in enumerate(xgrid)
            cs = Array{TF}(S)
            ns = Array{TF}(S)
            xprimes = Array{TF}(S)
            for j = 1:S
                cs[j], ns[j], xprimes[j] = cf[s, j](x), nf[s, j](x), xprimef[s, j](x)
            end
            z0[i_x, s] = vcat(cs, ns, xprimes, zeros(S))
        end
    end
    cFB, nFB, IFB, xFB, zFB = find_first_best(model, S, 2)
    return BellmanEquation_Recursive(model, S, xbar, time_0, z0, cFB, nFB, xFB, zFB)
end

"""
Finds the optimal policies
"""
function get_policies_time1(T::BellmanEquation_Recursive,
                            i_x::Integer,
                            x::Real,
                            s_::Integer,
                            Vf::AbstractArray{Function},
                            xbar::AbstractVector)
    model, S = T.model, T.S
    β, Θ, G, Π = model.β, model.Θ, model.G, model.Π
    U,Uc,Un = model.U, model.Uc, model.Un

    S_possible = sum(Π[s_, :].>0)
    sprimei_possible = find(Π[s_, :].>0)

    function objf(z, grad)
        c, xprime = z[1:S_possible], z[S_possible+1:2S_possible]
        n = (c + G[sprimei_possible]) ./ Θ[sprimei_possible]
        Vprime = [Vf[sprimei_possible[si]](xprime[si]) for si in 1:S_possible]
        return -dot(Π[s_, sprimei_possible], U.(c, n) + β * Vprime)
    end

    function cons(out, z, grad)
        c, xprime, TT =
            z[1:S_possible], z[S_possible + 1:2S_possible], z[2S_possible + 1:3S_possible]
        n = (c+G[sprimei_possible]) ./ Θ[sprimei_possible]
        u_c = Uc.(c, n)
        Eu_c = dot(Π[s_, sprimei_possible], u_c)
        out .= x * u_c/Eu_c - u_c .* (c - TT) - Un(c, n) .* n - β * xprime
    end
    function cons_no_trans(out, z, grad)
        c, xprime = z[1:S_possible], z[S_possible + 1:2S_possible]
        n = (c + G[sprimei_possible]) ./ Θ[sprimei_possible]
        u_c = Uc.(c, n)
        Eu_c = dot(Π[s_, sprimei_possible], u_c)
        out .= x * u_c / Eu_c - u_c .* c - Un(c, n) .* n - β * xprime
    end

    if model.transfers == true
        lb = vcat(zeros(S_possible), ones(S_possible)*xbar[1], zeros(S_possible))
        if model.n_less_than_one == true
            ub = vcat(ones(S_possible) - G[sprimei_possible],
                      ones(S_possible) * xbar[2], ones(S_possible))
        else
            ub = vcat(100 * ones(S_possible),
                      ones(S_possible) * xbar[2],
                      100 * ones(S_possible))
        end
        init = vcat(T.z0[i_x, s_][sprimei_possible],
                    T.z0[i_x, s_][2S + sprimei_possible],
                    T.z0[i_x, s_][3S + sprimei_possible])
        opt = Opt(:LN_COBYLA, 3S_possible)
        equality_constraint!(opt, cons, zeros(S_possible))
    else
        lb = vcat(zeros(S_possible), ones(S_possible)*xbar[1])
        if model.n_less_than_one == true
            ub = vcat(ones(S_possible)-G[sprimei_possible], ones(S_possible)*xbar[2])
        else
            ub = vcat(ones(S_possible), ones(S_possible) * xbar[2])
        end
        init = vcat(T.z0[i_x, s_][sprimei_possible],
                    T.z0[i_x, s_][2S+sprimei_possible])
        opt = Opt(:LN_COBYLA, 2S_possible)
        equality_constraint!(opt, cons_no_trans, zeros(S_possible))
    end
    init[init .> ub] = ub[init .> ub]
    init[init .< lb] = lb[init .< lb]

    min_objective!(opt, objf)
    lower_bounds!(opt, lb)
    upper_bounds!(opt, ub)
    maxeval!(opt, 10000000)
    maxtime!(opt, 10)
    ftol_rel!(opt, 1e-8)
    ftol_abs!(opt, 1e-8)

    (minf, minx, ret) = optimize(opt, init)

    if ret != :SUCCESS && ret != :ROUNDOFF_LIMITED && ret != :MAXEVAL_REACHED &&
         ret != :FTOL_REACHED && ret != :MAXTIME_REACHED
        error("optimization failed: ret = $ret")
    end

    T.z0[i_x, s_][sprimei_possible] = minx[1:S_possible]
    T.z0[i_x, s_][S + sprimei_possible] = minx[1:S_possible] + G[sprimei_possible]
    T.z0[i_x, s_][2S + sprimei_possible] = minx[S_possible + 1:2S_possible]
    if model.transfers == true
        T.z0[i_x, s_][3S + sprimei_possible] = minx[2S_possible + 1:3S_possible]
    else
        T.z0[i_x, s_][3S + sprimei_possible] = zeros(S)
    end

    return vcat(-minf, T.z0[i_x, s_])
end

"""
Finds the optimal policies
"""
function get_policies_time0(T::BellmanEquation_Recursive,
                            B_::Real,
                            s0::Integer,
                            Vf::AbstractArray{Function},
                            xbar::AbstractVector)
    model = T.model
    β, Θ, G = model.β, model.Θ, model.G
    U, Uc, Un = model.U, model.Uc, model.Un

    function objf(z, grad)
        c, xprime = z[1], z[2]
        n = (c + G[s0]) / Θ[s0]
        return -(U(c, n) + β * Vf[s0](xprime))
    end

    function cons(z,grad)
        c, xprime, TT = z[1], z[2], z[3]
        n = (c + G[s0]) / Θ[s0]
        return -Uc(c, n) * (c - B_ - TT) - Un(c, n) * n - β * xprime
    end
    cons_no_trans(z, grad) = cons(vcat(z, 0), grad)

    if model.transfers == true
        lb = [0.0, xbar[1], 0.0]
        if model.n_less_than_one == true
            ub = [1 - G[s0], xbar[2], 100]
        else
            ub = [100.0, xbar[2], 100.0]
        end
        init = vcat(T.zFB[s0][1], T.zFB[s0][3], T.zFB[s0][4])
        init = [0.95124922, -1.15926816,  0.0]
        opt = Opt(:LN_COBYLA, 3)
        equality_constraint!(opt, cons)
    else
        lb = [0.0, xbar[1]]
        if model.n_less_than_one == true
            ub = [1-G[s0], xbar[2]]
        else
            ub = [100, xbar[2]]
        end
        init = vcat(T.zFB[s0][1], T.zFB[s0][3])
        init = [0.95124922, -1.15926816]
        opt = Opt(:LN_COBYLA, 2)
        equality_constraint!(opt, cons_no_trans)
    end
    init[init .> ub] = ub[init .> ub]
    init[init .< lb] = lb[init .< lb]

    min_objective!(opt, objf)
    lower_bounds!(opt, lb)
    upper_bounds!(opt, ub)
    maxeval!(opt, 100000000)
    maxtime!(opt, 30)

    (minf, minx, ret) = optimize(opt, init)

    if ret != :SUCCESS && ret != :ROUNDOFF_LIMITED && ret != :MAXEVAL_REACHED && ret != :FTOL_REACHED
        error("optimization failed: ret = $ret")
    end

    if model.transfers == true
        return -minf, minx[1], minx[1]+G[s0], minx[2], minx[3]
    else
        return -minf, minx[1], minx[1]+G[s0], minx[2], 0
    end
end


text:
## FROM CLOUD POINT

using SparseArrays
@using LinearAlgebraicRepresentation as LAR
@using ViewerGL as GL

include("./../src/Lar.jl")
LarA = Lar.Arrangement

include("./models/casa.jl")

Lar.addModelCells!(casa, 2, LarA.biconnected_components(casa.T[1]))

bbox = LarA.getModelBoundingBoxes(casa, 2)
bboxXL = [box .+ [-0.5 +0.5] for box in bbox]

model = Lar.Model(zeros(3, 0));

for comp = 1 : length(bboxXL)
    pts = Lar.getModelCellVertices(casa, 2, comp);
    # plane : ax + by + cz = d
    d = 1
    a, b, c = [pts[1]'; pts[floor(Int, 1+end/2)]'; pts[end]']\ones(3, 1)*d

    x = bboxXL[comp][1, :]
    y = bboxXL[comp][2, :]
    z = bboxXL[comp][3, :]

    CM = Lar.Model(zeros(3, 0));
    V = [];

    px = zeros(2, 2)
    px[1, 1] = (d - b*y[1] - c*z[1]) / a
    px[1, 2] = (d - b*y[1] - c*z[2]) / a
    px[2, 1] = (d - b*y[2] - c*z[1]) / a
    px[2, 2] = (d - b*y[2] - c*z[2]) / a
    for j = 1 : 2  for i = 1 : 2  if x[1] < px[i, j] < x[2]
            push!(V, [px[i, j]; y[i]; z[j]])
            Lar.addModelVertex!(CM, [px[i, j]; y[i]; z[j]])
    end  end  end

    py = zeros(2, 2)
    py[1, 1] = (d - a*x[1] - c*z[1]) / b
    py[1, 2] = (d - a*x[1] - c*z[2]) / b
    py[2, 1] = (d - a*x[2] - c*z[1]) / b
    py[2, 2] = (d - a*x[2] - c*z[2]) / b
    for j = 1 : 2  for i = 1 : 2  if y[1] < py[i, j] < y[2]
            push!(V, [x[i]; py[i, j]; z[j]])
            Lar.addModelVertex!(CM, [x[i]; py[i, j]; z[j]])
    end  end  end

    pz = zeros(2, 2)
    pz[1, 1] = (d - a*x[1] - b*y[1]) / c
    pz[1, 2] = (d - a*x[1] - b*y[2]) / c
    pz[2, 1] = (d - a*x[2] - b*y[1]) / c
    pz[2, 2] = (d - a*x[2] - b*y[2]) / c
    for j = 1 : 2  for i = 1 : 2  if z[1] < pz[i, j] < z[2]
        push!(V, [x[i]; y[j]; pz[i, j]])
        Lar.addModelVertex!(CM, [x[i]; y[j]; pz[i, j]])
    end  end  end

    for i = 1 : size(CM, 0, 2)   for j = i + 1 : size(CM, 0, 2)
        if sum(CM.G[:, i] .== CM.G[:, j]) >= 1
            Lar.addModelCell!(CM, 1, [i, j])
        end
    end  end
    Lar.addModelCell!(CM, 2, collect(1:size(CM, 0, 2)))

    Lar.uniteModels!(model, CM)
end


text:
export UNet
# define a UNet class in Julia
@pydef mutable struct UNet <: nn.Module
    function __init__(self, num_channels=3, num_classes=3)
        pybuiltin(:super)(UNet, self).__init__()
        self.relu = nn.ReLU()
        self.maxpool = nn.MaxPool1d(kernel_size=2, stride=2)
        self.conv11 = nn.Conv1d(num_channels, 64, kernel_size=3, padding=1)
        self.conv12 = nn.Conv1d(64, 64, kernel_size=3, padding=1)

        self.conv21 = nn.Conv1d(64, 128, kernel_size=3, padding=1)
        self.conv22 = nn.Conv1d(128, 128, kernel_size=3, padding=1)

        self.conv31 = nn.Conv1d(128, 256, kernel_size=3, padding=1)
        self.conv32 = nn.Conv1d(256, 256, kernel_size=3, padding=1)

        self.conv41 = nn.Conv1d(256, 512, kernel_size=3, padding=1)
        self.conv42 = nn.Conv1d(512, 512, kernel_size=3, padding=1)

        self.conv51 = nn.Conv1d(512, 1024, kernel_size=3, padding=1)
        self.conv52 = nn.Conv1d(1024, 1024, kernel_size=3, padding=1)

        self.uconv6 = nn.ConvTranspose1d(1024, 512, kernel_size=2, stride=2)
        self.conv61 = nn.Conv1d(1024, 512, kernel_size=3, padding=1)
        self.conv62 = nn.Conv1d(512, 512, kernel_size=3, padding=1)

        self.uconv7 = nn.ConvTranspose1d(512, 256, kernel_size=2, stride=2)
        self.conv71 = nn.Conv1d(512, 256, kernel_size=3, padding=1)
        self.conv72 = nn.Conv1d(256, 256, kernel_size=3, padding=1)

        self.uconv8 = nn.ConvTranspose1d(256, 128, kernel_size=2, stride=2)
        self.conv81 = nn.Conv1d(256, 128, kernel_size=3, padding=1)
        self.conv82 = nn.Conv1d(128, 128, kernel_size=3, padding=1)

        self.uconv9 = nn.ConvTranspose1d(128, 64, kernel_size=2, stride=2)
        self.conv91 = nn.Conv1d(128, 64, kernel_size=3, padding=1)
        self.conv92 = nn.Conv1d(64, 64, kernel_size=3, padding=1)

        self.conv93 = nn.Conv1d(64, num_classes, kernel_size=1, padding=0)

        self.sigmoid = nn.Sigmoid()
    end

    function forward(self, x)
        x = self.conv11(x)
        x = self.relu(x)
        x = self.conv12(x)
        x1d = self.relu(x)
        x = self.maxpool(x1d)

        x = self.conv21(x)
        x = self.relu(x)
        x = self.conv22(x)
        x2d = self.relu(x)
        x = self.maxpool(x2d)

        x = self.conv31(x)
        x = self.relu(x)
        x = self.conv32(x)
        x3d = self.relu(x)
        x = self.maxpool(x3d)

        x = self.conv41(x)
        x = self.relu(x)
        x = self.conv42(x)
        x4d = self.relu(x)
        x = self.maxpool(x4d)

        x = self.conv51(x)
        x = self.relu(x)
        x = self.conv52(x)
        x5d = self.relu(x)

        x6u = self.uconv6(x5d)
        x = torch.cat((x4d, x6u), 1)
        x = self.conv61(x)
        x = self.relu(x)
        x = self.conv62(x)
        x = self.relu(x)

        x7u = self.uconv7(x)
        x = torch.cat((x3d, x7u), 1)
        x = self.conv71(x)
        x = self.relu(x)
        x = self.conv72(x)
        x = self.relu(x)

        x8u = self.uconv8(x)
        x = torch.cat((x2d, x8u), 1)
        x = self.conv81(x)
        x = self.relu(x)
        x = self.conv82(x)
        x = self.relu(x)

        x9u = self.uconv9(x)
        x = torch.cat((x1d, x9u), 1)
        x = self.conv91(x)
        x = self.relu(x)
        x = self.conv92(x)
        x = self.relu(x)

        x = self.conv93(x)

        x = self.sigmoid(x)
        #x = self.softmax(x)

        return x
    end
end


text:
#taken from github.com/jbierkens/ICML-boomerang
include("../../src/faberschauder.jl")
include("../../src/asvar.jl")
include("../../src/pdmp.jl")

max_stars_repo_name: SebaGraz/ZZDiffusionBridge

max_issues_repo_name: SebaGraz/ZZDiffusionBridge

max_forks_repo_name: SebaGraz/ZZDiffusionBridge


text:
function aalen(X, t, status; entry=[], id=[], beta=[], opt...)
end


text:
struct PropertySpec
    name::Symbol
    typ::Type
    mask::UInt16
end

const NON_PERSISTENT = TAMQPOctet(1)
const PERSISTENT = TAMQPOctet(2)

const PROPERTIES = Dict{Symbol, PropertySpec}(
    :content_type     => PropertySpec(:content_type,     TAMQPShortStr,   0x0001 << 15), # MIME content type (MIME typing)
    :content_encoding => PropertySpec(:content_encoding, TAMQPShortStr,   0x0001 << 14), # MIME content encoding (MIME typing)
    :headers          => PropertySpec(:headers,          TAMQPFieldTable, 0x0001 << 13), # message header field table (For applications, and for header exchange routing)
    :delivery_mode    => PropertySpec(:delivery_mode,    TAMQPOctet,      0x0001 << 12), # non-persistent (1) or persistent (2) (For queues that implement persistence)
    :priority         => PropertySpec(:priority,         TAMQPOctet,      0x0001 << 11), # message priority, 0 to 9 (For queues that implement priorities)
    :correlation_id   => PropertySpec(:correlation_id,   TAMQPShortStr,   0x0001 << 10), # application correlation identifier (For application use, no formal behaviour)
    :reply_to         => PropertySpec(:reply_to,         TAMQPShortStr,   0x0001 << 9),  # address to reply to (For application use, no formal behaviour)
    :expiration       => PropertySpec(:expiration,       TAMQPShortStr,   0x0001 << 8),  # message expiration specification (For application use, no formal behaviour)
    :message_id       => PropertySpec(:message_id,       TAMQPShortStr,   0x0001 << 7),  # application message identifier (For application use, no formal behaviour)
    :timestamp        => PropertySpec(:timestamp,        TAMQPTimeStamp,  0x0001 << 6),  # message timestamp (For application use, no formal behaviour)
    :message_type     => PropertySpec(:message_type,     TAMQPShortStr,   0x0001 << 5),  # message type name (For application use, no formal behaviour)
    :user_id          => PropertySpec(:user_id,          TAMQPShortStr,   0x0001 << 4),  # creating user id (For application use, no formal behaviour)
    :app_id           => PropertySpec(:app_id,           TAMQPShortStr,   0x0001 << 3),  # creating application id (For application use, no formal behaviour)
    :cluster_id       => PropertySpec(:cluster_id,       TAMQPShortStr,   0x0001 << 2)   # reserved, must be empty (Deprecated, was old cluster-id property)
)

const SORTED_PROPERTY_NAMES = [:content_type, :content_encoding, :headers, :delivery_mode, :priority, :correlation_id, :reply_to, :expiration, :message_id, :timestamp, :message_type, :user_id, :app_id, :cluster_id]
const SORTED_PROPERTIES = [PROPERTIES[k] for k in SORTED_PROPERTY_NAMES]

mutable struct Message
    data::Vector{UInt8}
    properties::Dict{Symbol,TAMQPField}

    filled::Int

    consumer_tag::String
    delivery_tag::TAMQPDeliveryTag
    redelivered::Bool
    exchange::String
    routing_key::String
    remaining::TAMQPMessageCount
end

function Message(data::Vector{UInt8}; kwargs...)
    msg = Message(data, Dict{Symbol,TAMQPField}(), length(data), "", TAMQPDeliveryTag(0), false, "", "", TAMQPMessageCount(0))
    set_properties(msg; kwargs...)
    msg
end

function set_properties(msg::Message; kwargs...)
    for (k,v) in kwargs
        if v === nothing
            delete!(msg.properties, k)
        else
            msg.properties[k] = convert(PROPERTIES[k].typ, v)
        end
    end
    nothing
end


text:
using HolyTraitPattern
using Test

@testset "HolyTraitPattern.jl" begin
    include("funcs.jl")
    @test trait_test_cash()
    @test trait_test_stock()
    @test trait_test_residence()
    @test trait_test_bond()
    @test trait_test_book()
end


text:
using Test, AutoARIMA

include("math.jl")
include("tsafac.jl")
include("forecast.jl")
include("simulated.jl")


text:
v = [9, 8, 7, 6, 5, 0, 1, 2, 3, 4]
@show v select(v, 1:10)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Mimi
using DataFrames
using CSVFiles

"""
    update_MimiFAIR_params!(m;
                            rcp_scenario::String="RCP85",
                            start_year::Int=1765,
                            end_year::Int=2500,
                            F2x::Float64=3.71,
                            TCR::Float64=1.6,
                            ECS::Float64=2.75,
                            d::Array{Float64,1}=[239.0, 4.1])

Update all parameter settings in the model `m` with default FAIR settings as pulled
from the MimiFAIR repository.  Optional arguments include `rcp_scenario` to indicate
what scenario to use, `start_year` and `end_year` to indicate the model time dimension,
and three parameter settings for the paramters `F2x`, `TCR`, and `ECS`.  See the
MimiFAIR model documentatoin for more.
"""
function update_MimiFAIR_params!(m; rcp_scenario::String="RCP85", start_year::Int=1765, end_year::Int=2500, F2x::Float64=3.71, TCR::Float64=1.6, ECS::Float64=2.75, d::Array{Float64,1}=[239.0, 4.1])

    # ---------------------------------------------
    # Set up some data.
    # ---------------------------------------------

    # Load RCP and other data needed to construct FAIR.
    rcp_emissions, volcano_forcing, solar_forcing, gas_data, gas_fractions, conversions = MimiFAIR.load_fair_data(start_year, end_year, rcp_scenario)

    # Names of minor greenhouse gases and ozone-depleting substances.
    other_ghg_names = ["CF4", "C2F6", "C6F14", "HFC23", "HFC32", "HFC43_10", "HFC125", "HFC134a", "HFC143a", "HFC227ea", "HFC245fa", "SF6", "CFC_11", "CFC_12", "CFC_113", "CFC_114", "CFC_115", "CARB_TET", "MCF", "HCFC_22", "HCFC_141B", "HCFC_142B", "HALON1211", "HALON1202", "HALON1301", "HALON2402", "CH3BR", "CH3CL"]
    ods_names = ["CFC_11", "CFC_12", "CFC_113", "CFC_114", "CFC_115", "CARB_TET", "MCF", "HCFC_22", "HCFC_141B", "HCFC_142B", "HALON1211", "HALON1202", "HALON1301", "HALON2402", "CH3BR", "CH3CL"]

    # Calculate CO₂ radiative forcing scaling term (to keep forcing from 2x CO₂ consistent).
    scale_co2_forcing = MimiFAIR.co2_rf_scale(F2x, gas_data[gas_data.gas .== "CO2", :pi_conc][1], gas_data[gas_data.gas .== "N2O", :pi_conc][1])

    # Calcualte coefficients of slow and fast temperature change in each timestep based on user=specified parameters.
    q = MimiFAIR.calculate_q(TCR, ECS, d, F2x)

    # Number of time periods to run model.
    n_steps = length(start_year:end_year)

    # ---------------------------------------------
    # Set component parameters
    # ---------------------------------------------

    # ---- Methane Cycle ---- #
    update_param!(m, :ch4_cycle, :fossil_emiss_CH₄, rcp_emissions.CH4)
    update_param!(m, :ch4_cycle, :natural_emiss_CH₄, rcp_emissions.NaturalCH4)
    update_param!(m, :ch4_cycle, :τ_CH₄, 9.3)
    update_param!(m, :ch4_cycle, :fossil_frac, gas_fractions.ch4_fossil)
    update_param!(m, :ch4_cycle, :oxidation_frac, 0.61)
    update_param!(m, :ch4_cycle, :mol_weight_CH₄, gas_data[gas_data.gas .== "CH4", :mol_weight][1])
    update_param!(m, :ch4_cycle, :mol_weight_C, gas_data[gas_data.gas .== "C", :mol_weight][1])
    update_param!(m, :ch4_cycle, :emiss2conc_ch4, conversions[conversions.gases .== "CH4", :emiss2conc][1])

    # ---- Nitrous Oxide Cycle ---- #
    update_param!(m, :n2o_cycle, :fossil_emiss_N₂O, rcp_emissions.N2O)
    update_param!(m, :n2o_cycle, :natural_emiss_N₂O, rcp_emissions.NaturalN2O)
    update_param!(m, :n2o_cycle, :τ_N₂O, 121.0)
    update_param!(m, :n2o_cycle, :emiss2conc_n2o, conversions[conversions.gases .== "N2O", :emiss2conc][1])

    # ---- Carbon Cycle ---- #
    update_param!(m, :co2_cycle, :CO2_0, gas_data[gas_data.gas .== "CO2", :pi_conc][1])
    update_param!(m, :co2_cycle, :r0, 35.0)
    update_param!(m, :co2_cycle, :rC, 0.019)
    update_param!(m, :co2_cycle, :rT, 4.165)
    update_param!(m, :co2_cycle, :iIRF_max, 97.0)
    update_param!(m, :co2_cycle, :a, [0.2173, 0.2240, 0.2824, 0.2763])
    update_param!(m, :co2_cycle, :τ_CO₂, [10.0^6, 394.4, 36.54, 4.304])
    update_param!(m, :co2_cycle, :E_CO₂, rcp_emissions.FossilCO2 .+ rcp_emissions.OtherCO2)

    # ---- Other Well-Mixed Greenhouse Gas Cycles ---- #
    update_param!(m, :other_ghg_cycles, :τ_other_ghg, gas_data[findall((in)(other_ghg_names), gas_data.gas), :lifetimes])
    update_param!(m, :other_ghg_cycles, :emiss_other_ghg, Matrix(rcp_emissions[!,Symbol.(other_ghg_names)]))
    update_param!(m, :other_ghg_cycles, :emiss2conc_other_ghg, conversions[findall((in)(other_ghg_names), conversions.gases), :emiss2conc])

    # ---- Methane Radiative Forcing ---- #
    update_param!(m, :ch4_rf, :H₂O_share, 0.12)
    update_param!(m, :ch4_rf, :scale_CH₄, 1.0)
    update_param!(m, :ch4_rf, :a₃, -1.3e-6)
    update_param!(m, :ch4_rf, :b₃, -8.2e-6)

    # ---- Nitrous Oxide Radiative Forcing ---- #
    update_param!(m, :n2o_rf, :a₂, -8.0e-6)
    update_param!(m, :n2o_rf, :b₂, 4.2e-6)
    update_param!(m, :n2o_rf, :c₂, -4.9e-6)

    # ---- Carbon Dioxide Radiative Forcing ---- #
    update_param!(m, :co2_rf, :a₁, -2.4e-7)
    update_param!(m, :co2_rf, :b₁, 7.2e-4)
    update_param!(m, :co2_rf, :c₁, -2.1e-4)
    update_param!(m, :co2_rf, :rf_scale_CO₂, scale_co2_forcing)

    # ---- Other Well-Mixed Greenhouse Gas Radiative Forcings ---- #
    update_param!(m, :other_ghg_rf, :radiative_efficiency, gas_data[findall((in)(other_ghg_names), gas_data.gas), :rad_eff])

    # ---- Tropospheric Ozone Radiative Forcing ---- #
    update_param!(m, :trop_o3_rf, :NOx_emissions, rcp_emissions.NOx)
    update_param!(m, :trop_o3_rf, :CO_emissions, rcp_emissions.CO)
    update_param!(m, :trop_o3_rf, :NMVOC_emissions, rcp_emissions.NMVOC)
    update_param!(m, :trop_o3_rf, :mol_weight_NO, gas_data[gas_data.gas .== "NO", :mol_weight][1])
    update_param!(m, :trop_o3_rf, :T0, 0.0)

    # ---- Stratospheric Ozone Radiative Forcing ---- #
    update_param!(m, :strat_o3_rf, :Br, gas_data[findall((in)(ods_names), gas_data.gas), :br_atoms])
    update_param!(m, :strat_o3_rf, :Cl, gas_data[findall((in)(ods_names), gas_data.gas), :cl_atoms])
    update_param!(m, :strat_o3_rf, :FC, gas_data[findall((in)(ods_names), gas_data.gas), :strat_frac])
    update_param!(m, :strat_o3_rf, :δ1, -1.46030698e-5)
    update_param!(m, :strat_o3_rf, :δ2, 2.05401270e-3)
    update_param!(m, :strat_o3_rf, :δ3, 1.03143308)
    update_param!(m, :strat_o3_rf, :ODS₀, gas_data[findall((in)(ods_names), gas_data.gas), :pi_conc])

    # ---- Aerosol Direct Radiative Forcing ---- #
    update_param!(m, :aerosol_direct_rf, :β_SOx, -6.2227e-3)
    update_param!(m, :aerosol_direct_rf, :β_CO, 0.0)
    update_param!(m, :aerosol_direct_rf, :β_NMVOC, -3.8392e-4)
    update_param!(m, :aerosol_direct_rf, :β_NOx, -1.16551e-3)
    update_param!(m, :aerosol_direct_rf, :β_BC, 1.601537e-2)
    update_param!(m, :aerosol_direct_rf, :β_OC, -1.45339e-3)
    update_param!(m, :aerosol_direct_rf, :β_NH3, -1.55605e-3)
    update_param!(m, :aerosol_direct_rf, :rf_scale_aero_direct, 0.0)
    update_param!(m, :aerosol_direct_rf, :CO_emiss, rcp_emissions.CO)
    update_param!(m, :aerosol_direct_rf, :NMVOC_emiss, rcp_emissions.NMVOC)
    update_param!(m, :aerosol_direct_rf, :NH3_emiss, rcp_emissions.NH3)

    # ---- Aerosol Indirect Radiative Forcing ---- #
    update_param!(m, :aerosol_indirect_rf, :ϕ, -1.95011431)
    update_param!(m, :aerosol_indirect_rf, :b_SOx, 0.01107147)
    update_param!(m, :aerosol_indirect_rf, :b_POM, 0.01387492)
    update_param!(m, :aerosol_indirect_rf, :rf_scale_aero_indirect, 0.0)
    update_param!(m, :aerosol_indirect_rf, :model_years, collect(start_year:end_year))
    update_param!(m, :aerosol_indirect_rf, :SOx_emiss_1765, 1.0)
    update_param!(m, :aerosol_indirect_rf, :BC_OC_emiss_1765, 11.2)
    update_param!(m, :aerosol_indirect_rf, :scale_AR5, true)
    update_param!(m, :aerosol_indirect_rf, :F_1765, -0.3002836449793625)
    update_param!(m, :aerosol_indirect_rf, :F_2011, -1.5236182344467388)

    # ---- Black Carbon on Snow Radiative Forcing ---- #

    # ---- Land Use Change Radiative Forcing ---- #
    update_param!(m, :landuse_rf, :landuse_emiss, rcp_emissions.OtherCO2)

    # ---- Contrails Radiative Forcing ---- #
    update_param!(m, :contrails_rf, :frac, gas_fractions.nox_aviation)
    update_param!(m, :contrails_rf, :E_ref, 2.946)
    update_param!(m, :contrails_rf, :F_ref, 0.0448)
    update_param!(m, :contrails_rf, :ref_is_NO2, true)
    update_param!(m, :contrails_rf, :mol_weight_NO₂, gas_data[gas_data.gas .== "NO2", :mol_weight][1])

    # ---- Total Radiative Forcing ---- #
    update_param!(m, :total_rf, :F_volcanic, volcano_forcing)
    update_param!(m, :total_rf, :F_solar, solar_forcing)
    update_param!(m, :total_rf, :F_exogenous, zeros(n_steps))
    update_param!(m, :total_rf, :efficacy_CO₂, 1.0)
    update_param!(m, :total_rf, :efficacy_CH₄, 1.0)
    update_param!(m, :total_rf, :efficacy_CH₄_H₂O, 1.0)
    update_param!(m, :total_rf, :efficacy_N₂O, 1.0)
    update_param!(m, :total_rf, :efficacy_other_ghg, ones(length(other_ghg_names)))
    update_param!(m, :total_rf, :efficacy_trop_O₃, 1.0)
    update_param!(m, :total_rf, :efficacy_strat_O₃, 1.0)
    update_param!(m, :total_rf, :efficacy_aerosol_direct, 1.0)
    update_param!(m, :total_rf, :efficacy_aerosol_indirect, 1.0)
    update_param!(m, :total_rf, :efficacy_bcsnow, 3.0)
    update_param!(m, :total_rf, :efficacy_landuse, 1.0)
    update_param!(m, :total_rf, :efficacy_contrails, 0.0) # Note: Efficacy set to 0.0 to match default settings in Python version of FAIR.

    # ---- Global Temperature Anomaly ---- #
    update_param!(m, :temperature, :d, d)
    update_param!(m, :temperature, :q, q)
    update_param!(m, :temperature, :F2x, F2x)

    # --- Set shared parameters ---- #

    # set_param!(m, :CO₂_0, gas_data[gas_data.gas .== "CO2", :pi_conc][1])
    add_shared_param!(m, :CO₂_0, gas_data[gas_data.gas .== "CO2", :pi_conc][1])
    connect_param!(m, :co2_rf, :CO₂_0, :CO₂_0)
    connect_param!(m, :n2o_rf, :CO₂_0, :CO₂_0)

    # set_param!(m, :CH₄_0, gas_data[gas_data.gas .== "CH4", :pi_conc][1])
    add_shared_param!(m, :CH₄_0, gas_data[gas_data.gas .== "CH4", :pi_conc][1])
    connect_param!(m, :ch4_cycle, :CH₄_0, :CH₄_0)
    connect_param!(m, :ch4_rf, :CH₄_0, :CH₄_0)
    connect_param!(m, :n2o_rf, :CH₄_0, :CH₄_0)
    connect_param!(m, :trop_o3_rf, :CH₄_0, :CH₄_0)

    # set_param!(m, :N₂O_0, gas_data[gas_data.gas .== "N2O", :pi_conc][1])
    add_shared_param!(m, :N₂O_0, gas_data[gas_data.gas .== "N2O", :pi_conc][1])
    connect_param!(m, :ch4_rf, :N₂O_0, :N₂O_0)
    connect_param!(m, :co2_rf, :N₂O_0, :N₂O_0)
    connect_param!(m, :n2o_rf, :N₂O_0, :N₂O_0)
    connect_param!(m, :n2o_cycle, :N₂O_0, :N₂O_0)

    # set_param!(m, :other_ghg_0, gas_data[findall((in)(other_ghg_names), gas_data.gas), :pi_conc])
    add_shared_param!(m, :other_ghg_0, gas_data[findall((in)(other_ghg_names), gas_data.gas), :pi_conc], dims = [:other_ghg])
    connect_param!(m, :other_ghg_cycles, :other_ghg_0, :other_ghg_0)
    connect_param!(m, :other_ghg_rf, :other_ghg_0, :other_ghg_0)

    # set_param!(m, :SOx_emiss, rcp_emissions.SOx)
    add_shared_param!(m, :SOx_emiss, rcp_emissions.SOx, dims = [:time])
    connect_param!(m, :aerosol_direct_rf, :SOx_emiss, :SOx_emiss)
    connect_param!(m, :aerosol_indirect_rf, :SOx_emiss, :SOx_emiss)

    # set_param!(m, :BC_emiss, rcp_emissions.BC)
    add_shared_param!(m, :BC_emiss, rcp_emissions.BC, dims = [:time])
    connect_param!(m, :aerosol_direct_rf, :BC_emiss, :BC_emiss)
    connect_param!(m, :aerosol_indirect_rf, :BC_emiss, :BC_emiss)
    connect_param!(m, :bc_snow_rf, :BC_emiss, :BC_emiss)

    # set_param!(m, :OC_emiss, rcp_emissions.OC)
    add_shared_param!(m, :OC_emiss, rcp_emissions.OC, dims = [:time])
    connect_param!(m, :aerosol_direct_rf, :OC_emiss, :OC_emiss)
    connect_param!(m, :aerosol_indirect_rf, :OC_emiss, :OC_emiss)

    # set_param!(m, :NOx_emiss, rcp_emissions.NOx)
    add_shared_param!(m, :NOx_emiss, rcp_emissions.NOx, dims = [:time])
    connect_param!(m, :aerosol_direct_rf, :NOx_emiss, :NOx_emiss)
    connect_param!(m, :contrails_rf, :NOx_emiss, :NOx_emiss)

    # set_param!(m, :fix_pre1850_RCP, true)
    add_shared_param!(m, :fix_pre1850_RCP, true; data_type = Bool)
    connect_param!(m, :aerosol_indirect_rf, :fix_pre1850_RCP, :fix_pre1850_RCP)
    connect_param!(m, :trop_o3_rf, :fix_pre1850_RCP, :fix_pre1850_RCP)

    # set_param!(m, :mol_weight_N, gas_data[gas_data.gas .== "N", :mol_weight][1])
    add_shared_param!(m, :mol_weight_N, gas_data[gas_data.gas .== "N", :mol_weight][1])
    connect_param!(m, :contrails_rf, :mol_weight_N, :mol_weight_N)
    connect_param!(m, :trop_o3_rf, :mol_weight_N, :mol_weight_N)

    # set_param!(m, :gtc2ppm, conversions[conversions.gases .== "CO2", :emiss2conc][1])
    add_shared_param!(m, :gtc2ppm, conversions[conversions.gases .== "CO2", :emiss2conc][1])
    connect_param!(m, :ch4_cycle, :gtc2ppm, :gtc2ppm)
    connect_param!(m, :co2_cycle, :gtc2ppm, :gtc2ppm)

    # ---------------------------------------------
    # Create connections between Mimi components.
    # ---------------------------------------------
    connect_param!(m, :co2_cycle, :T, :temperature, :T)
    connect_param!(m, :co2_cycle, :E_ox_CH₄, :ch4_cycle, :oxidised_CH₄_GtC)
    connect_param!(m, :trop_o3_rf, :CH₄, :ch4_cycle, :CH₄)
    connect_param!(m, :strat_o3_rf, :conc_ODS, :other_ghg_cycles, :conc_ods)
    connect_param!(m, :ch4_rf, :CH₄, :ch4_cycle, :CH₄)
    connect_param!(m, :ch4_rf, :N₂O, :n2o_cycle, :N₂O)
    connect_param!(m, :n2o_rf, :N₂O, :n2o_cycle, :N₂O)
    connect_param!(m, :n2o_rf, :CH₄, :ch4_cycle, :CH₄)
    connect_param!(m, :n2o_rf, :CO₂, :co2_cycle, :C)
    connect_param!(m, :co2_rf, :N₂O, :n2o_cycle, :N₂O)
    connect_param!(m, :co2_rf, :CO₂, :co2_cycle, :C)
    connect_param!(m, :other_ghg_rf, :conc_other_ghg, :other_ghg_cycles, :conc_other_ghg)
    connect_param!(m, :trop_o3_rf, :temperature, :temperature, :T)
    connect_param!(m, :total_rf, :F_CO₂, :co2_rf, :rf_co2)
    connect_param!(m, :total_rf, :F_CH₄, :ch4_rf, :forcing_CH₄)
    connect_param!(m, :total_rf, :F_CH₄_H₂O, :ch4_rf, :forcing_CH₄_H₂O)
    connect_param!(m, :total_rf, :F_N₂O, :n2o_rf, :forcing_N₂O)
    connect_param!(m, :total_rf, :F_other_ghg, :other_ghg_rf, :other_ghg_rf)
    connect_param!(m, :total_rf, :F_trop_O₃, :trop_o3_rf, :forcing_trop_O₃)
    connect_param!(m, :total_rf, :F_strat_O₃, :strat_o3_rf, :forcing_strat_O₃)
    connect_param!(m, :total_rf, :F_aerosol_direct, :aerosol_direct_rf, :F_aerosol_direct)
    connect_param!(m, :total_rf, :F_aerosol_indirect, :aerosol_indirect_rf, :ERF_aero_cloud)
    connect_param!(m, :total_rf, :F_bcsnow, :bc_snow_rf, :forcing_BC_snow)
    connect_param!(m, :total_rf, :F_landuse, :landuse_rf, :forcing_landuse)
    connect_param!(m, :total_rf, :F_contrails, :contrails_rf, :forcing_contrails)
    connect_param!(m, :temperature, :F, :total_rf, :total_forcing)
end

"""
    update_MimiFAIR162_params!(m; ar6_scenario::String="ssp245", start_year::Int=1750, end_year::Int=2300)

Update all parameter settings in the model `m` with default FAIR  v1.2.6 settings
as pulled from the MimiFAIRv1_6_2 repository.  Optional arguments include `ar6_scenario` to indicate
what scenario to use, `start_year` and `end_year` to indicate the model time dimension.
"""
function update_MimiFAIR162_params!(m; ar6_scenario::String="ssp245", start_year::Int=1750, end_year::Int=2300)

    # ---------------------------------------------
    # Set Up Data and Parameter Values
    # ---------------------------------------------

    # Load RCP and other data needed to construct FAIR (just hard-coded 1765 as start year because not using the RCP emissions from this function).
    rcp_emissions, volcano_forcing, solar_forcing, gas_data, gas_fractions, conversions = MimiFAIRv1_6_2.load_fair_data(1765, end_year, "RCP85")

    # Load IPCC AR6 emissions scenario used for FAIRv1.6.2 ensemble runs (options = "ssp119", "ssp126", "ssp245", "ssp370", "ssp460", "ssp585").
    ar6_emissions_raw = DataFrame(load(joinpath(@__DIR__, "..", "data", "model_data", "AR6_emissions_"*ar6_scenario*"_1750_2300.csv")))

    # Subset AR6 emissions to proper years.
    emission_indices = indexin(collect(start_year:end_year), ar6_emissions_raw.Year)
    ar6_emissions = ar6_emissions_raw[emission_indices, :]

    # Load IPCC AR6 natural CH₄ and N₂O emissions for FAIR (spans 1750-2500).
    ar6_natural_emissions_raw = DataFrame(load(joinpath(@__DIR__, "..", "data", "model_data", "fair_wg3_natural_ch4_n2o.csv")))

    # Subset natural CH₄ and N₂O emissions to proper years.
    natural_indices = indexin(collect(start_year:end_year), ar6_natural_emissions_raw.year)
    ar6_natural_emissions = ar6_natural_emissions_raw[natural_indices, [:CH4, :N2O]]

    # Load IPCC solar forcing scenario (note, dataset runs from -6755 to 2299).
    ar6_solar_forcing_raw = DataFrame(load(joinpath(@__DIR__, "..", "data", "model_data", "ar6_solar_erf.csv")))

    # For cases when you want to run the model out to 2300, assume solar forcing value in 2300 follows the trend from 2298-2299.
    push!(ar6_solar_forcing_raw, [2300 ar6_solar_forcing_raw[end,"solar_erf"] + (ar6_solar_forcing_raw[end,"solar_erf"] - ar6_solar_forcing_raw[(end-1),"solar_erf"])])

    # Subset solar forcing data to proper years.
    solar_indices = indexin(collect(start_year:end_year), ar6_solar_forcing_raw.year)
    ar6_solar_forcing = ar6_solar_forcing_raw[solar_indices, :solar_erf]

    # Load IPCC AR6 volcanic forcing scenario (note, dataset runs from -500 to 2019).
    ar6_volcanic_forcing_raw = DataFrame(load(joinpath(@__DIR__, "..", "data", "model_data", "ar6_volcanic_erf.csv")))

    # Extract indices for relevant years.
    if end_year > 2019
        volcanic_indices = indexin(collect(start_year:2019), ar6_volcanic_forcing_raw.year)
    else
        volcanic_indices = indexin(collect(start_year:end_year), ar6_volcanic_forcing_raw.year)
    end

    # Create an empty array and store subset of volcanic forcings.
    ar6_volcanic_forcing = zeros(length(start_year:end_year))
    ar6_volcanic_forcing[1:length(volcanic_indices)] = ar6_volcanic_forcing_raw[volcanic_indices, :volcanic_erf]

    # From FAIR AR6 code on volcanic forcing: "ramp down last 10 years to zero according to https://www.geosci-model-dev.net/9/3461/2016/gmd-9-3461-2016.html"
    # This copies that code exactly.
    index_2019 = findfirst(x->x==2019, start_year:end_year)
    ar6_volcanic_forcing[index_2019:(index_2019+10)] = ar6_volcanic_forcing[index_2019] .* collect(range(1,0,length=11))

    # Names of minor greenhouse gases and ozone-depleting substances (used or indexing).
    other_ghg_names = ["CF4", "C2F6", "C6F14", "HFC23", "HFC32", "HFC43_10", "HFC125", "HFC134a", "HFC143a", "HFC227ea", "HFC245fa", "SF6"]
    ods_names       = ["CFC_11", "CFC_12", "CFC_113", "CFC_114", "CFC_115", "CARB_TET", "MCF", "HCFC_22", "HCFC_141B", "HCFC_142B", "HALON1211", "HALON1202", "HALON1301", "HALON2402", "CH3BR", "CH3CL"]

    # ---------------------------------------------
    # Set component-specific parameters
    # ---------------------------------------------

    # ---- Carbon Cycle ---- #
    update_param!(m, :co2_cycle, :CO₂_0,  278.052989189439) # From FAIR model run.
    update_param!(m, :co2_cycle, :iirf_h,  100.0)
    update_param!(m, :co2_cycle, :r0_co2, 35.0)
    update_param!(m, :co2_cycle, :rT_co2, 4.165)
    update_param!(m, :co2_cycle, :rC_co2,  0.019)
    update_param!(m, :co2_cycle, :τ_co2, [1000000, 394.4, 36.54, 4.304])
    update_param!(m, :co2_cycle, :a_co2, [0.2173,0.2240,0.2824,0.2763])
    update_param!(m, :co2_cycle, :R0_co2, [0.0003062168651584551, 0.0003156584344017209, 0.0003979550976564552, 0.0003893590420767655]) # From FAIR model run.
    update_param!(m, :co2_cycle, :E_co2, ar6_emissions.FossilCO2 .+ ar6_emissions.OtherCO2)
    update_param!(m, :co2_cycle, :cumulative_emissions_CO2₀, 0.003)
    update_param!(m, :co2_cycle, :airborne_emissions_CO2₀, 0.0)
    update_param!(m, :co2_cycle, :iIRF_max, 97.0)
    connect_param!(m, :co2_cycle => :temperature, :temperature => :T)

    # ---- Methane Cycle ---- #
    update_param!(m, :ch4_cycle, :fossil_emiss_CH₄, ar6_emissions.CH4)
    update_param!(m, :ch4_cycle, :natural_emiss_CH₄, ar6_natural_emissions.CH4)
    update_param!(m, :ch4_cycle, :τ_CH₄, 9.3)
    update_param!(m, :ch4_cycle, :fossil_frac, ones(length(start_year:end_year)))
    update_param!(m, :ch4_cycle, :oxidation_frac, 0.61)
    update_param!(m, :ch4_cycle, :mol_weight_CH₄, gas_data[gas_data.gas .== "CH4", :mol_weight][1])
    update_param!(m, :ch4_cycle, :mol_weight_C, gas_data[gas_data.gas .== "C", :mol_weight][1])
    update_param!(m, :ch4_cycle, :emiss2conc_ch4, conversions[conversions.gases .== "CH4", :emiss2conc][1])
    update_param!(m, :ch4_cycle, :CH₄_0, gas_data[gas_data.gas .== "CH4", :pi_conc_ar6][1])

    # ---- Nitrous Oxide Cycle ---- #
    update_param!(m, :n2o_cycle, :fossil_emiss_N₂O, ar6_emissions.N2O)
    update_param!(m, :n2o_cycle, :natural_emiss_N₂O, ar6_natural_emissions.N2O)
    update_param!(m, :n2o_cycle, :τ_N₂O, 121.0)
    update_param!(m, :n2o_cycle, :emiss2conc_n2o, conversions[conversions.gases .== "N2O", :emiss2conc][1])
    update_param!(m, :n2o_cycle, :N₂O_0, gas_data[gas_data.gas .== "N2O", :pi_conc_ar6][1])

    # ---- Other Well-Mixed Greenhouse Gas Cycles ---- #
    update_param!(m, :other_ghg_cycles, :τ_other_ghg, gas_data[findall((in)(other_ghg_names), gas_data.gas), :lifetimes])
    update_param!(m, :other_ghg_cycles, :emiss_other_ghg, Matrix(ar6_emissions[!,Symbol.(other_ghg_names)]))
    update_param!(m, :other_ghg_cycles, :emiss2conc_other_ghg, conversions[findall((in)(other_ghg_names), conversions.gases), :emiss2conc])
    update_param!(m, :other_ghg_cycles, :other_ghg_0, gas_data[findall((in)(other_ghg_names), gas_data.gas), :pi_conc_ar6])

    # ---- Ozone-Depleting Substance Gas Cycles ---- #
    update_param!(m, :o3_depleting_substance_cycles, :τ_ods, gas_data[findall((in)(ods_names), gas_data.gas), :lifetimes])
    update_param!(m, :o3_depleting_substance_cycles, :emiss_ods, Matrix(ar6_emissions[!,Symbol.(ods_names)]))
    update_param!(m, :o3_depleting_substance_cycles, :emiss2conc_ods, conversions[findall((in)(ods_names), conversions.gases), :emiss2conc])
    update_param!(m, :o3_depleting_substance_cycles, :ods_0, gas_data[findall((in)(ods_names), gas_data.gas), :pi_conc_ar6])

    # ---- Carbon Dioxide Radiative Forcing ---- #
    update_param!(m, :co2_forcing, :F2x, 3.71)
    update_param!(m, :co2_forcing, :a₁, -2.4785e-07)
    update_param!(m, :co2_forcing, :b₁, 0.00075906)
    update_param!(m, :co2_forcing, :c₁, -0.0021492)
    update_param!(m, :co2_forcing, :d₁, 5.2488)
    update_param!(m, :co2_forcing, :adjust_F2x, true)

    connect_param!(m, :co2_forcing => :CO₂, :co2_cycle => :co2)
    connect_param!(m, :co2_forcing => :N₂O, :n2o_cycle => :N₂O)

    # ---- Methane Radiative Forcing ---- #
    update_param!(m, :ch4_forcing, :a₃, -8.9603e-05)
    update_param!(m, :ch4_forcing, :b₃, -0.00012462)
    update_param!(m, :ch4_forcing, :d₃, 0.045194)
    update_param!(m, :ch4_forcing, :h2o_from_ch4, 0.079047)

    connect_param!(m, :ch4_forcing => :N₂O, :n2o_cycle => :N₂O)
    connect_param!(m, :ch4_forcing => :CH₄, :ch4_cycle => :CH₄)

    # ---- Nitrous Oxide Radiative Forcing ---- #
    update_param!(m, :n2o_forcing, :a₂, -0.00034197)
    update_param!(m, :n2o_forcing, :b₂,  0.00025455)
    update_param!(m, :n2o_forcing, :c₂, -0.00024357)
    update_param!(m, :n2o_forcing, :d₂, 0.12173)

    connect_param!(m, :n2o_forcing => :CO₂, :co2_cycle => :co2)
    connect_param!(m, :n2o_forcing => :N₂O, :n2o_cycle => :N₂O)
    connect_param!(m, :n2o_forcing => :CH₄, :ch4_cycle => :CH₄)

    # ---- Ozone Radiative Forcing ---- #
    update_param!(m, :o3_forcing, :total_forcing_O₃_0, 0.0)
    update_param!(m, :o3_forcing, :Br, gas_data[findall((in)(ods_names), gas_data.gas), :br_atoms])
    update_param!(m, :o3_forcing, :Cl, gas_data[findall((in)(ods_names), gas_data.gas), :cl_atoms])
    update_param!(m, :o3_forcing, :FC, gas_data[findall((in)(ods_names), gas_data.gas), :strat_frac])
    update_param!(m, :o3_forcing, :feedback, -0.037)
    update_param!(m, :o3_forcing, :Ψ_CH₄, 2.33379720e-04)
    update_param!(m, :o3_forcing, :Ψ_N₂O, 1.27179106e-03)
    update_param!(m, :o3_forcing, :Ψ_ODS, -6.69347820e-05)
    update_param!(m, :o3_forcing, :Ψ_CO, 1.14647701e-04)
    update_param!(m, :o3_forcing, :Ψ_NMVOC, 5.14366051e-12)
    update_param!(m, :o3_forcing, :Ψ_NOx, 3.78354423e-03)
    update_param!(m, :o3_forcing, :CO_emiss_pi,    348.527359)
    update_param!(m, :o3_forcing, :NMVOC_emiss_pi, 60.0218262)
    update_param!(m, :o3_forcing, :NOx_emiss_pi,   3.87593407)

    connect_param!(m, :o3_forcing => :N₂O, :n2o_cycle => :N₂O)
    connect_param!(m, :o3_forcing => :CH₄, :ch4_cycle => :CH₄)
    connect_param!(m, :o3_forcing => :temperature, :temperature => :T)
    connect_param!(m, :o3_forcing => :conc_ODS, :o3_depleting_substance_cycles => :conc_ods)

    # ---- Aerosol Direct Radiative Forcing ---- #
    update_param!(m, :aerosol_direct_forcing, :β_SOx, -6.2227e-3)
    update_param!(m, :aerosol_direct_forcing, :β_CO, 0.0)
    update_param!(m, :aerosol_direct_forcing, :β_NMVOC, -3.8392e-4)
    update_param!(m, :aerosol_direct_forcing, :β_NOx, -1.16551e-3)
    update_param!(m, :aerosol_direct_forcing, :β_BC, 1.601537e-2)
    update_param!(m, :aerosol_direct_forcing, :β_OC, -1.45339e-3)
    update_param!(m, :aerosol_direct_forcing, :β_NH3, -1.55605e-3)
    update_param!(m, :aerosol_direct_forcing, :NH3_emiss, ar6_emissions.NH3)

    # ---- Aerosol Indirect Radiative Forcing ---- #
    update_param!(m, :aerosol_indirect_forcing, :ϕ, 0.07334277994353743)#-1.95011431)
    update_param!(m, :aerosol_indirect_forcing, :b_SOx, 3.452849302362568)#0.01107147)
    update_param!(m, :aerosol_indirect_forcing, :b_POM, 33.126485122209154)#0.01387492)
    update_param!(m, :aerosol_indirect_forcing, :rf_scale_aero_indirect, 1.0)
    update_param!(m, :aerosol_indirect_forcing, :SOx_emiss_pi,   1.22002422)
    update_param!(m, :aerosol_indirect_forcing, :BC_emiss_pi,    2.09777075)
    update_param!(m, :aerosol_indirect_forcing, :OC_emiss_pi,    15.4476682)

    # ---- Other Well-Mixed Greenhouse Gas Radiative Forcings ---- #
    update_param!(m, :other_ghg_forcing, :other_ghg_radiative_efficiency, gas_data[findall((in)(other_ghg_names), gas_data.gas), :rad_eff])
    update_param!(m, :other_ghg_forcing, :other_ghg_pi, gas_data[findall((in)(other_ghg_names), gas_data.gas), :pi_conc_ar6])

    connect_param!(m, :other_ghg_forcing => :conc_other_ghg, :other_ghg_cycles => :conc_other_ghg)

    # ---- Ozone-Depleting Substance Radiative Forcings ---- #
    update_param!(m, :o3_depleting_substance_forcing, :ods_radiative_efficiency, gas_data[findall((in)(ods_names), gas_data.gas), :rad_eff])

    connect_param!(m, :o3_depleting_substance_forcing => :conc_ods, :o3_depleting_substance_cycles => :conc_ods)

    # ---- Contrails Radiative Forcing ---- #
    update_param!(m, :contrails_forcing, :frac, zeros(length(start_year:end_year)))
    update_param!(m, :contrails_forcing, :E_ref_contrails, 2.946)
    update_param!(m, :contrails_forcing, :F_ref_contrails, 0.0448)
    update_param!(m, :contrails_forcing, :ref_is_NO2, true)
    update_param!(m, :contrails_forcing, :mol_weight_NO₂, gas_data[gas_data.gas .== "NO2", :mol_weight][1])
    update_param!(m, :contrails_forcing, :mol_weight_N, gas_data[gas_data.gas .== "N", :mol_weight][1])

    # ---- Black Carbon on Snow Radiative Forcing ---- #
    update_param!(m, :bc_snow_forcing, :E_ref_bc, 6.095)
    update_param!(m, :bc_snow_forcing, :F_ref_bc, 0.08)

    # ---- Land Use Change Radiative Forcing ---- #
    update_param!(m, :landuse_forcing, :α_CO₂_land, (-0.2/190))
    update_param!(m, :landuse_forcing, :landuse_emiss, ar6_emissions.OtherCO2)

    # ---- Total Radiative Forcing ---- #
    update_param!(m, :total_forcing, :scale_CO₂, 1.0)
    update_param!(m, :total_forcing, :scale_CH₄, 1.0)
    update_param!(m, :total_forcing, :scale_CH₄_H₂O, 1.0)
    update_param!(m, :total_forcing, :scale_N₂O, 1.0)
    update_param!(m, :total_forcing, :scale_O₃, 1.0)
    update_param!(m, :total_forcing, :scale_aerosol_indirect, 1.0)
    update_param!(m, :total_forcing, :scale_bcsnow, 1.0)
    update_param!(m, :total_forcing, :scale_landuse, 1.0)
    update_param!(m, :total_forcing, :scale_contrails, 0.0) # Default FAIR has contrail forcing switched off. Set scaling term to 0
    update_param!(m, :total_forcing, :scale_volcanic, 1.0)
    update_param!(m, :total_forcing, :scale_solar, 1.0)
    update_param!(m, :total_forcing, :scale_aerosol_direct_SOx, 1.0)
    update_param!(m, :total_forcing, :scale_aerosol_direct_CO_NMVOC, 1.0)
    update_param!(m, :total_forcing, :scale_aerosol_direct_NOx_NH3, 1.0)
    update_param!(m, :total_forcing, :scale_aerosol_direct_BC, 1.0)
    update_param!(m, :total_forcing, :scale_aerosol_direct_OC, 1.0)
    update_param!(m, :total_forcing, :scale_other_ghg, ones(length(other_ghg_names)))
    update_param!(m, :total_forcing, :scale_ods, ones(length(ods_names)))
    update_param!(m, :total_forcing, :F_volcanic, ar6_volcanic_forcing)
    update_param!(m, :total_forcing, :F_solar, ar6_solar_forcing)
    update_param!(m, :total_forcing, :F_exogenous, zeros(length(start_year:end_year)))

    connect_param!(m, :total_forcing => :F_CO₂, :co2_forcing => :rf_co2)
    connect_param!(m, :total_forcing => :F_CH₄, :ch4_forcing => :rf_ch4)
    connect_param!(m, :total_forcing => :F_CH₄_H₂O, :ch4_forcing => :rf_ch4_h2o)
    connect_param!(m, :total_forcing => :F_N₂O, :n2o_forcing => :rf_n2o)
    connect_param!(m, :total_forcing => :F_O₃, :o3_forcing => :total_forcing_O₃)
    connect_param!(m, :total_forcing => :F_aerosol_direct_SOx, :aerosol_direct_forcing => :F_SOx_aero)
    connect_param!(m, :total_forcing => :F_aerosol_direct_CO, :aerosol_direct_forcing => :F_CO_aero)
    connect_param!(m, :total_forcing => :F_aerosol_direct_NMVOC, :aerosol_direct_forcing => :F_NMVOC_aero)
    connect_param!(m, :total_forcing => :F_aerosol_direct_NOx, :aerosol_direct_forcing => :F_NOx_aero)
    connect_param!(m, :total_forcing => :F_aerosol_direct_BC, :aerosol_direct_forcing => :F_BC_aero)
    connect_param!(m, :total_forcing => :F_aerosol_direct_OC, :aerosol_direct_forcing => :F_OC_aero)
    connect_param!(m, :total_forcing => :F_aerosol_direct_NH3, :aerosol_direct_forcing => :F_NH3_aero)
    connect_param!(m, :total_forcing => :F_aerosol_indirect, :aerosol_indirect_forcing => :rf_aero_indirect)
    connect_param!(m, :total_forcing => :F_bcsnow, :bc_snow_forcing => :forcing_BC_snow)
    connect_param!(m, :total_forcing => :F_landuse, :landuse_forcing => :forcing_landuse)
    connect_param!(m, :total_forcing => :F_contrails, :contrails_forcing => :forcing_contrails)
    connect_param!(m, :total_forcing => :F_other_ghg, :other_ghg_forcing => :other_ghg_rf)
    connect_param!(m, :total_forcing => :F_ods, :o3_depleting_substance_forcing => :ods_rf)

    # ---- Temperature ---- #
    update_param!(m, :temperature, :earth_radius, 6371000)
    update_param!(m, :temperature, :seconds_per_year, (60*60*24*365.24219))
    update_param!(m, :temperature, :ocean_heat_exchange, 0.67)
    update_param!(m, :temperature, :deep_ocean_efficacy, 1.28)
    update_param!(m, :temperature, :lambda_global, 1.18)
    update_param!(m, :temperature, :T_mix₀, [5.30299681e-03, 6.41290105e-05]) # From Python-FAIR model run.
    update_param!(m, :temperature, :T_deep₀, [-1.33065374e-04,  1.50206328e-04]) # From Python-FAIR model run.
    update_param!(m, :temperature, :ocean_heat_capacity, [8.2, 109.0])
    connect_param!(m, :temperature => :forcing, :total_forcing => :total_forcing)

    # ---- Parameters Shared Across Multiple Components - Set shared parameters ---- #

    # set_param!(m, :dt, 1.0)
    add_shared_param!(m, :model_dt, 1.0)
    connect_param!(m, :co2_cycle,   :dt, :model_dt)
    connect_param!(m, :temperature, :dt, :model_dt)

    # set_param!(m, :emiss2conc_co2, conversions[conversions.gases .== "CO2", :emiss2conc][1])
    add_shared_param!(m, :model_emiss2conc_co2, conversions[conversions.gases .== "CO2", :emiss2conc][1])
    connect_param!(m, :co2_cycle, :emiss2conc_co2, :model_emiss2conc_co2)
    connect_param!(m, :ch4_cycle, :emiss2conc_co2, :model_emiss2conc_co2)

    # set_param!(m, :CH₄_pi, gas_data[gas_data.gas .== "CH4", :pi_conc_ar6][1])
    add_shared_param!(m, :model_CH₄_pi, gas_data[gas_data.gas .== "CH4", :pi_conc_ar6][1])
    connect_param!(m, :ch4_cycle,   :CH₄_pi, :model_CH₄_pi)
    connect_param!(m, :ch4_forcing, :CH₄_pi, :model_CH₄_pi)
    connect_param!(m, :o3_forcing,  :CH₄_pi, :model_CH₄_pi)

    # set_param!(m, :CO₂_pi, gas_data[gas_data.gas .== "CO2", :pi_conc_ar6][1])
    add_shared_param!(m, :model_CO₂_pi, gas_data[gas_data.gas .== "CO2", :pi_conc_ar6][1])
    connect_param!(m, :co2_cycle,   :CO₂_pi, :model_CO₂_pi)
    connect_param!(m, :co2_forcing, :CO₂_pi, :model_CO₂_pi)

    # set_param!(m, :N₂O_pi, gas_data[gas_data.gas .== "N2O", :pi_conc_ar6][1])
    add_shared_param!(m, :model_N₂O_pi, gas_data[gas_data.gas .== "N2O", :pi_conc_ar6][1])
    connect_param!(m, :co2_forcing, :N₂O_pi, :model_N₂O_pi)
    connect_param!(m, :n2o_forcing, :N₂O_pi, :model_N₂O_pi)
    connect_param!(m, :o3_forcing,  :N₂O_pi, :model_N₂O_pi)

    # set_param!(m, :ods_pi, gas_data[findall((in)(ods_names), gas_data.gas), :pi_conc_ar6])
    add_shared_param!(m, :model_ods_pi, gas_data[findall((in)(ods_names), gas_data.gas), :pi_conc_ar6], dims = [:ozone_depleting_substances])
    connect_param!(m, :o3_depleting_substance_forcing,  :ods_pi, :model_ods_pi)
    connect_param!(m, :o3_forcing,                      :ods_pi, :model_ods_pi)

    # set_param!(m, :SOx_emiss, ar6_emissions.SOx)
    add_shared_param!(m, :model_SOx_emiss, ar6_emissions.SOx, dims = [:time])
    connect_param!(m, :aerosol_direct_forcing,      :SOx_emiss, :model_SOx_emiss)
    connect_param!(m, :aerosol_indirect_forcing,    :SOx_emiss, :model_SOx_emiss)

    # set_param!(m, :BC_emiss, ar6_emissions.BC)
    add_shared_param!(m, :model_BC_emiss, ar6_emissions.BC, dims = [:time])
    connect_param!(m, :aerosol_direct_forcing,      :BC_emiss, :model_BC_emiss)
    connect_param!(m, :aerosol_indirect_forcing,    :BC_emiss, :model_BC_emiss)
    connect_param!(m, :bc_snow_forcing,             :BC_emiss, :model_BC_emiss)

    # set_param!(m, :OC_emiss, ar6_emissions.OC)
    add_shared_param!(m, :model_OC_emiss, ar6_emissions.OC, dims = [:time])
    connect_param!(m, :aerosol_direct_forcing,      :OC_emiss, :model_OC_emiss)
    connect_param!(m, :aerosol_indirect_forcing,    :OC_emiss, :model_OC_emiss)

    # set_param!(m, :CO_emiss, ar6_emissions.CO)
    add_shared_param!(m, :model_CO_emiss, ar6_emissions.CO, dims = [:time])
    connect_param!(m, :aerosol_direct_forcing,  :CO_emiss, :model_CO_emiss)
    connect_param!(m, :o3_forcing,              :CO_emiss, :model_CO_emiss)

    # set_param!(m, :NMVOC_emiss, ar6_emissions.NMVOC)
    add_shared_param!(m, :model_NMVOC_emiss, ar6_emissions.NMVOC, dims = [:time])
    connect_param!(m, :aerosol_direct_forcing,  :NMVOC_emiss, :model_NMVOC_emiss)
    connect_param!(m, :o3_forcing,              :NMVOC_emiss, :model_NMVOC_emiss)

    # set_param!(m, :NOx_emiss, ar6_emissions.NOx)
    add_shared_param!(m, :model_NOx_emiss, ar6_emissions.NOx, dims = [:time])
    connect_param!(m, :aerosol_direct_forcing,  :NOx_emiss, :model_NOx_emiss)
    connect_param!(m, :o3_forcing,              :NOx_emiss, :model_NOx_emiss)
    connect_param!(m, :contrails_forcing,       :NOx_emiss, :model_NOx_emiss)

end


text:
"""
    hermite(x, p::Integer)

Evaluate the Hermite polynomial of degree `p` at `x` using the three term recursion.
"""
function hermite(x, p::Integer)
    T = typeof( 2x )
    p₀ = one(T)
    p₁ = 2x

    if p <= 0
        return p₀
    elseif p == 1
        return p₁
    end

    for n in 2:p
        p₀, p₁ = p₁, 2x*p₁ - 2*(n-1)*p₀
    end

    p₁
end


text:
@testset "ProductArray" begin
    a = ProductArray(tuple, [1,2], [3.0,4.0])

    @test a isa ProductArray{Tuple{Int,Float64}, 2, typeof(tuple), Vector{Int}, Vector{Float64}}
    @test size(a) == (2,2)
    @test axes(a) === (Base.OneTo(2), Base.OneTo(2))
    @test a == [(1,3.0) (1,4.0); (2,3.0) (2,4.0)]
end


text:
using multimode
using QuantumOptics

E0 = 0.2
s = 0.25
Nparticles = 1
Nparticlemodes = 16

index1 = 12
cavitymodes1 = 14
delta1 = -4.
eta1 = 2.0
kappa1 = 1.
U0_1 = 0.

function QuantumOptics.qfunc(system)
    Ψ₀ = basis_ket(system.basis, 1)
    ρ₀ = Ψ₀⊗dagger(Ψ₀)
    Tss = [0,100]

    T, ρ_t = multimode.evolve_master(system, Tss, ρ₀)
    ρ_ss = ρ_t[end]

    rho_f = ptrace(ρ_ss, [1])
    rho_p = ptrace(ρ_ss, [2])

    x = linspace(-3,3,200)
    return qfunc(rho_f, x, x)
end

function n(system)
    Ψ₀ = basis_ket(system.basis, 1)
    ρ₀ = Ψ₀⊗dagger(Ψ₀)
    Tss = [0,20]

    T, ρ_t = multimode.evolve_master(system, Tss, ρ₀)
    ρ_ss = ρ_t[end]

    rho_f = ptrace(ρ_ss, [1])
    rho_p = ptrace(ρ_ss, [2])

    a = destroy(system.modes[1].basis)

    return expect(dagger(a)*a, rho_f)
end

#particle = Bosons(Nparticles, Nparticlemodes, E0, BoxPotential)
#eta1 = 3.0
#mode1 = CavityMode(index1, cavitymodes1, delta1, eta1, kappa1, U0_1)
#system = MultimodeSystem(s, particle, [mode1])
#Ψ₀ = basis_ket(system.basis, 1)
#println(n(system))
#println(multimode.evolve_mcwf(system, [0,1], Ψ₀))

# particle = Bosons(Nparticles, Nparticlemodes, E0, BoxPotential)
# eta1 = 3.0
# mode1 = CavityMode(index1, cavitymodes1, delta1, eta1, kappa1, U0_1)
# system = MultimodeSystem(s, particle, [mode1])
# M = qfunc(system)
# writecsv("qfunc_mode12_eta=3.0.csv", M)

# eta1 = 2.
#etas = [0:0.1:3.5]
deltas = [-3:0.1:3]
exp_n = Float64[]
particle = Bosons(Nparticles, Nparticlemodes, E0, BoxPotential)
for (i, delta1) in enumerate(deltas)
     println(delta1)
     mode1 = CavityMode(index1, cavitymodes1, delta1, eta1, kappa1, U0_1)
     system = MultimodeSystem(s, particle, [mode1])
     push!(exp_n, abs(n(system)))
end
#
# writecsv("expn_mode12.csv", exp_n)

using PyCall
@pyimport matplotlib.pyplot as plt

plt.plot(deltas, exp_n)
plt.show()