text:
import Pkg;
Pkg.activate(".")
path = pwd() * "/functions"
push!(LOAD_PATH,path)
# %%
using JLD2, FileIO
using Plots
using LaTeXStrings
using Optim
using StatsBase
pyplot()
# %%
using SharedArrays
using Distributed
addprocs()
# %%
@everywhere begin
import Pkg;
Pkg.activate(".")
path = pwd() * "/functions"
push!(LOAD_PATH,path)
end
# %%
@everywhere begin
using Optim
using SMM
using Macros
using CompStat
using CustomTypes
end
# %%
@load "output/calibrations.jld2" cal_y
cal_y.ϵ = 17/((2013-2005+1)*12)
ss = RegUncSS()
dual = DualSS()
# %%
targets = [
	4/3, 			# F/wp
	0.83,			# jcf / jc
	(1-0.12)*(1-0.26),	# np
	0.12*(1-0.26), 		# nf
	0.7, 			# s / ξ
	1-(1-0.7)^(1/3)		# q
]'
# %%
N = 1000
D = 1000
init_param = deepcopy([ cal_y.s, cal_y.α, cal_y.m, cal_y.b, cal_y.γ, cal_y.F ])
shocks = Vector{Vector{Float64}}(undef, N)
states = Vector{Vector{Int64}}(undef, N)
for k in 1:N
	shocks[k] = Float64[]
	states[k] = Int64[]
	simulate!(shocks[k], states[k], cal_y, D)
end
# %%
n_Δ = 20
grid_Δ = range(0.,stop=0.025,length=n_Δ)
grid_param = SharedArray{Float64,3}((n_Δ, N, 6))
grid_np_sim = SharedArray{Float64,2}((n_Δ, N))
grid_nf_sim = SharedArray{Float64,2}((n_Δ, N))
grid_flag = SharedArray{Bool,2}((n_Δ, N))
grid_np = zeros((n_Δ, 3))
grid_nf = zeros((n_Δ, 3))
init_param = deepcopy([ cal_y.s, cal_y.α, cal_y.m, cal_y.b, cal_y.γ, cal_y.F ])
# %%
@everywhere begin
	cal_y = $cal_y
	targets = $targets
	shocks = $shocks
	states = $states
	ss = $ss
	dual = $dual
	init_param = $init_param
	D = $D
end
# %%
# Non distributed version
# for i in 1:n_Δ
# 	println("$i")
# 	@time begin
# 		for j in 1:N
# 			res = optimize(	x-> loss!(x, cal_y, ss, shocks[j], states[j], grid_Δ[i], D, targets), init_param, NelderMead())
# 			grid_param[i,j,:] .= res.minimizer
# 			grid_flag[i,j] = res.x_converged
# 			@load_vec_into_struct res.minimizer cal_y s α m b γ F
# 			fill_DualSS!(dual, cal_y, [0.9, 1.3])
# 			grid_np_sim[i,j], grid_nf_sim[i,j] = (dual.np, dual.nf)
# 		end
# 		grid_np[i,1] = mean(grid_np_sim[i,:])
# 		grid_nf[i,1] = mean(grid_nf_sim[i,:])
# 		grid_np[i,2:3] .= quantile(grid_np_sim[i,:],[0.025,0.975])
# 		grid_nf[i,2:3] .= quantile(grid_nf_sim[i,:],[0.025,0.975])
# 		init_param .= reshape(mean(grid_param[i,:,:], dims = 1),6)
# 	end
# end
# %%
for i in 1:n_Δ
	println("$i")
	@time begin
		@sync @distributed for j in 1:N
			res = optimize(	x-> loss!(x, cal_y, ss, shocks[j], states[j], grid_Δ[i], D, targets), init_param, NelderMead())
			grid_param[i,j,:] .= res.minimizer
			grid_flag[i,j] = Optim.converged(res)
			@load_vec_into_struct res.minimizer cal_y s α m b γ F
			fill_DualSS!(dual, cal_y, [0.9, 1.3])
			grid_np_sim[i,j], grid_nf_sim[i,j] = (dual.np, dual.nf)
		end
		grid_np[i,1] = mean(grid_np_sim[i,:])
		grid_nf[i,1] = mean(grid_nf_sim[i,:])
		grid_np[i,2:3] .= quantile(grid_np_sim[i,:],[0.025,0.975])
		grid_nf[i,2:3] .= quantile(grid_nf_sim[i,:],[0.025,0.975])
		init_param .= reshape(mean(grid_param[i,:,:], dims = 1),6)
		@everywhere init_param = $init_param
	end
end
# %%
JLD2.@save "output/smm.jld" grid_Δ grid_flag grid_param grid_np_sim grid_nf_sim grid_np grid_nf D
# %%
@load "output/smm.jld" grid_Δ grid_param grid_np_sim grid_nf_sim grid_np grid_nf D
# %%
p1 = plot(100*grid_Δ,100*(grid_np[:,1] .- targets[3])./targets[3],lab=L"$n^p$", color = :black, line=(:solid, 2))
plot!(p1, 100*grid_Δ,100*(grid_np[:,2:3] .- targets[3])./targets[3],lab="", color = :black, line=(:dash, 2))
ylabel!("Deviation w.r.t to steady-state (%)")
p2 = plot(100*grid_Δ,100*(grid_nf[:,1] .- targets[4])./targets[4],lab=L"$n^f$", color = :black, line=(:solid, 2))
plot!(p2, 100*grid_Δ,100*(grid_nf[:,2:3] .- targets[4])./targets[4],lab="", color = :black, line=(:dash, 2))
plot(p1,p2,layout=(1,2))
xgrid!(:off)
ygrid!(:off)
xlabel!(L"100 \mid \frac{F_i - F}{F} \mid")
# %%
savefig("figures/risk_pol_smm.pdf")
# %%
@load "output/smm.jld" grid_Δ grid_param grid_np_sim grid_nf_sim grid_np grid_nf D
# %%
p1 = plot(100*grid_Δ,100*(grid_np[:,1] .- targets[3])./targets[3],lab="Open-ended", color = :red, line=(:solid, 1))
plot!(p1, 100*grid_Δ,100*(grid_np[:,2:3] .- targets[3])./targets[3],lab="", color = :red, line=(:dash, 1))
ylabel!("Deviation w.r.t the stochastic steady-state (%)")
p2 = plot(100*grid_Δ,100*(grid_nf[:,1] .- targets[4])./targets[4],lab="Fixed-term", color = :blue, line=(:solid, 1))
plot!(p2, 100*grid_Δ,100*(grid_nf[:,2:3] .- targets[4])./targets[4],lab="", color = :blue, line=(:dash, 1))
plot(p1,p2,layout=(1,2))
xgrid!(:off)
ygrid!(:off)
xlabel!(L"100 \mid \frac{F_i - F}{F} \mid")
# %%
savefig("figures/risk_pol_smm_c.pdf")


text:
# Check out mosfet.jl for more Julia-like usage

import NgSpice
n = NgSpice

netpath = joinpath(@__DIR__, "..", "inputs", "ac_ascii.cir") |> normpath

n.init()
n.cmd("source $netpath")
n.cmd("run")
n.cmd("display")
#n.cmd("<Standard Ngspice commands>")

ext: jl

max_stars_repo_name: JuliaComputing/Ngspice

max_issues_repo_name: JuliaComputing/Ngspice

max_forks_repo_name: JuliaComputing/Ngspice


text:
module Structs

using Mixtape
import Mixtape: CompilationContext, transform, allow, show_after_inference,
                show_after_optimization, debug, @load_call_interface
using MacroTools
using BenchmarkTools

struct Foo
    x::Float64
end

struct Bar
    x::Float64
end

function f(x::Foo, y::Float64)
    return Foo(x.x + y)
end

struct MyMix <: CompilationContext end

allow(ctx::MyMix, m::Module) = m == Structs
show_after_inference(ctx::MyMix) = false
show_after_optimization(ctx::MyMix) = false
debug(ctx::MyMix) = true

swap(e) = e
function swap(e::Expr)
    new = MacroTools.postwalk(e) do s
        isexpr(s, :call) || return s
        s.args[1] == Foo || return s
        return Expr(:call, Bar, s.args[2:end]...)
    end
    return new
end

function transform(::MyMix, b)
    for (v, st) in b
        replace!(b, v, swap(st))
    end
    display(b)
    return b
end

# Mixtape cached call.
Mixtape.@load_call_interface()
display(call(MyMix(), f, Foo(10.0), 5.0))
@btime call(MyMix(), f, Foo(10.0), 5.0)

end # module


text:
using Documenter
using Accounts
using Dates

makedocs(;
    modules=[Accounts],
    format=Documenter.HTML(),
    pages=[
        "Accounts" => "index.md",
        "1 - API" => "api.md",
        "2 - Infrastructure" => "infrastructure.md",
        "3 - Examples" => "example.md"
    ],
    sitename="Accounts.jl",
)


text:
import Statistics
import RestrictedBoltzmannMachines as RBMs
using Test: @test, @testset, @inferred
using Statistics: mean, var, cov
using LinearAlgebra: dot
using RestrictedBoltzmannMachines: repeat_size, sizedims, moving_average

@testset "two" begin
    @test RBMs.two(1) === RBMs.two(Int) === 2
    @test RBMs.two(Int8(1)) === RBMs.two(Int8) === Int8(2)
    @test RBMs.two(1f0) === RBMs.two(Float32) === 2f0
    @test RBMs.two(1.0) === RBMs.two(Float64) === 2.0
    @test_throws InexactError RBMs.two(Bool)
    @inferred RBMs.two(1)
end

@testset "inf" begin
    @test_throws InexactError RBMs.inf(1)
    @test Inf === @inferred RBMs.inf(1.0)
    @test Inf32 === @inferred RBMs.inf(1f0)
    @inferred RBMs.inf(1.0)
end

@testset "generate_sequences" begin
    @test collect(RBMs.generate_sequences(2, 1:3)) == reshape(
           [
               [1, 1], [2, 1], [3, 1],
               [1, 2], [2, 2], [3, 2],
               [1, 3], [2, 3], [3, 3]
           ],
           3, 3
        )
end

@testset "broadlike" begin
    A = randn(1,3)
    B = randn(2,1)
    @test RBMs.broadlike(A, B) ≈ A .+ B .- B
    @inferred RBMs.broadlike(A, B)
    @test RBMs.broadlike(A) == A
    @test RBMs.broadlike(A, 1) == A
end

@testset "wmean" begin
    A = randn(5)
    w = rand(5)
    @test dot(A, w) / sum(w) ≈ @inferred RBMs.wmean(A; wts=w)

    A = randn(4,3,5,2)
    @test mean(A) ≈ @inferred RBMs.wmean(A)
    @test mean(A; dims=(2,4)) ≈ @inferred RBMs.wmean(A; dims=(2,4))

    wts = rand(size(A)...)
    @test sum(A .* wts) ./ sum(wts) ≈ @inferred RBMs.wmean(A; wts)
    @test sum(A .* wts) ./ sum(wts) ≈ @inferred RBMs.wmean(A; wts, dims=:)

    wts = rand(3,2)
    @test sum(reshape(wts,1,3,1,2) .* A; dims=(2,4)) ./ sum(wts) ≈ @inferred RBMs.wmean(A; dims=(2,4), wts)

    wts = rand(2)
    @test sum(reshape(wts,1,1,1,2) .* A; dims=4) ./ sum(wts) ≈ @inferred RBMs.wmean(A; dims=4, wts)
end

@testset "wsum" begin
    A = randn(5)
    w = rand(5)
    @test dot(A, w) ≈ @inferred RBMs.wsum(A; wts=w)

    A = randn(4,3,5,2)
    @test sum(A) ≈ @inferred RBMs.wsum(A)
    @test sum(A; dims=(2,4)) ≈ @inferred RBMs.wsum(A; dims=(2,4))

    wts = rand(size(A)...)
    @test sum(A .* wts) ≈ @inferred RBMs.wsum(A; wts)
    @test sum(A .* wts) ≈ @inferred RBMs.wsum(A; wts, dims=:)

    wts = rand(3,2)
    @test sum(reshape(wts,1,3,1,2) .* A; dims=(2,4)) ≈ @inferred RBMs.wsum(A; dims=(2,4), wts)

    wts = rand(2)
    @test sum(reshape(wts,1,1,1,2) .* A; dims=4) ≈ @inferred RBMs.wsum(A; dims=4, wts)
end

@testset "reshape_maybe" begin
    @test RBMs.reshape_maybe(1, ()) == 1
    @test_throws Exception RBMs.reshape_maybe(1, 1)
    @test_throws Exception RBMs.reshape_maybe(1, (1,))

    @test RBMs.reshape_maybe(fill(1), ()) == 1
    @test RBMs.reshape_maybe(fill(1), (1,)) == [1]
    @test RBMs.reshape_maybe(fill(1), (1,1)) == hcat([1])
    @test_throws Exception RBMs.reshape_maybe(fill(1), (1,2))

    @test RBMs.reshape_maybe([1], ()) == 1
    @test RBMs.reshape_maybe([1], (1,)) == [1]
    @test RBMs.reshape_maybe([1], (1,1)) == hcat([1])
    @test_throws Exception RBMs.reshape_maybe([1], (1,2))

    A = randn(2,2)
    @test RBMs.reshape_maybe(A, 4) == reshape(A, 4)
end

@testset "repeat_size" begin
    ns = ((), (0,), (1,), (1,2), (2,1), (2,3))
    rs = ((), (0,), (1,), (1,2), (2,1), (2,3))
    for n in ns, r in rs
        @test repeat_size(n, r...) == size(repeat(trues(n...), r...))
    end
end

@testset "sizedims" begin
    A = randn(7,4,3,2)
    @test @inferred(sizedims(A)) == ()
    @test @inferred(sizedims(A, 2)) == (4,)
    @test @inferred(sizedims(A, :)) == (7,4,3,2)
    @test @inferred(sizedims(A, (2,))) == (4,)
    @test @inferred(sizedims(A, (2,3))) == (4,3)
    @test @inferred(sizedims(A, 2, 3)) == (4,3)
end

@testset "moving_average" begin
    for n in 1:10, m in 1:10
        @test length(@inferred moving_average(randn(n), m)) == n
    end
    x = randn(11)
    @test @inferred(moving_average(x, length(x)))[6] ≈ mean(x)
    @test @inferred(moving_average(x, 3))[5] ≈ mean([x[4], x[5], x[6]])
    @test @inferred(moving_average(x, 3))[1] ≈ mean([x[1], x[2]])
end


text:
facts("Outgoing events") do
    context("OutgoingMessageEvent") do
        s = DandelionSlack.serialize(OutgoingMessageEvent(utf8("Hello"), ChannelId("C0")))
        @fact s["type"] --> utf8("message")
        @fact s["text"] --> utf8("Hello")
        @fact s["channel"] --> ChannelId("C0")
        # Id is set when sending the message, not when serializing.
        @fact haskey(s, "id") --> false
    end

    context("OutgoingPingEvent") do
        s = DandelionSlack.serialize(DandelionSlack.OutgoingPingEvent())
        @fact s["type"] --> utf8("ping")
        @fact haskey(s, "id") --> false
    end
end


text:
"""
"""
abstract type SparseMatrixAssembler <: Assembler end

"""
"""
function get_matrix_type(a::SparseMatrixAssembler)
  @abstractmethod
end

"""
"""
function get_vector_type(a::SparseMatrixAssembler)
  @abstractmethod
end

"""
"""
function symbolic_loop_matrix!(A,a::SparseMatrixAssembler,matdata)
  @abstractmethod
end

"""
"""
function symbolic_loop_vector!(b,a::SparseMatrixAssembler,vecdata)
  @abstractmethod
end

"""
"""
function symbolic_loop_matrix_and_vector!(A,b,a::SparseMatrixAssembler,data)
  @abstractmethod
end

"""
"""
function numeric_loop_matrix!(A,a::SparseMatrixAssembler,matdata)
  @abstractmethod
end

"""
"""
function numeric_loop_vector!(b,a::SparseMatrixAssembler,vecdata)
  @abstractmethod
end

"""
"""
function numeric_loop_matrix_and_vector!(A,b,a::SparseMatrixAssembler,data)
  @abstractmethod
end

function allocate_vector(a::SparseMatrixAssembler,vecdata)
  v1 = nz_counter(get_vector_type(a),(get_rows(a),))
  symbolic_loop_vector!(v1,a,vecdata)
  v2 = nz_allocation(v1)
  symbolic_loop_vector!(v2,a,vecdata)
  v3 = create_from_nz(v2)
  v3
end

function assemble_vector!(b,a::SparseMatrixAssembler,vecdata)
  fill_entries!(b,zero(eltype(b)))
  assemble_vector_add!(b,a,vecdata)
end

function assemble_vector_add!(b,a::SparseMatrixAssembler,vecdata)
  numeric_loop_vector!(b,a,vecdata)
  b
end

function assemble_vector(a::SparseMatrixAssembler,vecdata)
  v1 = nz_counter(get_vector_type(a),(get_rows(a),))
  symbolic_loop_vector!(v1,a,vecdata)
  v2 = nz_allocation(v1)
  numeric_loop_vector!(v2,a,vecdata)
  v3 = create_from_nz(v2)
  v3
end

function allocate_matrix(a::SparseMatrixAssembler,matdata)
  m1 = nz_counter(get_matrix_type(a),(get_rows(a),get_cols(a)))
  symbolic_loop_matrix!(m1,a,matdata)
  m2 = nz_allocation(m1)
  symbolic_loop_matrix!(m2,a,matdata)
  m3 = create_from_nz(m2)
  m3
end

function assemble_matrix!(mat,a::SparseMatrixAssembler,matdata)
  fill_entries!(mat,zero(eltype(mat)))
  assemble_matrix_add!(mat,a,matdata)
end

function assemble_matrix_add!(mat,a::SparseMatrixAssembler,matdata)
  numeric_loop_matrix!(mat,a,matdata)
  mat
end

function assemble_matrix(a::SparseMatrixAssembler,matdata)
  m1 = nz_counter(get_matrix_type(a),(get_rows(a),get_cols(a)))
  symbolic_loop_matrix!(m1,a,matdata)
  m2 = nz_allocation(m1)
  numeric_loop_matrix!(m2,a,matdata)
  m3 = create_from_nz(m2)
  m3
end

function allocate_matrix_and_vector(a::SparseMatrixAssembler,data)
  m1 = nz_counter(get_matrix_type(a),(get_rows(a),get_cols(a)))
  v1 = nz_counter(get_vector_type(a),(get_rows(a),))
  symbolic_loop_matrix_and_vector!(m1,v1,a,data)
  m2 = nz_allocation(m1)
  v2 = nz_allocation(v1)
  symbolic_loop_matrix_and_vector!(m2,v2,a,data)
  m3 = create_from_nz(m2)
  v3 = create_from_nz(v2)
  m3,v3
end

function assemble_matrix_and_vector!(A,b,a::SparseMatrixAssembler, data)
  fill_entries!(A,zero(eltype(A)))
  fill_entries!(b,zero(eltype(b)))
  assemble_matrix_and_vector_add!(A,b,a,data)
  A, b
end

function assemble_matrix_and_vector_add!(A,b,a::SparseMatrixAssembler,data)
  numeric_loop_matrix_and_vector!(A,b,a,data)
  A, b
end

function assemble_matrix_and_vector(a::SparseMatrixAssembler, data)
  m1 = nz_counter(get_matrix_type(a),(get_rows(a),get_cols(a)))
  v1 = nz_counter(get_vector_type(a),(get_rows(a),))
  symbolic_loop_matrix_and_vector!(m1,v1,a,data)
  m2 = nz_allocation(m1)
  v2 = nz_allocation(v1)
  numeric_loop_matrix_and_vector!(m2,v2,a,data)
  m3 = create_from_nz(m2)
  v3 = create_from_nz(v2)
  m3,v3
end

function test_sparse_matrix_assembler(a::SparseMatrixAssembler,matdata,vecdata,data)
  test_assembler(a,matdata,vecdata,data)
  _ = get_matrix_type(a)
  _ = get_vector_type(a)
end

struct GenericSparseMatrixAssembler{M,V} <: SparseMatrixAssembler
  matrix_type::Type{M}
  vector_type::Type{V}
  rows::AbstractUnitRange
  cols::AbstractUnitRange
  strategy::AssemblyStrategy

  function GenericSparseMatrixAssembler(
    matrix_type::Type{M},
    vector_type::Type{V},
    rows::AbstractUnitRange,
    cols::AbstractUnitRange,
    strategy::AssemblyStrategy) where {M,V}
    new{M,V}(matrix_type,vector_type,rows,cols,strategy)
  end
end

function SparseMatrixAssembler(
  mat::Type,vec::Type,trial::FESpace,test::FESpace,strategy::AssemblyStrategy)
  rows = get_free_dof_ids(test)
  cols = get_free_dof_ids(trial)
  GenericSparseMatrixAssembler(mat,vec,rows,cols,strategy)
end

function SparseMatrixAssembler(mat::Type,vec::Type,trial::FESpace,test::FESpace)
  strategy = DefaultAssemblyStrategy()
  rows = get_free_dof_ids(test)
  cols = get_free_dof_ids(trial)
  GenericSparseMatrixAssembler(mat,vec,rows,cols,strategy)
end

function SparseMatrixAssembler(mat::Type,trial::FESpace,test::FESpace)
  strategy = DefaultAssemblyStrategy()
  rows = get_free_dof_ids(test)
  cols = get_free_dof_ids(trial)
  GenericSparseMatrixAssembler(mat,Vector{eltype(mat)},rows,cols,strategy)
end

"""
"""
function SparseMatrixAssembler(trial::FESpace,test::FESpace)
  T = get_dof_value_type(trial)
  matrix_type = SparseMatrixCSC{T,Int}
  vector_type = Vector{T}
  strategy = DefaultAssemblyStrategy()
  rows = get_free_dof_ids(test)
  cols = get_free_dof_ids(trial)
  GenericSparseMatrixAssembler(matrix_type,vector_type,rows,cols,strategy)
end

get_rows(a::GenericSparseMatrixAssembler) = a.rows

get_cols(a::GenericSparseMatrixAssembler) = a.cols

get_matrix_type(a::GenericSparseMatrixAssembler) = a.matrix_type

get_vector_type(a::GenericSparseMatrixAssembler) = a.vector_type

get_assembly_strategy(a::GenericSparseMatrixAssembler) = a.strategy

function symbolic_loop_matrix!(A,a::GenericSparseMatrixAssembler,matdata)
  if LoopStyle(A) == DoNotLoop()
    return A
  end
  for (cellmat,_cellidsrows,_cellidscols) in zip(matdata...)
    cellidsrows = map_cell_rows(a.strategy,_cellidsrows)
    cellidscols = map_cell_cols(a.strategy,_cellidscols)
    rows_cache = array_cache(cellidsrows)
    cols_cache = array_cache(cellidscols)
    @assert length(cellidscols) == length(cellidsrows)
    if length(cellidscols) > 0
      mat = first(cellmat)
      Is = _get_block_layout(mat)
      _symbolic_loop_matrix!(A,rows_cache,cols_cache,cellidsrows,cellidscols,Is)
    end
  end
  A
end

function _get_block_layout(a::Tuple)
  _get_block_layout(a[1])
end

function _get_block_layout(a::AbstractMatrix)
  nothing
end

function _get_block_layout(a::BlockArrayCoo)
  [(I,_get_block_layout(a[I])) for I in eachblockid(a) if is_nonzero_block(a,I) ]
end

@noinline function _symbolic_loop_matrix!(A,rows_cache,cols_cache,cell_rows,cell_cols,Is)
  for cell in 1:length(cell_cols)
    rows = getindex!(rows_cache,cell_rows,cell)
    cols = getindex!(cols_cache,cell_cols,cell)
    _symbolic_matrix_entries_at_cell!(A,rows,cols,Is)
  end
end

@inline function _symbolic_matrix_entries_at_cell!(A,rows,cols,Is)
  add_entries!(A,nothing,rows,cols)
end

@inline function _symbolic_matrix_entries_at_cell!(
  A,rows::BlockArrayCoo,cols::BlockArrayCoo,Is::AbstractArray)
  for (I,Is_next) in Is
    i,j = I.n
    _symbolic_matrix_entries_at_cell!(A,rows[Block(i)],cols[Block(j)],Is_next)
  end
end

function numeric_loop_matrix!(A,a::GenericSparseMatrixAssembler,matdata)
  for (cellmat,_cellidsrows,_cellidscols) in zip(matdata...)
    cellidsrows = map_cell_rows(a.strategy,_cellidsrows)
    cellidscols = map_cell_cols(a.strategy,_cellidscols)
    rows_cache = array_cache(cellidsrows)
    cols_cache = array_cache(cellidscols)
    vals_cache = array_cache(cellmat)
    @assert length(cellidscols) == length(cellidsrows)
    @assert length(cellmat) == length(cellidsrows)
    _numeric_loop_matrix!(
      A,vals_cache,rows_cache,cols_cache,cellmat,cellidsrows,cellidscols)
  end
  A
end

@noinline function _numeric_loop_matrix!(
  mat,vals_cache,rows_cache,cols_cache,cell_vals,cell_rows,cell_cols)
  for cell in 1:length(cell_cols)
    rows = getindex!(rows_cache,cell_rows,cell)
    cols = getindex!(cols_cache,cell_cols,cell)
    vals = getindex!(vals_cache,cell_vals,cell)
    _numeric_matrix_entries_at_cell!(mat,rows,cols,vals)
  end
end

@inline function _numeric_matrix_entries_at_cell!(mat,rows,cols,vals)
  add_entries!(mat,vals,rows,cols)
end

@inline function _numeric_matrix_entries_at_cell!(
  mat,rows::BlockArrayCoo,cols::BlockArrayCoo,vals::BlockArrayCoo)
  for I in eachblockid(vals)
    if is_nonzero_block(vals,I)
      i,j = I.n
      _numeric_matrix_entries_at_cell!(mat,rows[Block(i)],cols[Block(j)],vals[I])
    end
  end
end

function symbolic_loop_vector!(b,a::GenericSparseMatrixAssembler,vecdata)
  @notimplementedif LoopStyle(b) == Loop()
  b
end

function numeric_loop_vector!(b,a::GenericSparseMatrixAssembler,vecdata)
  for (cellvec, _cellids) in zip(vecdata...)
    cellids = map_cell_rows(a.strategy,_cellids)
    rows_cache = array_cache(cellids)
    vals_cache = array_cache(cellvec)
    _numeric_loop_vector!(b,vals_cache,rows_cache,cellvec,cellids)
  end
  b
end

@noinline function _numeric_loop_vector!(vec,vals_cache,rows_cache,cell_vals,cell_rows)
  @assert length(cell_vals) == length(cell_rows)
  for cell in 1:length(cell_rows)
    rows = getindex!(rows_cache,cell_rows,cell)
    vals = getindex!(vals_cache,cell_vals,cell)
    _numeric_vector_entries_at_cell!(vec,rows,vals)
  end
end

@inline function _numeric_vector_entries_at_cell!(vec,rows,vals)
  add_entries!(vec,vals,rows)
end

@inline function _numeric_vector_entries_at_cell!(vec,rows::BlockArrayCoo,vals::BlockArrayCoo)
  for I in eachblockid(vals)
    if is_nonzero_block(vals,I)
      _numeric_vector_entries_at_cell!(vec,rows[I],vals[I])
    end
  end
end

function symbolic_loop_matrix_and_vector!(A,b,a::GenericSparseMatrixAssembler,data)
  matvecdata, matdata, vecdata = data
  symbolic_loop_matrix!(A,a,matvecdata)
  symbolic_loop_matrix!(A,a,matdata)
  symbolic_loop_vector!(b,a,vecdata)
  A, b
end

function numeric_loop_matrix_and_vector!(A,b,a::GenericSparseMatrixAssembler,data)
  matvecdata, matdata, vecdata = data
  for (cellmatvec,_cellidsrows,_cellidscols) in zip(matvecdata...)
    cellidsrows = map_cell_rows(a.strategy,_cellidsrows)
    cellidscols = map_cell_cols(a.strategy,_cellidscols)
    rows_cache = array_cache(cellidsrows)
    cols_cache = array_cache(cellidscols)
    vals_cache = array_cache(cellmatvec)
    @assert length(cellidscols) == length(cellidsrows)
    @assert length(cellmatvec) == length(cellidsrows)
    _numeric_loop_mavec!(A,b,vals_cache,rows_cache,cols_cache,cellmatvec,cellidsrows,cellidscols)
  end
  numeric_loop_matrix!(A,a,matdata)
  numeric_loop_vector!(b,a,vecdata)
  A, b
end

@noinline function _numeric_loop_mavec!(
  A,b,vals_cache,rows_cache,cols_cache,cell_vals,cell_rows,cell_cols)
  for cell in 1:length(cell_cols)
    rows = getindex!(rows_cache,cell_rows,cell)
    cols = getindex!(cols_cache,cell_cols,cell)
    vals = getindex!(vals_cache,cell_vals,cell)
    matvals, vecvals = vals
    _numeric_matrix_entries_at_cell!(A,rows,cols,matvals)
    _numeric_vector_entries_at_cell!(b,rows,vecvals)
  end
end


text:
module TestDepthFirst

using TapirSchedulers
using TapirSchedulers:
    DEPTH_FIRST_SCHEDULER, Priority, rollback_priority, get_current_priority_range
using Test

current_priority() = first(get_current_priority_range())

function record_priority(n, sch = DEPTH_FIRST_SCHEDULER[])
    priorities = fill!(Vector{Priority}(undef, n), typemax(Priority))
    workers = zeros(Int, n)
    function recur(indices)
        if length(indices) == 0
        elseif length(indices) == 1
            @inbounds priorities[indices[1]] = current_priority()
            @inbounds workers[indices[1]] = Threads.threadid()
        else
            f = first(indices)
            l = last(indices)
            m = (l - f + 1) ÷ 2 + f - 1
            work = TapirSchedulers.spawn!(sch) do
                recur(m+1:l)
            end
            recur(f:m)
            fetch(work)::Nothing
        end
        return nothing
    end
    rollback_priority() do
        recur(eachindex(priorities))
    end
    return (; priorities, workers)
end

function test_priority()
    e = 5
    @assert e < sizeof(Priority) * 8
    (; priorities) = record_priority(2^e)
    @test all(<(typemax(Priority)), priorities)
    @test issorted(priorities)
end

end  # module


text:
function LinearInterpolator()
    return LinearInterpolator(())
end

function interpolate(obj::LinearInterpolator, arg0::Vector{RealFieldElement}, arg1::Vector{RealFieldElement})
    return jcall(obj, "interpolate", FieldPolynomialSplineFunction, (Vector{RealFieldElement}, Vector{RealFieldElement}), arg0, arg1)
end

function interpolate(obj::LinearInterpolator, arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return jcall(obj, "interpolate", PolynomialSplineFunction, (Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end


text:
using Dates: Millisecond
using TimeToLive: TTL, Node
using Test: @test, @testset, @test_throws

const p = Millisecond(250)

@testset "TimeToLive.jl" begin
    @testset "Constructors" begin
        t = TTL(p)
        @test t.ttl == p
        @test t.dict isa Dict{Any, Node{Any}}
        t = TTL{Int, String}(p)
        @test t.dict isa Dict{Int, Node{String}}
    end

    @testset "Basic expiry" begin
        t = TTL(p)
        t[0] = "!"
        @test get(t, 0, nothing) == "!"
        sleep(p)
        @test isempty(t)
    end

    @testset "Iteration" begin
        t = TTL(p)
        t[0] = "?"
        sleep(p)
        foreach(i -> t[i] = "!", 1:5)
        count = 0
        for pair in t
            @test pair isa Pair{Int, String}
            @test pair.second == "!"
            count += 1
        end
        @test count == 5
    end

    @testset "Refresh on access" begin
        t = TTL(p; refresh_on_access=true)
        t[0] = "!"
        sleep(3p/5)
        t[0]
        sleep(p/2)
        @test get(t, 0, nothing) == "!"
    end

    @testset "Disabled refresh on access" begin
        t = TTL(p; refresh_on_access=false)
        t[0] = "!"
        sleep(p/2)
        t[0]
        sleep(p)
        @test get(t, 0, nothing) === nothing
    end

    @testset "Base methods" begin
        @testset "delete!" begin
            t = TTL(p)
            t[0] = 1
            @test delete!(t, 0) === t
            @test isempty(t)
            @test delete!(t, 1) === t
        end

        @testset "empty!" begin
            t = TTL(p)
            t[0] = t[1] = 2
            @test empty!(t) === t
            @test isempty(t)
        end

        @testset "get" begin
            t = TTL(p)
            @test get(t, 0, nothing) === nothing
            t[0] = 1
            @test get(t, 0, nothing) == 1
            sleep(p)
            @test get(t, 0, nothing) === nothing

            @test get(() -> nothing, t, 0) === nothing
            t[0] = 1
            @test get(() -> nothing, t, 0) == 1
            sleep(p)
            @test get(() -> nothing, t, 0) === nothing
        end

        @testset "get!" begin
            t = TTL(p)
            @test get!(t, 0, 1) == 1
            @test t[0] == 1
            t[1] = 2
            @test get!(t, 1, 3) == 2
            @test t[1] == 2
            sleep(p)
            @test get!(t, 1, 3) == 3
            @test t[1] == 3
        end

        @testset "getkey" begin
            t = TTL(p)
            t[1] = 2
            @test getkey(t, 0, nothing) === nothing
            @test getkey(t, 1, nothing) == 1
            sleep(2p)
            @test getkey(t, 1, nothing) === nothing
        end

        @testset "length" begin
            t = TTL(p)
            @test length(t) == 0
            t[1] = t[2] = t[3] = 4
            @test length(t) == 3
            sleep(p/2)
            t[5] = 6
            sleep(p/2)
            @test length(t) == 1
        end

        @testset "pop!" begin
            t = TTL(p)
            t[1] = 2
            t[3] = 4
            t[5] = 6
            t[7] = 8
            @test pop!(t) == (7 => 8)
            @test length(t) == 3
            @test pop!(t, 1) == 2
            @test length(t) == 2
            sleep(p)
            @test_throws KeyError pop!(t, 5)
            @test_throws ArgumentError pop!(t)
        end

        @testset "push!" begin
            t = TTL(p)
            @test push!(t, 1 => 2) === t
            @test t[1] == 2
        end

        @testset "sizehint!" begin
            t = TTL(p)
            @test sizehint!(t, 10) === t
        end

        @testset "Others" begin
            t = TTL(p)
            t[1] = 2
            t[2] = 3
            t[3] = 4
            t[4] = 5
            t[5] = 6

            @test all(v -> v isa Int, values(t))
            @test length(filter!(p -> p.second > 2, t)) == 4
            @test !haskey(t, 1)
        end
    end
end


text:
module Sim

using DataFrames
using Distances
using CSV
using JuMP
using Gurobi
using Parameters
using DispatchOps.Nexus

import DispatchOps.Sim
import Base:
    show, isempty

include("interface.jl")

include("Events/plan.jl")
include("Events/transport.jl")
include("Events/fulfill.jl")

include("Utils/demand.jl")
include("Utils/flow.jl")
include("Utils/analysis.jl")

export
    # interface
    Libraries, Params, States, Accumulators, Simulation,

    # Events/plan
    buildGraph, hard_holdover_model, soft_holdover_model, optimizeModel, plan!,

    # Events/fulfill!
    fulfill!,

    # Events/transport!
    transport!,

    # Utils/flow
    reset!, initiate!, schedule!, run!,

    # Utils/demand
    noisify_fixed, noisify_varied,

    # Utils/analysis
    lost_sales, total_cost, fixed_cost, variable_cost

end


text:
# Autogenerated wrapper script for libiio_jll for x86_64-apple-darwin
export libiio

using libusb_jll
using XML2_jll
using libaio_jll
JLLWrappers.@generate_wrapper_header("libiio")
JLLWrappers.@declare_library_product(libiio, "@rpath/iio")
function __init__()
    JLLWrappers.@generate_init_header(libusb_jll, XML2_jll, libaio_jll)
    JLLWrappers.@init_library_product(
        libiio,
        "lib/iio.framework/Versions/0.23/iio",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using SatResample
using SatResample: resample
using Test

@testset "resample" begin

    @test isequal(resample([1,2,3],[1,2,3],[1,1,1],1:3,1:3,1,1),
        [1.0 NaN NaN
        NaN 1.0 NaN
        NaN NaN 1.0])

    @test isequal(resample([1,2,3],[1,2,3],[1,1,1],1:3,1:3,1+eps(),1),
    [1.0 1.0 NaN
     1.0 1.0 1.0
     NaN 1.0 1.0])

    @test isequal(resample([1,2,3],[1,2,3],[1,1,1],1:3,1:3,sqrt(2)+eps(),1), ones(3,3))

    @test isequal(resample([1,2,3],[1,2,3],[1,2,3],1:3,1:3,sqrt(2)+eps(),1),
      [1.0  2.0  2.0
       2.0  2.0  3.0
       2.0  3.0  3.0])

    @test isequal(resample([1,2,3],[1,2,3],[1,2,3],1:3,1:3,sqrt(2)+eps(),2),
    [1.5  1.5  2.0
     1.5  2.5  2.5
     2.0  2.5  2.5])

    @test isequal(resample([1,2,3],[1,2,3],[[1,1,1],[1,2,3]],1:3,1:3,sqrt(2)+eps(),2),
     [[1.0 1.0 1.0; 1.0 1.0 1.0; 1.0 1.0 1.0],
     [1.5 1.5 2.0; 1.5 2.5 2.5; 2.0 2.5 2.5]])

end

@testset "grids" begin
    #check NSIDC North and South grids for correct grid sizes
    #see https://nsidc.org/data/polar-stereo/ps_grids.html
    let genxy=SatResample.genxy
        let NSIDCN = SatResample.NSIDCN
            @test length.(genxy(NSIDCN(50000)))==(152,224)
            @test length.(genxy(NSIDCN(25000)))==(304,448)
            @test length.(genxy(NSIDCN(12500)))==(608,896)
            @test length.(genxy(NSIDCN( 6250)))==(1216,1792)
            @test length.(genxy(NSIDCN( 3125)))==(2432,3584)
        end

        let NSIDCS = SatResample.NSIDCS
            @test length.(genxy(NSIDCS(50000)))==(158,166)
            @test length.(genxy(NSIDCS(25000)))==(316,332)
            @test length.(genxy(NSIDCS(12500)))==(632,664)
            @test length.(genxy(NSIDCS( 6250)))==(1264,1328)
            @test length.(genxy(NSIDCS( 3125)))==(2528,2656)
        end
    end
end


text:
using TailRiskScenGen
using Cones
using LinearAlgebra: I
using Distributions

d = 30
β = 0.95
num_scen = 5000
A = rand(d,d)/sqrt(d)
dist = MvNormal(0.1*rand(d)-0.05, A'A)

K = PolyhedralCone(Array{Float64}(I, d, d))
Ω_gen = EllipticalRiskRegion(dist, K, β)
Ω_mono = MonotonicEllipticalRiskRegion(dist, K, β)
Ω_mono.max_frontier_points = d*5

rights = 0
wrongs = 0
for i in 1:200
    x = rand(dist)
    mon = x ∈ Ω_mono
    gen = x ∈ Ω_gen
    if mon == gen
        rights+=1
    else
        wrongs+=1
    end

    if i % 10 == 1
        println("i = ", i, ", num risk frontier points: ", length(Ω_mono.risk_frontier),
              ", num nonrisk frontier points: ", length(Ω_mono.nonrisk_frontier))
    end
end

println("Correct decisions: ", rights)
println("Wrong decisions: ", wrongs)

x = rand(dist)

x ∈ Ω_gen
x ∈ Ω_mono
y = rand(dist)
y ∈ Ω_gen
y ∈ Ω_mono

Ω_mono.risk_frontier
Ω_mono.nonrisk_frontier
Ω=Ω_mono

print("EllipticalRiskRegion...")
tic()
TailRiskScenGen._aggregation_sampling(dist, Ω_gen, num_scen)
toc()

print("MonotonicEllipticalRiskRegion...")
tic()
TailRiskScenGen._aggregation_sampling(dist, Ω_mono, num_scen)
toc()

println("Num non-risk frontier points: ", length(Ω_mono.risk_frontier))
println("Num risk frontier points: ", length(Ω_mono.nonrisk_frontier))


text:
# # Hydrostatic Equilibrium test for Richard's Equation
# Note that freezing and thawing are turned off in this example. That
# means that `θ_i`, initialized to zero by default, will remain zero.

# # Preliminary setup

# - load external packages

using MPI
using OrderedCollections
using StaticArrays
using Statistics
using Plots

# - load CLIMAParameters and ClimateMachine modules

using CLIMAParameters
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

using ClimateMachine
using ClimateMachine.Land
using ClimateMachine.Land.SoilWaterParameterizations
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Grids
using ClimateMachine.DGMethods
using ClimateMachine.DGMethods.NumericalFluxes
using ClimateMachine.DGMethods: BalanceLaw, LocalGeometry
using ClimateMachine.MPIStateArrays
using ClimateMachine.GenericCallbacks
using ClimateMachine.ODESolvers
using ClimateMachine.VariableTemplates
using ClimateMachine.SingleStackUtils
using ClimateMachine.BalanceLaws:
    BalanceLaw, Prognostic, Auxiliary, Gradient, GradientFlux, vars_state

# - define the float type desired (`Float64` or `Float32`)
const FT = Float64;

# - initialize ClimateMachine for CPU
ClimateMachine.init(; disable_gpu = true);

# # Set up the soil model

# We want to solve Richard's equation alone, i.e. do not solve
# the heat equation. The default is a constant temperature, but this
# only affects Richard's equation if one chooses a temperature dependent
# viscosity of water (see below).

soil_heat_model = PrescribedTemperatureModel()

# Define the porosity, Ksat, and specific storage values for the soil. Note
# that all values must be givne in mks.
soil_param_functions = SoilParamFunctions{FT}(
    porosity = 0.495,
    Ksat = 0.0443 / (3600 * 100),
    S_s = 1e-3,
)

# Define the boundary conditions. The user can specify two conditions,
# either at the top or at the bottom, and they can either be Dirichlet
# (on `ϑ_l`) or Neumann (on `-K∇h`). Note that fluxes are supplied as
# scalars, inside the code they are multiplied by ẑ. The two conditions
# not supplied must be set to `nothing`.

surface_flux = (aux, t) -> eltype(aux)(0.0)
bottom_flux = (aux, t) -> eltype(aux)(0.0)
surface_state = nothing
bottom_state = nothing

# Define the initial state function.
ϑ_l0 = (aux) -> eltype(aux)(0.494)

# Create the SoilWaterModel. The defaults are a temperature independent
# viscosity, and no impedance factor due to ice. We choose to make the
# hydraulic conductivity a function of the moisture content `ϑ_l`,
# and employ the vanGenuchten hydraulic model with `n` = 2.0. The van
# Genuchten parameter `m` is calculated form `n`, and we use the default
# value for `α`.
soil_water_model = SoilWaterModel(
    FT;
    moisture_factor = MoistureDependent{FT}(),
    hydraulics = vanGenuchten{FT}(n = 2.0),
    initialϑ_l = ϑ_l0,
    dirichlet_bc = Dirichlet(
        surface_state = surface_state,
        bottom_state = bottom_state,
    ),
    neumann_bc = Neumann(
        surface_flux = surface_flux,
        bottom_flux = bottom_flux,
    ),
)

# Create the soil model - the coupled soil water and soil heat models.
m_soil = SoilModel(soil_param_functions, soil_water_model, soil_heat_model)

# We are ignoring sources and sinks here.
sources = ()

# Define the function that initializes the prognostic variables. This
# in turn calls the functions supplied to soil_water_model. The default
# initialization for ice is zero volumetric ice fraction. This should not
# be altered unless freeze thaw is added.
function init_soil_water!(land, state, aux, coordinates, time)
    FT = eltype(state)
    state.soil.water.ϑ_l = FT(land.soil.water.initialϑ_l(aux))
    state.soil.water.θ_i = FT(land.soil.water.initialθ_i(aux))
end

# Create the land model - in this tutorial, it only includes the soil.
m = LandModel(
    param_set,
    m_soil;
    source = sources,
    init_state_prognostic = init_soil_water!,
)

# # Specify the numerical configuration and output data.

# Specify the polynomial order and vertical resolution.
N_poly = 5;
nelem_vert = 20;

# Specify the domain boundaries.
zmax = FT(0);
zmin = FT(-10)

# Create the driver configuration.
driver_config = ClimateMachine.SingleStackConfiguration(
    "LandModel",
    N_poly,
    nelem_vert,
    zmax,
    param_set,
    m;
    zmin = zmin,
    numerical_flux_first_order = CentralNumericalFluxFirstOrder(),
);

# Choose the initial and final times, as well as a timestep.
t0 = FT(0)
timeend = FT(60 * 60 * 24 * 4)
dt = FT(5)

# Create the solver configuration.
solver_config =
    ClimateMachine.SolverConfiguration(t0, timeend, driver_config, ode_dt = dt);

# Determine how often you want output.
const n_outputs = 5;

const every_x_simulation_time = ceil(Int, timeend / n_outputs);

# Create a place to store this output.
all_data = Dict([k => Dict() for k in 0:n_outputs]...)

# Obtain the initial state, and create the callback
# function to obtain the output from the simulation.
mygrid = solver_config.dg.grid;
Q = solver_config.Q;
aux = solver_config.dg.state_auxiliary;

t = ODESolvers.gettime(solver_config.solver)
state_vars = SingleStackUtils.get_vars_from_nodal_stack(
    mygrid,
    Q,
    vars_state(m, Prognostic(), FT),
)
grads = SingleStackUtils.get_vars_from_nodal_stack(
    mygrid,
    solver_config.dg.state_gradient_flux,
    vars_state(m, GradientFlux(), FT),
)
aux_vars = SingleStackUtils.get_vars_from_nodal_stack(
    mygrid,
    aux,
    vars_state(m, Auxiliary(), FT),
)
all_vars = OrderedDict(state_vars..., aux_vars..., grads...);
all_vars["t"] = [t]
all_data[0] = all_vars

iostep = [1];
callback = GenericCallbacks.EveryXSimulationTime(
    every_x_simulation_time,
) do (init = false)
    t = ODESolvers.gettime(solver_config.solver)
    grads = SingleStackUtils.get_vars_from_nodal_stack(
        mygrid,
        solver_config.dg.state_gradient_flux,
        vars_state(m, GradientFlux(), FT),
    )
    state_vars = SingleStackUtils.get_vars_from_nodal_stack(
        mygrid,
        Q,
        vars_state(m, Prognostic(), FT),
    )
    aux_vars = SingleStackUtils.get_vars_from_nodal_stack(
        mygrid,
        aux,
        vars_state(m, Auxiliary(), FT),
    )
    all_vars = OrderedDict(state_vars..., aux_vars..., grads...)
    all_vars["t"] = [t]
    all_data[iostep[1]] = all_vars

    iostep[1] += 1
    nothing
end;

# # Run the integration
ClimateMachine.invoke!(solver_config; user_callbacks = (callback,));

# # Get the final state back.
t = ODESolvers.gettime(solver_config.solver)
state_vars = SingleStackUtils.get_vars_from_nodal_stack(
    mygrid,
    Q,
    vars_state(m, Prognostic(), FT),
)
grads = SingleStackUtils.get_vars_from_nodal_stack(
    mygrid,
    solver_config.dg.state_gradient_flux,
    vars_state(m, GradientFlux(), FT),
)
aux_vars = SingleStackUtils.get_vars_from_nodal_stack(
    mygrid,
    aux,
    vars_state(m, Auxiliary(), FT),
)
all_vars = OrderedDict(state_vars..., aux_vars..., grads...);
all_vars["t"] = [t]
all_data[n_outputs] = all_vars

t = [all_data[k]["t"] for k in 0:n_outputs]

# # Create some plots.
slope = -1e-3

# The initial state.
plot(
    all_data[0]["soil.water.ϑ_l"],
    all_data[0]["z"],
    label = string("t = ", string(t[1][1])),
    xlim = [0.47, 0.501],
    ylabel = "z",
    xlabel = "ϑ_l",
    legend = :bottomleft,
    title = "Equilibrium test",
)

# A middle state.
plot!(
    all_data[3]["soil.water.ϑ_l"],
    all_data[3]["z"],
    label = string("t = ", string(t[4][1])),
)
# The final state
plot!(
    all_data[5]["soil.water.ϑ_l"],
    all_data[5]["z"],
    label = string("t = ", string(t[6][1])),
)
# The expected slope in hydrostatic equilibrium.
plot!(
    (all_data[0]["z"] .+ 10.0) .* slope .+ all_data[5]["soil.water.ϑ_l"][1],
    all_data[0]["z"],
    label = "expected",
)

# The porosity of the soil.
plot!(
    1e-3 .+ all_data[0]["soil.water.ϑ_l"],
    all_data[0]["z"],
    label = "porosity",
)
# save the output.
savefig("equilibrium_test_ϑ_l_vG.png")

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
@testset "frature" begin

    using Crayons

    include("shape.jl")
    include("display_obj.jl")
    include("display_stack.jl")

end


text:
"""
    Eratosthenes

This spacecraft simulation environment allows a user to quickly create simulations with different spacecraft, spacecraft components, and simulation options.

```julia
Pkg.clone("https://github.com/tuckermcclure/Eratosthenes.jl.git")
using Eratosthenes
scenario = setup(Scenario(), "my-scenario.json")
simulate(scenario)
```
"""
module Eratosthenes

import Base: convert, +, *, accumulate! # Required in order to extend these things.
import Base.Meta: parse
using Base.Iterators
using Random
using LinearAlgebra
import YAML
using HDF5Logger

# For now, our rotations are in another module included with this package. We
# might move them into their own package if people find them generally useful.
# To include this module from outside of Eratosthenes, use:
#
#    include(Pkg.dir("Eratosthenes") * "/src/EratosthenesRotations.jl")
#    using .EratosthenesRotations
#
include("EratosthenesRotations.jl")
using .EratosthenesRotations

export UtcDateTime
export RandSpec, RandSource, MultiRandSource, draw # Rand stuff
export ModelTiming, DynamicalModel # Top-level modeling types
export Computer, Vehicle # Containers
export LowEarthOrbit, TruthSensor, StarTracker, Gyro, IdealActuator # Models
export IdealActuatorCommand
export NoEffect, Gravity, BodyTorque, BodyForce
export find_effect, accumulate!
export setup, simulate, mc
export Scenario, SimParams

# Used to tell integrator how to add.
abstract type ModelState end

include("utilities.jl")
include("rand.jl")
include("physics.jl")
include("setup.jl")
include("effects.jl")

mutable struct ModelTiming
    dt::Float64
    t_start::Float64
    t_next::Float64
    count::Int64 # Total number of times this model has triggered
    ModelTiming(dt = 0., t_start = 0., t_next = 0., count = 0) = new(dt, t_start, t_next, count)
end

mutable struct DynamicalModel{F0, FI, FE, FC, FD, FU, FS, DC, DX, DU, DY, DV}

    name::String # Actually used as a key in a dictionary, as well as for logging

    # Functions
    #
    # NOTE: These can't be changed without changing the type of the object since
    # we've parameterized the type for speed. If we later make immutable things
    # out of the mutable things, then we'll remove this parameterization here
    # but keep it for the immutable version for speed there.
    #
    startup::F0 # Called to allow one to set up resources or whatever.
    init::FI # "Tell me what you produce."
    effects::FE # "Produce physical effects."
    constraints::FC # "Tell me how well a state satisfies your constraints."
    derivatives::FD # "Tell me how the state is changing."
    update::FU # "Update yourself for this time step."
    shutdown::FS # Called when the sim ends, even when there's an error.

    timing::ModelTiming

    # Data
    constants::DC # This will be set up during setup and will never change afterwards.
    state::DX # Bonus: return a ModelState to enable "struct addition and multiplication"
    inputs::DU # Inputs accepted by this model (e.g., commands); software can write to these.
    outputs::DY # Outputs produced by this model (e.g., measurements); software can consume these.
    implicit::DV # Implicit variables.

    rand::RandSpec # Provide the properties of the random number generator stream you'll need.

end

# Create a convenience constructor with lots of defaults.
DynamicalModel(name;
               startup = nothing,
               init = nothing,
               effects = nothing,
               constraints = nothing,
               derivatives = nothing,
               update = nothing,
               shutdown = nothing,
               timing = ModelTiming(),
               constants = nothing,
               state = nothing,
               inputs = nothing,
               outputs = nothing,
               implicit = Vector{Float64}(),
               rand = RandSpec()) =
    DynamicalModel(name, startup, init, effects, constraints, derivatives, update, shutdown, timing, constants, state, inputs, outputs, implicit, rand)
NullModel() = DynamicalModel("nothing")

# Include the other code in no particular order.
include("actuators.jl")
include("planets.jl")
include("sensors.jl")
include("software.jl")
include("computers.jl")
include("vehicles.jl")

# The things below set some defaults that depend on the above.

# Options for global simulation stuff.
mutable struct SimParams
    dt::Float64 # Time step (s)
    t_end::Float64 # End time of simulation (s)
    seed::Int64 # Seed for default random number generator
    date_0::UtcDateTime # UTC date at the beginning of the sim
    log_file::String # Name of HDF5 file to log to
    SimParams() = new(0.01, 10.0, 1, UtcDateTime(2017, 7, 28, 9, 27, 33.0), "")
    SimParams(a,b,c,d,e) = new(a,b,c,d,e)
end

# A place to keep track of dimensions
struct ScenarioDimensions
    nt::Int64 # Timesteps
    nx::Int64 # Primitive states
    nxd::Int64 # Primitive state dimensions
    nv::Int64 # Implicit variables
    x_to_xd::Vector{Int64}
end

# Top-level options describing what should happen in the sim.
mutable struct Scenario
    sim::SimParams
    environment::DynamicalModel
    vehicles::Array{Vehicle,1}
    dims::ScenarioDimensions
    Scenario() = new(SimParams(), LowEarthOrbit(), [Vehicle()], ScenarioDimensions(0,0,0,0,Vector{Int64}()))
    Scenario(s,e,va,d) = new(s,e,va,d)
end

include("simulate.jl")
include("mc.jl")

end


text:
function sawtooth(u0, s=0., n=1)
	m, = size(u0)
	n = n+1
	u_trj = zeros((m,n))
	u_trj[:,1] = u0
	for i = 2:n
		x = view(u_trj,:,i-1)
		u_trj[:,i] = (2*x + s*sin.(2*pi*x/4)) .% 1
	end
	return permutedims(u_trj,[2,1])
end
function dsawtooth(u, s=0)
	n, = size(u)
	return 2.0 .+ s*cos.(2*pi*u/4)*2*pi/4
end
function perturbation(u,s)
	n = size(u)[1]
	# the perturbation in row i in T_{u_(i+1)} M
	return sin.(2*pi*u/4) .% 1
end

max_stars_repo_name: nishaChandramoorthy/Lipschitz

max_issues_repo_name: nishaChandramoorthy/energies

max_forks_repo_name: nishaChandramoorthy/energies


text:
struct ConcaveTotalChunker{F}
    f::F
end

struct ConcaveTotalSplitter{F}
    f::F
end

function pack_stripe(A::SparseMatrixCSC{Tv, Ti}, method::ConcaveTotalChunker, args...) where {Tv, Ti}
    @inbounds begin
        (m, n) = size(A)

        f = oracle_stripe(RandomHint(), method.f, A, args...)

        ftr = CircularDeque{Tuple{Ti, Ti}}(n + 1)
        spl = zeros(Ti, n + 1)
        cst = fill(typemax(cost_type(f)), n + 1)
        cst[1] = zero(cost_type(f))
        f′(j, j′) = cst[j] + f(j, j′)
        chunk_concave!(cst, spl, f′, 1, n + 1, ftr)

        return unravel_chunks!(spl, n)
    end
end

function partition_stripe(A::SparseMatrixCSC{Tv, Ti}, K, method::ConcaveTotalSplitter, args...) where {Tv, Ti}
    @inbounds begin
        (m, n) = size(A)

        f = oracle_stripe(RandomHint(), method.f, A, args...)

        if K == 1
            return SplitPartition{Ti}(1, [Ti(1), Ti(n + 1)])
        end

        ftr = CircularDeque{Tuple{Ti, Ti}}(n + 1)

        ptr = zeros(Ti, n + 1, K)
        cst = fill(typemax(cost_type(f)), n + 1, K)

        for j′ = 1:n + 1
            cst[j′, 1] = f(1, j′, 1)
            ptr[j′, 1] = 1
        end
        for k = 2:K
            f′(j, j′) = cst[j, k - 1] + f(j, j′, k)
            for j′ = 1:n + 1
                cst[j′, k] = f′(j′, j′)
                ptr[j′, k] = j′
            end
            chunk_concave!((@view cst[:, k]), (@view ptr[:, k]), f′, 1, n + 1, ftr)
        end
        return unravel_splits(K, n, PermutedDimsArray(ptr, (2, 1)))
    end
end

function chunk_concave!(cst, ptr, f::F, j₀, j′₁, ftr) where {F}
    @inbounds begin
        empty!(ftr)
        push!(ftr, (j₀, j₀ + 1))
        for j′ = j₀ + 1:j′₁

            (j, h) = first(ftr)
            c = f(j, j′)
            c′ = f(j′ - 1, j′)
            if c′ <= c
                if c′ <= cst[j′]
                    cst[j′] = c′
                    ptr[j′] = j′ - 1
                end
                empty!(ftr)
                push!(ftr, (j′ - 1, j′ + 1))
            else
                if c <= cst[j′]
                    cst[j′] = c
                    ptr[j′] = j
                end
                while ((j, h) = last(ftr); (f(j′ - 1, h) <= f(j, h)))
                    pop!(ftr)
                end
                (j, h) = last(ftr)

                #=
                h_ref = h + 1
                while h_ref < j′₁ + 1 && f(j′ - 1, h_ref) > f(j, h_ref)
                    h_ref += 1
                end
                h = h_ref
                =#

                h_lo = h + 1
                h_hi = j′₁
                while h_lo <= h_hi
                    h = fld2(h_lo + h_hi)
                    if f(j′ - 1, h) > f(j, h)
                        h_lo = h + 1
                    else
                        h_hi = h - 1
                    end
                end
                h = h_lo

                if h != j′₁ + 1
                    push!(ftr, (j′ - 1, h))
                end

                (j, _) = popfirst!(ftr)
                if isempty(ftr) || ((_, h) = first(ftr); j′ + 1 != h)
                    pushfirst!(ftr, (j, j′ + 1))
                end
            end
        end
    end
end

#=
function chunk_concave!(cst, ptr, f, j₀, j′₁, ftr)
    @inbounds begin
        for a = j₀:j′₁
            for b = a:j′₁
                for c = b + 1:j′₁
                    for d = c:j′₁
                        @assert f(a, c) + f(b, d) <= f(a, d) + f(b, c)
                    end
                end
            end
        end
    end

    @inbounds begin
        for j′ = j₀ + 1:j′₁
            for j = j₀:j′-1
                if f(j, j′) <= cst[j′]
                    cst[j′] = f(j, j′)
                    ptr[j′] = j
                end
            end
        end
    end
end
=#

function partition_stripe(A::SparseMatrixCSC{Tv, Ti}, K, method::ConcaveTotalSplitter{<:ConstrainedCost}, args...) where {Tv, Ti}
    @inbounds begin
        (m, n) = size(A)

        f = oracle_stripe(RandomHint(), method.f.f, A, args...)
        w = oracle_stripe(StepHint(), method.f.w, A, args...)
        w_max = method.f.w_max

        ftr = CircularDeque{Tuple{Ti, Ti}}(n + 1)

        (j′_lo, j′_hi) = column_constraints(A, K, w, w_max)

        if j′_hi[K] < n + 1
            spl = ones(Ti, K + 1)
            spl[end] = n + 1
            #TODO throw(ArgumentError("infeasible"))
            return SplitPartition(K, spl)
        end

        ptr = WindowConstrainedMatrix{Ti, Ti}(zero(Ti), n + 1, K, j′_lo, j′_hi)
        cst = WindowConstrainedMatrix{extend(cost_type(f)), Ti}(infinity(cost_type(f)), n + 1, K, j′_lo, j′_hi, ptr.pos)

        for j′ = j′_lo[1] : j′_hi[1]
            cst[j′, 1] = extend(f(1, j′, 1))
            ptr[j′, 1] = 1
        end

        for k = 2:K
            f′(j, j′) = cst[j, k - 1] + extend(f(j, j′, k))
            for j′ = j′_lo[k]:j′_hi[k]
                cst[j′, k] = f′(j′, j′)
                ptr[j′, k] = j′
            end
            chunk_concave!((@view cst[:, k]), (@view ptr[:, k]), f′, j′_lo[k - 1], j′_hi[k], ftr)
        end

        return unravel_splits(K, n, PermutedDimsArray(ptr, (2, 1)))
    end
end


text:
# In this tutorial, we are exploring the application of Ridge and Lasso
# regression to the Hitters R dataset.
#
# ## Getting started

using MLJ
import RDatasets: dataset
using PrettyPrinting
MLJ.color_off() # hide
import Distributions
const D = Distributions

@load LinearRegressor pkg=MLJLinearModels
@load RidgeRegressor pkg=MLJLinearModels
@load LassoRegressor pkg=MLJLinearModels

# We load the dataset using the `dataset` function, which takes the Package and
# dataset names as arguments.

hitters = dataset("ISLR", "Hitters")
@show size(hitters)
names(hitters) |> pprint

# Let's unpack the dataset with the `unpack` function.
# In this case, the target is `Salary` (`==(:Salary)`) and all other columns are features (`col->true`).

y, X = unpack(hitters, ==(:Salary), col->true);

# The target has missing values which we will just ignore.
# We extract the row indices corresponding to non-missing values of the target.
# Note the use of the element-wise operator `.`.
no_miss = .!ismissing.(y);

# We collect the non missing values of the target in an Array.
# And keep only the corresponding features values.
y = collect(skipmissing(y))
X = X[no_miss, :]

# Let's now split our dataset into a train and test sets.
train, test = partition(eachindex(y), 0.5, shuffle=true, rng=424);

# Let's have a look at the target.

using PyPlot
ioff() # hide

figure(figsize=(8,6))
plot(y, ls="none", marker="o")

xticks(fontsize=12); yticks(fontsize=12)
xlabel("Index", fontsize=14), ylabel("Salary", fontsize=14)

savefig(joinpath(@OUTPUT, "ISL-lab-6-g1.svg")) # hide

# \figalt{Salary}{ISL-lab-6-g1.svg}

# That looks quite skewed, let's have a look at a histogram:

figure(figsize=(8,6))
hist(y, bins=50, density=true)

xticks(fontsize=12); yticks(fontsize=12)
xlabel("Salary", fontsize=14); ylabel("Density", fontsize=14)

edfit = D.fit_mle(D.Exponential, y)
xx = range(minimum(y), 2500, length=100)
yy = pdf.(edfit, xx)
plot(xx, yy, lw=3, label="Exponential distribution fit")

legend(fontsize=12)

savefig(joinpath(@OUTPUT, "ISL-lab-6-g2.svg")) # hide

# \figalt{Distribution of salary}{ISL-lab-6-g2.svg}
#
# ### Data preparation
#
# Most features are currently encoded as integers but we will consider them as continuous.
# To coerce `int` features to `Float`, we nest the `autotype` function in the `coerce` function.
# The `autotype` function returns a dictionary containing scientific types, which is then passed to the `coerce` function.
# For more details on the use of `autotype`, see the [Scientific Types](https://alan-turing-institute.github.io/DataScienceTutorials.jl/data/scitype/index.html#autotype)

Xc = coerce(X, autotype(X, rules=(:discrete_to_continuous,)))
scitype(Xc)

# There're a few features that are categorical which we'll one-hot-encode.

# ## Ridge pipeline
# ### Baseline
#
# Let's first fit a simple pipeline with a standardizer, a one-hot-encoder and a basic linear regression:

model = @pipeline(Standardizer(),
                     OneHotEncoder(),
                     LinearRegressor())

pipe  = machine(model, Xc, y)
fit!(pipe, rows=train)
ŷ = predict(pipe, rows=test)
round(rms(ŷ, y[test])^2, sigdigits=4)

# Let's get a feel for how we're doing

figure(figsize=(8,6))

res = ŷ .- y[test]
stem(res)

xticks(fontsize=12); yticks(fontsize=12)
xlabel("Index", fontsize=14); ylabel("Residual (ŷ - y)", fontsize=14)

ylim([-1300, 1000])

savefig(joinpath(@OUTPUT, "ISL-lab-6-g3.svg")) # hide

# \figalt{Residuals}{ISL-lab-6-g3.svg}

figure(figsize=(8,6))
hist(res, bins=30, density=true, color="green")

xx = range(-1100, 1100, length=100)
ndfit = D.fit_mle(D.Normal, res)
lfit  = D.fit_mle(D.Laplace, res)

plot(xx, pdf.(ndfit, xx), lw=3, color="orange", label="Normal fit")
plot(xx, pdf.(lfit, xx), lw=3, color="magenta", label="Laplace fit")

legend(fontsize=12)

xticks(fontsize=12); yticks(fontsize=12)
xlabel("Residual (ŷ - y)", fontsize=14); ylabel("Density", fontsize=14)
xlim([-1100, 1100])

savefig(joinpath(@OUTPUT, "ISL-lab-6-g4.svg")) # hide

# \figalt{Distribution of residuals}{ISL-lab-6-g4.svg}

# ### Basic Ridge
#
# Let's now swap the linear regressor for a Ridge one without specifying the penalty (`1` by default):
# We modify the supervised model in the pipeline directly.

pipe.model.linear_regressor = RidgeRegressor()
fit!(pipe, rows=train)
ŷ = predict(pipe, rows=test)
round(rms(ŷ, y[test])^2, sigdigits=4)

# Ok that's a bit better but surely we can do better with an appropriate selection of the hyperparameter.

# ### Cross validating

# What penalty should you use? Let's do a simple CV to try to find out:

r  = range(model, :(linear_regressor.lambda), lower=1e-2, upper=100_000, scale=:log10)
tm = TunedModel(model=model, ranges=r, tuning=Grid(resolution=50),
                resampling=CV(nfolds=3, rng=4141), measure=rms)
mtm = machine(tm, Xc, y)
fit!(mtm, rows=train)

best_mdl = fitted_params(mtm).best_model
round(best_mdl.linear_regressor.lambda, sigdigits=4)

# right, and  with that we get:

ŷ = predict(mtm, rows=test)
round(rms(ŷ, y[test])^2, sigdigits=4)

# Let's see:

figure(figsize=(8,6))

res = ŷ .- y[test]
stem(res)

xticks(fontsize=12); yticks(fontsize=12)
xlabel("Index", fontsize=14);
ylabel("Residual (ŷ - y)", fontsize=14)
xlim(1, length(res))

ylim([-1300, 1000])

savefig(joinpath(@OUTPUT, "ISL-lab-6-g5.svg")) # hide

# \figalt{Ridge residuals}{ISL-lab-6-g5.svg}
#
# You can compare that with the residuals obtained earlier.

# ## Lasso pipeline
#
# Let's do the same as above but using a Lasso model and adjusting the range a bit:

mtm.model.model.linear_regressor = LassoRegressor()
mtm.model.range = range(model, :(linear_regressor.lambda), lower=500, upper=100_000, scale=:log10)
fit!(mtm, rows=train)

best_mdl = fitted_params(mtm).best_model
round(best_mdl.linear_regressor.lambda, sigdigits=4)

# Ok and let's see how that does:

ŷ = predict(mtm, rows=test)
round(rms(ŷ, y[test])^2, sigdigits=4)

# Pretty good! and the parameters are reasonably sparse as expected:

coefs, intercept = fitted_params(mtm.fitresult).linear_regressor
@show coefs
@show intercept

# with around 50% sparsity:

coef_vals = [c[2] for c in coefs]
sum(coef_vals .≈ 0) / length(coefs)

# Let's visualise this:

figure(figsize=(8,6))
stem(coef_vals)

## name of the features including one-hot-encoded ones
all_names = [:AtBat, :Hits, :HmRun, :Runs, :RBI, :Walks, :Years,
             :CAtBat, :CHits, :CHmRun, :CRuns, :CRBI, :CWalks,
             :League__A, :League__N, :Div_E, :Div_W,
             :PutOuts, :Assists, :Errors, :NewLeague_A, :NewLeague_N]

idxshow = collect(1:length(coef_vals))[abs.(coef_vals) .> 10]
xticks(idxshow .- 1, all_names[idxshow], rotation=45, fontsize=12)
yticks(fontsize=12)
ylabel("Amplitude", fontsize=14)

savefig(joinpath(@OUTPUT, "ISL-lab-6-g6.svg")) # hide

# \figalt{Lasso coefficients}{ISL-lab-6-g6.svg}

# ## Elastic net pipeline

@load ElasticNetRegressor pkg=MLJLinearModels

mtm.model.model.linear_regressor = ElasticNetRegressor()
mtm.model.range = [range(model, :(linear_regressor.lambda), lower=0.1, upper=100, scale=:log10),
                    range(model, :(linear_regressor.gamma),  lower=500, upper=10_000, scale=:log10)]
mtm.model.tuning = Grid(resolution=10)
fit!(mtm, rows=train)

best_mdl = fitted_params(mtm).best_model
@show round(best_mdl.linear_regressor.lambda, sigdigits=4)
@show round(best_mdl.linear_regressor.gamma, sigdigits=4)

# And it's not too bad in terms of accuracy either

ŷ = predict(mtm, rows=test)
round(rms(ŷ, y[test])^2, sigdigits=4)

# But the simple ridge regression seems to work best here.
PyPlot.close_figs() # hide


text:
using Test

@testset "./test/integration/LocalPolytope.jl" begin

  println("running LocalPolytope.jl integration tests:")
  for test in readdir("./test/integration/LocalPolytope/")
      # run only julia files in test directory
      if occursin(r"^.*\.jl$", test)
          println("./integration/LocalPolytope/$test")
          @time include("./LocalPolytope/$test")
      end
  end

end


text:
import Base.@deprecate

@deprecate gmean geomean
@deprecate hmean harmmean

@deprecate acf autocor
@deprecate ccf crosscor

@deprecate indicators indicatormat


text:
function build_gram_matrix(q::Vector{MOI.VariableIndex},
                           basis::AbstractPolynomialBasis)
    return build_gram_matrix([MOI.SingleVariable(vi) for vi in q], basis)
end
function build_gram_matrix(q::Vector,
                           basis::AbstractPolynomialBasis)
    return GramMatrix(MultivariateMoments.SymMatrix(q, length(basis)),
                      basis)
end

function add_gram_matrix(model::MOI.ModelLike, matrix_cone_type::Type,
                         basis::AbstractPolynomialBasis)
    Q, cQ = MOI.add_constrained_variables(model, matrix_cone(matrix_cone_type, length(basis)))
    q = build_gram_matrix(Q, basis)
    return q, Q, cQ
end
function add_gram_matrix(model::MOI.ModelLike, matrix_cone_type::Type,
                         bases::Vector{<:AbstractPolynomialBasis})
    qQcQ = add_gram_matrix.(model, matrix_cone_type, bases)
    return SparseGramMatrix(getindex.(qQcQ, 1)), getindex.(qQcQ, 2), getindex.(qQcQ, 3)
end

function build_moment_matrix(q::Vector,
                             basis::AbstractPolynomialBasis)
    return MomentMatrix(MultivariateMoments.SymMatrix(q, length(basis)),
                        basis)
end

struct SOSPolynomialSet{DT <: AbstractSemialgebraicSet,
                        MT <: MP.AbstractMonomial,
                        MVT <: AbstractVector{MT},
                        CT <: Certificate.AbstractCertificate} <: MOI.AbstractVectorSet
    domain::DT
    monomials::MVT
    certificate::CT
end
MOI.dimension(set::SOSPolynomialSet) = length(set.monomials)


text:
# Need to create an instance of IVPTestSuite.Solver such that it can
# be called. This includes:
#
# - specify the call signature
# - transform fn!, jac! and mass! such that the solver can use it.

# should have it in a @require block but that doesn't currently work
# https://github.com/one-more-minute/Requires.jl/issues/1

function wrapped_solver(tr::TestRun)
    # Wraps the specific MyPkg.mysolver such that it works within
    # IVPTestSuite setup.
    tc = tr.tc
    so = tr.solver
    ###
    # 0) Wrap tc.fn!, tc.jac!, tc.mass! if necessary

    ###
    # 1) Make call signature

    args = ...
    kwargs = ...

    ###
    # 2) Call solver, if it does not succeed throw an error (if that
    # is not done anyway)
    out = so.solverfn(args...; kwargs...)
    # (probably no need to modify this section)

    ###
    # 3) Transform output to conform to standard:
    # tend -- end time reached
    # yend -- solution at tend
    # stats -- statistics, if available: (steps_total,steps_accepted, fn_evals, jac_evals, linear_solves)
    #                         otherwise  (-1, -1, -1, -1, -1)
    ...
    return tend, yend, stats
end

solverfn = ...::Function # the actual solver function, for example ODE.ode23s.
# choose from
typ =  [:ex, :im, :imex][2]                    # Typ: :ex (explicit method), :im (implicit method), :imex (IMEX method)
stiffness = [nonstiff, mildlystiff, stiff][1]  # stiff, mildly stiff or non-stiff solver
adaptive = [true, false][1]                    # is the solver adaptive
daeindex = 0                                   # maximum index of DAE it can solve, ==0 for ODE only
eq_type = [ explicit_eq, explicit_mass_eq, implicit_eq][1] # [dy/dt=F(t,y); M(t,y)*dy/dt=F(t,y); F(t,y,dy/dt)=0]

solverfn = MyPkg.myode
solverpkg = MyPkg

myodesolver = Solver{typ}(solverfn, solverpkg, wrapped_solver, stiffness, adaptive, daeindex, explicit_eq)

allsolvers[solverfn] = myodesolver


text:
using OrdinaryDiffEq: OrdinaryDiffEq

"""
    Classical <: DynamicsMethods.Method

Type for performing classical molecular dynamics.

```jldoctest
sim = Simulation{Classical}(Atoms(:H), Harmonic())

# output

Simulation{Classical}:
  Atoms{Float64}([:H], [1], [1837.4715941070515])
  Harmonic{Float64, Float64, Float64}
  m: Float64 1.0
  ω: Float64 1.0
  r₀: Float64 0.0
  dofs: Int64 1
```
"""
struct Classical <: DynamicsMethods.Method end

function NQCDynamics.Simulation{Classical}(atoms::Atoms, model::Model; kwargs...)
    NQCDynamics.Simulation(atoms, model, Classical(); kwargs...)
end
function NQCDynamics.RingPolymerSimulation{Classical}(atoms::Atoms, model::Model, n_beads::Integer; kwargs...)
    NQCDynamics.RingPolymerSimulation(atoms, model, Classical(), n_beads::Integer; kwargs...)
end

function NQCDynamics.Simulation(atoms::Atoms, model::Model; kwargs...)
    NQCDynamics.Simulation{Classical}(atoms, model; kwargs...)
end
function NQCDynamics.RingPolymerSimulation(atoms::Atoms, model::Model, n_beads::Integer; kwargs...)
    NQCDynamics.RingPolymerSimulation{Classical}(atoms, model, n_beads; kwargs...)
end

function DynamicsMethods.motion!(du, u, sim::AbstractSimulation{<:Classical}, t)
    dr = DynamicsUtils.get_positions(du)
    dv = DynamicsUtils.get_velocities(du)
    r = DynamicsUtils.get_positions(u)
    v = DynamicsUtils.get_velocities(u)
    DynamicsUtils.velocity!(dr, v, r, sim, t)
    acceleration!(dv, v, r, sim, t)
end

function DynamicsMethods.motion!(du, u, sim::RingPolymerSimulation{<:Classical}, t)
    dr = DynamicsUtils.get_positions(du)
    dv = DynamicsUtils.get_velocities(du)
    r = DynamicsUtils.get_positions(u)
    v = DynamicsUtils.get_velocities(u)
    DynamicsUtils.velocity!(dr, v, r, sim, t)
    ring_polymer_acceleration!(dv, v, r, sim, t)
end

"""
`f1` in `DifferentialEquations.jl` docs.
"""
function acceleration!(dv, v, r, sim::AbstractSimulation, t)
    Calculators.evaluate_derivative!(sim.calculator, r)
    dv .= -sim.calculator.derivative
    DynamicsUtils.divide_by_mass!(dv, sim.atoms.masses)
end

function ring_polymer_acceleration!(dv, v, r, sim::RingPolymerSimulation, t)
    Calculators.evaluate_derivative!(sim.calculator, r)
    dv .= -sim.calculator.derivative
    DynamicsUtils.divide_by_mass!(dv, sim.atoms.masses)
    DynamicsUtils.apply_interbead_coupling!(dv, r, sim)
end

function DynamicsMethods.create_problem(u0, tspan::Tuple, sim::Simulation{<:Classical})
    OrdinaryDiffEq.DynamicalODEProblem(acceleration!, DynamicsUtils.velocity!,
        DynamicsUtils.get_velocities(u0), DynamicsUtils.get_positions(u0), tspan, sim)
end

function DynamicsMethods.create_problem(u0, tspan::Tuple, sim::RingPolymerSimulation{<:Classical})
    OrdinaryDiffEq.DynamicalODEProblem(ring_polymer_acceleration!, DynamicsUtils.velocity!,
        DynamicsUtils.get_velocities(u0), DynamicsUtils.get_positions(u0), tspan, sim)
end

DynamicsMethods.select_algorithm(::AbstractSimulation{<:Classical}) = OrdinaryDiffEq.VelocityVerlet()


text:
abstract type AbstractHermiteBasis{P} <: AbstractMultipleOrthogonalBasis{P} end

polynomial_type(::Type{<:AbstractHermiteBasis}, V::Type) = MP.polynomialtype(V, Int)

even_odd_separated(::Type{<:AbstractHermiteBasis}) = true

reccurence_second_coef(::Type{<:AbstractHermiteBasis}, degree) = 0
reccurence_deno_coef(::Type{<:AbstractHermiteBasis}, degree) = 1

"""
    struct ProbabilistsHermiteBasis{P} <: AbstractHermiteBasis{P}
        polynomials::Vector{P}
    end

Orthogonal polynomial with respect to the univariate weight function ``w(x) = \\exp(-x^2/2)`` over the interval ``[-\\infty, \\infty]``.
"""
struct ProbabilistsHermiteBasis{P} <: AbstractHermiteBasis{P}
    polynomials::Vector{P}
end
reccurence_first_coef(::Type{<:ProbabilistsHermiteBasis}, degree) = 1
reccurence_third_coef(::Type{<:ProbabilistsHermiteBasis}, degree) = -(degree - 1)
degree_one_univariate_polynomial(::Type{<:ProbabilistsHermiteBasis}, variable::MP.AbstractVariable) = MA.@rewrite(1variable)

"""
    struct PhysicistsHermiteBasis{P} <: AbstractHermiteBasis{P}
        polynomials::Vector{P}
    end

Orthogonal polynomial with respect to the univariate weight function ``w(x) = \\exp(-x^2)`` over the interval ``[-\\infty, \\infty]``.
"""
struct PhysicistsHermiteBasis{P} <: AbstractHermiteBasis{P}
    polynomials::Vector{P}
end
reccurence_first_coef(::Type{<:PhysicistsHermiteBasis}, degree) = 2
reccurence_third_coef(::Type{<:PhysicistsHermiteBasis}, degree) = -2(degree - 1)
degree_one_univariate_polynomial(::Type{<:PhysicistsHermiteBasis}, variable::MP.AbstractVariable) = MA.@rewrite(2variable)

max_stars_repo_name: JuliaAlgebra/PolynomialBases

max_issues_repo_name: JuliaAlgebra/PolynomialBases

max_forks_repo_name: JuliaAlgebra/PolynomialBases


text:
### ChainStandard

# Stores the output of a Monte Carlo run
type ChainStandard{N<:Number,T<:AbstractFloat,V<:AbstractVector,A<:AbstractArray} <: AbstractChain
    values::A
    loglikelihood::V
    numburnin::Int
    accepted::Int
    proposed::Int
    runtime::Real
    ChainStandard(v::AbstractArray{N},ll::AbstractVector{T},b::Int,a::Int,p::Int,r::Real) = new(v,ll,b,a,p,r)
end

###Factory function
@inline function _chain{N<:Number,T<:AbstractFloat}(::Type{Val{:standard}},nparas::Int,nsamples::Int,::Type{N},::Type{T};numburnin::Int =0)
    v = zeros(N,_valuestuple(nparas,nsamples))
    ChainStandard{N,T,Vector,Array}(v,Vector{T}(nsamples),numburnin,0,0,0.0)
end

function store!(c::ChainStandard,s::AbstractSample,j::Int)
    nparas = numparas(s)
    @assert numparas(c) == nparas && j <= numsamples(s)
    if c.proposed < numsamples(c)
        c.proposed += 1
        copy!(c.values,(c.proposed-1)*nparas+1,s.values,(j-1)*nparas+1,nparas)
        c.loglikelihood[c.proposed] = s.loglikelihood[j]
    else
        warn("Chain is full. Additional results cannot be stored")
    end
end

function show(io::IO,c::ChainStandard)
  println("ChainStandard with numparas = $(numparas(c)) and numsamples = $(numsamples(c))")
  println("Samples proposed = $(c.proposed), samples accepted = $(c.accepted), acceptance rate = $(c.accepted/c.proposed)")
  println("Total runtime = $(c.runtime)")
  println("Additional fields: :values, :loglikelihood")
end


text:
# Record Helper
# =============
#
# Utilities to handle records in BioJulia packages.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/BioCore.jl/blob/master/LICENSE.md

module RecordHelper

function compare_memory(p1::Ptr, p2::Ptr, len::Integer)
    return ccall(:memcmp, Cint, (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t), p1, p2, len) % Int
end

function compare_memory(data1, offset1, data2, offset2, len)
    return compare_memory(pointer(data1, offset1), pointer(data2, offset2), len)
end

# r"[0-9]+" must match `data[range]`.
function unsafe_parse_decimal(::Type{T}, data::Vector{UInt8}, range::UnitRange{Int}) where {T<:Unsigned}
    x = zero(T)
    @inbounds for i in range
        x = Base.Checked.checked_mul(x, 10 % T)
        x = Base.Checked.checked_add(x, (data[i] - UInt8('0')) % T)
    end
    return x
end

# r"[-+]?[0-9]+" must match `data[range]`.
function unsafe_parse_decimal(::Type{T}, data::Vector{UInt8}, range::UnitRange{Int}) where {T<:Signed}
    lo = first(range)
    if data[lo] == UInt8('-')
        sign = T(-1)
        lo += 1
    elseif data[lo] == UInt8('+')
        sign = T(+1)
        lo += 1
    else
        sign = T(+1)
    end
    x = zero(T)
    @inbounds for i in lo:last(range)
        x = Base.Checked.checked_mul(x, 10 % T)
        x = Base.Checked.checked_add(x, (data[i] - UInt8('0')) % T)
    end
    return sign * x
end

end


text:
# ---
# title: 904. Fruit Into Baskets
# id: problem904
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Two Pointers
# link: <https://leetcode.com/problems/fruit-into-baskets/description/>
# hidden: true
# ---
#
# In a row of trees, the `i`-th tree produces fruit with type `tree[i]`.
#
# You **start at any tree  of your choice**, then repeatedly perform the
# following steps:
#
#   1. Add one piece of fruit from this tree to your baskets.  If you cannot, stop.
#   2. Move to the next tree to the right of the current tree.  If there is no tree to the right, stop.
#
# Note that you do not have any choice after the initial choice of starting
# tree: you must perform step 1, then step 2, then back to step 1, then step 2,
# and so on until you stop.
#
# You have two baskets, and each basket can carry any quantity of fruit, but you
# want each basket to only carry one type of fruit each.
#
# What is the total amount of fruit you can collect with this procedure?
#
#
#
# **Example 1:**
#
#
#
#     Input: [1,2,1]
#     Output: 3
#     Explanation: We can collect [1,2,1].
#
#
# **Example 2:**
#
#
#
#     Input: [0,1,2,2]
#     Output: 3
#     Explanation: We can collect [1,2,2].
#     If we started at the first tree, we would only collect [0, 1].
#
#
# **Example 3:**
#
#
#
#     Input: [1,2,3,2,2]
#     Output: 4
#     Explanation: We can collect [2,3,2,2].
#     If we started at the first tree, we would only collect [1, 2].
#
#
# **Example 4:**
#
#
#
#     Input: [3,3,3,1,2,1,1,2,3,3,4]
#     Output: 5
#     Explanation: We can collect [1,2,1,1,2].
#     If we started at the first tree or the eighth tree, we would only collect 4 fruits.
#
#
#
#
# **Note:**
#
#   1. `1 <= tree.length <= 40000`
#   2. `0 <= tree[i] < tree.length`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
"""
    LookupArray

Types defining the behaviour of a lookup index, how it is plotted
and how [`Selector`](@ref)s like [`Between`](@ref) work.

A `LookupArray` may be [`NoLookup`](@ref) indicating that the index is just the
underlying array axis, [`Categorical`](@ref) for ordered or unordered categories,
or a [`Sampled`](@ref) index for [`Points`](@ref) or [`Intervals`](@ref).
"""
abstract type LookupArray{T,N} <: AbstractArray{T,N} end

const LookupArrayTuple = Tuple{<:LookupArray,Vararg{<:LookupArray}}

span(lookup::LookupArray) = NoSpan()
sampling(lookup::LookupArray) = NoSampling()

dims(::LookupArray) = nothing
val(l::LookupArray) = parent(l)
index(l::LookupArray) = parent(l)
locus(l::LookupArray) = Center()

Base.parent(l::LookupArray) = l.data
Base.size(l::LookupArray) = size(parent(l))
Base.axes(l::LookupArray) = axes(parent(l))
Base.first(l::LookupArray) = first(parent(l))
Base.last(l::LookupArray) = last(parent(l))
Base.firstindex(l::LookupArray) = firstindex(parent(l))
Base.lastindex(l::LookupArray) = lastindex(parent(l))
function Base.:(==)(l1::LookupArray, l2::LookupArray)
    typeof(l1) == typeof(l2) && parent(l1) == parent(l2)
end

ordered_first(l::LookupArray) = l[ordered_firstindex(l)]
ordered_last(l::LookupArray) = l[ordered_lastindex(l)]

ordered_firstindex(l::LookupArray) = ordered_firstindex(order(l), l)
ordered_lastindex(l::LookupArray) = ordered_lastindex(order(l), l)
ordered_firstindex(o::ForwardOrdered, l::LookupArray) = firstindex(parent(l))
ordered_firstindex(o::ReverseOrdered, l::LookupArray) = lastindex(parent(l))
ordered_firstindex(o::Unordered, l::LookupArray) = firstindex(parent(l))
ordered_lastindex(o::ForwardOrdered, l::LookupArray) = lastindex(parent(l))
ordered_lastindex(o::ReverseOrdered, l::LookupArray) = firstindex(parent(l))
ordered_lastindex(o::Unordered, l::LookupArray) = lastindex(parent(l))

function Base.searchsortedfirst(lookup::LookupArray, val; lt=<)
    searchsortedfirst(parent(lookup), unwrap(val); order=ordering(order(lookup)), lt=lt)
end
function Base.searchsortedlast(lookup::LookupArray, val; lt=<)
    searchsortedlast(parent(lookup), unwrap(val); order=ordering(order(lookup)), lt=lt)
end

function Adapt.adapt_structure(to, l::LookupArray)
    rebuild(l; data=Adapt.adapt(to, parent(l)))
end

"""
    AutoLookup <: LookupArray

    AutoLookup()
    AutoLookup(index=AutoIndex(); kw...)

Automatic [`LookupArray`](@ref), the default lookup. It will be converted automatically
to another [`LookupArray`](@ref) when it is possible to detect it from the index.

Keywords will be used in the detected `LookupArray` constructor.
"""
struct AutoLookup{T,A<:AbstractVector{T},K} <: LookupArray{T,1}
    data::A
    kw::K
end
AutoLookup(index=AutoIndex(); kw...) = AutoLookup(index, kw)

order(lookup::AutoLookup) = hasproperty(lookup.kw, :order) ? lookup.kw.order : AutoOrder()
span(lookup::AutoLookup) = hasproperty(lookup.kw, :span) ? lookup.kw.span : AutoSpan()
sampling(lookup::AutoLookup) = hasproperty(lookup.kw, :sampling) ? lookup.kw.sampling : AutoSampling()
metadata(lookup::AutoLookup) = hasproperty(lookup.kw, :metadata) ? lookup.kw.metadata : NoMetadata()

Base.step(lookup::AutoLookup) = Base.step(parent(lookup))

bounds(lookup::LookupArray) = _bounds(order(lookup), lookup)

_bounds(::ForwardOrdered, l::LookupArray) = first(l), last(l)
_bounds(::ReverseOrdered, l::LookupArray) = last(l), first(l)
_bounds(::Unordered, l::LookupArray) = (nothing, nothing)

@noinline Base.step(lookup::T) where T <: LookupArray =
    error("No step provided by $T. Use a `Sampled` with `Regular`")

"""
    Aligned <: LookupArray

Abstract supertype for [`LookupArray`](@ref)s
where the index is aligned with the array axes.

This is by far the most common supertype for `LookupArray`.
"""
abstract type Aligned{T,O} <: LookupArray{T,1} end

order(lookup::Aligned) = lookup.order

"""
    NoLookup <: LookupArray

    NoLookup()

A [`LookupArray`](@ref) that is identical to the array axis.
[`Selector`](@ref)s can't be used on this lookup.

## Example

Defining a `DimArray` without passing an index
to the dimensions, it will be assigned `NoLookup`:

```jldoctest NoLookup
using DimensionalData

A = DimArray(rand(3, 3), (X, Y))
Dimensions.lookup(A)

# output

NoLookup, NoLookup
```

Which is identical to:

```jldoctest NoLookup
using .LookupArrays
A = DimArray(rand(3, 3), (X(NoLookup()), Y(NoLookup())))
Dimensions.lookup(A)

# output

NoLookup, NoLookup
```
"""
struct NoLookup{A<:AbstractVector{Int}} <: Aligned{Int,Order}
    data::A
end
NoLookup() = NoLookup(AutoIndex())

order(lookup::NoLookup) = ForwardOrdered()
span(lookup::NoLookup) = Regular(1)

Base.step(lookup::NoLookup) = 1

rebuild(l::NoLookup; data=parent(l), kw...) = NoLookup(data)

"""
    AbstractSampled <: Aligned

Abstract supertype for [`LookupArray`](@ref)s where the index is
aligned with the array, and is independent of other dimensions. [`Sampled`](@ref)
is provided by this package.

`AbstractSampled` must have  `order`, `span` and `sampling` fields,
or a `rebuild` method that accpts them as keyword arguments.
"""
abstract type AbstractSampled{T,O<:Order,Sp<:Span,Sa<:Sampling} <: Aligned{T,O} end

span(lookup::AbstractSampled) = lookup.span
sampling(lookup::AbstractSampled) = lookup.sampling
metadata(lookup::AbstractSampled) = lookup.metadata
locus(lookup::AbstractSampled) = locus(sampling(lookup))

Base.step(lookup::AbstractSampled) = step(span(lookup))

for f in (:getindex, :view, :dotview)
    @eval begin
        # Regular span may need its step size updated
        @propagate_inbounds Base.$f(l::AbstractSampled, i::AbstractRange) =
            rebuild(l; data=Base.$f(parent(l), i), span=slicespan(l, i))
    end
end

function Adapt.adapt_structure(to, l::AbstractSampled)
    rebuild(l; data=Adapt.adapt(to, parent(l)), metadata=NoMetadata())
end

# bounds
bounds(lookup::AbstractSampled) = _bounds(sampling(lookup), span(lookup), lookup)

_bounds(::Points, span, lookup::AbstractSampled) = _bounds(order(lookup), lookup)
_bounds(::Intervals, span::Irregular, lookup::AbstractSampled) = bounds(span)
_bounds(sampling::Intervals, span::Explicit, lookup::AbstractSampled) =
    _bounds(order(lookup), sampling, span, lookup)
_bounds(::ForwardOrdered, ::Intervals, span::Explicit, lookup::AbstractSampled) =
    (val(span)[1, 1], val(span)[2, end])
_bounds(::ReverseOrdered, ::Intervals, span::Explicit, lookup::AbstractSampled) =
    (val(span)[1, end], val(span)[2, 1])
_bounds(::Intervals, span::Regular, lookup::AbstractSampled) =
    _bounds(locus(lookup), order(lookup), span, lookup)
_bounds(::Start, ::ForwardOrdered, span, lookup) = first(lookup), last(lookup) + step(span)
_bounds(::Start, ::ReverseOrdered, span, lookup) = last(lookup), first(lookup) - step(span)
_bounds(::Center, ::ForwardOrdered, span, lookup) =
    first(lookup) - step(span) / 2, last(lookup) + step(span) / 2
_bounds(::Center, ::ReverseOrdered, span, lookup) =
    last(lookup) + step(span) / 2, first(lookup) - step(span) / 2
_bounds(::End, ::ForwardOrdered, span, lookup) = first(lookup) - step(span), last(lookup)
_bounds(::End, ::ReverseOrdered, span, lookup) = last(lookup) + step(span), first(lookup)

"""
    Sampled <: AbstractSampled

    Sampled(data::AbstractVector, order::Order, span::Span, sampling::Sampling, metadata)
    Sampled(; data=AutoIndex(), order=AutoOrder(), span=AutoSpan(), sampling=Points(), metadata=NoMetadata())

A concrete implementation of the [`LookupArray`](@ref)
[`AbstractSampled`](@ref). It can be used to represent
[`Points`](@ref) or [`Intervals`](@ref).

`Sampled` is capable of representing gridded data from a wide range of sources, allowing
correct `bounds` and [`Selector`](@ref)s for points or intervals of regular,
irregular, forward and reverse indexes.

On `AbstractDimArray` construction, `Sampled` lookup is assigned for all lookups of
`AbstractRange` not assigned to [`Categorical`](@ref).

## Arguments

- `data`: An `AbstractVector` of index values, matching the length of the curresponding
    array axis.
- `order`: [`Order`](@ref)) indicating the order of the index,
    [`AutoOrder`](@ref) by default, detected from the order of `data`
    to be [`ForwardOrdered`](@ref), [`ReverseOrdered`](@ref) or [`Unordered`](@ref).
    These can be provided explicitly if they are known and performance is important.
- `span`: indicates the size of intervals or distance between points, and will be set to
    [`Regular`](@ref) for `AbstractRange` and [`Irregular`](@ref) for `AbstractArray`,
    unless assigned manually.
- `sampling`: is assigned to [`Points`](@ref), unless set to [`Intervals`](@ref) manually.
    Using [`Intervals`](@ref) will change the behaviour of `bounds` and `Selectors`s
    to take account for the full size of the interval, rather than the point alone.
- `metadata`: a `Dict` or `Metadata` wrapper that holds any metadata object adding more
    information about the array axis - useful for extending DimensionalData for specific
    contexts, like geospatial data in GeoData.jl. By default it is `NoMetadata()`.

## Example

Create an array with [`Interval`] sampling, and `Regular` span for a vector with known spacing.

We set the [`Locus`](@ref) of the `Intervals` to `Start` specifying
that the index values are for the positions at the start of each interval.

```jldoctest Sampled
using DimensionalData, DimensionalData.LookupArrays

x = X(Sampled(100:-20:10; sampling=Intervals(Start())))
y = Y(Sampled([1, 4, 7, 10]; span=Regular(3), sampling=Intervals(Start())))
A = ones(x, y)

# output
5×4 DimArray{Float64,2} with dimensions:
  X Sampled 100:-20:20 ReverseOrdered Regular Intervals,
  Y Sampled Int64[1, 4, 7, 10] ForwardOrdered Regular Intervals
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
```
"""
struct Sampled{T,A<:AbstractVector{T},O,Sp,Sa,M} <: AbstractSampled{T,O,Sp,Sa}
    data::A
    order::O
    span::Sp
    sampling::Sa
    metadata::M
end
function Sampled(
    data=AutoIndex();
    order=AutoOrder(), span=AutoSpan(),
    sampling=AutoSampling(), metadata=NoMetadata()
)
    Sampled(data, order, span, sampling, metadata)
end

function rebuild(l::Sampled;
    data=parent(l), order=order(l), span=span(l), sampling=sampling(l), metadata=metadata(l), kw...
)
    Sampled(data, order, span, sampling, metadata)
end

"""
    AbstractCategorical <: Aligned

[`LookupArray`](@ref)s where the values are categories.

[`Categorical`](@ref) is the provided concrete implementation.
but this can easily be extended - all methods are defined for `AbstractCategorical`.

All `AbstractCategorical` must provide a `rebuild`
method with `data`, `order` and `metadata` keyword arguments.
"""
abstract type AbstractCategorical{T,O} <: Aligned{T,O} end

order(lookup::AbstractCategorical) = lookup.order
metadata(lookup::AbstractCategorical) = lookup.metadata

const CategoricalEltypes = Union{AbstractChar,Symbol,AbstractString}

function Adapt.adapt_structure(to, l::AbstractCategorical)
    rebuild(l; data=Adapt.adapt(to, parent(l)), metadata=NoMetadata())
end

"""
    Categorical <: AbstractCategorical

    Categorical(o::Order)
    Categorical(; order=Unordered())

An LookupArray where the values are categories.

This will be automatically assigned if the index contains `AbstractString`,
`Symbol` or `Char`. Otherwise it can be assigned manually.

[`Order`](@ref) will be determined automatically where possible.

## Arguments

- `data`: An `AbstractVector` of index values, matching the length of the curresponding
    array axis.
- `order`: [`Order`](@ref)) indicating the order of the index,
    [`AutoOrder`](@ref) by default, detected from the order of `data`
    to be `ForwardOrdered`, `ReverseOrdered` or `Unordered`.
    Can be provided if this is known and performance is important.
- `metadata`: a `Dict` or `Metadata` wrapper that holds any metadata object adding more
    information about the array axis - useful for extending DimensionalData for specific
    contexts, like geospatial data in GeoData.jl. By default it is `NoMetadata()`.

## Example

Create an array with [`Interval`] sampling.

```jldoctest Categorical
using DimensionalData

ds = X(["one", "two", "three"]), Y([:a, :b, :c, :d])
A = DimArray(rand(3, 4), ds)
Dimensions.lookup(A)

# output

Categorical String[one, two, three] Unordered,
Categorical Symbol[a, b, c, d] ForwardOrdered
```
"""
struct Categorical{T,A<:AbstractVector{T},O<:Order,M} <: AbstractCategorical{T,O}
    data::A
    order::O
    metadata::M
end
function Categorical(data=AutoIndex(); order=AutoOrder(), metadata=NoMetadata())
    Categorical(data, order, metadata)
end

function rebuild(l::Categorical;
    data=parent(l), order=order(l), metadata=metadata(l), kw...
)
    Categorical(data, order, metadata)
end

"""
    Unaligned <: LookupArray

Abstract supertype for [`LookupArray`](@ref) where the index is not aligned to the grid.

Indexing an [`Unaligned`](@ref) with [`Selector`](@ref)s must provide all
other [`Unaligned`](@ref) dimensions.
"""
abstract type Unaligned{T,N} <: LookupArray{T,N} end

"""
    Transformed <: Unaligned

    Transformed(f, dim::Dimension; metadata=NoMetadata())

[`LookupArray`](@ref) that uses an affine transformation to convert
dimensions from `dims(lookup)` to `dims(array)`. This can be useful
when the dimensions are e.g. rotated from a more commonly used axis.

Any function can be used to do the transformation, but transformations
from CoordinateTransformations.jl may be useful.

## Arguments

- `f`: transformation function
- `dim`: a dimension to transform to.

## Keyword Arguments

- `metdata`:

## Example

```jldoctest
using DimensionalData, DimensionalData.LookupArrays, CoordinateTransformations

m = LinearMap([0.5 0.0; 0.0 0.5])
A = [1 2  3  4
     5 6  7  8
     9 10 11 12];
da = DimArray(A, (t1=Transformed(m, X), t2=Transformed(m, Y)))

da[X(At(6)), Y(At(2))]

# output
9
```
"""
struct Transformed{T,A<:AbstractVector{T},F,D,M} <: Unaligned{T,1}
    data::A
    f::F
    dim::D
    metadata::M
end
function Transformed(f, dim; metadata=NoMetadata())
    Transformed(AutoIndex(), f, basetypeof(dim)(), metadata)
end

function rebuild(l::Transformed;
    data=data(l), f=f(l), dim=dim(l), metadata=metadata(l)
)
    Transformed(data, f, dim, metadata)
end

f(lookup::Transformed) = lookup.f
dim(lookup::Transformed) = lookup.dim

transformfunc(lookup::Transformed) = f(lookup)
transformdim(x) = nothing
transformdim(lookup::Transformed) = lookup.dim
transformdim(::Type{<:Transformed{<:Any,<:Any,<:Any,D}}) where D = D

Base.:(==)(l1::Transformed, l2::Transformed) = typeof(l1) == typeof(l2) && f(l1) == f(l2)

# TODO Transformed bounds

# Shared methods

slicespan(l::LookupArray, i::Colon) = span(l)
slicespan(l::LookupArray, i) = _slicespan(span(l), l, i)

_slicespan(span::Regular, l::LookupArray, i::Int) = span
_slicespan(span::Regular, l::LookupArray, i::AbstractRange) = Regular(step(l) * step(i))
_slicespan(span::Regular, l::LookupArray, i::AbstractArray) = _slicespan(Irregular(bounds(l)), l, i)
_slicespan(span::Explicit, l::LookupArray, i::Int) = Explicit(val(span)[:, i])
_slicespan(span::Explicit, l::LookupArray, i::AbstractArray) = Explicit(val(span)[:, i])
_slicespan(span::Irregular, l::LookupArray, i::StandardIndices) =
    _slicespan(sampling(l), span, l, i)
_slicespan(::Points, span::Irregular, l::LookupArray, i::StandardIndices) = span
_slicespan(::Intervals, span::Irregular, l::LookupArray, i::StandardIndices) =
    Irregular(_maybeflipbounds(l, _slicebounds(locus(l), span, l, i)))

function _slicebounds(locus::Start, span::Irregular, l::LookupArray, i::StandardIndices)
    l[first(i)], last(i) >= lastindex(l) ? _maybeflipbounds(l, bounds(span))[2] : l[last(i) + 1]
end
function _slicebounds(locus::End, span::Irregular, l::LookupArray, i::StandardIndices)
    first(i) <= firstindex(l) ? _maybeflipbounds(l, bounds(span))[1] : l[first(i) - 1], l[last(i)]
end
function _slicebounds(locus::Center, span::Irregular, l::LookupArray, i::StandardIndices)
    first(i) <= firstindex(l) ? _maybeflipbounds(l, bounds(span))[1] : (l[first(i) - 1] + l[first(i)]) / 2,
    last(i)  >= lastindex(l)  ? _maybeflipbounds(l, bounds(span))[2] : (l[last(i) + 1]  + l[last(i)]) / 2
end
# Have to special-case date/time so we work with seconds and add to the original
function _slicebounds(locus::Center, span::Irregular, l::LookupArray{T}, i::StandardIndices) where T<:Dates.AbstractTime
    op = T === Date ? div : /
    frst = if first(i) <= firstindex(l)
        _maybeflipbounds(l, bounds(span))[1]
    else
        if isrev(order(l))
            op(l[first(i)] - l[first(i) - 1], 2) + l[first(i) - 1]
        else
            op(l[first(i) - 1] - l[first(i)], 2) + l[first(i)]
        end
    end
    lst = if last(i) >= lastindex(l)
        _maybeflipbounds(l, bounds(span))[2]
    else
        if isrev(order(l))
            op(l[last(i)] - l[last(i) + 1], 2) + l[last(i) + 1]
        else
            op(l[last(i) + 1] - l[last(i)], 2) + l[last(i)]
        end
    end
    return (frst, lst)
end

# reducing methods
@inline reducelookup(lookup::NoLookup) = NoLookup(OneTo(1))
# TODO what should this do?
@inline reducelookup(lookup::Unaligned) = NoLookup(OneTo(1))
# Categories are combined.
@inline reducelookup(lookup::Categorical{<:AbstractString}) =
    rebuild(lookup; data=["combined"])
@inline reducelookup(lookup::Categorical) = rebuild(lookup; data=[:combined])
# Sampled is resampled
@inline reducelookup(lookup::AbstractSampled) = _reducelookup(span(lookup), lookup)

@inline _reducelookup(::Irregular, lookup::AbstractSampled) = begin
    rebuild(lookup; data=_reduceindex(lookup), order=ForwardOrdered())
end
@inline _reducelookup(span::Regular, lookup::AbstractSampled) = begin
    newstep = step(span) * length(lookup)
    newindex = _reduceindex(lookup, newstep)
    # Make sure the step type matches the new index eltype
    newstep = convert(promote_type(eltype(newindex), typeof(newstep)), newstep)
    newspan = Regular(newstep)
    rebuild(lookup; data=newindex, order=ForwardOrdered(), span=newspan)
end
@inline _reducelookup(
    span::Regular{<:Dates.CompoundPeriod}, lookup::AbstractSampled
) = begin
    newstep = Dates.CompoundPeriod(step(span).periods .* length(lookup))
    # We don't pass the step here - the range doesn't work with CompoundPeriod
    newindex = _reduceindex(lookup)
    # Make sure the step type matches the new index eltype
    newspan = Regular(newstep)
    rebuild(lookup; data=newindex, order=ForwardOrdered(), span=newspan)
end
@inline _reducelookup(span::Explicit, lookup::AbstractSampled) = begin
    bnds = val(span)
    newstep = bnds[2] - bnds[1]
    newindex = _reduceindex(lookup, newstep)
    # Make sure the step type matches the new index eltype
    newstep = convert(promote_type(eltype(newindex), typeof(newstep)), newstep)
    newspan = Explicit(reshape([bnds[1, 1]; bnds[2, end]], 2, 1))
    newlookup = rebuild(lookup; data=newindex, order=ForwardOrdered(), span=newspan)
end
# Get the index value at the reduced locus.
# This is the start, center or end point of the whole index.
@inline _reduceindex(lookup::LookupArray, step=nothing) = _reduceindex(locus(lookup), lookup, step)
@inline _reduceindex(locus::Start, lookup::LookupArray, step) = _mayberange(first(lookup), step)
@inline _reduceindex(locus::End, lookup::LookupArray, step) = _mayberange(last(lookup), step)
@inline _reduceindex(locus::Center, lookup::LookupArray, step) = begin
    index = parent(lookup)
    len = length(index)
    newval = centerval(index, len)
    _mayberange(newval, step)
end
# Ranges with a known step always return a range
_mayberange(x, step) = x:step:x
# Arrays return a vector
_mayberange(x, step::Nothing) = [x]

@inline centerval(index::AbstractArray{<:Number}, len) = (first(index) + last(index)) / 2
@inline function centerval(index::AbstractArray{<:DateTime}, len)
    f = first(index)
    l = last(index)
    if f <= l
        return (l - f) / 2 + first(index)
    else
        return (f - l) / 2 + last(index)
    end
end
@inline centerval(index::AbstractArray, len) = index[len ÷ 2 + 1]

ordering(::ForwardOrdered) = Base.Order.ForwardOrdering()
ordering(::ReverseOrdered) = Base.Order.ReverseOrdering()


text: @inline SIMD.Vec(g::Base.Generator{<:LengthRange}) = Vec(Tuple(g))


text:
import Neapolitan

function main()
  Neapolitan.render(STDIN)
end


text:
module FTUtils

export ftdiff, ftaxis, ft, ift, psd, ssb

function ftdiff(N::Int, dx::Real)
    return 1/(N*dx)
end

function ftaxis(N::Int, dx::Real, ssb::Bool=false)
    f = (collect(0:N-1)-fld(N,2)) * ftdiff(N,dx)
    if ssb
        f = f[Int(fld(N,2))+1:end]
    end
    return f
end

function ft(y::Vector{T}, dt::Real) where {T<:Number}
    return y*dt |> ifftshift |> fft |> fftshift
end

function ift(Y::Vector{T}, df::Real) where {T<:Number}
    return Y*df |> ifftshift |> bfft |> fftshift
end

function psd(y::Vector{T}, dt::Real, ssb=true) where {T<:Number}
    N = length(y)
    psd = abs2.(ft(y,dt)) * ftdiff(N-1,dt)
    if ssb
        dc = Int(fld(N,2)) + 1
        psd[dc+1:end] *= 2
        psd = psd[dc:end]
    end
    return psd
end

end # module


text:
#######
##FEM Heat Animation Test
#######

#Generates an animation for a solution of the heat equation
#Uses Plots.jl, requires matplotlib >=1.5
using FiniteElementDiffEq, DiffEqProblemLibrary
#using Plots, ImageMagick; gr()
prob = prob_femheat_moving

sol = solve(prob,FEMDiffEqHeatEuler(),save_everystep=true)
println("Generating Animation")
sol.tslocation = 1
#plot(sol)
#animate(sol,"test_animation.gif";zlims=(0,.1),cbar=false)

## Should have moved off the frame.
maximum(sol[end]) .< 1e-6


text:
"""
    Param

Struct containing interpolation parameters for a single point in a mesh.
* `lower_index  :: Int64`   : index of nearest neighbor in mesh with smaller value
* `upper_index  :: Int64`   : index of nearest neighbor in mesh with larger value
* `lower_weight :: Float64` : interpolation weight for mesh[lower_index]
* `upper_weight :: Float64` : interpolation weight for mesh[upper_index]
"""
struct Param
    lower_index  :: Int64
    upper_index  :: Int64
    lower_weight :: Float64
    upper_weight :: Float64
end

# generate param dummy
function get_param_empty() :: Param

    p = Param(0, 0, 0.0, 0.0)

    return p
end

# find nearest neighbor (lower and upper) indices in sorted list including zero
function get_indices(
    val  :: Float64,
    list :: Vector{Float64}
    )    :: NTuple{2, Int64}

    # init indices
    lower_index = 0
    upper_index = 0

    # check if in bounds otherwise search list
    if val >= list[end]
        lower_index, upper_index = length(list), length(list)
    else
        # iterate over list until upper_index is found (lower_index is then also determined)
        index_current = 1

        while val > list[index_current]
            index_current += 1
        end

        lower_index = index_current
        upper_index = index_current

        if val < list[index_current]
            lower_index -= 1
        end
    end

    return lower_index, upper_index
end

"""
    get_param(
        val  :: Float64,
        list :: Vector{Float64}
        )    :: Param

Compute interpolation parameters of val in a set of discrete points (list) and buffer result in Param struct.
"""
function get_param(
    val  :: Float64,
    list :: Vector{Float64}
    )    :: Param

    # get neighbors and init weights
    lower_index,  upper_index  = get_indices(val, list)
    lower_weight, upper_weight = 0.0, 0.0

    # compute weights
    if lower_index < upper_index
        d            = 1.0 / (list[upper_index] - list[lower_index])
        lower_weight = d * (list[upper_index] - val)
        upper_weight = d * (val - list[lower_index])
    else
        lower_weight = 1.0
    end

    # build param
    p = Param(lower_index, upper_index, lower_weight, upper_weight)

    return p
end


text:
using FourierGPE, Test

@testset "Parseval" begin include("test_parseval.jl") end
@testset "Ground state" begin include("test_groundstate.jl") end
@testset "Dynamics" begin include("test_dynamics.jl") end


text:
"Your optional docstring here"
function distance(a, b)

    if length(a) != length(b)
        throw(ArgumentError("differnt length"))
    end
    dist=0
    for (x,y) in zip(a,b)
        if x != y
            dist=dist+1
        end
    end

    return dist
end

max_stars_repo_name: stepinski/julia-in

max_issues_repo_name: stepinski/julia-init

max_forks_repo_name: stepinski/julia-init


text:
TOL = 0.5e-8

function falseposition(f::Function, a::Real, b::Real;
                       ε::Real = TOL, iₘₐₓ::Int = 100)
    i = 0
    xᵢ = 0
    while f(xᵢ) == 0 && i < iₘₐₓ
        fa, fb, fxᵢ = f(a), f(b), f(xᵢ)
        xᵢ = (b * fa - a * fb) / (fa - fb)

        if fa * fxᵢ < 0
            b = xᵢ
        else
            a = xᵢ
        end
        i += 1
    end
    xᵢ, i
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct UsagesApi <: SwaggerApi
    client::Swagger.Client
end

"""

List network usages for a subscription.
Param: location::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: UsagesListResult
"""
function _swaggerinternal_usagesList(_api::UsagesApi, location::String, api_version::String, subscriptionId::String; _mediaType=nothing)

    _ctx = Swagger.Ctx(_api.client, "GET", UsagesListResult, "/subscriptions/{subscriptionId}/providers/Microsoft.Network/locations/{location}/usages", ["azure_auth"])
    Swagger.set_param(_ctx.path, "location", location)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function usagesList(_api::UsagesApi, location::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_usagesList(_api, location, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function usagesList(_api::UsagesApi, response_stream::Channel, location::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_usagesList(_api, location, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export usagesList