text:
"""
Leaf energy balance module.

* [`bl_cond_heat`](@ref)
* [`bl_cond_vapor`](@ref)
* [`PAR_to_shortwave`](@ref)
* [`sensible_heat`](@ref)
* [`leaf_vapor_deficit`](@ref)
* [`water_flux`](@ref)
* [`latent_heat`](@ref)
* [`energy_imbalance`](@ref)

# References

## Leaf boundary layer conductances

* [CN98] Campbell, G. S. and Norman, J. M. (1998). *An Introduction to
  Environmental Biophysics* (2nd ed.), pp. 99--101. Springer Science+Business
  Media, New York, NY, USA.
* [B08] Bonan, G. (2008). *Ecological Climatology: Concepts and Applications*,
  pp. 229--232. Cambridge University Press, Cambridge, UK.

## Relationship between PAR and shortwave radiation

* [MHH84] Meek, D. W., Hatfield, J. L., Howell, T. A., Idso, S. B., and
  Reginato, R. J. (1984). A generalized relationship between photosynthetically
  active radiation and solar radiation. *Agron. J.*, 76(6), 939--945.
  <https://doi.org/10.2134/agronj1984.00021962007600060018x>
"""
module EBal

include("../docstring_style.jl")

using Canopy.Water: e_sat, latent_heat_vap
using Canopy.Air: air_molar, air_density
using Canopy.RadTrans: stefan_boltzmann
using Canopy.Transfer:
    dyn_visc_moistair,
    prandtl,
    therm_diff_moistair,
    diffus_air,
    heat_cap_moistair
using Canopy.Leaf.StomCond: total_cond_h2o  # TODO: TO BE IMPLEMENTED

export bl_cond_heat,
    bl_cond_vapor,
    PAR_to_shortwave,
    sensible_heat,
    leaf_vapor_deficit,
    water_flux,
    latent_heat,
    energy_imbalance

"""
Calculate leaf boundary layer conductance for heat transfer
(``g_\\mathrm{b,H}``) [mol m^-2 s^-1].

```math
g_\\mathrm{b,H} = \\frac{1.4\\,\\mathrm{Nu}\\cdot \\chi_\\mathrm{a}\\cdot
    D_\\mathrm{H}}{d}
```

where

* Nu = 0.664 Re``^{1/2}`` Pr``^{1/3}`` is the Nusselt number;
* ``\\chi_\\mathrm{a}`` [mol m^-3] is the molar concentration of air;
* ``D_\\mathrm{H}`` [m^2 s^-1] is the thermal diffusivity of air.

The multiplier 1.4 is a correction factor for leaves in field conditions.

# Arguments

* `temp`: Air temperature [K].
* `pressure`: Ambient pressure [Pa].
* `RH`: Relative humidity [0--1].
* `wind_speed`: Wind speed [m s^-1].
* `d_leaf`: Characteristic leaf size [m].

# See also

[`bl_cond_vapor`](@ref)
"""
function bl_cond_heat(temp, pressure, RH, wind_speed, d_leaf)
    # kinematic viscosity of moist air [m^2 s^-1]
    nu_a = dyn_visc_moistair(temp, pressure, RH) /
        air_density(temp, pressure, RH)
    Pr = prandtl(temp, pressure, RH)
    air_conc = air_molar(temp, pressure)  # [mol m^-3]

    return 1.4 * 0.664 * cbrt(Pr) * air_conc *
        therm_diff_moistair(temp, pressure, RH) *
        sqrt(wind_speed / (d_leaf * nu_a))
end

"""
Calculate leaf boundary layer conductance for water vapor transfer
(``g_\\mathrm{b,W}``).

```math
g_\\mathrm{b,W} = \\frac{1.4\\cdot 0.664\\ \\mathrm{Re}^{1/2}\\
    \\mathrm{Sc}^{1/3} \\cdot \\chi_\\mathrm{a} \\cdot D_\\mathrm{W}}{d}
```

where

* ``\\mathrm{Sc} = \\nu_\\mathrm{a} / D_\\mathrm{W}`` is the Schmidt number;
* ``\\nu_\\mathrm{a}`` [m^2 s^-1] is the kinematic viscosity of air;
* ``D_\\mathrm{W}`` [m^2 s^-1] is the diffusivity of water vapor in air;
* ``\\chi_\\mathrm{a}`` [mol m^-3] is the molar concentration of air.

The multiplier 1.4 is a correction factor for leaves in field conditions.

# Arguments

* `temp`: Air temperature [K].
* `pressure`: Ambient pressure [Pa].
* `RH`: Relative humidity [0--1].
* `wind_speed`: Wind speed [m s^-1].
* `d_leaf`: Characteristic leaf size [m].

# See also

[`bl_cond_heat`](@ref)
"""
function bl_cond_vapor(temp, pressure, RH, wind_speed, d_leaf)
    # kinematic viscosity of moist air [m^2 s^-1]
    nu_a = dyn_visc_moistair(temp, pressure, RH) /
        air_density(temp, pressure, RH)
    D_w = diffus_air("h2o", temp, pressure)
    Sc = nu_a / D_w  # Schmidt number
    air_conc = air_molar(temp, pressure)  # [mol m^-3]

    return 1.4 * 0.664 * cbrt(Sc) * air_conc * D_w *
        sqrt(wind_speed / (d_leaf * nu_a))
end

"""
Calculate the incoming shortwave radiation [W m^-2] from photosynthetically
active radiation (PAR) [µmol m^-2 s^-1]. The calculation is based on an
empirical equation fitted between PAR and shortwave radiation measured at
various sites in western US.
"""
PAR_to_shortwave(PAR) = 0.495785820525533 * PAR - 0.08081308874566188

"""
Calculate sensible heat transfer rate at the leaf surface [W m^-2].

```math
\\mathrm{SH} = 2 g_\\mathrm{b,H} \\cdot c_{p,\\mathrm{a}} \\cdot
    \\left( T_\\mathrm{leaf} - T_\\mathrm{air} \\right)
```

Note that both sides of a leaf are taken into account by the factor 2.

# Arguments

* `T_leaf`: Leaf temperature [K].
* `T_air`: Air temperature [K].
* `cp_a`: Molar heat capacity of the moist air [J K^-1 mol^-1].
* `g_bh`: Boundary layer conductance for heat transfer [mol m^-2 s^-1].
"""
sensible_heat(T_leaf, T_air, cp_a, g_bh) = 2. * g_bh * cp_a * (T_leaf - T_air)

"""
Calculate leaf-to-air vapor pressure deficit [Pa].

# Arguments

* `T_leaf`: Leaf temperature [K].
* `T_air`: Air temperature [K].
* `RH`: Relative humidity [0--1].
"""
leaf_vapor_deficit(T_leaf, T_air, RH) = e_sat(T_leaf) - e_sat(T_air) * RH

"""
Calculate leaf water flux [mol m^-2 s^-1].

```math
F_\\mathrm{H_2O} = g_\\mathrm{tot,W} \\cdot
    \\frac{\\mathrm{VPD}}{p_\\mathrm{atm}}
```

where

```math
g_\\mathrm{tot,W} = \\frac{1}{g_\\mathrm{b,W}^{-1} + g_\\mathrm{s,W}^{-1}}
```

is the total conductance of water vapor [mol m^-2 s^-1].

# Arguments

* `VPD_leaf`: Leaf-to-air vapor pressure deficit [Pa].
* `pressure`: Ambient pressure [Pa].
* `g_sw`: Stomatal conductance of water vapor [mol m^-2 s^-1].
* `g_bw`: Boundary layer conductance of water vapor [mol m^-2 s^-1].
"""
water_flux(VPD_leaf, pressure, g_sw, g_bw) = total_cond_h2o(g_bw, g_sw) *
    VPD_leaf / pressure

"""
Calculate latent heat transfer rate at the leaf surface [W m^-2].

```math
\\mathrm{LE} = L_\\mathrm{v} \\cdot g_\\mathrm{tot,W} \\cdot
\\frac{\\mathrm{VPD}}{p_\\mathrm{atm}}
```

where ``L_\\mathrm{v}`` is the latent heat of vaporization of water [J mol^-1],
and ``g_\\mathrm{tot,W}`` is the total conductance of water vapor.

# Arguments

* `T_leaf`: Leaf temperature [K].
* `VPD_leaf`: Leaf-to-air vapor pressure deficit [Pa].
* `pressure`: Ambient pressure [Pa].
* `g_sw`: Stomatal conductance of water vapor [mol m^-2 s^-1].
* `g_bw`: Boundary layer conductance of water vapor [mol m^-2 s^-1].
"""
latent_heat(T_leaf, VPD_leaf, pressure, g_sw, g_bw) =
    water_flux(VPD_leaf, pressure, g_sw, g_bw) *
    latent_heat_vap(T_leaf)

"""
Calculate the imbalance of energy transfer at the leaf surface [W m^-2]. This
function is used as a constraint for the coupled leaf temperature--flux solver.

# Arguments

* `T_leaf`: Leaf temperature [K].
* `T_air`: Air temperature [K].
* `pressure`: Ambient pressure [Pa].
* `RH`: Relative humidity [0--1].
* `R_sw`: Incoming shortwave radiation [W m^-2].
* `wind_speed`: Wind speed [m s^-1].
* `d_leaf`: Characteristic leaf size [m].
* `em_leaf`: Emissivity of leaf [0--1].
* `g_sw`: Stomatal conductance of water vapor [mol m^-2 s^-1].
"""
function energy_imbalance(T_leaf, T_air, pressure, RH, R_sw, wind_speed,
                          d_leaf, em_leaf, g_sw)
    R_net = R_sw - em_leaf * stefan_boltzmann(T_leaf)
    VPD_leaf = leaf_vapor_pressure_deficit(T_leaf, T_air, RH)
    cp_a = heat_cap_moistair(T_leaf, pressure, RH)
    g_bh = bl_cond_heat(T_air, pressure, RH, wind_speed, d_leaf)
    g_bw = bl_cond_vapor(T_air, pressure, RH, wind_speed, d_leaf)
    SH = sensible_heat(T_leaf, T_air, cp_a, g_bh)
    LE = latent_heat(T_leaf, VPD_leaf, pressure, g_sw, g_bw)

    return R_net - SH - LE
end

end  # module


text:
### A Pluto.jl notebook ###
# v0.14.6

using Markdown
using InteractiveUtils

# ╔═╡ 0279facc-3f77-401b-ab39-352a6ac9bc9f
begin
    import Pkg
    Pkg.activate(mktempdir())
    Pkg.add([
        Pkg.PackageSpec(name="Plots", version="1"),
        Pkg.PackageSpec(name="PlutoUI", version="0.7"),
	    Pkg.PackageSpec(name="DifferentialEquations", version="6"),
	    Pkg.PackageSpec(name="Parameters", version="0.12"),
        Pkg.PackageSpec(name="LabelledArrays", version="1"),
        Pkg.PackageSpec(name="Setfield", version="0.7"),
    ])
    using Plots, PlutoUI, LinearAlgebra, DifferentialEquations, Parameters, LabelledArrays, Setfield

	Plots.gr(fmt=:png, lw=2)

	PlutoUI.TableOfContents()
end

# ╔═╡ f0ab21d0-bdcb-11eb-39e8-59f7ff34b88c
md"""
# Fig 1.09 Hodgkin-Huxley model
"""

# ╔═╡ 51592b89-56fc-40b3-be86-8ae48241a80b
# Parameters for Hodgkin-Huxley model
    p = (E_N = 55.0,       # Reversal potential of Na (mV)
		 E_K = -72.0,      # Reversal potential of K (mV)
		 E_LEAK = -49.0,   # Reversal potential of leaky channels (mV)
		 G_N_BAR = 120.0,  # Max. Na channel conductance (mS/cm^2)
		 G_K_BAR = 36.0,   # Max. K channel conductance (mS/cm^2)
		 G_LEAK = 0.30,    # Max. leak channel conductance (mS/cm^2)
		 C_M = 1.0,        # membrane capacitance (uF/cm^2)
		 I_STIM = 0.0     # Stimulus current (uA/cm^2)
    )

# ╔═╡ 8018875c-258e-4516-a36c-2413fee97ccd
begin
	# Events and callbacks
	on_i1!(integrator) = integrator.p = @set p.I_STIM = -6.80
    off_i1!(integrator) = integrator.p = @set p.I_STIM = 0.0
    on_i2!(integrator) = integrator.p = @set p.I_STIM = -6.90

    cbs = CallbackSet(PresetTimeCallback([20.0], on_i1!),
                  PresetTimeCallback([21.0, 61.0], off_i1!),
                  PresetTimeCallback([60.0], on_i2!))
end

# ╔═╡ 999aae28-fe75-4a2f-a55d-12360b647f05
u0 = LVector(v=-59.8977, m=0.0536, h=0.5925, n=0.3192)

# ╔═╡ 23ec33a8-b6f3-4a89-bc24-95474ed5a9cb
tend = 100.0

# ╔═╡ 8e7633bb-1d79-43cc-81ec-8fc740ba13ff
md"""
## Appendix
Running environment and some auxillary functions.

**You can use [this helper tool](https://fonsp.com/article-test-3/pkghelper.html) to generate these commands!**
"""

# ╔═╡ 94adbc56-8f68-4048-bfbc-326acb359fde
begin
	# Convenience functions
    hill(x, k) = x / (x + k)
    hill(x, k, n) = hill(x^n, k^n)
    exprel(x) = ifelse(x≈zero(x), one(x), x / expm1(x))
end

# ╔═╡ 31549c4b-2aec-4c09-9888-a15b2b0ebbd6
"ODE system of HH model"
	function hh!(du, u, p, t)
		@unpack E_N, E_K, E_LEAK, G_N_BAR, G_K_BAR, G_LEAK, C_M, I_STIM = p
		@unpack v, m, h, n = u

		mαV = -0.10 * (v + 35)
		mα = exprel(mαV)
		mβ = 4.0 * exp(-(v + 60) / 18.0)

		hα = 0.07 * exp(-(v+60)/20)
		hβ = 1 / (exp(-(v+30)/10) + 1)

		nαV = -0.1 * (v+50)
		nα = 0.1 * exprel(nαV)
		nβ = 0.125 * exp( -(v+60) / 80)

		iNa = G_N_BAR * (v - E_N) * (m^3) * h
		iK  = G_K_BAR * (v - E_K) * (n^4)
		iLeak = G_LEAK * (v - E_LEAK)
		du.v = -(iNa + iK + iLeak + I_STIM) / C_M
		du.m = -(mα + mβ) * m + mα
		du.h = -(hα + hβ) * h + hα
		du.n = -(nα + nβ) * n + nα
		return du
	end

# ╔═╡ 2e369b99-0f70-4ff2-be75-d3eec6c46777
sol = solve(ODEProblem(hh!, u0, tend, p), callback = cbs)

# ╔═╡ 2cec94d6-fcc8-49b6-970c-b84f8a736f61
plot(sol, vars=(0, 1), ylabel="Membrane potential (mV)", xlabel="Time (ms)", legend=false, title="Figure 1.9")

# ╔═╡ Cell order:
# ╠═f0ab21d0-bdcb-11eb-39e8-59f7ff34b88c
# ╠═31549c4b-2aec-4c09-9888-a15b2b0ebbd6
# ╠═8018875c-258e-4516-a36c-2413fee97ccd
# ╠═51592b89-56fc-40b3-be86-8ae48241a80b
# ╠═999aae28-fe75-4a2f-a55d-12360b647f05
# ╠═23ec33a8-b6f3-4a89-bc24-95474ed5a9cb
# ╠═2e369b99-0f70-4ff2-be75-d3eec6c46777
# ╠═2cec94d6-fcc8-49b6-970c-b84f8a736f61
# ╠═8e7633bb-1d79-43cc-81ec-8fc740ba13ff
# ╠═94adbc56-8f68-4048-bfbc-326acb359fde
# ╠═0279facc-3f77-401b-ab39-352a6ac9bc9f

ext: jl

max_stars_repo_name: stevengogogo/pluto-notebooks

max_issues_repo_name: stevengogogo/pluto-notebooks

max_forks_repo_name: stevengogogo/pluto-notebooks


text:
mutable struct Stand
    treelist::Array{Tree,1}
    basal_area::Union{Real,Missing}
    trees_per_area::Union{Real,Missing}
    qmd::Union{Real,Missing}
end

Stand(tl::Array{Tree,1}) = Stand(tl,missing,missing,missing)


text:
using MetLP
using MetLP.MetNets
using MetLP.ProgressMeter
using Test
using Base.Threads

## ------------------------------------------------------
@testset "MetLP.jl" begin
    include("lp_model_tests.jl")
    include("fva_fba_tests.jl")
    # include("yLP_tests.jl")
end


text:
using CPLEX, Base.Test, MathOptInterface

const MOI  = MathOptInterface
const MOIT = MathOptInterface.Test
const MOIB = MathOptInterface.Bridges

@testset "Unit Tests" begin
    config = MOIT.TestConfig()
    solver = CPLEX.Optimizer()

    MOIT.unittest(solver, config, [
        "solve_affine_interval",     # not implemented
        "solve_qp_edge_cases",       # not implemented
        "solve_qcp_edge_cases",      # not implemented
        "solve_objbound_edge_cases", # Requires MOI.ObjectiveBound()
    ])
    @testset "solve_affine_interval" begin
        MOIT.solve_affine_interval(
            MOIB.SplitInterval{Float64}(CPLEX.Optimizer()),
            config
        )
    end
    MOIT.modificationtest(solver, config, [
        "solve_func_scalaraffine_lessthan"
    ])
end
@testset "Linear tests" begin
    linconfig = MOIT.TestConfig()
    @testset "Default Solver"  begin
        solver = CPLEX.Optimizer()
        MOIT.contlineartest(solver, linconfig, ["linear10","linear11",
                "linear12", "linear8a","linear8b","linear8c"])
    end
    @testset "linear10" begin
        MOIT.linear10test(
            MOIB.SplitInterval{Float64}(CPLEX.Optimizer()),
            MOIT.TestConfig()
        )
    end
    @testset "No certificate" begin
        MOIT.linear12test(
            CPLEX.Optimizer(),
            MOIT.TestConfig(infeas_certificates=false)
        )
    end
end
@testset "Integer Linear tests" begin
    intconfig = MOIT.TestConfig()
    solver = CPLEX.Optimizer()
    MOIT.intlineartest(solver, intconfig, ["int1", # requires objbound impl
                                           "int2", "int3"])
    # 3 is ranged, 2 has sos
    @testset "int3" begin
        MOIT.int3test(
            MOIB.SplitInterval{Float64}(CPLEX.Optimizer()),
            MOIT.TestConfig()
        )
    end
end
@testset "ModelLike tests" begin
    solver = CPLEX.Optimizer()
    @testset "nametest" begin
        MOIT.nametest(solver)
    end
    @testset "validtest" begin
        MOIT.validtest(solver)
    end
    @testset "emptytest" begin
        MOIT.emptytest(solver)
    end
    @testset "orderedindicestest" begin
        MOIT.orderedindicestest(solver)
    end
    @testset "copytest" begin
        MOIT.copytest(solver, CPLEX.Optimizer())
    end
end


text:
# The number, 1406357289, is a 0 to 9 pandigital number because it is made up
# of each of the digits 0 to 9 in some order, but it also has a rather
# interesting sub-string divisibility property.
#
# Let d1 be the 1st digit, d2 be the 2nd digit, and so on. In this way, we note
# the following:
#
# d2d3d4 = 406 is divisible by 2
# d3d4d5 = 063 is divisible by 3
# d4d5d6 = 635 is divisible by 5
# d5d6d7 = 357 is divisible by 7
# d6d7d8 = 572 is divisible by 11
# d7d8d9 = 728 is divisible by 13
# d8d9d10 = 289 is divisible by 17
#
# Find the sum of all 0 to 9 pandigital numbers with this property.

using ProjectEulerSolutions

# Returns whether an array of digits matches the criterion.
function legal_value(vals::Array)
    return (vals[6] == 0 || vals[6] == 5) &&
           mod(vals[4], 2) == 0 &&
           mod(vals[3] + vals[4] + vals[5], 3) == 0 &&
           mod(100 * vals[5] + 10 * vals[6] + vals[7], 7) == 0 &&
           mod(100 * vals[6] + 10 * vals[7] + vals[8], 11) == 0 &&
           mod(100 * vals[7] + 10 * vals[8] + vals[9], 13) == 0 &&
           mod(100 * vals[8] + 10 * vals[9] + vals[10], 17) == 0
end

# Uses a helper function to generate permutations that all obey the specified
# "legal_value" criterion.
function p043solution(n::Integer=10)::Integer

    matrix = zeros(Integer, 0)
    permutation_matrix!(collect(0:(n-1)), matrix, legal_value)
    matrix = reshape(matrix, (n, fld(length(matrix), n)))

    return sum(transpose(10 .^ collect((n-1):-1:0)) * matrix)
end

p043 = Problems.Problem(p043solution)

Problems.benchmark(p043, 10)

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
abstract type Converter <: InverterComponent end

#TODO: converter_ports for non-fixed DC converter will need:
## VDC as state of DC source
## IDC as algebraic inner_var output from cnv real output power

@def converter_ports begin
    state_input = Vector{Symbol}()
    inner_input = [md_var, mq_var, Vdc_var, Vdcnv_var,Vqcnv_var]
end

"""
Parameters of an average converter model
# Conmutable structor
```julia
AvgCnvFixedDC(v_rated, s_rated)
```

# Arguments
*  `v_rated`::Float64 : rated voltage
*  `s_rated`::Float64 : rated VA
"""
mutable struct AvgCnvFixedDC <: Converter
    v_rated::Float64
    s_rated::Float64
    n_states::Int64
    states::Vector{Symbol}
    ports::Ports

        function AvgCnvFixedDC(v_rated::Float64,
                                s_rated::Float64)

            @converter_ports

            new(v_rated,
                s_rated,
                0,
                Vector{Symbol}(),
                Ports(state_input, inner_input))
        end
end


text:
# # Extending `@set` and `@optic`
# This code demonstrates how to extend the `@set` and `@optic` mechanism with custom
# lenses.
# As a demo, we want to implement `@mylens!` and `@myreset`, which work much like
# `@optic` and `@set`, but mutate objects instead of returning modified copies.

using Accessors
using Accessors: IndexLens, PropertyLens, ComposedOptic

struct Lens!{L}
    pure::L
end

(l::Lens!)(o) = l.pure(o)
function Accessors.set(o, l::Lens!{<: ComposedOptic}, val)
    o_inner = l.pure.inner(o)
    set(o_inner, Lens!(l.pure.outer), val)
end
function Accessors.set(o, l::Lens!{PropertyLens{prop}}, val) where {prop}
    setproperty!(o, prop, val)
    o
end
function Accessors.set(o, l::Lens!{<:IndexLens}, val) where {prop}
    o[l.pure.indices...] = val
    o
end

# Now this implements the kind of `lens` the new macros should use.
# Of course there are more variants like `Lens!(<:DynamicIndexLens)`, for which we might
# want to overload `set`, but lets ignore that. Instead we want to check, that everything works so far:

using Test
mutable struct M
    a
    b
end

o = M(1,2)
l = Lens!(@optic _.b)
set(o, l, 20)
@test o.b == 20

l = Lens!(@optic _.foo[1])
o = (foo=[1,2,3], bar=:bar)
set(o, l, 100)
@test o == (foo=[100,2,3], bar=:bar)

# Now we can implement the syntax macros

using Accessors: setmacro, opticmacro, modifymacro

macro myreset(ex)
    setmacro(Lens!, ex)
end

macro mylens!(ex)
    opticmacro(Lens!, ex)
end

macro mymodify!(f, ex)
    modifymacro(Lens!, f, ex)
end

o = M(1,2)
@myreset o.a = :hi
@myreset o.b += 98
@test o.a == :hi
@test o.b == 100

o = M(1,3)
@mymodify!(x -> x+1, o.a)
@test o.a === 2
@test o.b === 3

deep = [[[[1]]]]
@myreset deep[1][1][1][1] = 2
@test deep[1][1][1][1] === 2

l = @mylens! _.foo[1]
o = (foo=[1,2,3], bar=:bar)
set(o, l, 100)
@test o == (foo=[100,2,3], bar=:bar)

# Everything works, we can do arbitrary nesting and also use `+=` syntax etc.


text:
function SimpleMotor(a::Body, b::Body, rate::Real)
    SimpleMotor(ccall(dlsym(libchipmunk, :cpSimpleMotorNew), Ptr{Void}, (Ptr{Void}, Ptr{Void}, Cdouble,), a.ptr, b.ptr, rate))
end

function get_rate(motor::SimpleMotor)
    ccall(dlsym(libchipmunk, :cpSimpleMotorGetRate), Cdouble, (Ptr{Void},), motor.ptr)
end

function set_rate(motor::SimpleMotor, rate::Real)
    ccall(dlsym(libchipmunk, :cpSimpleMotorSetRate), Void, (Ptr{Void}, Cdouble,), motor.ptr, rate)
end

export SimpleMotor, get_rate, set_rate


text:
using FFTW:fft
export MelFilters

mutable struct MelFilters
    winlen::Int   # 分帧参数-帧长
    stride::Int   # 分帧参数-帧移
    nbanks::Int   # 梅尔滤波器个数
    nffts ::Int   # 傅立叶变换点数
    alpha::Real   # 预加重系数，默认0.97
    fs::Int       # 采样率，一般16kHz
    maxfreq::Int  # 频域最大有效频率下标
    eps::Real     # 防止下溢系数
    winfunc::Function       # 窗计算函数，如 hanning
    window ::AbstractArray  # 窗函数数值表示
    melbank::AbstractArray  # mel滤波器矩阵

    function MelFilters(;
        winlen::Int = 256,
        stride::Int = winlen>>1,
        nffts ::Int = winlen,
        nbanks::Int = 32,
        alpha::Real = 0.97,
        fs::Int     = 16000,
        eps::Real   = 1e-6,
        winfunc::Function = hanning)
        @assert nffts>=winlen "nffts ≥ winlen, but got nffts=$nffts, winlen=$winlen"
        maxfreq = nffts>>1
        window  = winfunc(winlen)
        melbank = filterbanks(nbanks, nffts, fs)
        new(winlen,stride, nbanks,nffts,alpha,fs, maxfreq,eps, winfunc, window,melbank)
    end
end

function filterbanks(nbanks::Int, nffts::Int, fs::Int)
    MAX   = nffts>>1;                    # 正频率部分的下标最大值
    freq  = (0:(MAX-1))/MAX * fs/2;      # 下标映射到频率
    Fmel  = 2595*log10.(1 .+ freq/700);  # 频率映射到梅尔频率
    dFmel = Fmel[MAX]/(nbanks+1);        # 将Mel带平分
    bank  = zeros(nbanks, MAX);          # 滤波器的频域权重系数
    cFmel = 0.0;                         # 每个Mel频带的中心Mel频率
    for n = 1:nbanks
        cFmel = cFmel + dFmel
        for m = 1:MAX
            if ( Fmel[m] >= cFmel-dFmel ) && ( Fmel[m] <= cFmel+dFmel )
                bank[n,m] = 1.0 - abs( Fmel[m] - cFmel )/dFmel
            end
        end
    end
    return bank
end

function (filter::MelFilters)(wav, func::Union{Function,Nothing}=log)
    W = filter.melbank
    α = filter.alpha
    B = filter.eps
    F = filter.maxfreq
    winlen = filter.winlen
    stride = filter.stride
    window = filter.window
    nffts  = filter.nffts

    d = filterwav(wav, α)                            # 滤波
    X = spectrum(d, window, winlen, stride, nffts)   # 功率谱

    if func !== nothing
        return func.(W * X .+ B)
    else
        return W * X
    end
end

function Base.show(io::IO, f::MelFilters)
    winlen = trunc(f.winlen / f.fs * 1000, digits=3)
    stride = trunc(f.stride / f.fs * 1000, digits=3)
    println(io, "———————————————————————")
    println(io, " winfunc = $(f.winfunc)")
    println(io, "  winlen = $(f.winlen) ($winlen ms)")
    println(io, "  stride = $(f.stride) ($stride ms)")
    println(io, "   nffts = $(f.nffts)")
    println(io, " maxfreq = $(f.maxfreq)")
    println(io, "  nbanks = $(f.nbanks)")
    println(io, "   alpha = $(f.alpha)")
    println(io, "     eps = $(f.eps)")
    println(io, "      fs = $(f.fs) Hz")
    println(io, "———————————————————————")
end


text:
using QuantumWalk
using LightGraphs

importall QuantumWalk
##

abstract type AbstractStochastic <: QWModelDiscr end

struct UniformStochastic{G<:SimpleGraph} <: AbstractStochastic
  graph::G
end

UniformScaling(digraph::G) where G= UniformStochastic{G}(digraph)

function check_qwdynamics(::Type{QWSearch},
                          abs_stoch::UniformStochastic,
                          ::Vector{Int},
                          parameters::Dict{Symbol,Any})
  @assert :stochastic ∈ keys(parameters) "parameters needs to have key stochastic"
  n = nv(graph(abs_stoch))
  @assert isa(parameters[:stochastic], SparseMatrixCSC{<:Real}) "value for :stochastic needs to be sparse matrix with real numbers"
  @assert size(parameters[:stochastic], 1) == size(parameters[:stochastic], 2) "Stochastic matrix needs to be square stochastic matrix"
  @assert mapslices(sum, parameters[:stochastic], 1)[1,:] ≈ ones(n) "Stochastic matrix needs to be square stochastic matrix of order graph"
end

function QWSearch(stoch::AbstractStochastic,
                  marked::Vector{Int},
                  penalty::Real = 0.)
   parameters = Dict{Symbol,Any}(:stochastic => stochastic_matrix(graph(stoch)))

   QWSearch(stoch, marked, parameters, penalty)
end

function initial_state(qws::QWSearch{<:AbstractStochastic})
  n = nv(graph(qws))
  fill(1./n, n)
end

function stochastic_matrix(g::SimpleGraph)
  a = adjacency_matrix(g)
  a*spdiagm(mapslices(x->1/sum(x), a, 1)[1,:])
end

function stochastic_evolution(s::SparseMatrixCSC{T}, v::Vector{T}) where T<:Real
  s*v
end

function evolve(qws::QWSearch{<:AbstractStochastic}, state::Vector{<:Real})
  old_probability = measure(qws, state, marked(qws))
  state[marked(qws)] = zero(marked(qws))
  state = stochastic_evolution(parameters(qws)[:stochastic], state)
  state[marked(qws)] += old_probability
  state
end

function measure(::QWDynamics{<:AbstractStochastic}, state::Vector{<:Real})
   return state
end

function measure(::QWDynamics{<:AbstractStochastic},
                 state::Vector{<:Real},
                 vertices::Vector{Int})
   return state[vertices]
end

##

n = 100
penalty_stoch = log(n)
dynamic = QWSearch(UniformStochastic(CompleteGraph(n)), [1], penalty_stoch)

execute_single(dynamic, 40)


text:
using Test
using RoME

@testset "Test g2o file parsing." begin

  instructions = importG2o("octagon.g2o")
  # Choose random fields to test.
  @test instructions[1][1] == "EDGE_SE2"
  @test instructions[7][1] == "EDGE_SE2"
  @test instructions[6][12] == "6541.252776"
  @test instructions[3][7] == "1211.201664"
  # Input file has 8 lines, with 12 columns.
  @test length(instructions) == 8
  @test length(instructions[2]) == 12

  # build the factor graph from these instructions

  # generate g2o from factor graph and should match orginal file
end

@testset "Test g2o file parsing." begin

  reflines = ["LANDMARK 0 1 0.0 20.0 3.162277660168379 0.0 1.0";
              "EDGE_SE2 0 2 10.0 0.0 1.0471975511965976 10.0 10.0 0.0 0.0 10.0 0.0 0.0 0.0";
              "EDGE_SE2 2 3 10.0 0.0 1.0471975511965976 10.0 10.0 0.0 0.0 10.0 0.0 0.0 0.0";
              "EDGE_SE2 3 4 10.0 0.0 1.0471975511965976 10.0 10.0 0.0 0.0 10.0 0.0 0.0 0.0";
              "EDGE_SE2 4 5 10.0 0.0 1.0471975511965976 10.0 10.0 0.0 0.0 10.0 0.0 0.0 0.0";
              "EDGE_SE2 5 6 10.0 0.0 1.0471975511965976 10.0 10.0 0.0 0.0 10.0 0.0 0.0 0.0";
              "EDGE_SE2 6 7 10.0 0.0 1.0471975511965976 10.0 10.0 0.0 0.0 10.0 0.0 0.0 0.0";
              "LANDMARK 7 1 0.0 20.0 3.162277660168379 0.0 1.0"]

  # test file generation
  fg = generateCanonicalFG_Hexagonal(graphinit=false)

  filepath = exportG2o(fg)

  fid = open(filepath, "r")
  count = 0
  for line in eachline(fid)
    count += 1
    @test reflines[count] == line
  end
  close(fid)

end


text:
donothing(x...) = return
# nop(x) = x
returntrue(x) = true

# type ObjectWrapper
#   obj
# end
# wrap(obj...) = ObjectWrapper(obj)
# unwrap(wrapper::ObjectWrapper) = (obj = wrapper.obj; n = length(obj); (n == 1 ? obj[1] : (n == 0 ? nothing : obj)))


text:
using Flux, DiffEqFlux, OrdinaryDiffEq, Optim

#A desired MWE for now, not a test yet.

function f(du,u,p,t)
    y₁,y₂,y₃ = u
    k₁,k₂,k₃ = p
    du[1] = -k₁*y₁ + k₃*y₂*y₃
    du[2] =  k₁*y₁ - k₃*y₂*y₃ - k₂*y₂^2
    du[3] =  y₁ + y₂ + y₃ - 1
    nothing
end

u₀ = [1.0, 0, 0]
M = [1. 0  0
    0  1. 0
    0  0  0]
tspan = (0.0,10.0)
p = [0.04,3e7,1e4]
func = ODEFunction(f,mass_matrix=M)
prob = ODEProblem(f,u₀,tspan,(0.04,3e7,1e4))
sol = solve(prob,Rodas5())

dudt2 = Chain(x -> x.^3,Dense(3,50,tanh),Dense(50,2))

ndae = NeuralDAE(dudt2, (u,p,t) -> [u[1] + u[2] + u[3] - 1], tspan, M, Rodas5())

ndae(u₀)

function predict_n_dae(p)
    ndae(u₀,p)
end

function loss(p)
    pred = predict_n_dae(p)
    loss = sum(abs2,sol .- pred)
    loss,pred
end

p = p .+ rand(3) .* p

res = DiffEqFlux.sciml_train(loss, p, BFGS(initial_stepnorm = 0.0001))


text:
# random functions that dispatch either to CURAND or GPUArrays' generic RNG

using Random

export rand_logn!, rand_poisson!

# native RNG

"""
    CUDA.RNG()

A random number generator using `rand()` in a device kernel.

!!! warning

    This generator is experimental, and not used by default yet. Currently, it does not
    pass SmallCrush (but neither does the current GPUArrays.jl-based generator). If you
    can, use a CURAND-based generator for now.

See also: [SharedTauswortheGenerator](@ref)
"""
struct RNG <: AbstractRNG
    state::CuVector{UInt32}

    function RNG(seed)
        @assert length(seed) == 32
        new(seed)
    end
end

RNG() = RNG(rand(UInt32, 32))

function Random.seed!(rng::RNG, seed::AbstractVector{UInt32})
    @assert length(seed) == 32
    copyto!(rng.state, seed)
end

function Random.seed!(rng::RNG)
    Random.rand!(rng.state)
    return
end

function Random.rand!(rng::RNG, A::AnyCuArray)
    function kernel(A::AbstractArray{T}, state::AbstractVector{UInt32}) where {T}
        device_rng = Random.default_rng()

        # initialize the state
        tid = threadIdx().x
        if tid <= 32
            # we know for sure the seed will contain 32 values
            @inbounds Random.seed!(device_rng, state)
        end

        sync_threads()

        # grid-stride loop
        offset = (blockIdx().x - 1) * blockDim().x
        while offset < length(A)
            # generating random numbers synchronizes threads, so needs to happen uniformly.
            val = Random.rand(device_rng, T)

            i = tid + offset
            if i <= length(A)
                @inbounds A[i] = val
            end

            offset += (blockDim().x - 1) * gridDim().x
        end

        sync_threads()

        # save the device rng state of the first block (other blocks are derived from it)
        # so that subsequent launches generate different random numbers
        if blockIdx().x == 1 && tid <= 32
            @inbounds state[tid] = device_rng.state[tid]
        end

        return
    end

    kernel = @cuda launch=false name="rand!" kernel(A, rng.state)
    config = launch_configuration(kernel.fun; max_threads=64)
    threads = max(32, min(config.threads, length(A)))
    blocks = min(config.blocks, cld(length(A), threads))
    kernel(A, rng.state; threads=threads, blocks=blocks)

    # XXX: updating the state from within the kernel is racey,
    #      so we should have a RNG per stream

    A
end

# TODO: `randn!`; cannot reuse from Base or RandomNumbers, as those do scalar indexing

# generic functionality

function Random.rand!(rng::Union{RNG,CURAND.RNG,GPUArrays.RNG}, A::AbstractArray{T}) where {T}
    B = CuArray{T}(undef, size(A))
    Random.rand!(rng, B)
    copyto!(A, B)
end

function Random.rand(rng::Union{RNG,CURAND.RNG,GPUArrays.RNG}, T::Type)
    assertscalar("scalar rand")
    A = CuArray{T}(undef, 1)
    Random.rand!(rng, A)
    A[]
end

# RNG-less interface

# the interface is split in two levels:
# - functions that extend the Random standard library, and take an RNG as first argument,
#   will only ever dispatch to CURAND and as a result are limited in the types they support.
# - functions that take an array will dispatch to either CURAND or GPUArrays
# - non-exported functions are provided for constructing GPU arrays from only an eltype

const curand_rng = CURAND.default_rng
gpuarrays_rng() = GPUArrays.default_rng(CuArray)

function seed!(seed=Base.rand(UInt64))
    Random.seed!(curand_rng(), seed)
    Random.seed!(gpuarrays_rng(), seed)
end

# CURAND in-place
Random.rand!(A::CURAND.UniformArray) = Random.rand!(curand_rng(), A)
Random.randn!(A::CURAND.NormalArray; kwargs...) = Random.randn!(curand_rng(), A; kwargs...)
rand_logn!(A::CURAND.LognormalArray; kwargs...) = CURAND.rand_logn!(curand_rng(), A; kwargs...)
rand_poisson!(A::CURAND.PoissonArray; kwargs...) = CURAND.rand_poisson!(curand_rng(), A; kwargs...)

# CURAND out-of-place
rand(T::CURAND.UniformType, dims::Dims) = Random.rand(curand_rng(), T, dims)
randn(T::CURAND.NormalType, dims::Dims; kwargs...) = Random.randn(curand_rng(), T, dims; kwargs...)
rand_logn(T::CURAND.LognormalType, dims::Dims; kwargs...) = CURAND.rand_logn(curand_rng(), T, dims; kwargs...)
rand_poisson(T::CURAND.PoissonType, dims::Dims; kwargs...) = CURAND.rand_poisson(curand_rng(), T, dims; kwargs...)

# support all dimension specifications
rand(T::CURAND.UniformType, dim1::Integer, dims::Integer...) =
    Random.rand(curand_rng(), T, Dims((dim1, dims...)))
randn(T::CURAND.NormalType, dim1::Integer, dims::Integer...; kwargs...) =
    Random.randn(curand_rng(), T, Dims((dim1, dims...)); kwargs...)
rand_logn(T::CURAND.LognormalType, dim1::Integer, dims::Integer...; kwargs...) =
    CURAND.rand_logn(curand_rng(), T, Dims((dim1, dims...)); kwargs...)
rand_poisson(T::CURAND.PoissonType, dim1::Integer, dims::Integer...; kwargs...) =
    CURAND.rand_poisson(curand_rng(), T, Dims((dim1, dims...)); kwargs...)

# GPUArrays in-place
Random.rand!(A::AnyCuArray) = Random.rand!(gpuarrays_rng(), A)
Random.randn!(A::AnyCuArray; kwargs...) =
    error("CUDA.jl does not support generating normally-distributed random numbers of type $(eltype(A))")
# FIXME: GPUArrays.jl has a randn! nowadays, but it doesn't work with e.g. Cuint
rand_logn!(A::AnyCuArray; kwargs...) =
    error("CUDA.jl does not support generating lognormally-distributed random numbers of type $(eltype(A))")
rand_poisson!(A::AnyCuArray; kwargs...) =
    error("CUDA.jl does not support generating Poisson-distributed random numbers of type $(eltype(A))")

# GPUArrays out-of-place
rand(T::Type, dims::Dims) = Random.rand!(CuArray{T}(undef, dims...))
randn(T::Type, dims::Dims; kwargs...) = Random.randn!(CuArray{T}(undef, dims...); kwargs...)
rand_logn(T::Type, dims::Dims; kwargs...) = rand_logn!(CuArray{T}(undef, dims...); kwargs...)
rand_poisson(T::Type, dims::Dims; kwargs...) = rand_poisson!(CuArray{T}(undef, dims...); kwargs...)

# support all dimension specifications
rand(T::Type, dim1::Integer, dims::Integer...) =
    Random.rand!(CuArray{T}(undef, dim1, dims...))
randn(T::Type, dim1::Integer, dims::Integer...; kwargs...) =
    Random.randn!(CuArray{T}(undef, dim1, dims...); kwargs...)
rand_logn(T::Type, dim1::Integer, dims::Integer...; kwargs...) =
    rand_logn!(CuArray{T}(undef, dim1, dims...); kwargs...)
rand_poisson(T::Type, dim1::Integer, dims::Integer...; kwargs...) =
    rand_poisson!(CuArray{T}(undef, dim1, dims...); kwargs...)

# untyped out-of-place
rand(dim1::Integer, dims::Integer...) =
    Random.rand(curand_rng(), Dims((dim1, dims...)))
randn(dim1::Integer, dims::Integer...; kwargs...) =
    Random.randn(curand_rng(), Dims((dim1, dims...)); kwargs...)
rand_logn(dim1::Integer, dims::Integer...; kwargs...) =
    CURAND.rand_logn(curand_rng(), Dims((dim1, dims...)); kwargs...)
rand_poisson(dim1::Integer, dims::Integer...; kwargs...) =
    CURAND.rand_poisson(curand_rng(), Dims((dim1, dims...)); kwargs...)


text:
######################
### Mutable states ###
######################

mutable struct StanHMCAdaptorState
    i               ::  Int
    window_start    ::  Int         # start of mass matrix adaptation
    window_end      ::  Int         #   end of mass matrix adaptation
    window_splits   ::  Vector{Int} # iterations to update metric using mass matrix
end

StanHMCAdaptorState() = StanHMCAdaptorState(0, 0, 0, Vector{Int}(undef, 0))

# Ref: https://github.com/stan-dev/stan/blob/develop/src/stan/mcmc/windowed_adaptation.hpp
function initialize!(state::StanHMCAdaptorState, init_buffer::Int, term_buffer::Int, window_size::Int, n_adapts::Int)
    window_start = init_buffer + 1
    window_end = n_adapts - term_buffer

    # Update window-end points
    window_splits = Vector{Int}(undef, 0)
    next_window = init_buffer + window_size
    while next_window <= window_end
        # Extend the current window to the end of the full window
        # if the next window reaches the end of the full window
        next_window_boundary = next_window + 2 * window_size
        if next_window_boundary > window_end
            next_window = window_end
        end
        # Include the window split
        push!(window_splits, next_window)
        # Expand window and compute the next split
        window_size *= 2
        next_window += window_size
    end
    # Avoid updating in the end
    if window_splits[end] == n_adapts
        pop!(window_splits)
    end

    state.window_start  = window_start
    state.window_end    = window_end
    state.window_splits = window_splits
end

function Base.show(io::IO, state::StanHMCAdaptorState)
    print(io, "window($(state.window_start), $(state.window_end)), window_splits(" * string(join(state.window_splits, ", ")) * ")")
end

################
### Adaptors ###
################

# Acknowledgement: this adaption settings is mimicing Stan's 3-phase adaptation.
struct StanHMCAdaptor{M<:AbstractPreconditioner, Tssa<:StepSizeAdaptor} <: AbstractAdaptor
    pc          :: M
    ssa         :: Tssa
    init_buffer :: Int
    term_buffer :: Int
    window_size :: Int
    state       :: StanHMCAdaptorState
end
Base.show(io::IO, a::StanHMCAdaptor) =
    print(io, "StanHMCAdaptor(\n    pc=$(a.pc),\n    ssa=$(a.ssa),\n    init_buffer=$(a.init_buffer), term_buffer=$(a.term_buffer), window_size=$(a.window_size),\n    state=$(a.state)\n)")

function StanHMCAdaptor(
    pc::M,
    ssa::StepSizeAdaptor;
    init_buffer::Int=75,
    term_buffer::Int=50,
    window_size::Int=25
) where {M<:AbstractPreconditioner}
    return StanHMCAdaptor(pc, ssa, init_buffer, term_buffer, window_size, StanHMCAdaptorState())
end

getM⁻¹(adaptor::StanHMCAdaptor) = getM⁻¹(adaptor.pc)
getϵ(adaptor::StanHMCAdaptor)   = getϵ(adaptor.ssa)
function initialize!(adaptor::StanHMCAdaptor, n_adapts::Int)
    initialize!(adaptor.state, adaptor.init_buffer, adaptor.term_buffer, adaptor.window_size, n_adapts)
    return adaptor
end
finalize!(adaptor::StanHMCAdaptor) = finalize!(adaptor.ssa)

is_in_window(a::StanHMCAdaptor) = a.state.i >= a.state.window_start && a.state.i <= a.state.window_end
is_window_end(a::StanHMCAdaptor) = a.state.i in a.state.window_splits

function adapt!(
    tp::StanHMCAdaptor,
    θ::AbstractVecOrMat{<:AbstractFloat},
    α::AbstractScalarOrVec{<:AbstractFloat}
)
    tp.state.i += 1

    adapt!(tp.ssa, θ, α)

    resize!(tp.pc, θ) # Resize pre-conditioner if necessary.

    # Ref: https://github.com/stan-dev/stan/blob/develop/src/stan/mcmc/hmc/nuts/adapt_diag_e_nuts.hpp
    if is_in_window(tp)
        # We accumlate stats from θ online and only trigger the update of M⁻¹ in the end of window.
        is_update_M⁻¹ = is_window_end(tp)
        adapt!(tp.pc, θ, α, is_update_M⁻¹)
    end

    if is_window_end(tp)
        reset!(tp.ssa)
        reset!(tp.pc)
    end
end

# Deprecated functions

@deprecate StanHMCAdaptor(n_adapts, pc, ssa) initialize!(StanHMCAdaptor(pc, ssa), n_adapts)


text:
function loss_A(x, par)
    @unpack_Pars_Simple par
    k_I_b = x
    k_I, l = capital_and_labor_unc_A(k_I_b, par)
    # k_I = total_capital_A(k_I_b, k_I_P, par)
    # l = labor_demand_A(k_I, par)
    # Residual
    return loss_A_only(k_I, l, par)
end

"""
    A is selling, such that selling an additional unit yields p_I, but costs only
    (1-ϕ) units of the marginal product of capital, as ϕ units are kept.
"""
function loss_A_only(k_I, l, par)
    @unpack_Pars_Simple par
    p_I - (1 - ϕ) * γ * δ_hat * v * C^α_C * l^α_hat * k_I^(δ_hat - 1)
end

"""
    B is buying, so buying costs p_I/(1-τ) and yields the full marginal product.
"""
function loss_B_only(k_I, l, par)
    @unpack_Pars_Simple par
    p_I / (1 - τ) - (1 - γ) * δ_hat * v * C^α_C * l^α_hat * k_I^(δ_hat - 1)
end

function loss_A_con(k_I, par)
    @unpack_Pars_Simple par
    l = labor_demand_A(k_I, par)

    return mc_A(k_I, par) - p_I - ϕ^δ_hat * γ * δ_hat * v * C^α_C * l^α_hat * k_I^(δ_hat - 1)
end

function loss_B(x, par)
    @unpack_Pars_Simple par
    k_I_b = x
    # k_I = total_capital_B(k_I_b, k_I_P, par)
    # l = labor_demand_B(k_I, par)
    k_I, l = capital_and_labor_B(k_I_b, par)
    # Residual
    # very important: B firms face a higher marginal cost of investment!
    loss_B_only(k_I, l, par)
end

"""
    Find the equilibrium for a given C using Optim.jl
    Advantage: For some p_I, there might be no solution to the firm's problem.
    However, it might still be that there is a solution for a different p_I
    without resorting to the constrained solution. Only if this solution is
    unavailable, should I try the constrained solution.
"""
function find_price_info_optim(par)

    lower = [1.25, -2.1]
    upper = [Inf, Inf]
    initial_x = [2.0, 2.0]
    f(x) = residuals(x, par)
    inner_optimizer = GradientDescent()
    results = optimize(f, g!, lower, upper, initial_x, Fminbox(inner_optimizer))

    #! Problem, bounds depend on one parameter
end

function residuals( x, par)
        @unpack_Pars_Simple par
        p_I = x[1]
        @set! par.p_I = p_I
        k_I_b_A = x[2]
        k_I_b_B = x[3]
        # How much capital is produced
        k_I_P_A = k_I_prod_A(par)
        #k_I_P_B = k_I_prod_B(par)

        k_I_A, l_A = capital_and_labor_unc_A(k_I_b_A, par)
        k_I_B, l_B = capital_and_labor_unc_B(k_I_b_B, par)

        res_1 = loss_A_only(k_I_A, l_A, par)
        res_2 = loss_B_only(k_I_B, l_B, par)
        res_3 = market_clearing(k_I_b_B, k_I_b_A, par) # market clearing!
        res_4 = capicity_constraint(k_I_b_A, k_I_P_A) # This is the constraint that A cannot sell more than it produces!

        return res_1^2 + res_2^2 + res_3^2 + res_4^2
end

#! Try the following: Write your own solution algorithm.
# 1. Guess p_I
# 2. Check whether a sign change occurs
# 2.1. If yes, good, find solution.
# 2.2., If no, try lower p_I.
# 3. If there is no sign change even with some minimal p_I, switch to the constraint problem.

#! Can the constraint problem ever be a solution?

function capicity_constraint(k_I_b_A, k_I_P_A)
    # note that k_I_b_A is negative
    if k_I_P_A + k_I_b_A < 0
        1000*(k_I_P_A + k_I_b_A)^2
    else
        0
    end
end

function solve_agg_C(par)
    f(C) = begin
        @set! par.C = C^2
        return C - find_agg_C(par)
    end
    find_zero(f, 0.05)^2
end

function find_agg_C(par)
    @set! par.p_I = find_price_info(par)
    k_I_b_A, k_I_b_B = solve_firm_problems(par)

    if par.ϕ > 0 && loss_A(k_I_b_A, par) > 1e-8
        k_I_A, l_A = capital_and_labor_con_A(k_I_b_A, par)
    else
        k_I_A, l_A = capital_and_labor_unc_A(k_I_b_A, par)
    end
    k_I_B, l_B = capital_and_labor_B(k_I_b_B, par)

    C_A = firm_prod(l_A, k_I_A; par)
    C_B = firm_prod(l_B, k_I_B; par)

    return agg_prod(C_A, C_B, par), wage_rate(k_I_A, k_I_B; par)
end

function find_wage_rate(par)
    k_I_A, k_I_B = capital_prod(par)
    wage_rate(k_I_A, k_I_B; par)
end
function agg_prod(C_A, C_B, par)
    @unpack_Pars_Simple par
    (γ * C_A^((σ - 1) / σ) + (1 - γ) * C_B^((σ - 1) / σ))^(v * σ / (σ - 1))
end
function firm_prod(l, k_I; par)
    @unpack_Pars_Simple par
    l^α_hat * k_I^δ_hat
end

function solve_firm_problems(par)
    k_I_P_A = k_I_prod_A(par)
    if loss_A(-k_I_P_A, par) ≈ 0
        k_I_b_A = -k_I_P_A
    elseif par.ϕ > 0 && loss_A(-k_I_P_A, par) > 0
        #println("Constrained Optimization")
        # constrained, sell everything that is produced
        f_A_con(x) = loss_A_con(x^2, par)
        k_I_A = find_zero(f_A_con, 1)^2 # this is the production decision - always positive!
        k_I_P_A = k_I_A / (1 - par.ϕ) # more is produced than employed
        k_I_b_A = -k_I_P_A # production equals sharing
    else # unconstrained, sell less than what is produced
        k_I_P_A = k_I_prod_A(par)
        f_A(x) = loss_A(x, par)
        @show f_A(-k_I_P_A)
        @show f_A(100)
        @show par.p_I
        k_I_b_A = find_zero(f_A, [-k_I_P_A, 1000])
    end

    # B is always unconstrained, because they are buying capital
    k_I_P_B = k_I_prod_B(par)
    f_B(x) = loss_B(x, par)
    k_I_b_B = find_zero(f_B, [-k_I_P_B + 0.0001, 20])

    return k_I_b_A, k_I_b_B
end

function capital_prod(par)
    k_I_b_A, k_I_b_B = solve_firm_problems(par)
    if par.ϕ > 0 && loss_A(k_I_b_A, par) > 1e-8
        println("Constrained Optimization")
        k_I_A, l_A = capital_and_labor_con_A(k_I_b_A, par)
    else
        println("Unconstrained Optimization")
        k_I_A, l_A = capital_and_labor_unc_A(k_I_b_A, par)
    end
    k_I_B, l_B = capital_and_labor_B(k_I_b_B, par)

    return k_I_A, k_I_B
end

"""
    This is much more robust!
"""
function solve_agg_C_iterative(par; tol = 1e-10)
    diff = 10
    C_guess = par.C
    w_guess = par.w
    while diff > tol
        @set! par.C = C_guess
        #@show par.C
        C_sol, w_sol = find_agg_C(par)
        @set! par.w = w_sol
        diff = copy((C_sol - C_guess)^2 + (w_sol - w_guess)^2)
        if diff > tol
            C_guess = copy(C_sol)
            w_guess = copy(w_sol)
        else
            @set! par.w = w_sol
            @set! par.C = C_guess
            @set! par.p_I = find_price_info(par)
            return par
        end
    end
end

capital_of_tau(τ) = capital_prod(solve_agg_C_iterative(Pars_Simple(τ = τ)))

function capital_and_labor_unc_A(k_I_b_A, par)
    #! Cannot put this here, as I otherwise call loss_A recursively
    # if par.ϕ > 0 && loss_A(k_I_b_A, par) > 0
    #     # constrained, sell everything that is produced
    #     f_A_con(x) = loss_A_con(x^2, par)
    #     k_I_A = find_zero(f_A_con, 1)^2 # this is the production decision - always positive!
    # else
    # This is for an unconstrained firm
    k_I_P_A = k_I_prod_A(par)
    k_I_A = total_capital_A(k_I_b_A, k_I_P_A, par)
    # end
    l_A = labor_demand_A(k_I_A, par)
    return k_I_A, l_A
end

function capital_and_labor_con_A(k_I_b_A, par)
    #! Cannot put this here, as I otherwise call loss_A recursively
    @assert par.ϕ > 0
    # if par.ϕ > 0 && loss_A(k_I_b_A, par) > 0
    #     # constrained, sell everything that is produced
    f_A_con(x) = loss_A_con(x^2, par)
    k_I_A = find_zero(f_A_con, 1)^2 # this is the production decision - always positive!
    # else
    # This is for an unconstrained firm
    # k_I_P_A = k_I_prod_A(par)
    # k_I_A = total_capital_A(k_I_b_A, k_I_P_A, par)
    # # end
    l_A = labor_demand_A(k_I_A, par)
    return k_I_A, l_A
end

function capital_and_labor_B(k_I_b_B, par)
    k_I_P_B = k_I_prod_B(par)
    k_I_B = total_capital_B(k_I_b_B, k_I_P_B)
    l_B = labor_demand_B(k_I_B, par)
    return k_I_B, l_B
end

# function profit_no_effort(par)
#     @unpack_Pars_Simple par
#     sol_A, sol_B = solve_firm_problems(par)

#     l_B = sol_B[1]
#     k_I_b_B = sol_B[2]
#     k_I_P_B = k_I_prod_B(par)
#     k_I_B = k_I_b_B + k_I_P_B

#     l_A = sol_A[1]
#     k_I_b_A = sol_A[2]
#     k_I_P_A = k_I_prod_A(par)
#     k_I_A = k_I_b_A + k_I_P_A

#     profit_A = γ * v * C^α_C * l_A^α_hat * k_I_A^δ_hat - w * l_A - p_I * k_I_b_A - c_A(k_I_P_A, par)
#     profit_B = (1 - γ) * v * C^α_C * l_B^α_hat * k_I_B^δ_hat - w * l_B - p_I * k_I_b_B/(1-τ) - c_B(k_I_P_B, par)

#     return profit_A, profit_B
# end


text:
#

# The simulation gets some input from other simulations
# Every x step ?
#
# DataFrame
# Time
# susceptible
# Infected (?)
# Symptomatic
# Hospitalized
# Dead
# Recovered
#
# only SEII2R have compliance because HD are fully-compliant
# E is only knowable from "Nowcasting" or similar things
#

# Every agent has an agent-state from SEIIHDR, but not all are knowable for
# other agents

# From the delta between time stamps recreate the states in our agents.

using Agents

include("../agents/agent.jl")
include("../spaces/network.jl")
include("../dataimport/modelproperties.jl")
include("stepping.jl")
include("evaluation.jl")
include("../dataimport/loadmodel.jl")

# to allow a unique signature of the model configurtion
# the getAgentType function will determine the agent type
# from a configuration file.

mprops = getModelProperties()

function getAgentType()
    return SimplePolicyAgent
end

# same as above
function getSpace()
        props = Dict("type" => "barabasi",
                     "n" => 100,
                     "k" => 2,
                     "n0" => 10
        )
    generateSpace(props)
end

function getSpace()
        props = Dict("type" => "watts",
                     "n" => 100,
                     "k" => 2,
                     "beta" => 0.5
        )
    generateSpace(props)
end

# method for running
function runsimulation(abm)

 run!(abm, agents_step!, 100;
    adata =  a_data_props,
    mdata = m_data_props,
    replicates = 10,
    parallel = true)

end


text:
#!/usr/bin/env julia

using Distributions
using FileIO
using JLD2
using Random
using SpecialFunctions

"""
  gendat(nu::Int, t::Int)

Returns Matrix{Float64} - t realisations from t-student multivatiate distribution
with nu degress of freedom
"""

function gendat(nu::Int, t::Int = 150000000)
  cm = [[1. 0.7 0.7 0.7];[0.7 1. 0.7 0.7]; [0.7 0.7 1. 0.7]; [0.7 0.7 0.7 1]]
  p = MvTDist(nu, [0., 0., 0., 0.],cm)
  return Array(transpose(rand(p, t)))
end

"""
  tmom(nu::Int, k::Int)

Returns Float64, the k'th moment of standard t distribution with nu degreed of
freedom
"""
tmom(nu::Int, k::Int) = gamma((k+1)/2)*gamma((nu-k)/2)*nu^(k/2)/(sqrt(pi)*gamma(nu/2))
"""
  tcum(nu::Int, k::Int)

Returns Float64, the k'th cumulant of standard t distribution with nu degreed of
freedom
"""
function tcum(nu::Int, k::Int)
  if k in (1,3,5)
    return 0.
  elseif k == 2
    return tmom(nu, 2)
  elseif k == 4
    return tmom(nu, 4) - 3*tmom(nu, 2)^2
  elseif k == 6
    return tmom(nu, 6) - 15*tmom(nu, 4)*tmom(nu, 2) + 30*tmom(nu, 2)^3
  end
end

function main()
  nu = 14
  Random.seed!(42)
  data = gendat(nu::Int)
  d = Dict{String, Any}("theoretical diag" => Float64[tcum(14, k) for k in 1:6])
  push!(d, "data" => data)
  save("data/datafortests.jld2", d)
end

main()


text:
#=============================#
#== put needed modules here ==#
#=============================#
import JuliaChem

#================================#
#== JuliaChem execution script ==#
#================================#
function rhf_properties_xyz(input_file)
  try
    #== read in xyz input file ==#
    mol_charge = 0
    molecule = JuliaChem.JCInput.xyz_to_molecule(input_file, mol_charge)

    #== define other input parameters ==#
    driver = "energy"

    model = Dict(
      "method" => "RHF",
      "basis" => "6-31G*"
    )

    keywords = Dict(
      "scf" => Dict(),
      "prop" => Dict(
        "formation" => true,
        "mo energies" => true,
        "mulliken" => true,
        "multipole" => "dipole"
      )
    )

    #== generate basis set ==#
    mol, basis = JuliaChem.JCBasis.run(molecule, model;
      output=0)

    #JuliaChem.JCMolecule.run(mol)

    #== perform scf calculation ==#
    rhf_energy = JuliaChem.JCRHF.Energy.run(mol, basis, keywords["scf"];
      output=1)

    #== perform property calculations ==#
    rhf_properties = JuliaChem.JCRHF.Properties.run(mol, basis, rhf_energy,
        keywords["prop"]; output=2)

    return rhf_energy, rhf_properties
  catch e
    bt = catch_backtrace()
    msg = sprint(showerror, e, bt)
    println(msg)

    JuliaChem.finalize()
    exit()
  end
end


text:
export SCIPSolver

# dummy implementation to tell users about transition to MathOptInterface
function SCIPSolver(args...; kwargs...)
    error("Support for MathProgBase was dropped. " *
          "Please downgrade to v0.6.1, see README for details.")
end


text:
module CalculoNumerico

@doc raw"
	newton_raphson

Esta função implementa uma básica versão do método de Newton-Raphson. Essencialmente, esta implementação depende de um ponto inicial ``x_0``, uma função $f$ e sua derivada $f'$. Estabelecida uma precisão $ε$, o método de Newton-Raphson gera uma sequência iterativa
	$$x_{k+1} = x_{k} -frac{f(x_k)}{f'(x_k)},$$
que, sob certas condições converge para solução do problema.
	Esta função retorna a raiz, o valor da função na raiz e o número de iterações necessárias.

# Exemplo

colocar exemplo
"
function newton_raphson(x0,f,df,ε=10.0^(-4))
	x = x0 - f(x0)/df(x0)
	k = 1
	while abs(x-x0)>=ε && abs(f(x))>ε
		x0 = x
		x = x0 - f(x0)/df(x0)
		k += 1
	end
	return x,f(x),k
end

end # module


text:
abstract type ExactSolutionType end

export ExactSolutionType, ExactTrig, calcErrorL2Norm

mutable struct ExactTrig <: ExactSolutionType
end
function (obj::ExactTrig)(xy::AbstractArray{Tmsh}, q::AbstractArray{Tsol, 1}) where {Tmsh, Tsol}  # (numDofPerNode))
  k = 2.0
  q[:] = sin(2*k*pi*xy[1])*sin(2*k*pi*xy[2])
  return nothing
end

mutable struct ExactExpTrig <: ExactSolutionType
end
function (obj::ExactExpTrig)(xy::AbstractArray{Tmsh}, q::AbstractArray{Tsol, 1}) where {Tmsh, Tsol}  # (numDofPerNode))
  k = 2.0
  ex = exp(xy[1] + xy[2])
  q[:] = ex * sin(2*k*pi*xy[1])*sin(2*k*pi*xy[2])
  return nothing
end

mutable struct ExactPoly2nd <: ExactSolutionType
end
function (obj::ExactPoly2nd)(xy::AbstractArray{Tmsh}, q::AbstractArray{Tsol, 1}) where {Tmsh, Tsol}
  a = 1.0
  b = 1.0
  c = 1.0
  q[:] =  a*xy[1]*xy[1] + b*xy[1]*xy[2] + c*xy[2]*xy[2]
  return nothing
end

global const ExactDict = Dict{String, ExactSolutionType}(
 "ExactTrig" => ExactTrig(),
 "ExactPoly2nd" => ExactPoly2nd(),
 "ExactExpTrig" => ExactExpTrig()
)

function calcErrorL2Norm(mesh::AbstractMesh{Tmsh},
                         sbp::AbstractOperator,
                         eqn::AbstractEllipticData{Tsol, Tres},
                         opts) where {Tmsh, Tsol, Tres}
  l2norm::Float64 = 0.
  lInfnorm::Float64 = 0.
  qe = Array{Tsol}(mesh.numDofPerNode)
  exactFunc = ExactDict[opts["exactSolution"]]

  lInfnorm = -1.0
  elem_with_max_dq = 0
  for e = 1 : mesh.numEl
    for n = 1 : mesh.numNodesPerElement
      xy = sview(mesh.coords, :, n, e)
      exactFunc(xy, qe)
      q = sview(eqn.q, :, n, e)
      jac = mesh.jac[n, e]
      for v = 1:mesh.numDofPerNode
        # dq = Real(q[v] - qe[v])
        dq = Float64(q[v] - qe[v])
        l2norm += dq*dq*sbp.w[n]/jac
        if lInfnorm < abs(dq)
          lInfnorm = abs(dq)
        end
      end
    end
  end

  return sqrt(l2norm), lInfnorm
end


text:
# ---
# Test
# ---
using Test, StatsModels
using DataFrames

using unfold
include("test_utilities.jl")

data,evts = loadtestdata("testCase1") #
f  = @formula 0~1+continuousA+continuousB # 1

# prepare data
data_r = reshape(data,(1,:))
data_r = vcat(data_r,data_r)#add second channel
data_e,times = unfold.epoch(data=data_r,tbl=evts,τ=(-1.,1.9),sfreq=10) # cut the data into epochs
data_e_missing = data_e
data_e_missing[1,25,end-5:end] .= missing
data_missing = Array{Union{Missing,Number}}(undef,size(data_r))
data_missing .= data_r
data_missing[4500:4600] .= missing

## Mass Univariate Linear

m_mul,m_mul_results = fit(UnfoldLinearModel,f,evts,data_e,times)

@test all(m_mul_results[(m_mul_results.channel.==1).&(m_mul_results.colname_basis .==0.1),:estimate] .≈ [3.0 2.5 -1.5]')
# Timexpanded Univariate Linear
basisfunction = firbasis(τ=(-1,1),sfreq=10,name="A")
m_tul,m_tul_results = fit(UnfoldLinearModel,f,evts,data_r,basisfunction)
@test all(m_tul_results[(m_tul_results.channel.==1).&(m_tul_results.colname_basis .==0.1),:estimate] .≈ [3.0 2.5 -1.5]')

# Add Missing in Data
m_mul_missing,m_mul_missing_results = fit(UnfoldLinearModel,f,evts,data_e_missing,times)
@test m_mul_missing_results.estimate ≈ m_mul_results.estimate
# Timexpanded Univariate Linear
m_tul_missing,m_tul_missing_results = fit(UnfoldLinearModel,f,evts,data_missing,basisfunction)
@test  isapprox(m_tul_missing_results.estimate , m_tul_results.estimate,atol=1e-2)  # higher tol because we remove stuff

# runntime tests - does something explode?
for k in 1:3
    if k == 1
        f  = @formula 0~1
    elseif k == 2
        f  = @formula 0~1+continuousA
    elseif k == 3
        f  = @formula 0~1+continuousB
    end
    println("Testing Runtime $k with Formula:$f")
    Xs = designmatrix(UnfoldLinearModel,f,evts,basisfunction)

    # Fit the model
    df = unfoldfit(UnfoldLinearModel,Xs,data_e)
    c = condense_long(df,times)
    fit(UnfoldLinearModel,f,evts,data_e,times)
    fit(UnfoldLinearModel,f,evts,data,basisfunction)
end

##
data4,evts4 = loadtestdata("testCase4") #
data4 = reshape(data4,(1,:))

data4 = vcat(data4,data4)
f4  = @formula 0~1+conditionA+conditionB # 4
basisfunction4 = firbasis(τ=(-1,1),sfreq=1000,name="A")

@time designmatrix(UnfoldLinearModel,f4,evts4,basisfunction4)
# new version 7s-10s, dataset4, sfreq=1000, 1200stim,

###############################
##  Mixed Model tests
###############################
data,evts = loadtestdata("testCase3") #
append!(data,zeros(1000))
data = reshape(data,1,:)
data = vcat(data,data)
data = data.+ 1*randn(size(data)) # we have to add minimal noise, else mixed models crashes.
data_missing = Array{Union{Missing,Number}}(undef,size(data))
data_missing .= data
data_missing[4500:4600] .= missing

categorical!(evts,:subject)
f  = @formula 0~1+condA+condB + (1+condA+condB|subject)
#f  = @formula 0~1 + (1|subject)

# cut the data into epochs
# TODO This ignores subject bounds
data_e,times = unfold.epoch(data=data,tbl=evts,τ=(-1.,1.9),sfreq=10)
data_missing_e,times = unfold.epoch(data=data_missing,tbl=evts,τ=(-1.,1.9),sfreq=10)
evts_e,data_e = unfold.dropMissingEpochs(copy(evts),data_e)
evts_missing_e,data_missing_e = unfold.dropMissingEpochs(copy(evts),data_missing_e)

######################
##  Mass Univariate Mixed
@time m_mum = fit(UnfoldLinearMixedModel,f,evts_e,data_e    ,times,contrasts=Dict(:condA => EffectsCoding(), :condB => EffectsCoding()))
df = m_mum[2]
@test isapprox(df[(df.channel .== 1).&(df.term.=="condA: 1").&(df.colname_basis.==0.0),:estimate], [5.618,9.175],rtol=0.1)
@time m_mum = fit(UnfoldLinearMixedModel,f,evts_missing_e,data_missing_e    ,times,contrasts=Dict(:condA => EffectsCoding(), :condB => EffectsCoding()))
df = m_mum[2]
@test isapprox(df[(df.channel .== 1).&(df.term.=="condA: 1").&(df.colname_basis.==0.0),:estimate], [5.618,9.175],rtol=0.1)

# Timexpanded Univariate Mixed
f  = @formula 0~1+condA+condB + (1+condA+condB|subject)
basisfunction = firbasis(τ=(-0.2,0.3),sfreq=10,name="ABC")
@time m_tum = fit(UnfoldLinearMixedModel,f,evts,data,basisfunction, contrasts=Dict(:condA => EffectsCoding(), :condB => EffectsCoding()) )
df = m_tum[2]
@test isapprox(df[(df.channel .== 1).&(df.term.=="condA: 1").&(df.colname_basis.==0.0),:estimate], [5.618,9.175],rtol=0.1)

# missing data in LMMs
# not yet implemented
Test.@test_broken  m_tum = fit(UnfoldLinearMixedModel,f,evts,data_missing,basisfunction, contrasts=Dict(:condA => EffectsCoding(), :condB => EffectsCoding()) )

evts.subjectB = evts.subject;
evts1 = evts[evts.condA.==0,:]
evts2 = evts[evts.condA.==1,:]

f0_lmm  = @formula 0~1+condB+(1|subject) + (1|subjectB)
@time m_tum = fit(UnfoldLinearMixedModel,f0_lmm,evts,data,basisfunction)

f1_lmm  = @formula 0~1+condB+(1|subject)
f2_lmm  = @formula 0~1+condB+(1|subjectB)

b1 = firbasis(τ=(-0.2,0.3),sfreq=10,name="A")
b2 = firbasis(τ=(-0.1,0.3),sfreq=10,name="B")

X1_lmm  = designmatrix(UnfoldLinearMixedModel,f1_lmm,evts1,b1)
X2_lmm  = designmatrix(UnfoldLinearMixedModel,f2_lmm,evts2,b2)

r = unfold.unfoldfit(UnfoldLinearMixedModel,X1_lmm+X2_lmm,data);
df = condense_long(r)

@test isapprox(df[(df.channel .== 1).&(df.term.=="condB").&(df.colname_basis.==0.0),:estimate],[18.18,10.4],rtol=0.1)

if 1 == 0
    using WGLMakie,AlgebraOfGraphics
    m = mapping(:colname_basis,:estimate,color=:term,layout_x=:term,layout_y=:basisname)

    AlgebraOfGraphics.data(df[df.channel.==2,:]) * visual(Lines) * m  |> draw

end

##------

if 1 == 0
    # Fit mass-univariate 1st level for all subjects
    basisfunction = firbasis(τ=(-.1,.4),sfreq=10,name="A")
    resAll = DataFrame()
    f  = @formula 0~1+condA+condB
    for s in unique(evts.subject)
        ##
        from = minimum(evts[evts.subject.==s,:latency])-10
        to = maximum(evts[evts.subject.==s,:latency])+40
        to = min(to,size(data,1))
        evts_s = evts[evts.subject.==s,:]
        evts_s.latency .-= from
        m = fit(UnfoldLinearModel,f,evts_s,data[from:to],basisfunction,contrasts=Dict(:condA => EffectsCoding(), :condB => EffectsCoding()))
        m.results.subject = s
        append!(resAll,m.results)
    end

    results = resAll[resAll.term.=="(Intercept)",:]

    results[results.time .== .1,:]
end


text:
# TODO: when 0.5 support dropped, remove the isdefined check
const LINT_RESULT_COLORS = Dict(
    :INFO => Base.info_color(),
    :WARN => Base.warn_color(),
    :ERROR => isdefined(Base, :error_color) ? Base.error_color() :
                                              Base.warn_color(),
    :OK => :green)

"""
A collection of `LintMessage`s. This behaves similarly to a vector of
`LintMessage`s, but has different display behaviour.
"""
struct LintResult <: AbstractVector{LintMessage}
    messages::Array{LintMessage, 1}
end

function Base.show(io::IO, res::LintResult)
    print(io, "LintResult(")
    show(io, res.messages)
    print(io, ")")
end
function Base.show(io::IO, ::MIME"text/plain", res::LintResult)
    for m in res.messages
        Base.printstyled(io, string(m), color = LINT_RESULT_COLORS[level(m)])
    end
end

Base.size(r::LintResult) = size(r.messages)

# delegate getindex to parent collection
Base.getindex(r::LintResult, i::Int) = r.messages[i]

# this function is useful for filtering on severity level
Base.filter(p, r::LintResult) = LintResult(filter(p, r.messages))


text:
# Autogenerated wrapper script for SSGraphBLAS_jll for i686-linux-musl
export libgraphblas

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("SSGraphBLAS")
JLLWrappers.@declare_library_product(libgraphblas, "libgraphblas.so.6")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libgraphblas,
        "lib/libgraphblas.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
    CHELSA{BioClim} <: RasterDataSource

Data from CHELSA, currently only the `BioClim` layer is implemented.

See: [chelsa-climate.org](https://chelsa-climate.org/)
"""
struct CHELSA{X} <: RasterDataSource end

layers(::Type{CHELSA{BioClim}}) = 1:19

"""
    getraster(T::Type{CHELSA{BioClim}}, [layer::Integer]) => String

Download CHELSA BioClim data, choosing layers from: `$(layers(CHELSA{BioClim}))`.

Without a layer argument, all layers will be downloaded, and a tuple of paths is returned.
If the data is already downloaded the path will be returned.
"""
function getraster(T::Type{CHELSA{BioClim}}, layer::Integer)
    _check_layer(T, layer)
    path = rasterpath(T, layer)
    url = rasterurl(T, layer)
    return _maybe_download(url, path)
end

rastername(::Type{CHELSA{BioClim}}, layer::Integer) = "CHELSA_bio10_$(lpad(layer, 2, "0")).tif"

rasterpath(::Type{CHELSA}) = joinpath(rasterpath(), "CHELSA")
rasterpath(::Type{CHELSA{BioClim}}) = joinpath(rasterpath(CHELSA), "BioClim")
rasterpath(T::Type{CHELSA{BioClim}}, layer::Integer) = joinpath(rasterpath(T), rastername(T, layer))

# https://os.zhdk.cloud.switch.ch/envicloud/chelsa/chelsa_V1/climatologies/bio/CHELSA_prec_01_V1.2_land.tif
rasterurl(::Type{CHELSA}) = URI(scheme="https", host="os.zhdk.cloud.switch.ch", path="/envicloud/chelsa/chelsa_V1/")
rasterurl(::Type{CHELSA{BioClim}}) = joinpath(rasterurl(CHELSA), "climatologies/bio/")
rasterurl(T::Type{CHELSA{BioClim}}, layer::Integer) = joinpath(rasterurl(T), rastername(T, layer))


text:
function onedQuadratureRule(rule::Int)
#= """
  onedQuadratureRule - Return Gauss integration points over (-1,1)

  Author: Jeff Borggaard, Virginia Tech
          part of FEMfunctions.jl

  Licensing:
     This code is distributed under the MIT license.

  Usage:
  ```julia
    r,w = onedQuadratureRule(rule)
  ```

  Arguments:
  - `rule`: Number of Gauss points:

  Output arguments:
  - `r`: Gauss points located between (-1,1)
  - `w`: Gauss weights corresponding to r
""" =#

  r = zeros(Float64,rule);
  w = zeros(Float64,rule);

  if rule == 1         # up to order 1 polynomials integrated exactly
    r[1] = 0.0
    w[1] = 2.0

  elseif rule == 2     # up to order 3 polynomials integrated exactly
    r[1] =-1.0 / sqrt(3.0);
    r[2] =-r[1]
    w[1] = 1.0
    w[2] = 1.0

  elseif rule == 3     # up to order 5 polynomials integrated exactly
    r[1] =-sqrt(3.0/5.0)
    r[2] = 0.0
    r[3] =-r[1]
    w[1] = 5.0 / 9.0
    w[2] = 8.0 / 9.0
    w[3] = w[1]

  elseif rule == 4    # up to order 7 polynomials integrated exactly
    r[1] =-sqrt((3.0+2.0*sqrt(6.0/5.0))/7.0)
    r[2] =-sqrt((3.0-2.0*sqrt(6.0/5.0))/7.0)
    r[3] =-r[2]
    r[4] =-r[1]
    w[1] = 0.5 - 1.0 / ( 6.0 * sqrt(6.0/5.0) )
    w[2] = 0.5 + 1.0 / ( 6.0 * sqrt(6.0/5.0) )
    w[3] = w[2]
    w[4] = w[1]

  elseif rule == 5    # up to order 9 polynomials integrated exactly
    r[1] =-sqrt(5.0+4.0*sqrt(5.0/14.0)) / 3.0
    r[2] =-sqrt(5.0-4.0*sqrt(5.0/14.0)) / 3.0
    r[3] = 0.0
    r[4] =-r[2]
    r[5] =-r[1]
    w[1] = 161.0/450.0-13.0/(180.0*sqrt(5.0/14.0))
    w[2] = 161.0/450.0+13.0/(180.0*sqrt(5.0/14.0))
    w[3] = 128.0/225.0
    w[4] = w[2]
    w[5] = w[1]

  elseif rule == 6    # up to order 11 polynomials integrated exactly
    r[1] = -0.2386191860831969
    r[2] = -0.6612093864662645
    r[3] = -0.9324695142031521
    r[4] = - r[1]
    r[5] = - r[2]
    r[6] = - r[3]
    w[1] = 0.4679139345726910
    w[2] = 0.3607615730481386
    w[3] = 0.1713244923791704
    w[4] = w[1]
    w[5] = w[2]
    w[6] = w[3]

  elseif (rule == 7)    # up to order 13 polynomials integrated exactly
    r[1] = -0.9491079123427585
    r[2] = -0.7415311855993945
    r[3] = -0.4058451513773972
    r[4] =  0.0000000000000000
    r[5] = - r[3]
    r[6] = - r[2]
    r[7] = - r[1]
    w[1] = 0.1294849661688697
    w[2] = 0.2797053914892766
    w[3] = 0.3818300505051189
    w[4] = 0.4179591836734694
    w[5] = w[3]
    w[6] = w[2]
    w[7] = w[1]

  elseif (rule == 8)    # up to order 15 polynomials integrated exactly
    r[1] = -0.9602898564975363
    r[2] = -0.7966664774136267
    r[3] = -0.5255324099163290
    r[4] = -0.1834346424956498
    r[5] = - r[4]
    r[6] = - r[3]
    r[7] = - r[2]
    r[8] = - r[1]
    w[1] = 0.1012285362903763
    w[2] = 0.2223810344533745
    w[3] = 0.3137066458778873
    w[4] = 0.3626837833783620
    w[5] = w[4]
    w[6] = w[3]
    w[7] = w[2]
    w[8] = w[1]

  elseif (rule == 9)     #  up to order 17 polynomials integrated exactly
    r[1] = -0.9681602395076261
    r[2] = -0.8360311073266358
    r[3] = -0.6133714327005904
    r[4] = -0.3242534234038089
    r[5] =  0.0000000000000000
    r[6] = - r[4]
    r[7] = - r[3]
    r[8] = - r[2]
    r[9] = - r[1]
    w[1] = 0.0812743883615744
    w[2] = 0.1806481606948574
    w[3] = 0.2606106964029354
    w[4] = 0.3123470770400029
    w[5] = 0.3302393550012598
    w[6] = w[4]
    w[7] = w[3]
    w[8] = w[2]
    w[9] = w[1]

  elseif (rule == 10)    #  up to order 19 polynomials integrated exactly
    r[ 1] = -0.9739065285171717
    r[ 2] = -0.8650633666889845
    r[ 3] = -0.6794095682990244
    r[ 4] = -0.4333953941292472
    r[ 5] = -0.1488743389816312
    r[ 6] = - r[5]
    r[ 7] = - r[4]
    r[ 8] = - r[3]
    r[ 9] = - r[2]
    r[10] = - r[1]
    w[ 1] = 0.0666713443086881
    w[ 2] = 0.1494513491505806
    w[ 3] = 0.2190863625159820
    w[ 4] = 0.2692667193099963
    w[ 5] = 0.2955242247147529
    w[ 6] = w[5]
    w[ 7] = w[4]
    w[ 8] = w[3]
    w[ 9] = w[2]
    w[10] = w[1]

  elseif (rule == 11)     #  up to order 21 integrated exactly
    r[ 1] = -0.9782286581460570
    r[ 2] = -0.8870625997680953
    r[ 3] = -0.7301520055740494
    r[ 4] = -0.5190961292068118
    r[ 5] = -0.2695431559523450
    r[ 6] =  0.0000000000000000
    r[ 7] = - r[5]
    r[ 8] = - r[4]
    r[ 9] = - r[3]
    r[10] = - r[2]
    r[11] = - r[1]
    w[ 1] = 0.0556685671161737
    w[ 2] = 0.1255803694649046
    w[ 3] = 0.1862902109277343
    w[ 4] = 0.2331937645919905
    w[ 5] = 0.2628045445102467
    w[ 6] = 0.2729250867779006
    w[ 7] = w[5]
    w[ 8] = w[4]
    w[ 9] = w[3]
    w[10] = w[2]
    w[11] = w[1]

  else
    error("onedQuadratureRule, rule not supported")

  end

  return r, w
end


text:
export densitymatrix, ket

"""
    densitymatrix(net, prob, norm=true)

Returns the Density matrix encoded by the neural network `net`, and normalizes
it if `norm==true`.
"""
function densitymatrix(net, prob, norm=true)
    ρ = DenseOperator(basis(prob))
    v = state(prob, net)
    if v isa DiagonalStateWrapper
        v = v.parent
    end
    for i=1:spacedimension(row(v))
        set_index!(row(v), i)
        for j=1:spacedimension(col(v))
            set_index!(col(v),j)
            ρ.data[i,j] = exp(net(v))
        end
    end
    if norm
        ρ.data ./= tr(ρ)
    end
    ρ
end

QuantumOptics.dm(net::NeuralNetwork, prob::AbstractProblem, norm=false) =
    densitymatrix(net, prob, norm)

"""
    ket(net, prob, norm=true)

Returns the state (ket) encoded by the neural network `net`, and normalizes
it if `norm==true`.
"""
function ket(net, prob, norm=true)
    psi = Ket(basis(prob))
    v = state(prob, net)
    for i=1:spacedimension(v)
        set_index!(v, i)
        psi.data[i] = exp(net(v))
    end
    if norm
        normalize!(psi)
    end
    return psi
end


text:
@testset "aggregation" begin

    function test_mean(x)
        z = mean(x)
    end

    function test_mean2(x)
        z = mean(x, 1)
    end

    test_compare(test_mean; x=rand(3))
    test_compare(test_mean; x=rand(3,4))

    function test_sum_1(x)
        y = sum(x, dims=1)
        y2 = y .+ 2
        z = sum(y2)
    end

    function test_sum_2(x)
        y = sum(x, dims=2)
        y2 = y .+ 2
        z = sum(y2)
    end

    test_compare(test_sum_1; x=rand(3,4))
    test_compare(test_sum_2; x=rand(3,4))
end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

using Base.Test
# Check that the fallback test set throws immediately
@test 1 == 2


text:
module DiskCaches

using Serialization
using Pidfile

export DiskCache

"""
    DiskCache(path)

`DiskCache{K,V}(path)` constructs a cache with keys of type `K` and values of type `V`, mapped to the filesystem `path` using `serialize` and `deserialize`. `DiskCache` behaves like `Dict` by default, but `DiskCache{K,V,D}` behaves like `D`, an `AbstractDict{K, V}`.

A `DiskCache` supports fast reads of keys already on chip and slow writes of new keys to disk, but does not support modification of existing keys. While writes to existing keys will raise an error, it is always safe to call `get!`.
"""
struct DiskCache{K,V,D<:AbstractDict{K, V}} <: AbstractDict{K,V}
    path::String
    data::D
    function DiskCache{K,V,D}(path::AbstractString) where {K,V,D}
        path, ext = splitext(path)
        if !(ext == ".jls" || ext == "")
            throw(ArgumentError("Extension for DiskCache files is .jls"))
        end
        path = "$path.jls"

        dir = splitdir(path)[1]
        mkpath(dir)
        data = mkpidlock("$(path).lock") do
            if !(ispath(path))
                data = D()
                serialize(path, data)
            end
            deserialize(path)
        end
        return new{K,V,D}(String(path), data)
    end
end

DiskCache{K, V}(path) where {K, V} = DiskCache{K, V, Dict{K, V}}(path)
DiskCache(path) = DiskCache{Any, Any}(path)

function _load(c)
    _data = deserialize(c.path)
    for (key, value) in pairs(_data)
        if !haskey(c.data, key)
            c.data[key] = value
        end
    end
end

function _save(c)
    serialize(c.path, c.data)
end

function Base.length(c::DiskCache)
    mkpidlock("$(c.path).lock") do
        _load(c)
    end
    length(c.data)
end

function Base.isempty(c::DiskCache)
    if !isempty(c.data)
        return false
    else
        mkpidlock("$(c.path).lock") do
            _load(c)
        end
        return isempty(c.data)
    end
end

function Base.haskey(c::DiskCache, key)
    mkpidlock("$(c.path).lock") do
        _load(c)
    end
    return haskey(c.data, key)
end

function Base.get(c::DiskCache, key, default)
    if haskey(c.data, key)
        return c.data[key]
    else
        mkpidlock("$(c.path).lock") do
            _load(c)
        end
        return get(c.data, key, default)
    end
end

function Base.get(f::Union{Function, Type}, c::DiskCache, key)
    if haskey(c.data, key)
        return c.data[key]
    else
        mkpidlock("$(c.path).lock") do
            _load(c)
        end
        return get(f, c.data, key)
    end
end

function Base.getindex(c::DiskCache, key)
    if haskey(c.data, key)
        return c.data[key]
    else
        mkpidlock("$(c.path).lock") do
            _load(c)
        end
        if haskey(c.data, key)
            return c.data[key]
        else
            throw(KeyError(key))
        end
    end
end

function Base.keys(c::DiskCache)
    mkpidlock("$(c.path).lock") do
        _load(c)
    end
    keys(c.data)
end
function Base.values(c::DiskCache)
    mkpidlock("$(c.path).lock") do
        _load(c)
    end
    values(c.data)
end
Base.pairs(c::DiskCache) = c
function Base.iterate(c::DiskCache)
    mkpidlock("$(c.path).lock") do
        _load(c)
    end
    iterate(c.data)
end
Base.iterate(c::DiskCache, state) = iterate(c.data, state)

function Base.get!(c::DiskCache, key, default)
    if haskey(c.data, key)
        return c.data[key]
    else
        mkpidlock("$(c.path).lock") do
            _load(c)
            if !haskey(c.data, key)
                c.data[key] = default
                _save(c)
            end
        end
        return c.data[key]
    end
end

function Base.get!(f::Union{Function, Type}, c::DiskCache, key)
    if haskey(c.data, key)
        return c.data[key]
    else
        mkpidlock("$(c.path).lock") do
            _load(c)
            if !haskey(c.data, key)
                c.data[key] = f()
                _save(c)
            end
        end
        return c.data[key]
    end
end

function Base.setindex!(c::DiskCache, value, key)
    mkpidlock("$(c.path).lock") do
        _load(c)
        if haskey(c.data, key)
            throw(ErrorException("To ensure the validity of the on-chip cache, DiskCaches do not support value modification."))
        else
            c.data[key] = value
            _save(c)
        end
    end
    return c.data[key]
end

function Base.delete!(c::DiskCache, key)
    throw(ErrorException("To ensure the validity of the on-chip cache, DiskCaches do not support value modification."))
end

function Base.pop!(c::DiskCache)
    throw(ErrorException("To ensure the validity of the on-chip cache, DiskCaches do not support value modification."))
end

function Base.pop!(c::DiskCache, key)
    throw(ErrorException("To ensure the validity of the on-chip cache, DiskCaches do not support value modification."))
end

function Base.pop!(c::DiskCache, key, default)
    throw(ErrorException("To ensure the validity of the on-chip cache, DiskCaches do not support value modification."))
end

function Base.empty!(c::DiskCache) where {K, V}
    throw(ErrorException("To ensure the validity of the on-chip cache, DiskCaches do not support value modification."))
end

end


text:
@testset "FastNPU" begin

    npu = NPU(2,3)
    x = rand(Float32,2,10)
    z = npu(x)

    @test size(z) == (3,10)
    ps = params(npu)
    @test length(ps) == 3

    loss1() = sum(abs2, npu(x) .- reshape(x[1,:],1,:))
    gs = Flux.gradient(loss1, ps)
    gs = vcat([vec(gs[p]) for p in ps]...)

    fastnpu = FastNPU(2,3)
    loss2(p) = sum(abs2, fastnpu(x,p) .- reshape(x[1,:],1,:))
    @test DiffEqFlux.paramlength(fastnpu) == 14
    @test length(DiffEqFlux.initial_params(fastnpu)) == 14

    p = Flux.destructure(npu)[1]
    gf = Flux.gradient(loss2, p)[1]
    @test all(gs .== gf)

    loss3(x) = sum(abs2, npu(x) .- reshape(x[1,:],1,:))
    gs = Flux.gradient(loss3, x)

    loss4(x) = sum(abs2, fastnpu(x,p) .- reshape(x[1,:],1,:))
    gf = Flux.gradient(loss4, x)

    @test all(gs .== gf)
end

@testset "FastRealNPU" begin

    npu = RealNPU(2,3)
    x = rand(Float32,2,10)
    z = npu(x)

    @test size(z) == (3,10)
    ps = params(npu)
    @test length(ps) == 2

    loss1() = sum(abs2, npu(x) .- reshape(x[1,:],1,:))
    gs = Flux.gradient(loss1, ps)
    gs = vcat([vec(gs[p]) for p in ps]...)

    fastnpu = FastRealNPU(2,3)
    loss2(p) = sum(abs2, fastnpu(x,p) .- reshape(x[1,:],1,:))
    @test DiffEqFlux.paramlength(fastnpu) == 8
    @test length(DiffEqFlux.initial_params(fastnpu)) == 8

    p = Flux.destructure(npu)[1]
    gf = Flux.gradient(loss2, p)[1]
    @test all(gs .== gf)

    loss3(x) = sum(abs2, npu(x) .- reshape(x[1,:],1,:))
    gs = Flux.gradient(loss3, x)

    loss4(x) = sum(abs2, fastnpu(x,p) .- reshape(x[1,:],1,:))
    gf = Flux.gradient(loss4, x)

    @test all(gs .== gf)
end


text:
function action_matrix(a::Vector{Int}, num_actions::Int)
    n = length(a)
    A = zeros(num_actions, n)
    for i in 1:n
        A[a[i], i] = 1.0
    end
    return A
end

"""
    NeuralEncoder(d_int::Int64, d_out::Int64, layer_sizes::Vector{Int64})

Construct a callable object that returns the final layer activations of a
neural network **(which can be trained with a contextual multi armed bandit trajectory)**.
"""
mutable struct NeuralEncoder{T1<:Chain,T2<:Function}
    nn::T1
    loss::T2
    d_in::Int64
    d_out::Int64
    layer_sizes::Vector{Int64}

    function NeuralEncoder(d_in::Int64, d_out::Int64, layer_sizes::Vector{Int64})
        if length(layer_sizes) < 1
            throw(ArgumentError("must provide at least 1 layer"))
        end
        enc = [Dense(d_in, layer_sizes[1], relu)]
        if length(layer_sizes) > 1
            for i in 2:length(layer_sizes)
                push!(enc, Dense(layer_sizes[i - 1], layer_sizes[i], relu))
            end
        end
        nn = Chain(; enc=gpu(Chain(enc...)), dec=gpu(Dense(layer_sizes[end], d_out)))

        # Create partial loss function
        function loss(X::AbstractMatrix, A::AbstractMatrix, r::AbstractMatrix)
            Ypred = nn(X)
            mse = sum(((Ypred .- r) .* A) .^ 2) / size(r, 2)
            return mse
        end

        return new{typeof(nn),typeof(loss)}(nn, loss, d_in, d_out, layer_sizes)
    end
end

function (enc::NeuralEncoder)(X::AbstractArray)
    X = gpu(X)
    return cpu(enc.nn[:enc](X))
end

"""
    fit!(enc::NeuralEncoder, X::AbstractMatrix, a::AbstractVector{<:Int}, r::AbstractMatrix,
         epochs::Int64)

Update a the neural network parameters using a contextual multi-armed bandit trajectory.

# Keyword Arguments

- `batch_size::Int64=32`: The batch size to use while training the neural network.
- `opt=ADAM()`: The optimizer to use for updating the parameters.
- `verbose::Bool=true`: Print details of the fitting procedure.
"""
function fit!(
    enc::NeuralEncoder,
    X::AbstractMatrix,
    a::AbstractVector{<:Int},
    r::AbstractMatrix,
    epochs::Int64;
    batch_size::Int64=32,
    opt=ADAM(),
    verbose::Bool=true,
)
    check_regression_data(X, r)
    if length(a) != size(X, 2)
        throw(ArgumentError("length of a must equal the second dimension of X/r"))
    end
    X = gpu(X)
    A = gpu(action_matrix(a, enc.d_out))
    r = gpu(r)
    data = DataLoader((X, A, r); batchsize=batch_size, shuffle=true)
    for i in 1:epochs
        train!(enc.loss, params(enc.nn), data, opt)
        if verbose
            println("\rEpoch: ", i, "; Training loss: ", enc.loss(X, A, r))
        end
    end
end

"""
    NeuralLinear(d::Int64, num_arms::Int64, layer_sizes::Vector{Int64}, inital_batches::Int64,
                 retrain_freq::Int64, epochs::Int64)

NeuralLinear policy introduced in the paper [Deep Bayesian Bandits Showdown: An Empirical
Comparison of Bayesian Deep Networks for Thompson Sampling](https://arxiv.org/abs/1802.09127).

# Keyword Arguments

- `opt=ADAM()`: Optimizer used to update the neural network parameters.
- `λ::Float64=1.0`: Prior scaling.
- `shape0::Float64=1e-3`: Inverse-gamma prior shape hyperparameter.
- `scale0::Float64=1e-3`: Inverse-gamma prior scale hyperparameter.
- `verbose_retrain::Bool=false`: Print the details of the neural network training
    procedure.
"""
mutable struct NeuralLinear{T1<:NeuralEncoder,T2} <: AbstractPolicy
    t::Int64
    batches::Int64
    data::BanditDataset
    arms::Vector{BayesLM}
    layer_sizes::Vector{Int64}
    enc::T1
    initial_batches::Int64
    retrain_freq::Int64
    epochs::Int64
    opt::T2
    λ::Float64
    shape0::Float64
    scale0::Float64
    verbose_retrain::Bool

    function NeuralLinear(
        d::Int64,
        num_arms::Int64,
        layer_sizes::Vector{Int64},
        inital_batches::Int64,
        retrain_freq::Int64,
        epochs::Int64;
        opt=ADAM(),
        λ::Float64=1.0,
        shape0::Float64=1e-3,
        scale0::Float64=1e-3,
        verbose_retrain::Bool=false,
    )
        arms = Vector{BayesLM}(undef, num_arms)
        enc = NeuralEncoder(d, num_arms, layer_sizes)
        data = BanditDataset(d)
        return new{typeof(enc),typeof(opt)}(
            0,
            0,
            data,
            arms,
            layer_sizes,
            enc,
            inital_batches,
            retrain_freq,
            epochs,
            opt,
            λ,
            shape0,
            scale0,
            verbose_retrain,
        )
    end
end

function update!(
    pol::NeuralLinear, X::AbstractMatrix, a::AbstractVector{<:Int}, r::AbstractMatrix
)
    check_regression_data(X, r)
    pol.t += size(X, 2)
    pol.batches += 1
    append_data!(pol.data, X, a, r)
    if pol.batches < pol.initial_batches
        return nothing
    end
    if pol.batches % pol.retrain_freq == 0 || pol.batches == pol.initial_batches
        pol.enc = NeuralEncoder(size(X, 1), length(pol.arms), pol.layer_sizes)
        fit!(
            pol.enc,
            pol.data.X,
            pol.data.a,
            pol.data.r,
            pol.epochs;
            opt=pol.opt,
            verbose=pol.verbose_retrain,
        )
        for a in 1:length(pol.arms)
            Xa, ra = arm_data(pol.data, a)
            Za = pol.enc(Xa)
            pol.arms[a] = BayesLM(
                size(Za, 1); λ=pol.λ, shape0=pol.shape0, scale0=pol.scale0
            )
            fit!(pol.arms[a], Za, ra)
        end
    else
        Z = pol.enc(X)
        for i in unique(a)
            fit!(pol.arms[i], Z[:, a .== i], r[:, a .== i])
        end
    end
end

function (pol::NeuralLinear)(X::AbstractMatrix)
    n = size(X, 2)
    num_arms = length(pol.arms)
    actions = zeros(Int64, n)

    # Check if inital batches have been completed
    if pol.batches < pol.initial_batches
        for i in 1:n
            actions[i] = (pol.t + i) % num_arms + 1
        end
        return actions
    end

    Z = pol.enc(X)
    thompson_samples = zeros(num_arms)
    for i in 1:n
        for a in 1:num_arms
            shape, scale = shape_scale(pol.arms[a])
            z = Z[:, i]
            varsim = rand(InverseGamma(shape, scale))
            β = pol.arms[a].β
            Σ = varsim * pol.arms[a].Σ
            βsim = rand(MvNormal(β[:, 1], Σ))
            thompson_samples[a] = βsim' * z
        end
        actions[i] = argmax(thompson_samples)
    end
    return actions
end


text:
@testset "$TEST $G" begin

if G <: ANetwork
## EDGE
g = G(10, 20)
m = eprop!(g, "label", Int)

@test valtype(m) == Int
@test typeof(m) <: AEdgeMap
e = first(edges(g))
m[e] = 2
@test m[e] == 2

@test eprop(g, "label") === m
@test_throws ErrorException eprop!(g, "label", Int)

rem_eprop!(g, "label")
@test_throws ErrorException rem_eprop!(g, "label")
@test_throws KeyError eprop(g, "label")

m = eprop!(g, "hi", String)
@test valtype(m) == String
@test typeof(m) <: AEdgeMap
e = first(edges(g))
m[e] = "ciao"
@test m[e] == "ciao"

m = EdgeMap(g, String)
m2 = eprop!(g, "bye", m)
@test m === m2
@test eprop_names(g) == ["bye","hi"]

edg = collect(edges(g))
@test eprop(g, edg[1]) == Dict("hi"=>"ciao")
@test eprop(g, edg[2]) == Dict()
eprop(g, "bye")[edg[1]] = "e1"
eprop(g, "bye")[edg[2]] = "e2"
@test eprop(g, edg[1]) == Dict("hi"=>"ciao", "bye"=>"e1")
@test eprop(g, edg[2]) == Dict("bye"=>"e2")

@test has_eprop(g, "bye", edg[1])
@test has_eprop(g, "hi", edg[1])
@test has_eprop(g, "bye", edg[2])
@test !has_eprop(g, "hi", edg[2])

## VERTEX
g = G(3)
add_edge!(g,1,2)
add_edge!(g,2,3)
m = vprop!(g, "label", [1,2,3])
@test valtype(m) == Int
@test typeof(m) <: AVertexMap
@test vprop(g, "label") === m

rem_vertex!(g, 1)
@test m[1] == 3
@test m[2] == 2

swap_vertices!(g, 1, 2)
@test m[1] == 2
@test m[2] == 3

m = vprop!(g, "size", Int)
m[1] =  4
@test vprop(g, 1) == Dict("size"=>4, "label"=>2)
@test vprop(g, 2) == Dict("label"=>3)

@test has_vprop(g, "size", 1)
@test has_vprop(g, "label", 1)
@test !has_vprop(g, "size", 2)
@test has_vprop(g, "label", 2)

## GRAPH
g = DG(3, 5)
gprop!(g, "lab", "mygraph")
@test gprop(g, "lab") == "mygraph"
@test gprop_names(g) == ["lab"]
rem_gprop!(g, "lab")
@test gprop_names(g) == []

end

end #testset


text:
import SpecialFunctions: trigamma

"""
    OhmicBath

Ohmic bath object to hold a particular parameter set.

**Fields**
- `η` -- strength.
- `ωc` -- cutoff frequence.
- `β` -- inverse temperature.
"""
struct OhmicBath <: AbstractBath
    η::Float64
    ωc::Float64
    β::Float64
end

"""
$(SIGNATURES)

Construct OhmicBath from parameters with physical unit: `η`--unitless interaction strength; `fc`--cutoff frequency in GHz; `T`--temperature in mK.
"""
function Ohmic(η, fc, T)
    ωc = 2 * π * fc
    β = temperature_2_β(T)
    OhmicBath(η, ωc, β)
end

"""
$(SIGNATURES)

Calculate spectral density ``γ(ω)`` of `bath`.
"""
function spectrum(ω, bath::OhmicBath)
    if isapprox(ω, 0.0, atol=1e-9)
        return 2 * pi * bath.η / bath.β
    else
        return 2 * pi * bath.η * ω * exp(-abs(ω) / bath.ωc) / (1 - exp(-bath.β * ω))
    end
end

"""
$(SIGNATURES)

Calculate the two point correlation function ``C(τ)`` of `bath`.
"""
function correlation(τ, bath::OhmicBath)
    x2 = 1 / bath.β / bath.ωc
    x1 = 1.0im * τ / bath.β
    bath.η * (trigamma(-x1 + 1 + x2) + trigamma(x1 + x2)) / bath.β^2
end

"""
$(SIGNATURES)

Calculate the polaron transformed correlation function of `bath`.
"""
function polaron_correlation(τ, bath::OhmicBath)
    res = (1 + 1.0im * bath.ωc * τ)^(-4 * bath.η)
    if !isapprox(τ, 0, atol=1e-9)
        x = π * τ / bath.β
        res *= (x / sinh(x))^(4 * bath.η)
    end
    res
end

@inline polaron_correlation(t1, t2, bath::OhmicBath) = polaron_correlation(t1 - t2, bath)

function info_freq(bath::OhmicBath)
    println("ωc (GHz): ", bath.ωc / pi / 2)
    println("T (GHz): ", temperature_2_freq(β_2_temperature(bath.β)))
end

function Base.show(io::IO, ::MIME"text/plain", m::OhmicBath)
    print(
        io,
        "Ohmic bath instance:\n",
        "η (unitless): ",
        m.η,
        "\n",
        "ωc (GHz): ",
        m.ωc / pi / 2,
        "\n",
        "T (mK): ",
        β_2_temperature(m.β),
    )
end


text:
const REQUEST = 0
const RESPONSE = 1
const NOTIFICATION = 2

type NvimClient{S} <: NvimObject
    input::S #input to nvim
    output::S

    channel_id::Int
    next_reqid::Int
    waiting::Dict{Int,RemoteRef}
    reader::Task
    NvimClient(a,b,c,d,e) = new(a,b,c,d,e)
end

function NvimClient{S}(input::S, output::S, handler=DummyHandler())
    c = NvimClient{S}(input, output, -1, 0, Dict{Int,RemoteRef}())
    c.reader = @async readloop(c,handler)
    c.channel_id, metadata = send_request(c, "vim_get_api_info", [])
    #println("CONNECTED $(c.channel_id)"); flush(STDOUT)
    if symbolize(metadata) != _metadata
        println("warning: possibly incompatible api metadata")
    end
    c
end

# this is probably not most efficient in the common case (no contention)
# but it's the simplest way to assure task-safety of reading from the stream
function readloop(c::NvimClient, handler)
    while !eof(c.output)
        msg = unpack(c.output)
        kind = msg[1]::Int
        if kind == RESPONSE
            serial = msg[2]::Int
            ref = pop!(c.waiting, serial)
            put!(ref, (msg[3], msg[4]))
        elseif kind == NOTIFICATION
            name = bytestring(msg[2])
            args = retconvert(Any, c, msg[3])
            try
                on_notify(handler, c, name, args)
            catch err
                logerr(err, catch_backtrace(), "handler", "notification", name, args)
            end
        elseif kind == REQUEST
            serial = msg[2]::Int
            method = bytestring(msg[3])
            args = retconvert(Any, c, msg[4])
            try
                on_request(handler, c, serial, method, args)
            catch err
                reply_error(c, serial, "Caught Exception in request handler.")
                logerr(err, catch_backtrace(), "handler", "request", method, args)
            end
        end
        flush(STDERR)
    end
end

function logerr(e::Exception, bt, where, what, name, args)
    println(STDERR, "Caught Exception in $where for $what \"$name\" with arguments $args:")
    showerror(STDERR, e, bt)
    println(STDERR, "")
    flush(STDERR)
end

Base.wait(c::NvimClient) = wait(c.reader)

# we cannot use pack(stream, msg) as it's not synchronous
_send(c, msg) = (write(c.input, pack(msg)), flush(c.input))

function send_request(c::NvimClient, meth, args)
    reqid = c.next_reqid
    c.next_reqid += 1
    # TODO: are these things cheap to alloc or should they be reused
    res = RemoteRef()
    c.waiting[reqid] = res
    meth = string(meth)

    _send(c, Any[REQUEST, reqid, meth, args])
    (err, res) = take!(res) #blocking
    # TODO: make these recoverable
    if err !== nothing
        error(string(meth, ": ", bytestring(err[2])))
    end
    res
end

function reply_error(c, serial, err)
    _send(c, Any[RESPONSE, serial, err, nothing])
end

function reply_result(c, serial, res)
    _send(c, Any[RESPONSE, serial, nothing, res])
end

# when overriding these, note that this runs in the reader task,
# use @async/@spawn when doing anything long-running or blocking,
# including method calls to nvim
immutable DummyHandler; end
function on_notify(::DummyHandler, c, name, args)
    println("notification: $name $args")
end

function on_request(::DummyHandler, c, serial, name, args)
    println("WARNING: ignoring request $name $args, please override `on_request`")
    reply_error(c, serial, "Client cannot handle request, please override `on_request`")
end


text:
using Test
using DynamicPolynomials
# using Revise
using PolynomialBasis

function allequal(v1,v2)
    return all(v1 .≈ v2)
end

coeffs = PolynomialBasis.hermite_polynomial_coefficients(-1.0,1.0)
testcoeffs = 0.25*[2    1  2  -1
                   -3  -1  3  -1
                   0   -1  0  +1
                   1    1 -1  +1]
@test allequal(coeffs,testcoeffs)

@polyvar x

hermitepolynomials = PolynomialBasis.hermite_polynomials(x,-1.0,1.0)
derivatives = differentiate.(hermitepolynomials,x)
v1 = [f(-1.0) for f in hermitepolynomials]
v2 = [f(-1.0) for f in derivatives]
v3 = [f(1.0) for f in hermitepolynomials]
v4 = [f(1.0) for f in derivatives]
@test allequal(v1,[1.0,0.,0.,0.])
@test allequal(v2,[0.,1.,0.,0.])
@test allequal(v3,[0.,0.,1.,0.])
@test allequal(v4,[0.,0.,0.,1.])

hermitepolynomials = PolynomialBasis.hermite_polynomials(x,0.0,1.0)
derivatives = differentiate.(hermitepolynomials,x)
v1 = [f(0.0) for f in hermitepolynomials]
v2 = [f(0.0) for f in derivatives]
v3 = [f(1.0) for f in hermitepolynomials]
v4 = [f(1.0) for f in derivatives]
@test allequal(v1,[1.0,0.,0.,0.])
@test allequal(v2,[0.,1.,0.,0.])
@test allequal(v3,[0.,0.,1.,0.])
@test allequal(v4,[0.,0.,0.,1.])


text:
#
# JSON Server with julia
#
# author: Atsushi Sakai
#

using HTTP
using JSON
using Dates

function show_http_request(http)
	@show http.message
	@show HTTP.header(http, "Content-Type")
	while !eof(http)
		println("body data: ", String(readavailable(http)))
	end
end

function get_config()
    ip = "0.0.0.0" #default
    if haskey(ENV,"SERVER_IP")
        ip = ENV["SERVER_IP"]
    end

	port = 8000 #default
    if haskey(ENV,"SERVER_PORT")
        port =parse(Int64, ENV["SERVER_PORT"])
    end

    println("serving at ip:", ip)
    println("serving at port:", port)

    return ip, port
end

function main()
    println(PROGRAM_FILE," start!!")

    ip, port = get_config()

	HTTP.listen(ip, port) do http::HTTP.Stream

		show_http_request(http)

		#response json
	    jdic = Dict("title"=>"Julia JSON server", "PI" => pi)
    	jdic["time"] = Dates.format(Dates.now(), "yyyymmddHHMMSS")

		write(http, JSON.json(jdic))
	end

    println(PROGRAM_FILE," Done!!")
end

main()

max_stars_repo_name: AtsushiSakai/JSONServerAndClient

max_issues_repo_name: AtsushiSakai/JSONServerAndClient

max_forks_repo_name: AtsushiSakai/JSONServerAndClient


text:
export concretize

untype(::Type{T}) where {T} = T
parameters(::Type{T}) where {T} = T.parameters
parameters(t) = t.parameters

function concretize(ntt::Core.SimpleVector, gid) # typeof(ntt)==Core.SimpleVector
    k_lev = ntt[1]
    v_lev = tuple((
    begin
        if v <: NamedTuple
            val, gid = concretize(parameters(v), gid)
        elseif v <: NTuple
            v_ntuple = parameters(v)
            N = length(v_ntuple)
            val = ntuple(i->begin
                nt, gid = concretize(parameters(v_ntuple[i]), gid)
                nt
                end
                , N)
        else
            val = (type=Val(v), gid=Val(gid))
            gid+=get_len(v)
        end
        val
    end
    for v in parameters(untype(ntt[2])))...)
    nt = (; zip(k_lev, v_lev)...)
    return (nt, gid)
end

"""
    concretize(ntt::Type{NamedTuple}, gid=1)

A nested NamedTuple given a nested
NamedTuple type `ntt`, whose values
are types. The concretized nested
NamedTuple is a NamedTuple with
`type` (the original type values)
and `gid`, a global index, which
is incremented by `get_len` of the
type.
"""
function concretize(ntt::Type{T}, gid=1) where {T} # typeof(ntt)==DataType
    nt,gid = concretize(parameters(ntt), gid)
    return nt
end


text:
# heuristics.jl
# heuristic policies

mutable struct Simple{M} <: Policy #
  mdp::M
  sweeping_up::Bool
end
mutable struct SimpleSolver <: Solver end

Simple(mdp) = Simple(mdp,true)
solve(solver::SimpleSolver, problem::MDP) = Simple(problem)
solve(solver::SimpleSolver, problem::EmbeddedBehaviorMDP) = Simple(problem.base)
solve(solver::SimpleSolver, problem::POMDP) = Simple(problem)
solve(solver::SimpleSolver, problem::AggressivenessBeliefMDP) = Simple(get(problem.up.problem))
solve(solver::SimpleSolver, problem::BehaviorBeliefMDP) = Simple(get(problem.up.problem))
solve(solver::SimpleSolver, problem::QMDPWrapper) = Simple(problem.mdp)
solve(solver::SimpleSolver, problem::OutcomeMDP) = Simple(problem.mdp)
POMDPs.updater(::Simple) = POMDPToolbox.FastPreviousObservationUpdater{MLObs}()
create_policy(s::SimpleSolver, problem::MDP) = Simple(problem)
create_policy(s::SimpleSolver, problem::POMDP) = Simple(problem)

set_rng!(solver::SimpleSolver, rng::AbstractRNG) = nothing
Base.srand(p::Simple, s) = p

function action(p::Simple,s::Union{MLState,MLObs})
# lane changes if there is an opportunity
    a = MLAction()
    goal_lane = p.mdp.rmodel.target_lane
    y_desired = goal_lane
    dmodel = p.mdp.dmodel
    lc = sign(y_desired-s.cars[1].y) * dmodel.lane_change_rate
    acc = dmodel.adjustment_acceleration

    #if can't move towards desired lane sweep through accelerating and decelerating

    if is_safe(p.mdp, s, MLAction(0.,lc))
        return MLAction(0.,lc)
    end

    # maintain distance
    nbhd = get_neighborhood(dmodel.phys_param,s,1)

    if nbhd[2] == 0 && nbhd[5] == 0
        return MLAction(0.,0.)
    end

    dist_ahead = nbhd[2] != 0 ? s.cars[nbhd[2]].x - s.cars[1].x : Inf
    dist_behind = nbhd[5] != 0 ? s.cars[nbhd[5]].x - s.cars[1].x : Inf

    sgn = abs(dist_ahead) <= abs(dist_behind) ? -1 : 1

    accel = sgn * acc

    max_accel = max_safe_acc(p.mdp, s, 0.0)

    return MLAction(min(accel, max_accel),0.)
end
action(p::Simple, b::BehaviorBelief) = action(p, b.physical)
function action(pol::Simple, s::QMDPState)
    if s.isstate
        return action(pol, get(s.s))
    else
        return action(pol, get(s.b))
    end
end
action(p::Simple, b::AbstractParticleBelief) = action(p, MLPhysicalState(first(particles(b))))

mutable struct BehaviorSolver <: Solver
    b::BehaviorModel
    keep_lane::Bool
    rng::AbstractRNG
end
mutable struct BehaviorPolicy <: Policy
    problem::NoCrashProblem
    b::BehaviorModel
    keep_lane::Bool
    rng::AbstractRNG
end
solve(s::BehaviorSolver, p::NoCrashProblem) = BehaviorPolicy(p, s.b, s.keep_lane, s.rng)

function action(p::BehaviorPolicy, s::MLState, a::MLAction=MLAction(0.0,0.0))
    nbhd = get_neighborhood(p.problem.dmodel.phys_param, s, 1)
    acc = gen_accel(p.b, p.problem.dmodel, s, nbhd, 1, p.rng)
    if p.keep_lane
        lc = 0.0
    else
        lc = gen_lane_change(p.b, p.problem.dmodel, s, nbhd, 1, p.rng)
    end
    return MLAction(acc, lc)
end
action(p::BehaviorPolicy, b::AggressivenessBelief, a::MLAction=MLAction(0.0,0.0)) = action(p, most_likely_state(b))
action(p::BehaviorPolicy, b::BehaviorParticleBelief, a::MLAction=MLAction(0.0,0.0)) = action(p, most_likely_state(b))

mutable struct IDMLaneSeekingSolver <: Solver
    b::BehaviorModel
    rng::AbstractRNG
end

mutable struct IDMLaneSeekingPolicy <: Policy
    problem::NoCrashProblem
    b::BehaviorModel
    rng::AbstractRNG
end
solve(s::IDMLaneSeekingSolver, p::NoCrashProblem) = IDMLaneSeekingPolicy(p, s.b, s.rng)

function action(p::IDMLaneSeekingPolicy, s::MLState, a::MLAction=MLAction(0.0,0.0))
    nbhd = get_neighborhood(p.problem.dmodel.phys_param, s, 1)
    acc = gen_accel(p.b, p.problem.dmodel, s, nbhd, 1, p.rng)
    # try to positive lanechange
    # lc = problem.dmodel.lane_change_rate * !is_lanechange_dangerous(pp,s,nbhd,1,1)
    lc = p.problem.dmodel.lane_change_rate
    if is_safe(p.problem, s, MLAction(acc, lc))
        return MLAction(acc, lc)
    end
    return MLAction(acc, 0.0)
end

struct OptimisticValue end

function MCTS.estimate_value(v::OptimisticValue, m::NoCrashProblem, s::Union{MLState,MLPhysicalState}, steps::Int)
    rm = m.rmodel
    if rm isa NoCrashRewardModel
        rw = rm.reward_in_target_lane
    elseif rm isa SuccessReward
        rw = 1.0
    else
        error("Unrecognized reward model.")
    end
    dm = m.dmodel
    dt = dm.phys_param.dt
    steps_away = abs(rm.target_lane - first(s.cars).y) / (dm.lane_change_rate * dt)
    if s.x + first(s.cars).vel*steps_away*dt > dm.max_dist
        return 0.0
    elseif dm.lane_terminate
        return discount(m)^steps_away*rw
    else
        gamma = discount(m)
        return gamma^steps_away*rw/(1-gamma)
    end
end

MCTS.estimate_value(v::OptimisticValue, m::MLPOMDP, s::Union{MLState,MLPhysicalState}, h, steps::Int) = estimate_value(v, m, s, steps)