text:
#CSV, Serialize, Dates, Satistics, and Dataframes are essential for working with data.
#GZip is for gzip files. The other packages relate to drawing graphs using Gadfly.
using Revise
using  GZip, CSV, DataFrames, Serialization, Dates,Statistics, Cairo, Fontconfig, Gadfly

#define path names as constants
const DATAPATH = "data"
const OUTPUTPATH = "output"
const CRSPNAME = "CRSPData"
const CRSPEVENTNAME = "CRSPDataEvent"

#make sure we can load things using relative paths
if  pwd() ∉ LOAD_PATH
  push!(LOAD_PATH,pwd())
end

function runSolution(;refresh=true, buildsolution=true)::Nothing

  #load the data
  crsp::DataFrame = loadCRSPData("CRSPData", refresh=refresh)
  crspevent::DataFrame = loadCRSPData("CRSPDataEvent", refresh=refresh)

  #examine the data
#  show(describe(crsp), allrows=true, allcols=true)

  buildsolution && buildStudyPanel(crsp, crspevent)
  makegraph()

  return nothing
end

function loadCRSPData(crspname;datapath=DATAPATH, refresh=true)::DataFrame
  local p::String = "$datapath\\$crspname"
  local df::DataFrame

  #makes a serialization object so we don't have to keep parsing the CSV
  if refresh || (!isfile("$p.jls"))
    df = CSV.read("$p.csv", allowmissing=:all, silencewarnings=true)
    serialize("$p.jls", df)
  end

  #read in the dataframe
  df = deserialize("$p.jls")
  sort!(df, [:PERMNO])
  return df
end

function buildStudyPanel(crsp, crspevent)::Nothing

  local dateFormat::Dates.DateFormat = Dates.DateFormat("yyyymmdd")

  #get the total returns of each security
  crsp[:ret] = ((ret,dlret)->
    ismissing(ret) && ismissing(dlret) ? missing :
      ismissing(ret) ? dlret :
      ismissing(dlret) ? ret :
      (1+ret)*(1+dlret) - 1.
    ).(crsp[:RET],crsp[:DLRET])

  #drop rows with missing data
  crsp = crsp[completecases(crsp[[:ret, :SHRCD, :EXCHCD, :sprtrn]]),:]
  crspevent = crspevent[completecases(crspevent),:]

  #uncomment major exchanges only
  #crsp = crsp[((exchcd,shrcd)->
  #  (exchcd ∈ [1,2,3]) && (shrcd ∈ [10,11])).(crsp[:EXCHCD], crsp[:SHRCD]),:]

  #make actual dates
  crsp[:DATE] = (s::String->Dates.Date(s,dateFormat)).((string).(crsp[:date]))
  crspevent[:DATE] = (s::String->Dates.Date(s,dateFormat)).((string).(crspevent[:EXDT]))

  #form excess returns
  crsp[:excess] = crsp[:ret] .- crsp[:sprtrn]

  sort!(crsp,[:PERMNO,:DATE])
  sort!(crspevent, [:PERMNO, :DATE])

  #make fields to contain the event data
  local dayfields::Vector{Symbol} = (d->d<0 ? Symbol("TM",abs(d)) : Symbol("D",d)).(collect(-20:20))
  for f::Symbol ∈ dayfields
    crspevent[f] = Vector{Union{Float64,Missing}}(missing, size(crspevent,1))
  end

  #We want to be careful with the dates- crsp will not adjust the index
  #return for gaps in the holding period returns! (See R solution)
  crspdates::Vector{Date} = unique(crsp[:,:DATE])
  sort!(crspdates)
  validdateranges::Dict =
    Dict(crspdates[i]=>(crspdates[i-20], crspdates[i+20]) for i ∈ 21:(length(crspdates)-20))

  #we will use this to set the boundaries for the dates
  crspevent[:tokeep]=trues(size(crspevent,1))
  crspevent[:begindate] = Vector{Union{Missing,Date}}(missing, size(crspevent,1))
  crspevent[:enddate] = Vector{Union{Missing,Date}}(missing, size(crspevent,1))

  print("\nStarting # of events: $(size(crspevent,1))\n")

  #build a set of crsp events
  subevents::GroupedDataFrame = groupby(crspevent, :PERMNO)

  #create a lookup file mapping permnos to events
  subeventsindex::Dict = Dict((subevents[i])[1,:PERMNO] => i for i ∈ 1:length(subevents))

  permctr::Int = 0
  for subcrsp::SubDataFrame ∈ groupby(crsp, [:PERMNO])
    permctr += 1

    #this could take a while looking at so many sotcks
    if permctr % 1000 == 0
      println("$permctr permnos complete")
    end

    #make a lookup table mapping dates to crsp rows
    dateindex::Dict = Dict(subcrsp[i,:DATE]=>i for i ∈ 1:(size(subcrsp,1)))
    permno::Int = subcrsp[1,:PERMNO]

    #partition out the events for a single stock
    if haskey(subeventsindex, permno)
      subevent::SubDataFrame = subevents[subeventsindex[permno]]

      #make vectors for fast inner-loop access
      excesses::Vector{SubArray{Union{Float64,Missing},1}} =
        (i::Int->subevent[dayfields[i]]).(1:length(dayfields))
      tokeeps::SubArray{Bool,1} = subevent[:tokeep]
      excess::SubArray{Float64,1} = subcrsp[:excess]

      for i::Int ∈ 1:(size(subevent,1))

        local ind::Int #this willbe used as a row pointer later
        #make sure we are not too close to the ends, and the ex-div date is a trading day
        #strictly this may not be necessary, but it does keep a consistent timeline
        eventdate::Date = subevent[i,:DATE]
        if haskey(dateindex, eventdate) && haskey(validdateranges,eventdate)
          ind = dateindex[eventdate]
        else
          tokeeps[i] = false
        end

        #check the endpoints and make sure we capture the entire range
        #also look up the associated ate range against the valid date ranges
        if ((tokeeps[i]) && (ind > 20) && (ind + 20 ≤ size(subcrsp,1)) &&
          (validdateranges[eventdate] == (subcrsp[ind-20, :DATE], subcrsp[ind+20, :DATE])))

          subevent[i, :begindate] = subcrsp[ind-20, :DATE] #set the date boundaries to ensure no overlap
          subevent[i, :enddate] = subcrsp[ind+20, :DATE]

          #innermost loop- tell julia to fully optimize
          @fastmath @inbounds @simd for j ∈ 1:(length(dayfields))
            excesses[j][i] = excess[ind-21+j]
          end

          #now we check the boundaries
          if i > 1
            nocollision::Bool = (subevent[i-1,:DATE]<subevent[i, :begindate])
            tokeeps[i] =tokeeps[i] && nocollision
            tokeeps[i-1] = tokeeps[i-1] && nocollision
          end

          if i < length(tokeeps)
            nocollision = (subevent[i+1,:DATE]>subevent[i, :enddate])
            tokeeps[i] = tokeeps[i] && nocollision
            tokeeps[i+1] = tokeeps[i+1] && nocollision
          end
        else
          tokeeps[i] = false
        end
      end
    end
  end

  crspevent = crspevent[crspevent[:tokeep],:]
  crspevent = crspevent[completecases(crspevent[dayfields]),:]

  serialize("$DATAPATH\\results.jls", crspevent)

  print("Ending # of events: $(size(crspevent,1))\n")
  return nothing
end

function makegraph()::Nothing
  local dayfields::Vector{Symbol} = (d->d<0 ? Symbol("TM",abs(d)) : Symbol("D",d)).(collect(-20:20))

  #----
  crspevent::DataFrame = deserialize("$DATAPATH\\results.jls")

  meanexcess::DataFrame = DataFrame(day=collect(-20:20),
    excess=vec(Statistics.mean(Matrix{Float64}(crspevent[dayfields]),dims=1)).*10_000)

  p = plot(meanexcess, x=:day, y=:excess, Guide.xlabel("Day"), Guide.ylabel("Excess Returns (bp)"),
    Geom.line, Geom.point, Guide.title("Mean Total Returns Over S&P500 Near Ex-Div Dates"))

  draw(PNG("ExcessReturns.png",9inch, 7inch),p)
  return nothing
end

#to run the full code, set both of these to true
@time runSolution(refresh=false, buildsolution=false)

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
# deep copying

# Note: deepcopy_internal(::Any, ::ObjectIdDict) is
#       only exposed for specialization by libraries

deepcopy(x) = deepcopy_internal(x, ObjectIdDict())

deepcopy_internal(x::Union(Symbol,LambdaStaticData,TopNode,QuoteNode,
                           DataType,UnionType,Task),
                  stackdict::ObjectIdDict) = x
deepcopy_internal(x::Tuple, stackdict::ObjectIdDict) =
    ntuple(length(x), i->deepcopy_internal(x[i], stackdict))
deepcopy_internal(x::Module, stackdict::ObjectIdDict) = error("deepcopy of Modules not supported")

function deepcopy_internal(x::Function, stackdict::ObjectIdDict)
    if isa(x.env, Union(MethodTable, Symbol)) || x.env === ()
        return x
    end
    invoke(deepcopy_internal, (Any, ObjectIdDict), x, stackdict)
end

function deepcopy_internal(x, stackdict::ObjectIdDict)
    if haskey(stackdict, x)
        return stackdict[x]
    end
    _deepcopy_t(x, typeof(x), stackdict)
end

function _deepcopy_t(x, T::DataType, stackdict::ObjectIdDict)
    if T.names===() || !T.mutable
        return x
    end
    ret = ccall(:jl_new_struct_uninit, Any, (Any,), T)
    stackdict[x] = ret
    for i in 1:length(T.names)
        if isdefined(x,i)
            ret.(i) = deepcopy_internal(x.(i), stackdict)
        end
    end
    return ret
end

function deepcopy_internal(x::Array, stackdict::ObjectIdDict)
    if haskey(stackdict, x)
        return stackdict[x]
    end
    _deepcopy_array_t(x, eltype(x), stackdict)
end

function _deepcopy_array_t(x, T, stackdict::ObjectIdDict)
    if isbits(T)
        return copy(x)
    end
    dest = similar(x)
    stackdict[x] = dest
    for i=1:length(x)
        if isdefined(x,i)
            arrayset(dest, deepcopy_internal(x[i], stackdict), i)
        end
    end
    return dest
end

ext: jl

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
using RCall, JLD

R"
library(gridExtra)
library(cowplot)
library(tidyverse)
library(latex2exp)
library(viridis)
"

DF = load("/Users/dariabystrova/Documents/GitHub/approximatingBNPpriors/saves_for_figures/figure_6_5.jld")
EV__NGG_SB = DF["NGG_FSB"]

## SB truncation was appled on reduced grid due to computational issues

DF_NGG_SB= DataFrame(Exp = first.(EV__NGG_SB),Var =last.(EV__NGG_SB), beta = first.(grid_vec[11:100]), sigma = last.(grid_vec[11:100]))
DF_NGG_SB.std = sqrt.(DF_NGG_SB.Var)
DF_NGG_SB.beta_log = log.(DF_NGG_SB.beta)
DF_NGG_SB

R"
range01 <- function(x){(x-min(x))/(max(x)-min(x))}"
R"
df_r= $DF_NGG_SB
df_r$beta_scaled= range01(df_r$beta_log)
df_r$color = rgb(df_r$beta_scaled, df_r$sigma,0.5,1)
df_r$color_sigma = rgb(0, df_r$sigma,0.5,1)
df_r$color_beta = rgb(df_r$beta_scaled, 0,0.5,1)
"
R"
p <- df_r %>% ggplot(aes(x=Exp, y = std, group=sigma))  + geom_line( alpha = 0.8,linetype = 'longdash') + geom_point(aes(colour=color), size=4)
p_ngg_sb <- p  + geom_line(data =df_r, aes(x=Exp, y = std, group=beta), alpha = 0.9,linetype = 'dotted') +
         xlim(0, 100)+ ylim(0,15)+labs(y='Std', x='Expectation')+scale_colour_identity()+theme_classic()+ theme(plot.title = element_text(hjust = 0.5,size = 15), axis.text.x = element_text(size=15), axis.text.y = element_text(size=15),legend.position='none')
 "

R"p_ngg_sb
pdf('/Users/dariabystrova/Documents/GitHub/approximatingBNPpriors/figures/figure_6_5.pdf', width= 4, height = 4)
plot(p_ngg_sb)
dev.off()"

max_stars_repo_name: konkam/approximatingBNPpriors

max_issues_repo_name: konkam/approximatingBNPpriors

max_forks_repo_name: konkam/approximatingBNPpriors


text:
struct SquaresVector <: AbstractArray{Int, 1}
    count::Int
end

Base.size(S::SquaresVector) = (S.count, )

Base.IndexStyle(::Type{<:SquaresVector}) = IndexLinear()

Base.getindex(S::SquaresVector, i::Int) = i*i

struct SparseArray{T,N} <: AbstractArray{T,N}
        data::Dict{NTuple{N,Int}, T}
        dims::NTuple{N,Int}
    end

SparseArray(::Type{T}, dims::Int...) where {T} = SparseArray(T, dims)

SparseArray(::Type{T}, dims::NTuple{N,Int}) where {T,N} = SparseArray{T,N}(Dict{NTuple{N,Int}, T}(), dims)

Base.size(A::SparseArray) = A.dims

Base.similar(A::SparseArray, ::Type{T}, dims::Dims) where {T} = SparseArray(T, dims)

Base.getindex(A::SparseArray{T,N}, I::Vararg{Int,N}) where {T,N} = get(A.data, I, zero(T))

Base.setindex!(A::SparseArray{T,N}, v, I::Vararg{Int,N}) where {T,N} = (A.data[I] = v)

A = SparseArray(Float64, 3, 3)
println(A)

fill!(A, 2)
println(A)

A[:] = 1:length(A)
println(A)

println(A[1:2, :])

B = copy(A)
println(B)

struct ArrayAndChar{T,N} <: AbstractArray{T,N}
    data::Array{T,N}
    char::Char
end
Base.size(A::ArrayAndChar) = size(A.data)
Base.getindex(A::ArrayAndChar{T,N}, inds::Vararg{Int,N}) where {T,N} = A.data[inds...]
Base.setindex!(A::ArrayAndChar{T,N}, val, inds::Vararg{Int,N}) where {T,N} = A.data[inds...] = val
Base.showarg(io::IO, A::ArrayAndChar, toplevel) = print(io, typeof(A), " with char '", A.char, "'")

function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{ArrayAndChar}}, ::Type{ElType}) where ElType
    A = find_aac(bc)
    ArrayAndChar(similar(Array{ElType}, axes(bc)), A.char)
end

find_aac(bc::Base.Broadcast.Broadcasted) = find_aac(bc.args)
find_aac(args::Tuple) = find_aac(find_aac(args[1]), Base.tail(args))
find_aac(x) = x
find_aac(a::ArrayAndChar, rest) = a
find_aac(::Any, rest) = find_aac(rest)

a = ArrayAndChar([1 2; 3 4], 'x')
println(a)
println(a .+ 1)
println(a .+ [5,10])

max_stars_repo_name: JulialangOrgCN/howtojulia

max_issues_repo_name: JulialangOrgCN/howtojulia

max_forks_repo_name: JulialangOrgCN/howtojulia


text:
# StdVector

eval(cxxnewfile("StdVector.cxx", """
    #include <string>
    #include <memory>
    #include <vector>

    static_assert(sizeof(bool) == 1, "");
    """))

function generate(::Type{StdVector{T}}) where {T}
    CT = T == Bool ? "bool" :
         T == RefStdString ? "std::string" :
         T == SharedStdString ? "std::shared_ptr<std::string>" :
         T == RefStdVector{RefStdString} ? "std::vector<std::string>" :
         T == RefStdVector{SharedStdString} ? "std::vector<std::shared_ptr<std::string>>" :
         T == SharedStdVector{RefStdString} ? "std::shared_ptr<std::vector<std::string>>" :
         T == SharedStdVector{SharedStdString} ? "std::shared_ptr<std::vector<std::shared_ptr<std::string>>>" : cxxtype[T]
    NT = cxxname(CT)

    # RefStdVector

    eval(cxxfunction(FnName(:RefStdVector_new, "std_vector_$(NT)_new", libSTL),
                     FnResult(Ptr{StdVector{T}}, "std::vector<$CT> *", RefStdVector{T}, expr -> :(RefStdVector{$T}($expr))),
                     [FnArg(:type, Nothing, "type", "void", Type{T}, identity; skip=true)], "return new std::vector<$CT>;"))
    eval(cxxfunction(FnName(:RefStdVector_new, "std_vector_$(NT)_new_std_size_t", libSTL),
                     FnResult(Ptr{StdVector{T}}, "std::vector<$CT> *", RefStdVector{T}, expr -> :(RefStdVector{$T}($expr))),
                     [FnArg(:type, Nothing, "type", "void", Type{T}, identity; skip=true),
                      FnArg(:size, Csize_t, "size", "std::size_t", Integer, identity)], "return new std::vector<$CT>(size);"))

    eval(cxxfunction(FnName(:RefStdVector_delete, "std_vector_$(NT)_delete", libSTL), FnResult(Nothing, "void"),
                     [FnArg(:vec, Ptr{StdVector{T}}, "vec", "std::vector<$CT> * restrict", RefStdVector{T}, identity)],
                     "delete vec;"))

    eval(cxxfunction(FnName(:(Base.copy), "std_vector_$(NT)_copy", libSTL),
                     FnResult(Ptr{StdVector{T}}, "std::vector<$CT> *", RefStdVector{T}, expr -> :(RefStdVector{$T}($expr))),
                     [FnArg(:vec, Ptr{StdVector{T}}, "vec", "const std::vector<$CT> * restrict", RefStdVector{T}, identity)],
                     "return new std::vector<$CT>(*vec);"))

    # GCStdVector

    eval(cxxfunction(FnName(:GCStdVector_sizeof, "std_vector_$(NT)_sizeof", libSTL), FnResult(Csize_t, "std::size_t"),
                     [FnArg(:valuetype, Nothing, "valuetype", "void", Type{T}, identity; skip=true)],
                     "return sizeof(std::vector<$CT>);"))
    push!(selftests, () -> @assert GCStdVector_sizeof(T) <= GCStdVector_size)

    eval(cxxfunction(FnName(:GCStdVector_construct, "std_vector_$(NT)_construct", libSTL), FnResult(Nothing, "void"),
                     [FnArg(:ptr, Ptr{StdVector{T}}, "ptr", "void *", GCStdVector{T}, identity)], "new(ptr) std::vector<$CT>;"))
    eval(cxxfunction(FnName(:GCStdVector_construct, "std_vector_$(NT)_construct_std_size_t", libSTL), FnResult(Nothing, "void"),
                     [FnArg(:ptr, Ptr{StdVector{T}}, "ptr", "void *", GCStdVector{T}, identity),
                      FnArg(:size, Csize_t, "size", "std::size_t", Integer, identity)], "new(ptr) std::vector<$CT>(size);"))

    eval(cxxfunction(FnName(:GCStdVector_destruct, "std_vector_$(NT)_destruct", libSTL), FnResult(Nothing, "void"),
                     [FnArg(:ptr, Ptr{StdVector{T}}, "ptr", "std::vector<$CT> * restrict", GCStdVector{T}, identity)],
                     "ptr->~vector();"))

    eval(cxxfunction(FnName(:GCStdVector_copy_construct, "std_vector_$(NT)_copy_construct", libSTL), FnResult(Nothing, "void"),
                     [FnArg(:ptr, Ptr{StdVector{T}}, "ptr", "std::vector<$CT> * restrict", GCStdVector{T}, identity),
                      FnArg(:vec, Ptr{StdVector{T}}, "vec", "const std::vector<$CT> * restrict", GCStdVector{T}, identity)],
                     "new(ptr) std::vector<$CT>(*vec);"))

    # SharedStdVector

    eval(cxxfunction(FnName(:SharedStdVector_sizeof, "std_shared_ptr_std_vector_$(NT)_sizeof", libSTL),
                     FnResult(Csize_t, "std::size_t"),
                     [FnArg(:valuetype, Nothing, "valuetype", "void", Type{T}, identity; skip=true)],
                     "return sizeof(std::shared_ptr<std::vector< $CT>>);"))
    push!(selftests, () -> @assert SharedStdVector_sizeof(T) <= SharedStdVector_size)

    eval(cxxfunction(FnName(:SharedStdVector_construct, "std_shared_ptr_std_vector_$(NT)_placement_new", libSTL),
                     FnResult(Nothing, "void"),
                     [FnArg(:ptr, Ptr{Cvoid}, "ptr", "void *", SharedStdVector{T}, expr -> :(pointer_from_objref($expr)))],
                     """
                     auto res = new(ptr) std::shared_ptr<std::vector<$CT>>;
                     *res = std::make_shared<std::vector<$CT>>();
                     """))
    eval(cxxfunction(FnName(:SharedStdVector_construct, "std_shared_ptr_std_vector_$(NT)_placement_new_std_size_t", libSTL),
                     FnResult(Nothing, "void"),
                     [FnArg(:ptr, Ptr{Cvoid}, "ptr", "void *", SharedStdVector{T}, expr -> :(pointer_from_objref($expr))),
                      FnArg(:size, Csize_t, "size", "std::size_t", Integer, identity)],
                     """
                     auto res = new(ptr) std::shared_ptr<std::vector<$CT>>;
                     *res = std::make_shared<std::vector<$CT>>(size);
                     """))

    eval(cxxfunction(FnName(:SharedStdVector_destruct, "std_shared_ptr_std_vector_$(NT)_placement_delete", libSTL),
                     FnResult(Nothing, "void"),
                     [FnArg(:ptr, Ptr{Cvoid}, "ptr", "std::shared_ptr<std::vector<$CT>> * restrict", SharedStdVector{T},
                            expr -> :(pointer_from_objref($expr)))], "ptr->~shared_ptr();"))

    eval(cxxfunction(FnName(:SharedStdVector_copy_construct, "std_shared_ptr_std_vector_$(NT)_placement_copy", libSTL),
                     FnResult(Nothing, "void"),
                     [FnArg(:ptr, Ptr{Cvoid}, "ptr", "void *", SharedStdVector{T}, expr -> :(pointer_from_objref($expr))),
                      FnArg(:vec, Ptr{StdVector{T}}, "vec", "const std::vector<$CT> * restrict", SharedStdVector{T}, identity)],
                     """
                     auto res = new(ptr) std::shared_ptr<std::vector<$CT>>;
                     *res = std::make_shared<std::vector<$CT>>(*vec);
                     """))

    eval(cxxfunction(FnName(:SharedStdVector_get, "std_shared_ptr_std_vector_$(NT)_get", libSTL),
                     FnResult(Ptr{StdVector}, "std::vector<$CT> *"),
                     [FnArg(:ptr, Ptr{Cvoid}, "ptr", "std::shared_ptr<std::vector<$CT>> * restrict", SharedStdVector{T},
                            expr -> :(pointer_from_objref($expr)))], "return ptr->get();"))

    # StdVector

    eval(cxxfunction(FnName(:(Base.resize!), "std_vector_$(NT)_resize_", libSTL), FnResult(Nothing, "void"),
                     [FnArg(:vec, Ptr{StdVector{T}}, "vec", "std::vector<$CT> * restrict", StdVector{T}, identity),
                      FnArg(:size, Csize_t, "size", "std::size_t", Integer, identity)], "vec->resize(size);"))

    eval(cxxfunction(FnName(:(Base.length), "std_vector_$(NT)_length", libSTL),
                     FnResult(Csize_t, "std::size_t", Int, expr -> :(convert(Int, $expr))),
                     [FnArg(:vec, Ptr{StdVector{T}}, "vec", "const std::vector<$CT> * restrict", StdVector{T}, identity)],
                     "return vec->size();"))

    if T == Bool
        eval(cxxfunction(FnName(:(Base.getindex), "std_vector_$(NT)_getindex", libSTL), FnResult(T, CT),
                         [FnArg(:vec, Ptr{StdVector{T}}, "vec", "std::vector<$CT> * restrict", StdVector{T}, identity),
                          FnArg(:idx, Csize_t, "idx", "std::size_t", Integer, identity)], "return (*vec)[idx];"))
    else
        eval(cxxfunction(FnName(:(Base.getindex), "std_vector_$(NT)_getindex", libSTL),
                         FnResult(Ptr{T}, "$CT *", T, expr -> :(convert_result($T, $expr))),
                         [FnArg(:vec, Ptr{StdVector{T}}, "vec", "std::vector<$CT> * restrict", StdVector{T}, identity),
                          FnArg(:idx, Csize_t, "idx", "std::size_t", Integer, identity)], "return &(*vec)[idx];"))
    end

    eval(cxxfunction(FnName(:(Base.setindex!), "std_vector_$(NT)_setindex_", libSTL), FnResult(Nothing, "void"),
                     [FnArg(:vec, Ptr{StdVector{T}}, "vec", "std::vector<$CT> * restrict", StdVector{T}, identity),
                      FnArg(:elt, Ptr{T}, "elt", "$CT const *", Any, expr -> :(convert_arg(Ptr{$T}, convert($T, $expr)))),
                      FnArg(:idx, Csize_t, "idx", "std::size_t", Integer, identity)], "(*vec)[idx] = *elt;"))

    eval(cxxfunction(FnName(:(Base.:(==)), "std_vector_$(NT)_equals", libSTL), FnResult(Bool, "bool"),
                     [FnArg(:vec1, Ptr{StdVector{T}}, "vec1", "const std::vector<$CT> * restrict", StdVector{T}, identity),
                      FnArg(:vec2, Ptr{StdVector{T}}, "vec2", "const std::vector<$CT> * restrict", StdVector{T}, identity)],
                     "return *vec1 == *vec2;"))
    # eval(cxxfunction(FnName(:(Base.:(<)), "std_vector_$(NT)_less", libSTL), FnResult(Bool, "bool"),
    #                  [FnArg(:vec1, Ptr{StdVector{T}}, "vec1", "const std::vector<$CT> * restrict", StdVector{T}, identity),
    #                   FnArg(:vec2, Ptr{StdVector{T}}, "vec2", "const std::vector<$CT> * restrict", StdVector{T}, identity)],
    #                  "return *vec1 < *vec2;"))

    eval(cxxfunction(FnName(:(Base.push!), "std_vector_$(NT)_push_", libSTL), FnResult(Nothing, "void"),
                     [FnArg(:vec, Ptr{StdVector{T}}, "vec", "std::vector<$CT> * restrict", StdVector{T}, identity),
                      FnArg(:elt, Ptr{T}, "elt", "$CT const *", Any, expr -> :(convert_arg(Ptr{$T}, convert($T, $expr))))],
                     "vec->push_back(*elt);"))
    eval(cxxfunction(FnName(:pop_back!, "std_vector_$(NT)_pop_back_", libSTL), FnResult(Nothing, "void"),
                     [FnArg(:vec, Ptr{StdVector{T}}, "vec", "std::vector<$CT> * restrict", StdVector{T}, identity)],
                     "vec->pop_back();"))

    return nothing
end

const StdVector_types = value_types ∪ Set([RefStdString, SharedStdString, RefStdVector{RefStdString}, RefStdVector{SharedStdString},
                                           SharedStdVector{RefStdString}, SharedStdVector{SharedStdString}])
for T in sort!(collect(StdVector_types); by=string)
    generate(StdVector{T})
end

Base.copy(vec::GCStdVector{T}) where {T} = GCStdVector{T}(Base.copy, vec)
Base.copy(vec::SharedStdVector{T}) where {T} = SharedStdVector{T}(Base.copy, vec)
Base.pop!(vec::StdVector) = (elt = vec[end]; pop_back!(vec); elt)

free(vec::RefStdVector) = RefStdVector_delete(vec)

################################################################################

Base.convert(::Type{Vector{T}}, vec::StdVector{T}) where {T} = T[elt for elt in vec]
Base.convert(::Type{Vector}, vec::StdVector{T}) where {T} = convert(Vector{T}, vec)

Base.empty!(vec::StdVector) = resize!(vec, 0)
Base.isempty(vec::StdVector) = length(vec) == 0
Base.size(vec::StdVector) = (length(vec),)
Base.firstindex(vec::StdVector) = 0
Base.lastindex(vec::StdVector) = length(vec) - 1

Base.eltype(::StdVector{T}) where {T} = T

function Base.iterate(vec::StdVector, pos::Int=0)
    pos ≥ length(vec) && return nothing
    return vec[pos], pos + 1
end


text:
using ParallelAccelerator

immutable newType
    w::Vector{Float64}
    b::Float64
end
@acc function evaluate(f::newType, X::AbstractMatrix)
    if f.b != 0
	X = X + f.b
    end
    return X
end

function main(m::Int, n::Int)
    X = Array(Float64, m, n)
    W = Array(Float64, m)
    fill!(X, 5)
    fill!(W, 3)
    b = 1.0
    f = newType(W,b)
    evaluate(f, X)
    println("done")

end

main(100, 200)


text:
using BenchmarkTools, MixedModels, StatsModels
using MixedModels: dataset

const SUITE = BenchmarkGroup()

const global contrasts = Dict(
    :mrk17_exp1 => merge(Dict(n => HelmertCoding() for n in (:F, :P, :Q, :lQ, :lT)), Dict(n => Grouping() for n in (:item, :subj))),
    :kb07 => merge(Dict(n => HelmertCoding() for n in (:spkr, :prec, :load)), Dict(n => Grouping() for n in (:item, :subj))),
    :insteval => Dict(n => Grouping() for n in (:s, :d)),
)
const global fms = Dict(
    :dyestuff => [
        @formula(yield ~ 1 + (1|batch)),
        ],
    :dyestuff2 => [
        @formula(yield ~ 1 + (1|batch)),
        ],
    :d3 => [
        @formula(y ~ 1 + u + (1+u|g) + (1+u|h) + (1+u|i)),
        ],
    :insteval => [
        @formula(y ~ 1 + service + (1|s) + (1|d) + (1|dept)),
        @formula(y ~ 1 + service*dept + (1|s) + (1|d)),
        ],
    :kb07 => [
        @formula(rt_trunc ~ 1+spkr+prec+load+(1|subj)+(1|item)),
        @formula(rt_trunc ~ 1+spkr*prec*load+(1|subj)+(1+prec|item)),
        @formula(rt_trunc ~ 1+spkr*prec*load+(1+spkr+prec+load|subj)+(1+spkr+prec+load|item)),
        ],
    :machines => [
        @formula(score ~ 1 + (1|Worker) + (1|Machine)),
        ],
    :ml1m => [
        @formula(y ~ 1 + (1|g) + (1|h)),
        ],
    :mrk17_exp1 => [
        @formula(1000/rt ~ 1+F*P*Q*lQ*lT + (1|item) + (1|subj)),
        @formula(1000/rt ~ 1+F*P*Q*lQ*lT + (1+P+Q+lQ+lT|item) + (1+F+P+Q+lQ+lT|subj)),
        ],
    :pastes => [
        @formula(strength ~ 1 + (1|batch&cask)),
        @formula(strength ~ 1 + (1|batch/cask)),
        ],
    :penicillin => [
        @formula(diameter ~ 1 + (1|plate) + (1|sample)),
        ],
    :sleepstudy => [
        @formula(reaction ~ 1 + days + (1|subj)),
        @formula(reaction ~ 1 + days + zerocorr(1+days|subj)),
        @formula(reaction ~ 1 + days + (1|subj) + (0+days|subj)),
        @formula(reaction ~ 1 + days + (1+days|subj)),
        ],
)

function fitbobyqa(dsname::Symbol, index::Integer)
    fit(
        MixedModel,
        fms[dsname][index],
        dataset(dsname),
        contrasts=get!(contrasts, dsname, Dict{Symbol,StatsModels.AbstractContrasts}()),
        )
end

SUITE["simplescalar"] = BenchmarkGroup(["single", "simple", "scalar"])
for (ds, i) in [
    (:dyestuff, 1,),
    (:dyestuff2, 1,),
    (:pastes, 1),
    (:sleepstudy, 1,),
]
    SUITE["simplescalar"][string(ds, ':', i)] = @benchmarkable fitbobyqa($(QuoteNode(ds)), $(QuoteNode(i)))
end

SUITE["singlevector"] = BenchmarkGroup(["single", "vector"])
for (ds, i) in [
    (:sleepstudy, 2,),
    (:sleepstudy, 3,),
    (:sleepstudy, 4,),
]
    SUITE["singlevector"][string(ds, ':', i)] = @benchmarkable fitbobyqa($(QuoteNode(ds)), $(QuoteNode(i)))
end

SUITE["nested"] = BenchmarkGroup(["multiple", "nested", "scalar"])
for (ds, i) in [
    (:pastes, 2,),
]
    SUITE["nested"][string(ds, ':', i)] = @benchmarkable fitbobyqa($(QuoteNode(ds)), $(QuoteNode(i)))
end

SUITE["crossed"] = BenchmarkGroup(["multiple", "crossed", "scalar"])
for (ds, i) in [
    (:insteval, 1),
    (:insteval, 2),
    (:kb07, 1),
    (:machines, 1),
    (:ml1m, 1),
    (:mrk17_exp1, 1),
    (:penicillin, 1),
]
    SUITE["crossed"][string(ds, ':', i)] = @benchmarkable fitbobyqa($(QuoteNode(ds)), $(QuoteNode(i)))
end

SUITE["crossedvector"] = BenchmarkGroup(["multiple", "crossed", "vector"])
for (ds, i) in [
    (:d3, 1),
    (:kb07, 2),
    (:kb07, 3),
    (:mrk17_exp1, 2),
]
    SUITE["crossedvector"][string(ds, ':', i)] = @benchmarkable fitbobyqa($(QuoteNode(ds)), $(QuoteNode(i)))
end


text:
# statistics
cd("G:\\My Drive\\18. Github Repo\\DataScience");
# load the libraries
using Statistics, StatsBase, Plots, StatsPlots, Distributions, LinearAlgebra, MLBase, HypothesisTests, PyCall, DataFrames, RDatasets

# Load data and clean
D = dataset("datasets", "faithful");
@show names(D);
describe(D)

# plot
eruptions = D[!, :Eruptions];
scatter(eruptions, lab="Eruptions");
waittime = D[!, :Waiting];
scatter!(waittime, lab="waittime")

# boxplot
boxplot(["eruption length"], eruptions, legend=false, size=(200,400), whisker_width=1, ylabel="time in minutes")

# histogram
histogram(eruptions, lab="eruptions", bins=:sqrt)

# fit using kernel density
using KernelDensity
p = kde(eruptions);

# fits
histogram(eruptions,label="eruptions", bins=:sqrt)
plot!(p.x,p.density .* length(eruptions)*0.4, linewidth=3,color=2,label="kde fit")
# nb of elements*bin width

histogram(eruptions,bins=:sqrt,label="eruptions")
plot!(p.x,p.density .* length(eruptions) .*0.2, linewidth=3,color=2,label="kde fit")

#
myrandomvector = randn(100_000)
histogram(myrandomvector)
p=kde(myrandomvector)
plot!(p.x,p.density .* length(myrandomvector) .*0.1, linewidth=3,color=2,label="kde fit") # nb of elements*bin width

# Probability distributions
d = Normal()
myrandomvector = rand(d,100000)
histogram(myrandomvector)
p=kde(myrandomvector)
plot!(p.x,p.density .* length(myrandomvector) .*0.1, linewidth=3,color=2,label="kde fit") # nb of elements*bin width

b = Binomial(40)
myrandomvector = rand(b,1000000)
histogram(myrandomvector)
p=kde(myrandomvector)
plot!(p.x,p.density .* length(myrandomvector) .*0.5,color=2,label="kde fit") # nb of elements*bin width

# fit distributions
x = rand(1000)
d = fit(Normal, x)
myrandomvector = rand(d,1000)
histogram(myrandomvector,nbins=20,fillalpha=0.3,label="fit")
histogram!(x,nbins=20,linecolor = :red,fillalpha=0.3,label="myvector")

x = eruptions
d = fit(Normal, x)
myrandomvector = rand(d,1000)
histogram(myrandomvector,nbins=20,fillalpha=0.3)
histogram!(x,nbins=20,linecolor = :red,fillalpha=0.3)

# T-Test
ttest = OneSampleTTest(myrandomvector)

# MLBase
gt = [1,1,1,1,1,1,1,2];
pred = [1,1,2,2,1,1,1,1];
C = confusmat(2, gt, pred);
C = C ./ sum(C, dims=2)
correctrate(gt, pred)

#
scatter(eruptions, waittime)

max_stars_repo_path: 03. Statistics.jl

max_issues_repo_path: 03. Statistics.jl

max_forks_repo_path: 03. Statistics.jl


text:
#=

@author : Spencer Lyon <spencer.lyon@nyu.edu>

=#

using QuantEcon
using Interpolations
using Distributions

"""
Unemployment/search problem where offer distribution is unknown

##### Fields

- `β::Real` : Discount factor on (0, 1)
- `c::Real` : Unemployment compensation
- `F::Distribution` : Offer distribution `F`
- `G::Distribution` : Offer distribution `G`
- `f::Function` : The pdf of `F`
- `g::Function` : The pdf of `G`
- `n_w::Int` : Number of points on the grid for w
- `w_max::Real` : Maximum wage offer
- `w_grid::AbstractVector` : Grid of wage offers w
- `n_π::Int` : Number of points on grid for π
- `π_min::Real` : Minimum of π grid
- `π_max::Real` : Maximum of π grid
- `π_grid::AbstractVector` : Grid of probabilities π
- `quad_nodes::Vector` : Notes for quadrature ofer offers
- `quad_weights::Vector` : Weights for quadrature ofer offers

"""
struct SearchProblem{TR<:Real, TI<:Integer, TF<:AbstractFloat,
                    TAVw<:AbstractVector{TF}, TAVpi<:AbstractVector{TF}}
    β::TR
    c::TR
    F::Distribution
    G::Distribution
    f::Function
    g::Function
    n_w::TI
    w_max::TR
    w_grid::TAVw
    n_π::TI
    π_min::TR
    π_max::TR
    π_grid::TAVpi
    quad_nodes::Vector{TF}
    quad_weights::Vector{TF}
end

"""
Constructor for `SearchProblem` with default values

##### Arguments

- `β::Real(0.95)` : Discount factor in (0, 1)
- `c::Real(0.6)` : Unemployment compensation
- `F_a::Real(1), F_b::Real(1)` : Parameters of `Beta` distribution for `F`
- `G_a::Real(3), G_b::Real(1.2)` : Parameters of `Beta` distribution for `G`
- `w_max::Real(2)` : Maximum of wage offer grid
- `w_grid_size::Integer(40)` : Number of points in wage offer grid
- `π_grid_size::Integer(40)` : Number of points in probability grid

##### Notes

There is also a version of this function that accepts keyword arguments for
each parameter

"""

# use key word argment
function SearchProblem(;β::Real=0.95, c::Real=0.6, F_a::Real=1, F_b::Real=1,
                       G_a::Real=3, G_b::Real=1.2, w_max::Real=2.0,
                       w_grid_size::Integer=40, π_grid_size::Integer=40)

    F = Beta(F_a, F_b)
    G = Beta(G_a, G_b)

    # NOTE: the x./w_max)./w_max in these functions makes our dist match
    #       the scipy one with scale=w_max given
    f(x) = pdf.(F, x/w_max)/w_max
    g(x) = pdf.(G, x/w_max)/w_max

    π_min = 1e-3  # avoids instability
    π_max = 1 - π_min

    w_grid = linspace(0, w_max, w_grid_size)
    π_grid = linspace(π_min, π_max, π_grid_size)

    nodes, weights = qnwlege(21, 0.0, w_max)

    SearchProblem(β, c, F, G, f, g,
                  w_grid_size, w_max, w_grid,
                  π_grid_size, π_min, π_max, π_grid, nodes, weights)
end

function q(sp::SearchProblem, w, π_val)
    new_π = 1.0 / (1 + ((1 - π_val) * sp.g(w)) / (π_val * sp.f(w)))

    # Return new_π when in [π_min, π_max] and else end points
    return clamp(new_π, sp.π_min, sp.π_max)
end

"""
Apply the Bellman operator for a given model and initial value.

##### Arguments

- `sp::SearchProblem` : Instance of `SearchProblem`
- `v::Matrix`: Current guess for the value function
- `out::Matrix` : Storage for output.
- `;ret_policy::Bool(false)`: Toggles return of value or policy functions

##### Returns

None, `out` is updated in place. If `ret_policy == true` out is filled with the
policy function, otherwise the value function is stored in `out`.

"""
function bellman_operator!(sp::SearchProblem, v::Matrix, out::Matrix;
                          ret_policy::Bool=false)
    # Simplify names
    f, g, β, c = sp.f, sp.g, sp.β, sp.c
    nodes, weights = sp.quad_nodes, sp.quad_weights

    vf = extrapolate(interpolate((sp.w_grid, sp.π_grid), v,
                     Gridded(Linear())), Flat())

    # set up quadrature nodes/weights
    # q_nodes, q_weights = qnwlege(21, 0.0, sp.w_max)

    for (w_i, w) in enumerate(sp.w_grid)
        # calculate v1
        v1 = w / (1 - β)

        for (π_j, _π) in enumerate(sp.π_grid)
            # calculate v2
            integrand(m) = [vf[m[i], q.(sp, m[i], _π)] *
                          (_π * f(m[i]) + (1 - _π) * g(m[i])) for i in 1:length(m)]
            integral = do_quad(integrand, nodes, weights)
            # integral = do_quad(integrand, q_nodes, q_weights)
            v2 = c + β * integral

            # return policy if asked for, otherwise return max of values
            out[w_i, π_j] = ret_policy ? v1 > v2 : max(v1, v2)
        end
    end
    return out
end

function bellman_operator(sp::SearchProblem, v::Matrix;
                          ret_policy::Bool=false)
    out_type = ret_policy ? Bool : Float64
    out = Array{out_type}(sp.n_w, sp.n_π)
    bellman_operator!(sp, v, out, ret_policy=ret_policy)
end

"""
Extract the greedy policy (policy function) of the model.

##### Arguments

- `sp::SearchProblem` : Instance of `SearchProblem`
- `v::Matrix`: Current guess for the value function
- `out::Matrix` : Storage for output

##### Returns

None, `out` is updated in place to hold the policy function

"""
get_greedy!(sp::SearchProblem, v::Matrix, out::Matrix) =
    bellman_operator!(sp, v, out, ret_policy=true)

get_greedy(sp::SearchProblem, v::Matrix) =
    bellman_operator(sp, v, ret_policy=true)

"""
Updates the reservation wage function guess ϕ via the operator Q.

##### Arguments

- `sp::SearchProblem` : Instance of `SearchProblem`
- `ϕ::Vector`: Current guess for ϕ
- `out::Vector` : Storage for output

##### Returns

None, `out` is updated in place to hold the updated levels of ϕ
"""
function res_wage_operator!(sp::SearchProblem, ϕ::Vector, out::Vector)
    # Simplify name
    f, g, β, c = sp.f, sp.g, sp.β, sp.c

    # Construct interpolator over π_grid, given ϕ
    ϕ_f = LinInterp(sp.π_grid, ϕ)

    # set up quadrature nodes/weights
    q_nodes, q_weights = qnwlege(7, 0.0, sp.w_max)

    for (i, _π) in enumerate(sp.π_grid)
        integrand(x) = max.(x, ϕ_f.(q.(sp, x, _π))) .* (_π * f(x) + (1 - _π) * g(x))
        integral = do_quad(integrand, q_nodes, q_weights)
        out[i] = (1 - β) * c + β * integral
    end
end

"""
Updates the reservation wage function guess ϕ via the operator Q.

See the documentation for the mutating method of this function for more details
on arguments
"""
function res_wage_operator(sp::SearchProblem, ϕ::Vector)
    out = similar(ϕ)
    res_wage_operator!(sp, ϕ, out)
    return out
end


text:
module PoissonRandom

using Random

export pois_rand

count_rand(λ) = count_rand(Random.GLOBAL_RNG, λ)
function count_rand(rng::AbstractRNG, λ)
    n = 0
    c = randexp(rng)
    while c < λ
        n += 1
        c += randexp(rng)
    end
    return n
end

# Algorithm from:
#
#   J.H. Ahrens, U. Dieter (1982)
#   "Computer Generation of Poisson Deviates from Modified Normal Distributions"
#   ACM Transactions on Mathematical Software, 8(2):163-179
#
#   For μ sufficiently large, (i.e. >= 10.0)
#
ad_rand(λ) = ad_rand(Random.GLOBAL_RNG, λ)
function ad_rand(rng::AbstractRNG, λ)
  s = sqrt(λ)
  d = 6.0*λ^2
  L = floor(Int,λ-1.1484)
  # Step N
  G = λ + s*randn(rng)

  if G >= 0.0
      K = floor(Int,G)
      # Step I
      if K >= L
          return K
      end

      # Step S
      U = rand(rng)
      if d*U >= (λ-K)^3
          return K
      end

      # Step P
      px,py,fx,fy = procf(λ,K,s)

      # Step Q
      if fy*(1-U) <= py*exp(px-fx)
          return K
      end
  end

  while true
      # Step E
      E = randexp(rng)
      U = 2.0*rand(rng)-1.0
      T = 1.8+copysign(E,U)
      if T <= -0.6744
          continue
      end

      K = floor(Int,λ + s*T)
      px,py,fx,fy = procf(λ,K,s)
      c = 0.1069/λ

      # Step H
      @fastmath if c*abs(U) <= py*exp(px+E)-fy*exp(fx+E)
          return K
      end
  end
end

# log(1+x)-x
# accurate ~2ulps for -0.227 < x < 0.315
function log1pmx_kernel(x::Float64)
    r = x/(x+2.0)
    t = r*r
    w = @evalpoly(t,
                6.66666666666666667e-1, # 2/3
                4.00000000000000000e-1, # 2/5
                2.85714285714285714e-1, # 2/7
                2.22222222222222222e-1, # 2/9
                1.81818181818181818e-1, # 2/11
                1.53846153846153846e-1, # 2/13
                1.33333333333333333e-1, # 2/15
                1.17647058823529412e-1) # 2/17
    hxsq = 0.5*x*x
    r*(hxsq+w*t)-hxsq
end

# use naive calculation or range reduction outside kernel range.
# accurate ~2ulps for all x
function log1pmx(x::Float64)
    if !(-0.7 < x < 0.9)
        return log1p(x) - x
    elseif x > 0.315
        u = (x-0.5)/1.5
        return log1pmx_kernel(u) - 9.45348918918356180e-2 - 0.5*u
    elseif x > -0.227
        return log1pmx_kernel(x)
    elseif x > -0.4
        u = (x+0.25)/0.75
        return log1pmx_kernel(u) - 3.76820724517809274e-2 + 0.25*u
    elseif x > -0.6
        u = (x+0.5)*2.0
        return log1pmx_kernel(u) - 1.93147180559945309e-1 + 0.5*u
    else
        u = (x+0.625)/0.375
        return log1pmx_kernel(u) - 3.55829253011726237e-1 + 0.625*u
    end
end

# Procedure F
function procf(λ, K::Int, s::Float64)
    # can be pre-computed, but does not seem to affect performance
    ω = 0.3989422804014327/s
    b1 = 0.041666666666666664/λ
    b2 = 0.3*b1*b1
    c3 = 0.14285714285714285*b1*b2
    c2 = b2 - 15.0*c3
    c1 = b1 - 6.0*b2 + 45.0*c3
    c0 = 1.0 - b1 + 3.0*b2 - 15.0*c3

    if K < 10
        px = -float(λ)
        py = λ^K/factorial(K)
    else
        δ = 0.08333333333333333/K
        δ -= 4.8*δ^3
        V = (λ-K)/K
        px = K*log1pmx(V) - δ # avoids need for table
        py = 0.3989422804014327/sqrt(K)
    end
    X = (K-λ+0.5)/s
    X2 = X^2
    fx = -0.5*X2 # missing negation in pseudo-algorithm, but appears in fortran code.
    fy = ω*(((c3*X2+c2)*X2+c1)*X2+c0)
    return px,py,fx,fy
end

pois_rand(λ) = pois_rand(Random.GLOBAL_RNG, λ)
pois_rand(rng::AbstractRNG, λ) = λ < 6 ? count_rand(rng, λ) : ad_rand(rng, λ)

end # module


text:
# Authors: Ali Siahkoohi, alisk@gatech.edu
# Date: March 2022

using DrWatson
@quickactivate :ConditionalNFs4Imaging

using InvertibleNetworks
using PyPlot
using Seaborn

font_prop, sfmt = sef_plot_configs()
args = read_config("amortized_MAP_visualize.json")
args = parse_input_args(args)

max_epoch = args["max_epoch"]
epoch = args["epoch"]
lr = args["lr"]
idx = args["idx"]
sigma = args["sigma"]
nsrc = args["nsrc"]
nrec = args["nrec"]
sim_name = args["sim_name"]

# Loading the experiment—only network weights and training loss
loaded_keys = load_experiment(args, ["fval", "mval", "X_hat", "X_OOD", "Y_OOD", "X0"])
fval = loaded_keys["fval"]
mval = loaded_keys["mval"]
X_hat = loaded_keys["X_hat"]
X_OOD = loaded_keys["X_OOD"]
Y_OOD = loaded_keys["Y_OOD"]
X0 = loaded_keys["X0"]

save_dict = @strdict sim_name idx sigma lr max_epoch nsrc nrec epoch
save_path = plotsdir(sim_name, savename(save_dict; digits = 6))

spacing = [20.0, 12.5]
extent = [0.0, size(X_OOD, 1) * spacing[1], size(X_OOD, 2) * spacing[2], 0.0] / 1e3

# Training loss
fig = figure("training logs", figsize = (7, 4))
plot(range(0, max_epoch, length = length(fval)), fval, color = "#d48955")
title("Amortized MAP objective")
ylabel("Objective value")
xlabel("Numbr of passes over shots")
wsave(joinpath(save_path, "log.png"), fig)
close(fig)

fig = figure(figsize = (7, 4))
plot(range(0, max_epoch, length = length(mval)), mval, color = "#4fc9af")
title("Average L2-squared error")
xlabel("Numbr of passes over shots")
wsave(joinpath(save_path, "err_log.png"), fig)
close(fig)

# Plot the true model
fig = figure("x", figsize = (7.68, 4.8))
imshow(
    X_OOD[:, :, 1, 1],
    vmin = -1.5e3,
    vmax = 1.5e3,
    aspect = 1,
    cmap = "Greys",
    resample = true,
    interpolation = "lanczos",
    filterrad = 1,
    extent = extent,
)
title(L"High-fidelity OOD image, $\mathbf{x}$")
colorbar(fraction = 0.03, pad = 0.01, format = sfmt)
grid(false)
xlabel("Horizontal distance (km)")
ylabel("Depth (km)")
safesave(joinpath(save_path, string(idx), "true_model.png"), fig)
close(fig)

# Plot the true model
fig = figure("x", figsize = (7.68, 4.8))
imshow(
    X0[:, :, 1, 1],
    vmin = -1.5e3,
    vmax = 1.5e3,
    aspect = 1,
    cmap = "Greys",
    resample = true,
    interpolation = "lanczos",
    filterrad = 1,
    extent = extent,
)
title("Initial guess")
colorbar(fraction = 0.03, pad = 0.01, format = sfmt)
grid(false)
xlabel("Horizontal distance (km)")
ylabel("Depth (km)")
safesave(joinpath(save_path, string(idx), "X0.png"), fig)
close(fig)

# Plot the observed data
fig = figure("y", figsize = (7.68, 4.8))
imshow(
    Y_OOD[:, :, 1, 1],
    vmin = -1.5e6,
    vmax = 1.5e6,
    aspect = 1,
    cmap = "Greys",
    resample = true,
    interpolation = "lanczos",
    filterrad = 1,
    extent = extent,
)
title(L"Low-fidelity OOD reverse-time migrated image, $\mathbf{y}$")
colorbar(fraction = 0.03, pad = 0.01, format = sfmt)
grid(false)
xlabel("Horizontal distance (km)")
ylabel("Depth (km)")
safesave(joinpath(save_path, string(idx), "observed_data.png"), fig)
close(fig)

# Plot MAP
fig = figure("x_cm", figsize = (7.68, 4.8))
imshow(
    X_hat[:, :, 1, 1],
    vmin = -1.5e3,
    vmax = 1.5e3,
    aspect = 1,
    cmap = "Greys",
    resample = true,
    interpolation = "lanczos",
    filterrad = 1,
    extent = extent,
)
title("MAP estimate")
colorbar(fraction = 0.03, pad = 0.01, format = sfmt)
grid(false)
xlabel("Horizontal distance (km)")
ylabel("Depth (km)")
safesave(joinpath(save_path, string(idx), "MAP.png"), fig)
close(fig)


text:
module Broken

greet() = print("Hello World!")

foo() = 1

end # module


text:
import Pkg

modelproto(;kwargs...) = ModelProto(;
    ir_version=7,
    opset_import=[OperatorSetIdProto(version=14)],
    producer_name="ONNX.jl",
    producer_version=string(Pkg.Types.Context().env.project.version), # TODO: Ugh....
    kwargs...)

"""
    graphproto()
Return an [`ONNX.GraphProto`](@ref) with all fields initialized to empty arrays.
"""
graphproto(;kwargs...) = GraphProto(;
    node = NodeProto[],
    initializer = TensorProto[],
    input = ValueInfoProto[],
    output = ValueInfoProto[],
    value_info = ValueInfoProto[],
    kwargs...
)

add!(gp::GraphProto, np::NodeProto) = push!(gp.node, np)

add!(gp::GraphProto, tp::TensorProto) = push!(gp.initializer, tp)

##############################################################################
#                                 Utils                                    #
##############################################################################

# can we make it more robust?
iskwfunc(f) = endswith(string(f), "##kw")

function kwargs2dict(op::Ghost.Call)
    kw = iskwfunc(op.fn) ? op.args[1] : (;)
    return Dict(zip(keys(kw), values(kw)))
end

macro opconfig_kw(backend, fn)
    return quote
        $OpConfig{$backend, <:Union{typeof($fn), typeof(Core.kwfunc($fn))}}
    end
end

function NodeProto(op_type::String, op::Ghost.Call, attrs::Dict=Dict())
    args = iskwfunc(op.fn) ? op.args[3:end] : op.args
    return NodeProto(
        input=[onnx_name(v) for v in args],
        output=[onnx_name(op)],
        name=onnx_name(op),
        attribute=AttributeProto.(keys(attrs), values(attrs)),
        op_type=op_type
    )
end

ValueInfoProto(op::Ghost.AbstractOp) = ValueInfoProto(
    onnx_name(op),
    # utils in write.jl reverse the shape, so we don't do it here
    # try the following for example:
    #     TypeProto_Tensor((4, 3), Float64).shape.dim[1].dim_value
    # which gives 3 instead of 4
    size(op.val),
    eltype(op.val)
)

##############################################################################
#                                 Methods                                    #
##############################################################################

onnx_name(v::Variable) = "x$(v.id)"
onnx_name(op::Ghost.AbstractOp) = "x$(op.id)"

"""
    save_node!(g::GraphProto, op::Ghost.Call)
    save_node!(g::GraphProto, ::OpConfig{:Backend, Fn}, op::Ghost.Call)

Serialize a single operation from a tape to graph.
"""
function save_node!(g::GraphProto, op::Ghost.Call)
    save_node!(g, OpConfig{:ONNX, typeof(op.fn)}(), op)
end

function save_node!(g::GraphProto, ::OpConfig{:ONNX, typeof(getfield)}, op::Ghost.Call)
    # Do nothing: getfield is only used to destructure multi-ouput nodes
    # and doesn't need to be written to ONNX graph.
    # Using getfield() for anything other then destructuring is thus a mistake.
end

function save_node!(g::GraphProto, ::OpConfig{:ONNX, typeof(*)}, op::Ghost.Call)
    nd = NodeProto(
        input=[onnx_name(v) for v in reverse(op.args)],
        output=[onnx_name(op)],
        name=onnx_name(op),
        attribute=AttributeProto[],
        op_type="Gemm"
    )
    push!(g.node, nd)
end

function save_node!(g::GraphProto, ::@opconfig_kw(:ONNX, onnx_gemm), op::Ghost.Call)
    kw_dict = kwargs2dict(op)
    attrs = rename_keys(kw_dict, Dict(
        :tA => :transA,
        :tB => :transB,
        :α => :alpha,
        :β => :beta
    ))
    nd = NodeProto("Gemm", op, attrs)
    push!(g.node, nd)
end

function save_node!(g::GraphProto, ::@opconfig_kw(:ONNX, conv), op::Ghost.Call)
    args = iskwfunc(op.fn) ? op.args[3:end] : op.args
    w = args[2]._op.val
    # ONNXRuntime gives the following error for Float64:
    # NOT_IMPLEMENTED : Could not find an implementation for the node x3:Conv(11)')
    eltype(w) == Float64 && @warn "Not all ONNX runtimes support input & weights as Float64"
    attrs = from_nnlib_conv(kwargs2dict(op), ndims(w) - 2)
    nd = NodeProto("Conv", op, attrs)
    push!(g.node, nd)
end

function save_node!(g::GraphProto, ::@opconfig_kw(:ONNX, maxpool), op::Ghost.Call)
    args = iskwfunc(op.fn) ? op.args[3:end] : op.args
    x = args[1]._op.val
    attrs = from_nnlib_conv(kwargs2dict(op), ndims(x) - 2)
    nd = NodeProto("MaxPool", op, attrs)
    push!(g.node, nd)
end

function save_node!(g::GraphProto, ::OpConfig{:ONNX, typeof(global_average_pool)}, op::Ghost.Call)
    nd = NodeProto("GlobalAveragePool", op)
    push!(g.node, nd)
end

function save_node!(g::GraphProto, ::@opconfig_kw(:ONNX, onnx_flatten), op::Ghost.Call)
    nd = NodeProto("Flatten", op)
    push!(g.node, nd)
end

function save_node!(g::GraphProto, ::OpConfig{:ONNX, typeof(add)}, op::Ghost.Call)
    nd = NodeProto("Add", op)
    push!(g.node, nd)
end

function save_node!(g::GraphProto, ::OpConfig{:ONNX, typeof(mul)}, op::Ghost.Call)
    nd = NodeProto("Mul", op)
    push!(g.node, nd)
end

function save_node!(g::GraphProto, ::OpConfig{:ONNX, typeof(relu)}, op::Ghost.Call)
    nd = NodeProto("Relu", op)
    push!(g.node, nd)
end

function save_node!(g::GraphProto, ::OpConfig{:ONNX, typeof(NNlib.batched_mul)}, op::Ghost.Call)
    nd = NodeProto(
        input=[onnx_name(v) for v in reverse(op.args)],
        output=[onnx_name(op)],
        name=onnx_name(op),
        attribute=AttributeProto[],
        op_type="MatMul"
    )
    push!(g.node, nd)
end

function save_node!(g::GraphProto, ::@opconfig_kw(:ONNX, batch_norm), op::Ghost.Call)
    kw_dict = kwargs2dict(op)
    attrs = from_nnlib_norm(kw_dict)
    args = iskwfunc(op.fn) ? op.args[3:end] : op.args
    output = if Bool(get(attrs, :training_mode, 0))
        vars = unpacked_vars(op)
        @assert(all([v isa V for v in vars]),
            "Not all output vars of batch_norm are unpacked to the tape")
        [onnx_name(v) for v in vars]
    else
        [onnx_name(op)]
    end
    nd = NodeProto(
        input=[onnx_name(v) for v in args],
        output=output,
        name=onnx_name(op),
        attribute=AttributeProto.(keys(attrs), values(attrs)),
        op_type="BatchNormalization"
    )
    push!(g.node, nd)
end

function save_node!(g::GraphProto, ::OpConfig{:ONNX, typeof(size)}, op::Ghost.Call)
    nd = NodeProto("Shape", op)
    push!(g.node, nd)
end

function save_node!(g::GraphProto, ::OpConfig{:ONNX, <:Any}, op::Ghost.Constant)
    @assert(
        op.val isa AbstractArray,
        "ONNX.jl currently doesn't support saving constants of type $(typeof(op.val))"
    )
    attr_name = :value
    attr_value = from_nnlib(op.val)
    nd = NodeProto(
        input=[],
        output=[onnx_name(op)],
        name=onnx_name(op),
        attribute=AttributeProto.([attr_name], [attr_value]),
        op_type="Constant"
    )
    push!(g.node, nd)
end

function save_node!(g::GraphProto, ::@opconfig_kw(:ONNX, onnx_gather), op::Ghost.Call)
    data = iskwfunc(op.fn) ? op.args[3]._op.val : op.args[1]._op.val
    kw_dict = kwargs2dict(op)
    dim = get(kw_dict, :dim, ndims(data))
    axis = ndims(data) - dim
    nd = NodeProto("Gather", op, Dict(:axis => axis))
    push!(g.node, nd)
end

function save_node!(g::GraphProto, ::@opconfig_kw(:ONNX, onnx_unsqueeze), op::Ghost.Call)
    nd = NodeProto("Unsqueeze", op)
    push!(g.node, nd)
end

function save_node!(g::GraphProto, ::OpConfig{:ONNX, typeof(onnx_slice)}, op::Ghost.Call)
    nd = NodeProto("Slice", op)
    push!(g.node, nd)
end

##############################################################################
#                                    API                                     #
##############################################################################

"""
    save(io::IO, tape::Ghost.Tape{ONNXCtx})
    save(filename::String, tape::Ghost.Tape{ONNXCtx})

Save tape as an ONNX model. The way a particular operation is serialized is
controlled by methods of [save_node!](@ref).

See also: [`load!`](@ref)
"""
function save(io::IO, tape::Tape{ONNXCtx})
    g = graphproto()
    g.name = "generated_model"
    for (i, op) in enumerate(tape)
        if op isa Ghost.Input
            # add input to g.input, but not to g.initializer
            push!(g.input, ValueInfoProto(op))
        elseif op isa Ghost.Constant
            # add constant to g.initializer, but not to g.input
            # some models out there also put constants & parameters
            # to g.init, but it seems to be an outdated practise
            push!(g.initializer, TensorProto(op.val, onnx_name(op)))
        elseif op isa Ghost.Call
            save_node!(g, op)
        else
            error("$(typeof(op)) is not yet supported in model export")
        end
    end
    res = tape[tape.result]
    if res.val isa Tuple
        # if the last operation in the graph is multi-output, there must be
        # unpacked elements of that var
        vars = unpacked_vars(res)
        @assert(all(v isa V for v in vars), "Cannot save the tape because the result " *
            "is multi-output, but its elements aren't destructured to the tape")
        for v in vars
            push!(g.output, ValueInfoProto(tape[v]))
        end
    else
        push!(g.output, ValueInfoProto(tape[tape.result]))
    end
    m = modelproto();
    m.graph = g;
    writeproto(io, m)
end

function save(filename::String, tape::Tape{ONNXCtx})
    open(filename, "w") do io
        save(io, tape)
    end
end


text:
@testset "utils" begin
    @test @inferred(TreeArrays.nfill(2, Val(2))) === (2,2)
    @test @inferred(TreeArrays.nfill(2, Val(3))) === (2,2,2)
end

@testset "Power2" begin
    @test TreeArrays.Power2(0) == 2^0
    @test TreeArrays.Power2(1) == 2^1
    @test TreeArrays.Power2(2) == 2^2

    x = TreeArrays.Power2(3)
    x′ = 2^3
    @test @inferred(2 * x) === 2 * x′
    @test @inferred(x * 2) === x′ * 2
    @test @inferred(x * x) === TreeArrays.Power2(6) == x′ * x′

    @test @inferred(x + 1) === x′ + 1
    @test @inferred(x - 1) === x′ - 1
    @test @inferred(1 + x) === 1 + x′
    @test @inferred(1 - x) === 1 - x′

    @test @inferred(div(20, x)) === div(20, x′)
    @test @inferred(rem(20, x)) === rem(20, x′)
    @test @inferred(divrem(20, x)) === divrem(20, x′)

    @test @inferred(zero(x)) === TreeArrays.Power2(0)
end


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

"""
    create_endpoint(outpost_id, security_group_id, subnet_id)
    create_endpoint(outpost_id, security_group_id, subnet_id, params::Dict{String,<:Any})

S3 on Outposts access points simplify managing data access at scale for shared datasets in
Amazon S3 on Outposts. S3 on Outposts uses endpoints to connect to Outposts buckets so that
you can perform actions within your virtual private cloud (VPC).  This action creates an
endpoint and associates it with the specified Outpost.   Related actions include:
DeleteEndpoint     ListEndpoints

# Arguments
- `outpost_id`: The ID of the AWS Outpost.
- `security_group_id`: The ID of the security group to use with the endpoint.
- `subnet_id`: The ID of the subnet in the selected VPC.

"""
create_endpoint(OutpostId, SecurityGroupId, SubnetId; aws_config::AbstractAWSConfig=global_aws_config()) = s3outposts("POST", "/S3Outposts/CreateEndpoint", Dict{String, Any}("OutpostId"=>OutpostId, "SecurityGroupId"=>SecurityGroupId, "SubnetId"=>SubnetId); aws_config=aws_config)
create_endpoint(OutpostId, SecurityGroupId, SubnetId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = s3outposts("POST", "/S3Outposts/CreateEndpoint", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("OutpostId"=>OutpostId, "SecurityGroupId"=>SecurityGroupId, "SubnetId"=>SubnetId), params)); aws_config=aws_config)

"""
    delete_endpoint(endpoint_id, outpost_id)
    delete_endpoint(endpoint_id, outpost_id, params::Dict{String,<:Any})

S3 on Outposts access points simplify managing data access at scale for shared datasets in
Amazon S3 on Outposts. S3 on Outposts uses endpoints to connect to Outposts buckets so that
you can perform actions within your virtual private cloud (VPC).  This action deletes an
endpoint.   Related actions include:    CreateEndpoint     ListEndpoints

# Arguments
- `endpoint_id`: The ID of the end point.
- `outpost_id`: The ID of the AWS Outpost.

"""
delete_endpoint(endpointId, outpostId; aws_config::AbstractAWSConfig=global_aws_config()) = s3outposts("DELETE", "/S3Outposts/DeleteEndpoint", Dict{String, Any}("endpointId"=>endpointId, "outpostId"=>outpostId); aws_config=aws_config)
delete_endpoint(endpointId, outpostId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = s3outposts("DELETE", "/S3Outposts/DeleteEndpoint", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("endpointId"=>endpointId, "outpostId"=>outpostId), params)); aws_config=aws_config)

"""
    list_endpoints()
    list_endpoints(params::Dict{String,<:Any})

S3 on Outposts access points simplify managing data access at scale for shared datasets in
Amazon S3 on Outposts. S3 on Outposts uses endpoints to connect to Outposts buckets so that
you can perform actions within your virtual private cloud (VPC).  This action lists
endpoints associated with the Outpost.   Related actions include:    CreateEndpoint
DeleteEndpoint

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The max number of endpoints that can be returned on the request.
- `"nextToken"`: The next endpoint requested in the list.
"""
list_endpoints(; aws_config::AbstractAWSConfig=global_aws_config()) = s3outposts("GET", "/S3Outposts/ListEndpoints"; aws_config=aws_config)
list_endpoints(params::AbstractDict{String, Any}; aws_config::AbstractAWSConfig=global_aws_config()) = s3outposts("GET", "/S3Outposts/ListEndpoints", params; aws_config=aws_config)


text:
module NonlinearWrapper

using JavaCall

include("ScalarWrapper/ScalarWrapper.jl")
include("VectorWrapper/VectorWrapper.jl")

end


text:
include("core.jl")
include("algorithms.jl")
include("io.jl")
include("readnrrd.jl")
include("readremote.jl")


text:
macro Ax_min()
	code = quote
		# minimum longitudinal acceleration for given speed
		Ax_min = zeros(length(U),1)
		for i in eachindex(U)
			Ax_min[i,1] = AXC[5]*U[i]^3 + AXC[6]*U[i]^2 + AXC[7]*U[i] + AXC[8]
		end
		Ax_min
	end
	return esc(code)
end


text:
export GeneticMap

using Interpolations

const DECODE_GENETIC_MAP = "/odinn/data/extdata/deCODE/GenMap/20150202_hg38/decode.genmap_hg38.gor"

abstract type AbstractGeneticMap end

struct GeneticMap <: AbstractGeneticMap
    path::String
    interp::Dict
end

function GeneticMap(path = DECODE_GENETIC_MAP)
    chromDict = Dict()
    fh = open(path)
    for line in eachline(fh)
        chrom, pos, name, cM = split(chomp(line), "\t")

        if !haskey(chromDict, chrom)
            chromDict[chrom] = (Float64[], Float64[])
        end

        push!(chromDict[chrom][1], parse(Float64, pos))
        push!(chromDict[chrom][2], parse(Float64, cM))
    end

    interp = Dict()
    for (chrom, pos_cM) in chromDict
        interp[chrom] = LinearInterpolation(pos_cM[1], pos_cM[2]; extrapolation_bc = Flat())
    end
    close(fh)

    GeneticMap(path, interp)
end

Base.show(io::IO, gm::GeneticMap) = print(io, "(GeneticMap: ", gm.path, ")")

(gm::GeneticMap)(chrom, pos) = gm.interp[chrom](pos)


text:
using Adapt: adapt

default_collect(A::AbstractArray{T,N}) where {T,N} =
    invoke(collect, Tuple{AbstractArray{T,N}}, A)

Base.collect(A::UnionArray) = default_collect(adapt(Array, A))


text:
#=export A_CODEC,
    V_CODEC,
    CODECS_EACH_RES,
    isaudio,
    name,
    bitrate,
    resolution,
    fps,
    sellect_format,
    sellect_ext=#

const A_CODEC = Dict(
    251 => ("opus", 160),
    250 => ("opus", 70),
    249 => ("opus", 50),
    140 => ("aac", 128),
    139 => ("aac", 48)
)
const V_CODEC = Dict(
    # VP9
    272 => (4320, "vp9", "60HDR"),
    337 => (2160, "vp9", "60HDR"),
    315 => (2160, "vp9", "60"),
    313 => (2160, "vp9", "30"),
    336 => (1440, "vp9", "60HDR"),
    308 => (1440, "vp9", "60"),
    271 => (1440, "vp9", "30"),
    335 => (1080, "vp9", "60HDR"),
    303 => (1080, "vp9", "60"),
    248 => (1080, "vp9", "30"),
    334 => (720, "vp9", "60HDR"),
    302 => (720, "vp9", "60"),
    247 => (720, "vp9", "30"),
    333 => (480, "vp9", "60HDR"),
    244 => (480, "vp9", "30"),
    332 => (360, "vp9", "60HDR"),
    243 => (360, "vp9", "30"),
    331 => (240, "vp9", "60HDR"),
    242 => (240, "vp9", "30"),
    330 => (144, "vp9", "60HDR"),
    278 => (144, "vp9", "30"),
    # AV1
    571 => (4320, "av1", "60"), # higher bitrate than 402
    402 => (4320, "av1", "60"),
    401 => (2160, "av1", "60"),
    400 => (1440, "av1", "60"),
    399 => (1080, "av1", "60"),
    398 => (720, "av1", "60"),
    397 => (480, "av1", "30"),
    396 => (360, "av1", "30"),
    395 => (240, "av1", "30"),
    394 => (144, "av1", "30"),
    # AVC (H.264)
    299 => (1080, "avc", "60"),
    137 => (1080, "avc", "30"),
    298 => (720, "avc", "60"),
    136 => (720, "avc", "30"),
    135 => (480, "avc", "30"),
    134 => (360, "avc", "30"),
    133 => (240, "avc", "30"),
    160 => (144, "avc", "30"),
)
const CODECS_EACH_RES = Dict(
    4320 => (272, 571, 402),
    2160 => (337, 315, 401, 313),
    1440 => (336, 308, 400, 271),
    1080 => (335, 303, 399, 299, 248, 137),
    720  => (334, 302, 398, 298, 247, 136),
    480  => (333, 244, 397, 135),
    360  => (332, 243, 396, 134),
    240  => (331, 242, 395, 133),
    144  => (330, 278, 394, 160)
)

# Codec Utils
isaudio(c) = haskey(A_CODEC, c)

name(c) = isaudio(c) ? A_CODEC[c][begin] : V_CODEC[c][2]

bitrate(c) = A_CODEC[c][end]

resolution(c) = V_CODEC[c][begin]

fps(c) = V_CODEC[c][end]

# ====================================================
# ===================    Codec    ====================
# ====================================================
# codec sellection for native extracted codecs
function sellect_format(codecs_dict::AbstractDict, opts)
    sellect_format(keys(codecs_dict), opts)
end

# codec sellection (compatatible with backend mode)
function sellect_format(cs, opts)
    p = get(opts, :preset, "bestaudio") |>
        x -> replace(x, '-'=>'_')
    if p == "bestaudio"
        a = preset_a_best(cs)
        v = nothing
    elseif p == "best"
        a = preset_a_best(cs)
        v = preset_limited_best(cs)
    else
        a,v = try_preset_func(p, cs)
    end
    return (a=a, v=v)
end

# call user-defined preset function
function try_preset_func(preset, cs)
    f = Symbol(preset)
    if isdefined(@__MODULE__, f)
        Expr(:call, f, cs) |> eval
    else
        return error("""

            Preset function 【$(preset)】 is not defined.
            Call the predefined function below instead.
                best      : best_audio + best_video
                bestaudio : best_audio
            Or, define your preset following the rules:
                - Define in  ~/.julia/config/ytdl.jl
                - Return integer tuple (audio_id, video_id)
                - Have valid name
                    can include a-z, 0-9 and "_"
                    not starts with numbers
                    in lowercase
        """)
    end
end

# ====================================================
# =================    Extension    ==================
# ====================================================

# sellect audio ext
sellect_ext(a) = a=="aac" ? "m4a" : "opus"

# extension for audio codec
function sellect_ext(a_id, v_id::Nothing)
    a = name(a_id)
    ext = sellect_ext(a)
    temp = sellect_temp(a)
    (a=ext, v=nothing, atemp=temp, vtemp=nothing)
end

# extension for audio and video codec
function sellect_ext(a_id, v_id)
    a = name(a_id)
    v = name(v_id)
    aext = sellect_ext(a)
    if a=="aac" && (v=="av1" || v=="avc")
        vext = "mp4"
    elseif a=="opus" && v=="vp9"
        vext = "webm"
    else
        vext = "mkv"
    end
    atemp = sellect_temp(a)
    vtemp = sellect_temp(v)
    (a=aext, v=vext, atemp=atemp, vtemp=atemp)
end

#
function sellect_temp(c)
    if c=="opus" || c=="vp9"
        "webm"
    elseif c=="aac" || c=="av1" || c=="avc"
        "mp4"
    end
end


text:
# branch_basics.jl
# ======================
#
# Types and methods for manipulating the branches of phylogenies and networks.
#
#
# This file is a part of BioJulia. License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

"""
    branchdata{C,B}(tree::Phylogeny{C,B}, edge::Edge)

Getter for metadata associated with branch represented by `edge`.
"""
function branchdata{C,B}(tree::Phylogeny{C,B}, edge::Edge)
    return get(tree.edgedata, edge, B())
end

"""
    branchdata!{C,B}(tree::Phylogeny{C,B}, edge::Edge, data::B)

Setter for metadata associated with branch represented by `edge`.
"""
function branchdata!{C,B}(tree::Phylogeny{C,B}, edge::Edge, data::B)
    if has_edge(tree.graph, edge)
        tree.edgedata[edge] = data
    else
        warn("You tried to assign branch data to a non-existant branch. Consequently, nothing was done.")
    end
    return tree
end

"""
    create_branch!{C,B}(tree::Phylogeny{C,B}, branch::Edge, branchdata::B)

Create a branch between two nodes in a phylogenetic tree.
Metadata `branchdata` will be associated with the branch.
"""
function create_branch!{C,B}(tree::Phylogeny{C,B}, branch::Edge, branchdata::B)
    add_edge!(tree.graph, branch)
    branchdata!(tree, branch, branchdata)
    return tree
end

"""
    create_branch!{C,B}(tree::Phylogeny{C,B}, branch::Edge)

Create a branch between two nodes in a phylogenetic tree.
"""
function create_branch!{C,B}(tree::Phylogeny{C,B}, branch::Edge)
    return add_branch!(tree, branch, B())
end

"""
    destroy_branch!(tree::Phylogeny, branch::Edge)

Destroy a branch between two nodes in a phylogenetic tree.
When the branch is destroyed so will any associated metadata object.
"""
function destroy_branch!(tree::Phylogeny, branch::Edge)
    rem_edge!(tree.graph, src(branch), dst(branch))
    delete!(tree.edgedata, branch)
    return tree
end

function parent_branch(tree::Phylogeny, vertex::Int)
    @assert hasparent(tree, vertex) error("Vertex is not connected to a parent.")
    return in_edges(tree.graph, vertex)[1]
end

function child_branches(tree::Phylogeny, vertex::Int)
    return out_edges(tree.graph, vertex)
end

"""
    branchlength(tree::Phylogeny, edge::Edge)

Get the branchlength of branch defined by `edge` in a phylogeny.
"""
function branchlength(tree::Phylogeny, edge::Edge)
    return branchlength(branchdata(tree, edge))
end

"""
    branchlength!{C,B}(tree::Phylogeny{C,B}, edge::Edge, value::Nullable{Float64})

Set the branchlength of branch defined by `edge` in a phylogeny.

The assumption is metadata values on branches are immutables or value types,
hence the reassignment using the branchdata! method.
"""
@generated function branchlength!{C,B}(tree::Phylogeny{C,B}, edge::Edge, value::Nullable{Float64})
    if !B.mutable
        setting = :(branchdata!(tree, edge, new_branchlength(branchdata(tree, edge), value)))
    else
        setting = :(branchlength!(branchdata(tree, edge), value))
    end

    quote
        $setting
        return tree
    end
end

function branchlength!(tree::Phylogeny, edge::Edge, value::Float64)
    bl = ifelse(value <= 0.0, Nullable{Float64}(), Nullable(value))
    return branchlength!(tree, edge, bl)
end


text:
struct Solution
    ret_val::Int
    primal::Vector{Float64}
    dual::Vector{Float64}
    slack::Vector{Float64}
    objval::Float64
end
Solution() = Solution(0, # SCS_UNFINISHED
                      Float64[], Float64[], Float64[], NaN)

mutable struct Data
    m::Int
    n::Int
    I::Vector{Int}
    J::Vector{Int}
    V::Vector{Float64}
    b::Vector{Float64}
    objconstant::Float64
    c::Vector{Float64}
end

mutable struct Cone
    f::Int # number of linear equality constraints
    l::Int # length of LP cone
    q::Int # length of SOC cone
    qa::Vector{Int} # array of second-order cone constraints
    s::Int # length of SD cone
    sa::Vector{Int} # array of semi-definite constraints
    ep::Int # number of primal exponential cone triples
    ed::Int # number of dual exponential cone triples
    p::Vector{Float64} # array of power cone params
    setconstant::Dict{Int, Float64} # For the constant of EqualTo, LessThan and GreaterThan
    nrows::Dict{Int, Int} # The number of rows of each vector sets
    function Cone()
        new(0, 0,
            0, Int[],
            0, Int[],
            0, 0, Float64[],
            Dict{Int, Float64}(),
            Dict{Int, UnitRange{Int}}())
    end
end


text:
using Mimi

@defcomp grosseconomy begin
    YGROSS  = Variable(index=[time])    # Gross output
    K       = Variable(index=[time])    # Capital
    l       = Parameter(index=[time])   # Labor
    tfp     = Parameter(index=[time])   # Total factor productivity
    s       = Parameter(index=[time])   # Savings rate
    depk    = Parameter()               # Depreciation rate on capital - Note that it has no time index
    k0      = Parameter()               # Initial level of capital
    share   = Parameter()               # Capital share

    function run_timestep(p, v, d, t)
        # Define an equation for K
        if is_first(t)
            v.K[t]  = p.k0  # Note the use of v. and p. to distinguish between variables and parameters
        else
            v.K[t]  = (1 - p.depk)^5 * v.K[t-1] + v.YGROSS[t-1] * p.s[t-1] * 5
        end

        # Define an equation for YGROSS
        v.YGROSS[t] = p.tfp[t] * v.K[t]^p.share * p.l[t]^(1-p.share)
    end
end

@defcomp emissions begin
    E       = Variable(index=[time])    # Total greenhouse gas emissions
    sigma   = Parameter(index=[time])   # Emissions output ratio
    YGROSS  = Parameter(index=[time])   # Gross output - Note that YGROSS is now a parameter

    function run_timestep(p, v, d, t)
        # Define an eqation for E
        v.E[t] = p.YGROSS[t] * p.sigma[t]   # Note the p. in front of YGROSS
    end
end

@Mimi.defmodel model begin

    index[time] = 2015:5:2110

    # Order matters here. If the emissions component were defined first, the model would not run.
    component(grosseconomy)
    component(emissions)

    # Set parameters for the grosseconomy component
    grosseconomy.l = [(1. + 0.015)^t * 6404 for t in 1:20]
    grosseconomy.tfp = [(1 + 0.065)^t * 3.57 for t in 1:20]
    grosseconomy.s = ones(20) * 0.22
    grosseconomy.depk = 0.1
    grosseconomy.k0 = 130.
    grosseconomy.share = 0.3

    # Set parameters for the emissions component
    emissions.sigma = [(1. - 0.05)^t * 0.58 for t in 1:20]

    # Connect parameters
    grosseconomy.YGROSS => emissions.YGROSS
end

# Above macro yields this:
# quote
#     tutorial = (Mimi.Model)()
#     (Mimi.set_dimension!)(tutorial, :time, 2015:5:2110)
#     (Mimi.add_comp!)(tutorial, Main.grosseconomy, :grosseconomy)
#     (Mimi.add_comp!)(tutorial, Main.emissions, :emissions)
#     (Mimi.set_param!)(tutorial, :grosseconomy, :l, [(1.0 + 0.015) ^ t * 6404 for t = 1:20])
#     (Mimi.set_param!)(tutorial, :grosseconomy, :tfp, [(1 + 0.065) ^ t * 3.57 for t = 1:20])
#     (Mimi.set_param!)(tutorial, :grosseconomy, :s, ones(20) * 0.22)
#     (Mimi.set_param!)(tutorial, :grosseconomy, :depk, 0.1)
#     (Mimi.set_param!)(tutorial, :grosseconomy, :k0, 130.0)
#     (Mimi.set_param!)(tutorial, :grosseconomy, :share, 0.3)
#     (Mimi.set_param!)(tutorial, :emissions, :sigma, [(1.0 - 0.05) ^ t * 0.58 for t = 1:20])
#     (Mimi.connect_param!)(tutorial, :emissions, :YGROSS, :grosseconomy, :YGROSS)
# end

run(model)

# Show model results
model[:emissions, :E]

# Or, as a DataFrame
getdataframe(model, :emissions, :E)


text:
function octree(charts::Vector{S} where {S <: CompScienceMeshes.Simplex})

    ncells = length(charts)

    T = coordtype(charts[1])
    P = eltype(charts[1].vertices)

    points = zeros(P, ncells)
    radii = zeros(T, ncells)

    for (i,ch) in enumerate(charts)
        points[i] = cartesian(center(ch))
        radii[i] = maximum(norm(v-points[i]) for v in ch.vertices)
    end

    return Octree(points, radii)
end

"""
    Hi
"""
function grideval(points, coeffs, basis; type=nothing)

    # charts: active charts
    # ad: assembly data (active_cell_idx, local_shape_idx) -> [dof1, dfo2, ...]
    # ag: active_cell_idx -> global_cell_idx
    charts, ad, ag = assemblydata(basis)
    refs = refspace(basis)

    V = valuetype(refs, eltype(charts))
    T = promote_type(eltype(coeffs), eltype(V))
    P = similar_type(V, T)

    type != nothing && (P = type)

    values = zeros(P, size(points))

    chart_tree = BEAST.octree(charts)
    for (j,point) in enumerate(points)
        i = CompScienceMeshes.findchart(charts, chart_tree, point)
        if i != nothing
            # @show i
            chart = charts[i]
            u = carttobary(chart, point)
            vals = refs(neighborhood(chart,u))
            for r in 1 : numfunctions(refs)
                for (m,w) in ad[i, r]
                    values[j] += w * coeffs[m] * vals[r][1]
                end
            end
        end
    end
    return values
end


text:
include("../examples/baker.jl")

include("get_objective.jl")
using PyPlot
using JLD
function plot_obj_erg_avg(ind)
	X = load("../data/obj_erg_avg/cos4y_s$(ind)_s3.jld")
	s4 = X["s$ind"]
	J = X["J"]
	fig, ax = subplots(1,1)
	ax.plot(s4, J, ".", ms=4.0)
	ax.grid(true)
	ax.set_xlim([minimum(s4), maximum(s4)])
	ax.xaxis.set_tick_params(labelsize=28)
	ax.yaxis.set_tick_params(labelsize=28)
	ax.set_xlabel("\$ s_$(ind), s_3 \$", fontsize=28)
	ax.set_ylabel(L"\langle \cos(4 x_2)\rangle", fontsize=28)
end
function plot_obj_fun()

	n_xpts, n_ypts = 100, 100
	eps = 1.e-8
	x_pts = LinRange(0.,2*pi-eps, n_xpts)
	y_pts = LinRange(0.,2*pi-eps, n_ypts)
	x_g = repeat(x_pts, 1, n_ypts)
	y_g = repeat(y_pts', n_xpts, 1)
	J_g = zeros(n_xpts, n_ypts)

	nSteps = 10000
	s = zeros(4)
	s[2] = 1.0
	for n = 1:nSteps
		J_g .+= obj_fun(x_g, y_g)/nSteps
		for j = 1:n_ypts, i = 1:n_xpts
			x_g[i,j], y_g[i,j] = step(
						[x_g[i,j], y_g[i,j]],
						s, 1)[:,end]
		end
	end
	@show minimum(J_g), maximum(J_g)
	fig, ax = subplots(1,1)
	cm = plt.get_cmap("coolwarm")
	J_g = J_g[:]
	clrs = cm(abs.(J_g)/maximum(abs.(J_g)))
	cplot = ax.scatter(x_g[:], y_g[:],
					   c=clrs, cmap=cm)
	ax.xaxis.set_tick_params(labelsize=28)
	ax.yaxis.set_tick_params(labelsize=28)
	#cbar = fig.colorbar(cplot, ax=ax)
	#cbar.ax.tick_params(labelsize=28)
	ax.axis("scaled")

end

max_stars_repo_name: nishaChandramoorthy/linearResponse

max_issues_repo_name: nishaChandramoorthy/linearResponse

max_forks_repo_name: nishaChandramoorthy/linearResponse


text:
#unum-multiplication.jl
#does multiplication for unums.

doc"""
  `Unums.frac_mul!(carry, ::Unum, fraction)`
  multiplies fraction into the the fraction value of unum.
"""
function frac_mul!{ESS,FSS}(a::UnumSmall{ESS,FSS}, multiplier::UInt64)
  (carry, a.fraction, ubit) = i64mul(a.fraction, multiplier, Val{FSS})
  a.flags |= ubit
  return carry + o64
end

frac_mul!{ESS,FSS}(a::UnumLarge{ESS,FSS}, multiplier::ArrayNum{FSS}) = frac_mul!(a, multiplier, Val{__cell_length(FSS)}, Val{true})
function frac_mul!{ESS,FSS, cells, addin}(a::UnumLarge{ESS,FSS}, multiplier::ArrayNum{FSS}, ::Type{Val{cells}}, ::Type{Val{addin}})
  (carry, ubit) = i64mul!(a.fraction, multiplier, Val{cells}, Val{addin})
  a.flags |= ubit
  return carry + o64
end

doc"""
  `Unums.mul(::Unum, ::Unum)` outputs a Unum OR Ubound corresponding to the product
  of two unums.  This is bound to the (\*) operation if options[:usegnum] is not
  set.  Note that in the case of degenerate unums, mul may change the bit values
  of the individual unums, but the values will not be altered.
"""
@universal function mul(a::Unum, b::Unum)
  #some basic checks out of the gate.

  (is_nan(a) || is_nan(b)) && return nan(U)
  is_zero(a) && return is_inf(b) ? nan(U) : zero(b)
  is_zero(b) && return is_inf(a) ? nan(U) : zero(a)

  result_sign = (a.flags $ b.flags) & UNUM_SIGN_MASK

  is_unit(a) && return coerce_sign!(copy(b), result_sign)
  is_unit(b) && return coerce_sign!(copy(a), result_sign)

  is_inf(a) && return inf(U, result_sign)
  is_inf(b) && return inf(U, result_sign)

  #resolve degenerate conditions in both A and B before calculating the exponents.
  resolve_degenerates!(a)
  resolve_degenerates!(b)

  if (is_exact(a) && is_exact(b))
    mul_exact(a, b, result_sign)
  else
    mul_inexact(a, b, result_sign)
  end
end

@universal function mul_exact(a::Unum, b::Unum, result_sign::UInt16)
  #first, calculate the exponents.
  _asubnormal = is_subnormal(a)
  _bsubnormal = is_subnormal(b)

  _aexp = decode_exp(a)
  _bexp = decode_exp(b)

  #solve the exponent.
  exponent = _aexp + _bexp

  (exponent > max_exponent(ESS)) && return mmr(U, result_sign)
  (exponent < min_exponent(ESS,FSS) && return sss(U, result_sign))

  exponent += _asubnormal * 1 + _bsubnormal * 1

  #TODO: convert this into "supranormal form"

  #solve the fractional product.
  if _asubnormal
    result = copy(a)
    multiplicand = b
    is_ulp(b) && make_ulp!(result)

    exponent -= normalize!(result)

    #this code is only necessary when ESS == 0, because you can have subnormal 1's
    #which don't exponent-degrade.
    if (ESS == 0)
      if _bsubnormal
        multiplicand = copy(b)
        exponent -= normalize!(multiplicand)
      end
    end
  else
    result = copy(b)
    multiplicand = a
    is_ulp(a) && make_ulp!(result)

    if _bsubnormal
      #normalize b (regardless of exponent)
      exponent -= normalize!(result)
    end
  end

  #check for a really small exponent (again!)
  (exponent < min_exponent(ESS,FSS) && return sss(U, result_sign))

  carry = frac_mul!(result, multiplicand.fraction)

  if (FSS < 6)
    result.fraction &= mask_top(FSS)
  end

  exponent = resolve_carry!(carry, result, exponent)

  if (exponent < min_exponent(ESS))
    #in the case that we need to make this subnormal.
    right_shift = to16(min_exponent(ESS) - exponent)
    frac_rsh_underflow_check!(result, right_shift)
    frac_set_bit!(result, right_shift)
    result.esize = max_esize(ESS)
    result.exponent = z64
  elseif (exponent > max_exponent(ESS))
    #possible to exceed max_exponent via the additive portios fo the multiplication.
    return mmr(U)
  end

  is_exact(result) && exact_trim!(result)
  trim_and_set_ubit!(result)

  return coerce_sign!(result, result_sign)
end

@universal function mul_inexact(a::Unum, b::Unum, result_sign::UInt16)
  is_inf_ulp(a) && (return inf_ulp_mult(a, b, result_sign))
  is_inf_ulp(b) && (return inf_ulp_mult(b, a, result_sign))

  is_zero_ulp(a) && return zero_ulp_mult(a, b, result_sign)
  is_zero_ulp(b) && return zero_ulp_mult(b, a, result_sign)

  inner_result = mul_exact(a, b, result_sign)
  inner_result.fsize = max_fsize(FSS)

  is_exact(inner_result) && outer_ulp!(inner_result)

  result_exponent = decode_exp(inner_result) + is_subnormal(inner_result) * 1

  outer_a = is_ulp(a) ? outer_exact(a) : a
  outer_b = is_ulp(b) ? outer_exact(b) : b

  outer_result = mul_exact(outer_a, outer_b, result_sign)
  is_exact(outer_result) && inner_ulp!(outer_result)

  if is_ulp(inner_result) && is_ulp(outer_result)
    return (result_sign == z16) ? resolve_as_utype!(inner_result, outer_result) : resolve_as_utype!(outer_result, inner_result)
  else
    return (result_sign == z16) ? B(inner_result, outer_result) : B(outer_result, inner_result)
  end
end

@universal function allsignedreals(result_sign::UInt16, T::Type{Unum})
  (result_sign == 0) ? B(cheapest_zero_ulp(U, z16), cheapest_inf_ulp(U, z16)) : B(cheapest_inf_ulp(U, UNUM_SIGN_MASK), cheapest_zero_ulp(U, UNUM_SIGN_MASK))
end

@universal function inf_ulp_mult(inf_ulp::Unum, b::Unum, result_sign::UInt16)
  #it is possible to multiply by a number bigger than but very close to one.
  (decode_exp(b) > 0) && return mmr(U, result_sign) #mmr is lucky enough to have this cheat.

  is_zero_ulp(b) && return allsignedreals(result_sign, U)

  inner_result = mul_exact(inner_exact(inf_ulp), b, result_sign)
  is_exact(inner_result) && outward_ulp!(inner_result)
  is_inf_ulp(inner_result) && return inner_result

  return (result_sign == 0) ? resolve_as_utype!(inner_result, mmr(U, result_sign)) :
                              resolve_as_utype!(mmr(U, result_sign), inner_result)
end

@universal function zero_ulp_mult(zero_ulp::Unum, b::Unum, result_sign::UInt16)

  _ulp_exact = outer_exact(zero_ulp)
  _val_exact = is_exact(b) ? b : outer_exact(b)
  res_exp = decode_exp(_ulp_exact) + 1 * is_subnormal(_ulp_exact) + decode_exp(_val_exact) + 1 * is_subnormal(_val_exact)

  (res_exp < min_exponent(ESS, FSS)) && return sss(U, result_sign)
  #cannot be inf_ulp because that case is handled above.

  outer_value = mul_exact(_ulp_exact, _val_exact, result_sign)

  is_exact(outer_value) && inner_ulp!(outer_value)
  is_zero_ulp(outer_value) && return outer_value
  (result_sign == z16) ? resolve_as_utype!(sss(U, result_sign), outer_value) : resolve_as_utype!(outer_value, sss(U, result_sign))
end

#import the Base add operation and bind it to the add and add! functions
import Base.*
@bind_operation(*, mul)

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
abstract KalmanFilter

abstract LinearKalmanFilter <: KalmanFilter

abstract Model

abstract ObservationModel

abstract AbstractState

type State{T} <: AbstractState
    x::Vector{T}
    p::Matrix
end

Base.(:(==))(x1::State,x2::State) = x1.x==x2.x && x1.p == x2.p

type LinearModel <: Model
    a::Matrix
    g::Matrix
    q::Matrix
end

## This is the first function a Kalman filter needs to
# implement: Apply its model to its state to generate a
# new state

function ap(f::LinearModel,x::State)
    x1 = f.a*x.x
    p1 = f.a*x.p*f.a' + f.g*f.q*f.g'
    State(x1,p1)
end

type Observation{T}
    y::Vector{T}
end

Base.convert(::Type{Observation},y) = Observation([y])

type LinearObservationModel <: ObservationModel
    h::Matrix
    r::Matrix
end

type BasicKalmanFilter <: LinearKalmanFilter
    x::State
    f::LinearModel
    z::LinearObservationModel
end

Base.copy(kf::KalmanFilter) = deepcopy(kf)

## This is the second of two functions a Kalman filter needs to implement
# Return three matrices:
# 1. The residual (y-Hx) (`res`)
# 2. The cross-covariance between state and measurement (`ph`)
# 3. The innovation covariance matrix
function covs(kf::BasicKalmanFilter,y::Observation)
    res = y.y - kf.z.h * kf.x.x
    ph = kf.x.p * kf.z.h'
    s = kf.z.h * ph + kf.z.r
    (res,ph,s)
end


text:
export min_time,
    appropriate_vegetation,
    min_acceleration,
    min_cost,
    clear_state,
    clear_state_region,
    managed_clear_eutrophic,
    make_decision!

##############################################################
# Predefined objective functions
##############################################################

"""
    min_time(model)

Objective function that returns the time of the model at the end of a run. If a target
function is interested in moving from one state to the next in the quickest amount of
time, this is a useful objective.
"""
min_time(model::ABM) = model.lake.t / 365.0

"""
    min_acceleration(model)

Objective function that returns the sum of the absolute value of the second derivative of
all lake variables. Span is from start of the optimisation to the final `model.year`
with monthly increments.

Helps to mitigate large spikes in transitions.
"""
min_acceleration(model::ABM) = sum(abs.(model.lake.sol(model.lake.sol.t, Val{2})))

"""
    min_cost(model)

Objective function that returns a "cost" of future `Planting` and `Trawling`
interventions, which ultimately is just a sum of all proposed rates.

!!! note

    When using `nutrient_series = Noise(...)` and only the `min_cost` objective,
    `opt_replicates` *must* be large (`≳ 10`), since the noise process may cause the
    optimiser to identify a solution in theory but fail in practice. The solution to
    failing meet the target scenario in this case is increase `opt_replicates` and retry.
"""
function min_cost(model::ABM)
    # By the time we get here, we're looking at interventions still to do (some may have already been completed, but they have already been budgeted.
    # This decision is only for however much needs to be done in the future.

    # Our budget is not pinned to money directly, but there is an assumed cost based on years to implement at the given rate.
    budget = 0.0

    # Get all Trawling and Planting actions. We assume there is no municipal cost for WastewaterTreatment -> that cost is shifted to the Household for now
    # We also assume no (monetary) cost for Angling, which neglects revenue raised by licensing for the moment.
    for municipality in municipalities(model)
        interventions = vcat(values(municipality.interventions)...)
        plant = Iterators.filter(i -> i isa Planting, interventions)
        trawl = Iterators.filter(i -> i isa Trawling, interventions)

        isempty(plant) || (budget += sum(p -> p.rate * p.cost, plant))
        isempty(trawl) || (budget += sum(t -> t.rate * t.cost, trawl))
    end
    # Above assumes weighting is different for municipalities. If that ends up not to be the case, we can grab all active interventions regardless of municipality with
    # `interventions = merge(vcat, (m.interventions for m in municipalities(model))...)`
    # To simplify the above loop.

    return budget
end

"""
    appropriate_vegetation(model)

Objective function that returns a penalty if vegetation is higher than an operational
density of 60. Higher densities cause recreational issues that are considered
unacceptable.
"""
appropriate_vegetation(model::ABM) = sum(model.lake.sol[3, model.lake.sol[3, :].>60.0])

##############################################################
# Predefined target functions
##############################################################

"""
    clear_state(model, s)

Targets the clear lake steady state, with an extra stopping condition if that state was
not reached within 100 years.

Calculated via instantaneous comparisons at latest model time of all lake variables, with
an additional check to verify a near-zero first derivative.

**NOTE:** This target is hard coded to the default `Martin` parameters.
"""
function clear_state(model, s)
    s >= 100 && return true # Escape if we dont converge after 100 years

    # We also want the system to stabilise a bit, so we wait until the derivatives calm down too.
    # That cutoff can be a complication with noisy nutrients, so we relax the criteria in that case.
    if model.nutrient_series isa Noise &&
       hasfield(typeof(model.nutrient_series.process.dist), :σ)
        cutoff = model.nutrient_series.process.dist.σ
    else
        cutoff = 5e-4
    end
    sum(abs.(interpolated_position(model.lake.p.nutrients, Clear) .- model.lake.u,)) <
    cutoff && sum(abs.(model.lake.sol(model.lake.t, Val{1}))) < cutoff
end

"""
    clear_state_region(model, s)

Targets the region around the clear lake steady state, with an extra stopping condition if that state was
not reached within 100 years.

This method does not expect the solution to reach the exact value of the stable bifurcation like
[`clear_state`](@ref), instead it will exit once various conditions are met that define
the clear state region, namely: `V >= 33.4, B < 31.6 and if n > 1 then P > 0.5`.

Calculated via instantaneous comparisons at latest model time of all lake variables, with
an additional check to verify a near-zero first derivative.
"""
function clear_state_region(model, s)
    s >= 100 && return true # Escape if we dont converge after 100 years

    B, P, V = model.lake.u

    # Based on bifurcation analysis: for a clear state, we need
    # V >= 33.4, B < 31.6 and if n > 1 then P > 0.5
    # We also want the system to stabilise a bit, so we wait until the derivatives calm down too.
    # That cutoff can be a complication with noisy nutrients, so we relax the criteria in that case.
    if model.nutrient_series isa Noise &&
       hasfield(typeof(model.nutrient_series.process.dist), :σ)
        cutoff = model.nutrient_series.process.dist.σ
    else
        cutoff = 5e-4
    end
    B < 31.6 &&
        V >= 33.4 &&
        (model.lake.p.nutrients > 1.0 ? P > 0.5 : true) &&
        sum(abs.(model.lake.sol(model.lake.t, Val{1}))) < cutoff
end

"""
    managed_clear_eutrophic(model, s)

Targets the `T3` state, which is a high nutrient (`N>=3`), unstable state with a high
pike population. Will stop at 100 years if not successful.

**Note:** For the moment this targets the region of T3, not the explicit starting point.
"""
function managed_clear_eutrophic(model, s)
    s >= 100 && return true # Escape if we dont converge after 100 years

    B, P, V = model.lake.u

    P >= 3.5 && V > B
end

##############################################################
# Optimisation capacity
##############################################################

"""
    create_test_model(model)

Creates a complete copy of the current model, with a modified set of interventions.
All "completed" interventions (i.e. ones that have happened in the models' past) are
ignored. This test model is then used in the optimisation procedure.
"""
function create_test_model(model::ABM)
    m = deepcopy(model)

    push!(m.properties, :test => true)

    # Drop previous (completed) interventions.
    for municipality in municipalities(m)
        newi = Dict{Int,Vector{Intervention}}()
        for (k, v) in municipality.interventions
            if k == -1
                push!(newi, k => v)
            elseif k - model.year >= 0 && (
                (
                    model.policy.current_term_only &&
                    k <= model.year + model.policy.every - 1
                ) || !model.policy.current_term_only
            )
                # We also drop future interventions that are outside
                # the policy window if such a flag is active.
                push!(newi, k - model.year => v)
            end
        end
        municipality.interventions = newi
    end
    return m
end

"""
    apply_policies!(x, test::ABM)

Applies new test values to all active intervention properties. Is used before each
optimisation call, but also to finalise the decision process via a `best_candidate`
call.
"""
function apply_policies!(x, test)
    xs = Iterators.Stateful(x)

    for municipality in municipalities(test)
        # Update intervention values with new x value
        for intervention in Iterators.flatten(values(municipality.interventions))
            # We must do this in the order of policy range
            if haskey(municipality.policies, typeof(intervention))
                targets = keys(municipality.policies[typeof(intervention)])
                assign = collect(zip(targets, Iterators.take(xs, length(targets))))
                map(t -> setproperty!(intervention, t[1], t[2]), assign)
            end
        end
    end
    return nothing
end

"""
    update_true_model!(test, model)

Reconstructs each municipalities `interventions` property in the model after a successful
optimisation run. Be careful with the model order here, as there's no simple way to
differentiate the two of them via types.
"""
function update_true_model!(test::ABM, model::ABM)
    # Both models have the same municipality ids
    for mid in (m.id for m in municipalities(model))
        for (k, v) in test[mid].interventions
            if k == -1
                # "Always on" values can be replaced directly
                model[mid].interventions[k] = test[mid].interventions[k]
            else
                model[mid].interventions[k+model.year] = test[mid].interventions[k]
            end
        end
    end
    return nothing
end

"""
    cost(x, u0, p, test::ABM)

Optimization function. To be used in in conjunction with the `bboptimize` call only.
Needs to be overloaded slightly, since we want to reset the lake dynamics at each call,
so `model.lake.u` and `model.lake.p` are expected to be passed into the second and third
variables respectively, with the cut down "test" version of the model being the last
value.
"""
function cost(x, u0::Vector{Float64}, p::L, test::ABM) where {L<:LakeParameters}
    # Re-initialise certain parameters we need to satisfy in the optimisation
    prob = OrdinaryDiffEq.ODEProblem(lake_dynamics!, u0, (0.0, Inf), p)

    test.lake = OrdinaryDiffEq.init(prob, OrdinaryDiffEq.Tsit5())
    if haskey(test.properties, :nutrient_stabilise)
        test.nutrient_stabilise -= test.year
    end
    test.year = 0
    test.init_nutrients = p.nutrients
    test.init_pike_mortality = p.mp

    # Search range has everything from the policies section now in search.
    # This means that x is the values associated with all of those results.
    # `test` gives us a copy of model, which we can now edit with the new values in x for
    # each parameter
    apply_policies!(x, test)

    if test.policy.opt_replicates > 0
        results = NTuple{length(test.policy.objectives),Float64}[]
        for _ in 1:test.policy.opt_replicates
            push!(results, calculate_objectives(deepcopy(test)))
        end
        #results = Agents.Distributed.pmap(
        #j -> calculate_objectives(deepcopy(test)),
        #1:test.policy.opt_replicates,
        #)
        return Tuple(mean.(Iterators.zip(results...)))
    else
        return calculate_objectives(test)
    end
end

function calculate_objectives(test)
    if test.nutrient_series isa Noise
        # Don't assume we use the same seed
        if typeof(test.rng) isa Random.RandomDevice
            @set test.rng = Random.RandomDevice()
        else
            seed!(test)
        end
        Random.seed!(test.nutrient_series.process.rng)
    end

    Agents.step!(test, agent_step!, model_step!, test.policy.target)
    objectives = first.(test.policy.objectives)
    if test.policy.target(test, 1)
        return map(o -> o(test), objectives)
    else
        # Penalise this result, as it failed to
        # reach the target before cutoff.
        return map(o -> (o(test) * 10.0) + 10.0, objectives)
    end
end

"""
    make_decision!(model)

Runs the optimisation routine, calling on policy ranges set via [`policy`](@ref).
Decisions are made only from the year of the call onwards.
"""
function make_decision!(model::ABM)
    # Create our test model outside of the loop, it will therefore be cut down to
    # appropriate points already.
    test = create_test_model(model)

    objective_dimension = length(test.policy.objectives)
    w = last.(test.policy.objectives)

    # We only need the tuples for the search range here. So long as we use the same
    # scheduler, there's no need to worry about anything else on the other side
    # of the function guard. The only nuance is to iterate through each municipality
    # in order so that we replicate the ranges of everything correctly.
    search = Tuple{Float64,Float64}[]
    x0 = Float64[]
    for municipality in municipalities(test)
        for intervention in Iterators.flatten(values(municipality.interventions))
            if haskey(municipality.policies, typeof(intervention))
                append!(search, collect(municipality.policies[typeof(intervention)]))
                push!(x0, intervention.rate)
            end
        end
    end

    # No need to optimise if there are no more interventions
    isempty(search) && return nothing

    if model.policy.trace_mode != :silent
        word = model.year == model.policy.start ? "Starting" : "Adjusting"
        println("$(word) policy decisions in year $(model.year)")
    end

    result = bboptimize(
        x -> cost(x, model.lake.u, model.lake.p, test),
        x0;
        Method = :borg_moea,
        FitnessScheme = ParetoFitnessScheme{objective_dimension}(
            is_minimizing = true,
            aggregator = f -> weightedfitness(f, w),
        ),
        SearchRange = search,
        MaxSteps = model.policy.max_steps,
        MaxTime = model.policy.max_time,
        TraceMode = model.policy.trace_mode,
        TraceInterval = 10,
    )

    x = best_candidate(result)

    apply_policies!(x, test)
    update_true_model!(test, model)

    return nothing
end

"""
    weightedfitness(f, weights)

Handle objective weights. Should be used as an `aggregator` in `bboptimize`.
"""
weightedfitness(f, w) = sum(f .+ w)