text:
"""
Basic chemical kinetics.

Types

* [`TempDep`](@ref)
* [`TempDepQ10`](@ref)
* [`TempDepArrhenius`](@ref)
* [`TempDepEnzyme`](@ref)

Functions

* [`arrhenius`](@ref)
* [`q10_temp_dep`](@ref)
* [`enzyme_temp_dep`](@ref)
* [`enzyme_temp_opt`](@ref)
* [`eval_temp_dep`](@ref)
"""
module ChemKinet

include("../docstring_style.jl")

using Canopy.Constants: R

export TempDep,
    TempDepQ10,
    TempDepArrhenius,
    TempDepEnzyme,
    arrhenius,
    q10_temp_dep,
    enzyme_temp_dep,
    enzyme_temp_opt,
    eval_temp_dep

"""A abstract type for structs of temperature dependence parameters."""
abstract type TempDep end

"""
Exponential (\$Q_{10}\$) temperature dependence parameters.

# Fields

* `q10`: The \$Q_{10}\$ parameter in [`q10_temp_dep`](@ref), dimensionless.
* `temp_ref`: Reference temperature [K].
* `v_max_ref`: Reaction rate at the reference temperature, arbitrary unit.
"""
struct TempDepQ10
    q10::Number
    temp_ref::Number
    v_max_ref::Number
end

"""
Arrhenius temperature dependence parameters.

# Fields

* `e_act`: The activation energy [J mol^-1] in [`arrhenius`](@ref).
* `temp_ref`: Reference temperature [K].
* `v_max_ref`: Reaction rate at the reference temperature, arbitrary unit.
"""
struct TempDepArrhenius
    e_act::Number
    temp_ref::Number
    v_max_ref::Number
end

"""
Enzyme optimum temperature dependence parameters.

# Fields

* `delta_G_a`, `delta_H_d`, `delta_S_d`: Enzyme energetic parameters in
  [`enzyme_temp_dep`](@ref).
* `temp_ref`: Reference temperature [K].
* `v_max_ref`: Reaction rate at the reference temperature, arbitrary unit.
"""
struct TempDepEnzyme
    delta_G_a::Number
    delta_H_d::Number
    delta_S_d::Number
    temp_ref::Number
    v_max_ref::Number
end

"""
Calculate the ratio of reaction rate at a given temperature (`temp`) to that at
a reference temperature (`temp_ref`) using the Arrhenius function. The
activation energy is given by `e_act`.

```math
\\dfrac{k}{k_\\mathrm{ref}} = \\exp\\left[ -\\dfrac{E_\\mathrm{act}}{R} \\cdot
\\left( \\dfrac{1}{T} - \\dfrac{1}{T_\\mathrm{ref}}\\right) \\right]
```

# Examples

```jldoctest
julia> arrhenius(5e4, 298.15, 283.15)
0.3435224255653669
```
"""
function arrhenius(e_act, temp_ref, temp)
    exp(e_act * (temp - temp_ref) / (R * temp_ref * temp))
end

"""
Calculate the ratio of reaction rate at a given temperature (`temp`) to that at
a reference temperature (`temp_ref`) using an exponential function defined by a
`q10` value.

```math
\\dfrac{k}{k_\\mathrm{ref}} = \\exp\\left( Q_{10} \\cdot
\\dfrac{T - T_\\mathrm{ref}}{10\\ \\mathrm{K}} \\right)
```

# Examples

```jldoctest
julia> q10_temp_dep(2.0, 298.15, 283.15)
0.3535533905932738
```
"""
function q10_temp_dep(q10, temp_ref, temp)
    q10^(0.1 * (temp - temp_ref))
end

"""
A helper function to calculate the enzyme reaction rate at a given temperature
(`temp`) in an arbitrary unit.

```math
f(T) = \\dfrac{\\exp\\left(-\\dfrac{\\Delta G_\\mathrm{a}}{R T}\\right)}{
    1 + \\exp\\left(
        -\\dfrac{\\Delta H_\\mathrm{d} - T\\Delta S_\\mathrm{d}}{RT}\\right)}
```

# Arguments

* `delta_G_a`: The standard Gibbs free energy of activation of the active state
  of the enzyme [J mol^-1].
* `delta_H_d`: The standard enthalpy change when the enzyme switches from the
  active to the deactivated state [J mol^-1].
* `delta_S_d`: The standard entropy change when the enzyme switches from the
  active to the deactivated state [J mol^-1 K^-1].
* `temp`: Current temperature [K].

!!! note

    This function is reserved for internal calls. Use
    [`enzyme_temp_dep`](@ref) for calculation.
"""
function f_enzyme_temp_dep(delta_G_a, delta_H_d, delta_S_d, temp)
    exp(-delta_G_a / (R * temp)) /
    (1.0 + exp((delta_S_d - delta_H_d / temp) / R))
end

"""
Calculate the ratio of reaction rate at a given temperature (`temp`) to that at
a reference temperature (`temp_ref`) for an enzyme reaction with a temperature
optimum.

```math
\\dfrac{V_\\mathrm{max}(T)}{V_\\mathrm{max}(T_\\mathrm{ref})} =
\\dfrac{f(T)}{f(T_\\mathrm{ref})} \\text{, where}\\quad f(T) =
\\dfrac{\\exp\\left(-\\dfrac{\\Delta G_\\mathrm{a}}{RT}\\right)}{1 + \\exp
\\left(-\\dfrac{\\Delta H_\\mathrm{d} - T\\Delta S_\\mathrm{d}}{RT}\\right)}
```

# Arguments

* `delta_G_a`: The standard Gibbs free energy of activation of the active state
  of the enzyme [J mol^-1].
* `delta_H_d`: The standard enthalpy change when the enzyme switches from the
  active to the deactivated state [J mol^-1].
* `delta_S_d`: The standard entropy change when the enzyme switches from the
  active to the deactivated state [J mol^-1 K^-1].
* `temp_ref`: Reference temperature [K].
* `temp`: Current temperature [K].

# References

* Johnson, F. H., Eyring, H., & Williams, R. W. (1942). The nature of enzyme
  inhibitions in bacterial luminescence: sulfanilamide, urethane, temperature
  and pressure. _Journal of Cellular and Comparative Physiology_, _20_(3),
  247--268. <https://doi.org/10.1002/jcp.1030200302>
* Peterson, M. E., Eisenthal, R., Danson, M. J., Spence, A., & Daniel, R. M.
  (2004). A new intrinsic thermal parameter for enzymes reveals true
  temperature optima. _Journal of Biological Chemistry_, _279_(20),
  20,717--20,722. <https://doi.org/10.1074/jbc.M309143200>
* Sharpe, P. J. H., & DeMichele, D. W. (1977). Reaction kinetics of
  poikilotherm development. _Journal of Theoretical Biology_, _64_(4),
  649--670. <https://doi.org/10.1016/0022-5193(77)90265-X>

!!! note

    Interpretations of the thermodynamic parameters in this equation may differ
    in the literature, but the mathematical form remains largely the same.

# Examples

```jldoctest
julia> enzyme_temp_dep(4e4, 2e5, 660.0, 298.15, 283.15)
0.5393220862804745
```

# See also

[`enzyme_temp_opt`](@ref)
"""
function enzyme_temp_dep(delta_G_a, delta_H_d, delta_S_d, temp_ref, temp)
    f_enzyme_temp_dep(delta_G_a, delta_H_d, delta_S_d, temp) /
    f_enzyme_temp_dep(delta_G_a, delta_H_d, delta_S_d, temp_ref)
end

"""
Calculate the temperature optimum of an enzyme reaction [K].

```math
T_\\mathrm{opt} = \\dfrac{\\Delta H_\\mathrm{d}}{
    \\Delta S_\\mathrm{d} +
    R \\ln \\left(\\Delta H_\\mathrm{d} / \\Delta G_\\mathrm{a} - 1 \\right)}
```

# Arguments

* `delta_G_a`: The standard Gibbs free energy of activation of the active state
  of the enzyme [J mol^-1].
* `delta_H_d`: The standard enthalpy change when the enzyme switches from the
  active to the deactivated state [J mol^-1].
* `delta_S_d`: The standard entropy change when the enzyme switches from the
  active to the deactivated state [J mol^-1 K^-1].

# Examples

```jldoctest
julia> enzyme_temp_opt(4e4, 2e5, 660.0)
297.8289937283252
```

# See also

[`enzyme_temp_dep`](@ref)
"""
function enzyme_temp_opt(delta_G_a, delta_H_d, delta_S_d)
    delta_H_d / (delta_S_d + R * log(delta_H_d / delta_G_a - 1.0))
end

"""
Evaluate the reaction rate at a given temperature (`temp`). The temperature
dependence of the reaction is characterized by a set of parameters (`p`), which
can be one of the three types: [`TempDepQ10`](@ref),
[`TempDepArrhenius`](@ref), and [`TempDepEnzyme`](@ref). This function is a
unified interface to different temperature dependence functions.
"""
function eval_temp_dep(p::TempDepQ10, temp)
    p.v_max_ref * q10_temp_dep(p.q10, p.temp_ref, temp)
end

function eval_temp_dep(p::TempDepArrhenius, temp)
    p.v_max_ref * arrhenius(p.e_act, p.temp_ref, temp)
end

function eval_temp_dep(p::TempDepEnzyme, temp)
    p.v_max_ref *
    enzyme_temp_dep(p.delta_G_a, p.delta_H_d, p.delta_S_d, p.temp_ref, temp)
end

end  # module


text:
using Test

using StaticArrays

using Polyhedra

include("solvers.jl")

@testset "Redundancy removal" begin
    @testset "LP-based" begin
        @testset "HRepresentation with $(typeof(x))" for (x, y, z, w) in [
            ([1, 0], [0, 1], [1, 1], [0, 0]),
            ([1.0, 0.0], [0.0, 1.0], [1.0, 1.0], [0.0, 0.0]),
            ((@SVector [1, 0]), (@SVector [0, 1]), (@SVector [1, 1]), (@SVector [0, 0]))]
            T = eltype(x)
            AT = typeof(x)
            d = Polyhedra.FullDim(x)
            @testset "HAffineSpace" begin
                for hr in [hrep(typeof(HyperPlane(x, one(T)))[]; d=d),
                           intersect(HyperPlane(w, zero(T)))]
                    @test dim(hr) == 2
                    for h in [detecthlinearity(hr, lp_solver),
                              removehredundancy(hr, lp_solver)]
                        @test dim(h) == 2
                        @test h isa HyperPlanesIntersection{T, AT}
                        @test !hashyperplanes(h)
                        @test !hashalfspaces(h)
                    end
                end
                hr = intersect(HyperPlane(x, one(T)))
                @test dim(hr) == 1
                for h in [detecthlinearity(hr, lp_solver),
                          removehredundancy(hr, lp_solver)]
                    @test dim(h) == 1
                    @test h isa HyperPlanesIntersection{T, AT}
                    @test nhyperplanes(h) == 1
                    @test !hashalfspaces(h)
                end
                hr = intersect(HyperPlane(w, one(T)))
                @test dim(hr) == -1
                for h in [detecthlinearity(hr, lp_solver),
                          removehredundancy(hr, lp_solver)]
                    @test dim(h) == -1
                    @test h isa HyperPlanesIntersection{T, AT}
                    @test nhyperplanes(h) == 3
                    @test !hashalfspaces(h)
                end
            end
            @testset "Intersection" begin
                for hr in [HalfSpace(z + x, 3) ∩ HalfSpace(-x, -one(T)) ∩ HalfSpace(-2y, -2),
                           HalfSpace(x, one(T)) ∩ HalfSpace(y, one(T)) ∩ HalfSpace(-z, -2)]
                    h = detecthlinearity(hr, lp_solver, verbose=2)
                    @test h isa Polyhedra.Intersection{T, AT}
                    @test nhyperplanes(h) == 2
                    @test !hashalfspaces(h)
                end
                for hr in [HalfSpace(z, zero(T)) ∩ HalfSpace(-z, -zero(T)),
                           HalfSpace(x, one(T)) ∩ HalfSpace(-x, -one(T))]
                    h = detecthlinearity(hr, lp_solver, verbose=2)
                    @test h isa Polyhedra.Intersection{T, AT}
                    @test nhyperplanes(h) == 1
                    @test !hashalfspaces(h)
                end
                for hr in [HalfSpace(-x, one(T)) ∩ HalfSpace(-y, one(T)) ∩ HalfSpace(x, one(T)) ∩ HalfSpace(y, one(T)) ∩ HalfSpace(-x, -9),
                           HalfSpace(x, zero(T)) ∩ HalfSpace(-x, -one(T))]
                    h = detecthlinearity(hr, lp_solver, verbose=2)
                    @test h isa Polyhedra.Intersection{T, AT}
                    @test nhyperplanes(h) == 3
                    @test !hashalfspaces(h)
                end
            end
        end
        @testset "VRepresentation with $(typeof(x))" for (x, y, z, w) in [
            ([1, 0], [0, 1], [1, 1], [0, 0]),
            ([1.0, 0.0], [0.0, 1.0], [1.0, 1.0], [0.0, 0.0]),
            ((@SVector [1, 0]), (@SVector [0, 1]), (@SVector [1, 1]), (@SVector [0, 0]))]
            T = eltype(x)
            AT = typeof(x)
            d = Polyhedra.FullDim(x)
            @testset "VLinearSpace" begin
                for vr in [vrep(typeof(Line(x))[]; d=d),
                           convexhull(Line(w))]
                    for v in [detectvlinearity(vr, lp_solver),
                              removevredundancy(vr, lp_solver)]
                        @test v isa LinesHull{T, AT}
                        @test !haspoints(v)
                        @test !hasrays(v)
                        @test !haslines(v)
                    end
                end
                vr = convexhull(Line(x))
                for v in [detectvlinearity(vr, lp_solver),
                          removevredundancy(vr, lp_solver)]
                    @test v isa LinesHull{T, AT}
                    @test npoints(v) == 1
                    @test !hasrays(v)
                    @test nlines(v) == 1
                end
                for vr in [convexhull(Line(x), Line(y)),
                           convexhull(Line(x), Line(y), Line(z)),
                           convexhull(Line(x), Line(z), Line(y)),
                           convexhull(Line(x), Line(z), Line(z))]
                    for v in [detectvlinearity(vr, lp_solver),
                              removevredundancy(vr, lp_solver)]
                        @test v isa LinesHull{T, AT}
                        @test npoints(v) == 1
                        @test !hasrays(v)
                        @test nlines(v) == 2
                    end
                end
                vrepm1 = Polyhedra.emptyspace(vr)
                @test vrepm1 === @inferred removevredundancy(vrepm1, lp_solver)
            end
            @testset "Polyhedra.RaysHull" begin
                vr = vrep(typeof(Line(x))[], typeof(Ray(x))[]; d=d)
                rm = removevredundancy(vr, lp_solver)
                @test rm isa Polyhedra.RaysHull{T, AT}
                @test !haspoints(rm)
                @test !hasrays(rm)
                @test !haslines(rm)
                vr = convexhull(Ray(x))
                rm = removevredundancy(vr, lp_solver)
                @test rm isa Polyhedra.RaysHull{T, AT}
                @test npoints(rm) == 1
                @test nrays(rm) == 1
                @test !haslines(rm)
                vr = convexhull(Ray(x), Ray(y))
                rm = removevredundancy(vr, lp_solver)
                @test rm isa Polyhedra.RaysHull{T, AT}
                @test npoints(rm) == 1
                @test nrays(rm) == 2
                @test !haslines(rm)
                for vr in [convexhull(Ray(x), Ray(y), Ray(z)),
                           convexhull(Ray(x), Ray(z), Ray(y)),
                           convexhull(Ray(z), Ray(x), Ray(y))]
                    rm = removevredundancy(vr, lp_solver)
                    @test rm isa Polyhedra.RaysHull{T, AT}
                    @test npoints(rm) == 1
                    @test nrays(rm) == 2
                    rs = collect(rays(rm))
                    @test rs[1] == Ray(x)
                    @test rs[2] == Ray(y)
                    @test !haslines(rm)
                end
                vr = convexhull(Ray(x), Ray(y), Ray(-z))
                rm = removevredundancy(vr, lp_solver)
                @test rm isa Polyhedra.RaysHull{T, AT}
                @test npoints(rm) == 1
                @test nrays(rm) == 3
                @test !haslines(rm)
                vr = convexhull(Ray(x), Ray(y), Line(z))
                rm = removevredundancy(vr, lp_solver)
                @test rm isa Polyhedra.RaysHull{T, AT}
                @test npoints(rm) == 1
                @test nrays(rm) == 2
                @test nlines(rm) == 1
                vr = convexhull(Ray(x), Ray(-y), Line(z))
                rm = removevredundancy(vr, lp_solver)
                @test rm isa Polyhedra.RaysHull{T, AT}
                @test npoints(rm) == 1
                @test nrays(rm) == 1
                @test nlines(rm) == 1
                vr = conichull(Ray(x), Ray(y), Ray(-z))
                rm = detectvlinearity(vr, lp_solver)
                @test rm isa Polyhedra.RaysHull{T, AT}
                @test npoints(rm) == 1
                @test nrays(rm) == 0
                @test nlines(rm) == 2
            end
            @testset "Polyhedra.PointsHull" begin
                vr = vrep(typeof(x)[]; d=d)
                rm = removevredundancy(vr, lp_solver)
                @test rm isa Polyhedra.PointsHull{T, AT}
                @test !haspoints(rm)
                @test !hasrays(rm)
                @test !haslines(rm)
                vr = convexhull(x)
                rm = removevredundancy(vr, lp_solver)
                @test rm isa Polyhedra.PointsHull{T, AT}
                @test npoints(rm) == 1
                @test !hasrays(rm)
                @test !haslines(rm)
                vr = convexhull(x, y)
                rm = removevredundancy(vr, lp_solver)
                @test rm isa Polyhedra.PointsHull{T, AT}
                @test npoints(rm) == 2
                @test !hasrays(rm)
                @test !haslines(rm)
                for vr in [convexhull(2x, 2y, z),
                           convexhull(2x, z, 2y),
                           convexhull(z, 2x, 2y)]
                    rm = removevredundancy(vr, lp_solver)
                    @test rm isa Polyhedra.PointsHull{T, AT}
                    @test npoints(rm) == 2
                    ps = collect(points(rm))
                    @test ps[1] == 2x
                    @test ps[2] == 2y
                    @test !hasrays(rm)
                    @test !haslines(rm)
                end
                vr = convexhull(x, y, z)
                rm = removevredundancy(vr, lp_solver)
                @test rm isa Polyhedra.PointsHull{T, AT}
                @test npoints(rm) == 3
                @test !hasrays(rm)
                @test !haslines(rm)
                vr = convexhull(x, y, -z)
                rm = removevredundancy(vr, lp_solver)
                @test rm isa Polyhedra.PointsHull{T, AT}
                @test npoints(rm) == 3
                @test !hasrays(rm)
                @test !haslines(rm)
                vr = convexhull(x, y) + Ray(z)
                rm = removevredundancy(vr, lp_solver)
                @test rm isa Polyhedra.Hull{T, AT}
                @test npoints(rm) == 2
                @test nrays(rm) == 1
                @test !haslines(rm)
                vr = convexhull(x, y) + Line(z)
                rm = removevredundancy(vr, lp_solver)
                @test rm isa Polyhedra.Hull{T, AT}
                @test npoints(rm) == 2
                @test !hasrays(rm)
                @test nlines(rm) == 1
                vr = convexhull(x, -y) + Ray(z)
                rm = removevredundancy(vr, lp_solver)
                @test rm isa Polyhedra.Hull{T, AT}
                @test npoints(rm) == 1
                @test first(points(rm)) == -y
                @test nrays(rm) == 1
                @test !haslines(rm)
                vr = convexhull(x, -y) + Line(z)
                rm = removevredundancy(vr, lp_solver)
                @test rm isa Polyhedra.Hull{T, AT}
                @test npoints(rm) == 1
                @test !hasrays(rm)
                @test nlines(rm) == 1
            end
        end
    end
    @testset "Representation-based" begin
        @testset "isredundant" begin
            @test isredundant(HyperPlane([0, 1], 0) ∩ HalfSpace([1, 1], 1), [0, 0], nl=0)
            @test !isredundant(HalfSpace([1, 1], 1) ∩ HyperPlane([0, 1], 0) ∩ HalfSpace([-1, -1], 1), [1, 0], nl=0)
            #@test !ishredundant(convexhull(SymPoint([1, 0])), HyperPlane([0, 1], 0), d=1)
            #@test isredundant(convexhull(SymPoint([1, 0])), HalfSpace([0, 1], 1), d=1)
            @test isredundant(HalfSpace([1, 1], 2) ∩ HalfSpace([-1, 0], 0) ∩ HalfSpace([0, -1], 0), [1, 1], nl=0)
            @test !isredundant(HalfSpace([1, 1], 2) ∩ HalfSpace([-1, 0], 0) ∩ HalfSpace([0, -1], 0), [1, 1], nl=0, strongly=true)
        end
#    @testset "Non-redundant SymPoint" begin
#        vr = convexhull(SymPoint([-1, 0]), SymPoint([0, 1]))
#        hr = HalfSpace([-1, -1], 1) ∩ HalfSpace([1, -1], 1) ∩ HalfSpace([1, 1], 1) ∩ HalfSpace([-1, 1], 1)
#        vr = @inferred Polyhedra.removevredundancy(vr, hr)
#        @test vr isa Polyhedra.PointsHull{2,Int,Vector{Int}}
#        @test collect(sympoints(vr)) == [SymPoint([-1, 0]), SymPoint([0, 1])]
#        @test !haspoints(vr)
#        @test !hasallrays(vr)
#    end
        @testset "Point" begin
            # x - y ≤  1
            # x + y ≥ -1
            hr = HalfSpace([-1, -1], 1) ∩ HalfSpace([1, -1], 1)
            # [-1, 0] is weakly redundant as it belongs to the facet x + y = -1
            @test isredundant(hr, [-1, 0], nl=0)
            @test !isredundant(hr, [-1, 0], strongly=true, nl=0)
            # [0, -1] is not redundant as it is an extreme point
            @test !isredundant(hr, [0, -1], nl=0)
            @test !isredundant(hr, [0, -1], strongly=true, nl=0)
            # [0, 0] is strongly redundant as it is in the relative interior
            @test isredundant(hr, [0, 0], nl=0)
            @test isredundant(hr, [0, 0], strongly=true, nl=0)
            vr = convexhull([-1, 0], [0, -1], [0, 0]) + conichull([1, 1], [-1, 1])
            @test collect(removevredundancy(points(vr), hr, nl=0)) == [[0, -1]]
            @test collect(removevredundancy(points(vr), hr, strongly=true, nl=0)) == [[-1, 0], [0, -1]]
            vrr = Polyhedra.removevredundancy(vr, hr)
            @test collect(points(vrr)) == [[0, -1]]
            @test !haslines(vrr)
            @test collect(rays(vrr)) == [Ray([1, 1]), Ray([-1, 1])]
            vrr = Polyhedra.removevredundancy(vr, hr, strongly=true)
            @test collect(points(vrr)) == [[-1, 0], [0, -1]]
            @test !haslines(vrr)
            @test collect(rays(vrr)) == [Ray([1, 1]), Ray([-1, 1])]
            p = polyhedron(vr, DefaultLibrary{Rational{BigInt}}(Polyhedra.OppositeMockOptimizer))
            Polyhedra.computehrep!(p)
            removevredundancy!(p, strongly=true)
            @test collect(points(p)) == [[-1, 0], [0, -1]]
            @test !haslines(p)
            @test collect(rays(p)) == [Ray([1, 1]), Ray([-1, 1])]
            removevredundancy!(p, planar=false)
            @test collect(points(p)) == [[0, -1]]
            @test !haslines(p)
            @test collect(rays(p)) == [Ray([1, 1]), Ray([-1, 1])]
        end
#    @testset "Split SymPoint" begin
#        vr = convexhull(SymPoint([-1, 0]), SymPoint([0, 1])) + conichull([1, 1], [-1, 1])
#        hr = HalfSpace([-1, -1], 1) ∩ HalfSpace([1, -1], 1)
#        vr = Polyhedra.removevredundancy(vr, hr)
#        @test !hassympoints(vr)
#        @test collect(points(vr)) == [[0, -1]]
#        @test !haslines(vr)
#        @test collect(rays(vr)) == [Ray([1, 1]), Ray([-1, 1])]
#
#        vr = convexhull(SymPoint([1, 0]), SymPoint([0, 1])) + conichull([-1, 1])
#        hr = HalfSpace([-1, -1], 1) ∩ HalfSpace([1, -1], 1) ∩ HalfSpace([1, 1], 1)
#        vr = Polyhedra.removevredundancy(vr, hr)
#        @test !hassympoints(vr)
#        @test collect(points(vr)) == [[1, 0], [0, -1]]
#        @test !haslines(vr)
#        @test collect(rays(vr)) == [Ray([-1, 1])]
#    end
    end
end
@testset "Duplicate removal" begin
    h = removeduplicates(hrep([1 1; -1 -1; 1 0; 0 -1], [1, -1, 1, 0]))
    @test nhyperplanes(h) == 1
    @test first(hyperplanes(h)) == HyperPlane([1, 1], 1)
    @test first(hyperplanes(h)) == HyperPlane([-1, -1], -1)
    @test nhalfspaces(h) == 1
    h = removeduplicates(hrep([1 1; -1 -1], [2, -1]))
    @test !hashyperplanes(h)
    @test nhalfspaces(h) == 2
#    v = removeduplicates(convexhull([1, 2], SymPoint([1, 2]), [-1, 1]))
#    @test typeof(v) == Polyhedra.PointsHull{2,Int,Vector{Int}}
#    @test collect(sympoints(v)) == [SymPoint([1, 2])]
#    @test collect(points(v)) == [[-1, 1]]
#    @test !hasallrays(v)
    for v in (#removeduplicates(convexhull([1, 2], [2, 1], SymPoint([1, 2]), [-1, 1]) + Line([0, 1])),
              #removeduplicates(Line([0, 1]) + convexhull([1, 2], [2, 1], SymPoint([1, 2]), [-1, 1])),
              removeduplicates(Line([0, 1]) + convexhull([-1, -2], [2, 1], [1, 2], [-1, 1], [-1, -2])),)
        @test typeof(v) == Polyhedra.Hull{Int, Vector{Int}, Int}
        #@test collect(sympoints(v)) == [SymPoint([1, 2])]
        @test collect(points(v)) == [[-1, -2], [2, 1], [1, 2]]
        @test !hasrays(v)
        @test collect(lines(v)) == [Line([0, 1])]
    end
end

ext: jl


text:
"""
    ReMat

A representation of the model matrix for a random-effects term
"""
abstract ReMat

"""
    ScalarReMat

The representation of the model matrix for a scalar random-effects term

# Members
* `f`: the grouping factor as a `PooledDataVector`
* `z`: the raw random-effects model matrix as a `Vector`
* `fnm`: the name of the grouping factor as a `Symbol`
* `cnms`: a `Vector` of column names
"""
immutable ScalarReMat{T <: AbstractFloat, S, R <: Integer} <: ReMat
    f::PooledDataVector{S,R}
    z::Vector{T}
    fnm::Symbol
    cnms::Vector
end

"""
    VectorReMat

The representation of the model matrix for a vector-valued random-effects term

# Members
* `f`: the grouping factor as a `PooledDataVector`
* `z`: the transposed raw random-effects model matrix
* `fnm`: the name of the grouping factor as a `Symbol`
* `cnms`: a `Vector` of column names (row names after transposition) of `z`
"""
immutable VectorReMat{T <: AbstractFloat, S, R <: Integer} <: ReMat
    f::PooledDataVector{S,R}
    z::Matrix{T}
    fnm::Symbol
    cnms::Vector
end

"""
     remat(e::Expr, df::DataFrames.DataFrame)

A factory for `ReMat` objects.

`e` should be of the form `:(e1 | e2)` where `e1` is a valid rhs of a `Formula` and
`pool(e2)` can be evaluated within `df`.  The result is a
[`ScalarReMat`](@ref) or a [`VectorReMat`](@ref), as appropriate.
"""
function remat(e::Expr, df::DataFrame)
    e.args[1] == :| || throw(ArgumentError("$e is not a call to '|'"))
    fnm = e.args[3]
    gr = getindex(df, fnm)
    gr = isa(gr,PooledDataArray) ? gr : pool(gr)
    if e.args[2] == 1
        return ScalarReMat(gr, ones(length(gr)), fnm, ["(Intercept)"])
    end
    mf = ModelFrame(Formula(nothing, e.args[2]), df)
    z = ModelMatrix(mf).m
    cnms = coefnames(mf)
    size(z,2) == 1 ? ScalarReMat(gr, vec(z), fnm, cnms) : VectorReMat(gr, z', fnm, cnms)
end

Base.eltype(R::ReMat) = eltype(R.z)

function Base.copy!{S,R}(d::PooledDataVector{S,R}, s::PooledDataVector{S,R})
    copy!(d.pool, s.pool)
    copy!(d.refs, s.refs)
    d
end

function Base.copy!{T,S,R}(d::ScalarReMat{T,S,R}, s::ScalarReMat{T,S,R})
    copy!(d.f, s.f)
    copy!(d.z, s.z)
    d
end

function Base.copy!{T,S,R}(d::VectorReMat{T,S,R}, s::VectorReMat{T,S,R})
    copy!(d.f, s.f)
    copy!(d.z, s.z)
    d
end

"""
    vsize(A::ReMat)

Return the size of vector-valued random effects (i.e. 1 for a [`ScalarReMat`](@ref)).
"""
vsize(A::ReMat) = isa(A,ScalarReMat) ? 1 : size(A.z, 1)

"""
    nlevs(A::ReMat)

Return the number of levels in the grouping factor of `A`.
"""
nlevs(A::ReMat) = length(A.f.pool)

Base.size(A::ReMat) = (length(A.f), vsize(A) * nlevs(A))

Base.size(A::ReMat, i::Integer) =
    i < 1 ? throw(BoundsError()) : i == 1 ? length(A.f) :  i == 2 ? vsize(A)*nlevs(A) : 1

==(A::ReMat,B::ReMat) = (A.f == B.f) && (A.z == B.z)

function Base.A_mul_B!{T}(α::Real, A::ReMat, B::StridedVecOrMat{T}, β::Real, R::StridedVecOrMat{T})
    n,q = size(A)
    k = size(B, 2)
    if size(R, 1) ≠ n || size(B, 1) ≠ q || size(R, 2) ≠ k
        throw(DimensionMismatch())
    end
    if β ≠ 1
        β == 0 ? fill!(R, 0) : scale!(β, R)
    end
    rr, zz = A.f.refs, A.z
    if isa(A, ScalarReMat)
        for j in 1 : k, i in 1 : n
            R[i, j] += α * zz[i] * B[rr[i],j]
        end
    else
        l = size(zz,1)
        Bt = reshape(B, (l, div(q,l), k))
        for j in 1:k, i in 1:n
            R[i, j] += α * dot(Compat.view(Bt, :, Int(rr[i]), j), Compat.view(zz, :, i))
        end
    end
    R
end

Base.A_mul_B!{T}(A::ReMat, B::StridedVecOrMat{T}, R::StridedVecOrMat{T}) = A_mul_B!(one(T), A, B, zero(T), R)

function Base.Ac_mul_B!{T}(α::Real, A::ReMat, B::StridedVecOrMat{T}, β::Real, R::StridedVecOrMat{T})
    n, q = size(A)
    k = size(B, 2)
    if size(R, 1) ≠ q || size(B, 1) ≠ n || size(R, 2) ≠ k
        throw(DimensionMismatch())
    end
    if β ≠ 1
        β == 0 ? fill!(R, 0) : scale!(β, R)
    end
    rr, zz = A.f.refs, A.z
    if isa(A, ScalarReMat)
        for j in 1 : k, i in 1 : n
            R[rr[i], j] += α * zz[i] * B[i, j]
        end
    else
        l = size(zz, 1)
        for j in 1 : k, i in 1 : n
            roffset = (rr[i] - 1) * l
            mul = α * B[i, j]
            for ii in 1 : l
                R[roffset + ii, j] += mul * zz[ii, i]
            end
        end
    end
    R
end

Base.Ac_mul_B!{T}(R::StridedVecOrMat{T}, A::ReMat, B::StridedVecOrMat{T}) = Ac_mul_B!(one(T), A, B, zero(T), R)

function Base.Ac_mul_B(A::ReMat, B::DenseVecOrMat)
    k = size(A, 2)
    Ac_mul_B!(zeros(eltype(B), isa(B, Vector) ? (k,) : (k, size(B, 2))), A, B)
end

function Base.Ac_mul_B(A::ScalarReMat, B::ScalarReMat)
    Az = A.z
    Ar = A.f.refs
    if is(A, B)
        v = zeros(eltype(Az), nlevs(A))
        for i in eachindex(Ar)
            v[Ar[i]] += abs2(Az[i])
        end
        return Diagonal(v)
    end
    densify(sparse(convert(Vector{Int32}, Ar), convert(Vector{Int32}, B.f.refs), Az .* B.z))
end

function Base.Ac_mul_B!{T}(C::Diagonal{T}, A::ScalarReMat{T}, B::ScalarReMat{T})
    c, a, r, b = C.diag, A.z, A.f.refs, B.z
    if r ≠ B.f.refs
        throw(ArgumentError("A'B is not diagonal"))
    end
    fill!(c, 0)
    for i in eachindex(a)
        c[r[i]] += a[i] * b[i]
    end
    C
end

function Base.Ac_mul_B!{Tv, Ti}(C::SparseMatrixCSC{Tv, Ti}, A::ScalarReMat{Tv}, B::ScalarReMat{Tv})
    m, n = size(A)
    p, q = size(B)
    if m ≠ p || size(C, 1) ≠ n || size(C, 2) ≠ q
        throw(DimensionMismatch("size(A) = $(size(A)), size(B) = $(size(B)), size(C) = $(size(C))"))
    end
    SparseArrays.sparse!(convert(Vector{Ti}, A.f.refs), convert(Vector{Ti}, B.f.refs), A.z .* B.z,
        n, q, +, Array(Ti, q), Array(Ti, n + 1), Array(Ti, m), Array(Tv, m), C.colptr, C.rowval, C.nzval)
    C
end

function Base.Ac_mul_B!{T}(C::HBlkDiag{T}, A::VectorReMat{T}, B::VectorReMat{T})
    c, a, r = C.arr, A.z, A.f.refs
    _, m, n = size(c)
    fill!(c, 0)
    if !is(A, B)
        throw(ArgumentError("Currently defined only for A === B"))
    end
    for k in eachindex(r)
        ri = Int(r[k])
        for j in 1 : m
            aj = a[j, k]
            c[j, j, ri] += abs2(aj)
            for i in 1 : j - 1
                aij = a[i, k] * aj
                c[i, j, ri] += aij
                c[j, i, ri] += aij
            end
        end
    end
    C
end

function Base.Ac_mul_B!{T}(C::Matrix{T}, A::ScalarReMat{T}, B::ScalarReMat{T})
    m, n = size(C)
    ma, na = size(A)
    mb, nb = size(B)
    if m ≠ na || n ≠ nb || ma ≠ mb
        throw(DimensionMismatch())
    end
    a, b, ra, rb = A.z, B.z, A.f.refs, B.f.refs
    fill!(C, 0)
    for i in eachindex(a)
        C[ra[i], rb[i]] += a[i] * b[i]
    end
    C
end

function Base.Ac_mul_B{T}(A::VectorReMat{T}, B::VectorReMat{T})
    if is(A, B)
        l = size(A.z, 1)
        return Ac_mul_B!(HBlkDiag(Array(T, (l, l, length(A.f.pool)))), A, B)
    end
    Az = A.z
    Bz = B.z
    if (m = size(Az, 2)) ≠ size(Bz, 2)
        throw(DimensionMismatch("$m = size(Az,2) ≠ size(Bz,2) = $(size(Bz, 2))"))
    end
    a, b = size(Az,1), size(Bz, 1)
    ab = a * b
    nz = ab * m
    I, J, V = sizehint!(Int[], nz), sizehint!(Int[], nz), sizehint!(T[], nz)
    Ar, Br = A.f.refs, B.f.refs
    Ipat = Compat.repeat(1 : a, outer = b)
    Jpat = Compat.repeat(1 : b, inner = a)
    for i in 1 : m
        append!(I, Ipat + (Ar[i] - 1) * a)
        append!(J, Jpat + (Br[i] - 1) * b)
        append!(V, vec(Compat.view(Az, :, i) * Compat.view(Bz, :, i)'))
    end
    sparse(I, J, V)
end

function Base.Ac_mul_B!{T}(R::DenseVecOrMat{T}, A::DenseVecOrMat{T}, B::ReMat)
    m = size(A, 1)
    n = size(A, 2)
    p, q = size(B)
    if m ≠ p || size(R, 1) ≠ n || size(R, 2) ≠ q
        throw(DimensionMismatch(""))
    end
    fill!(R, 0)
    r, z = B.f.refs, B.z
    if isa(B, ScalarReMat)
        for j in 1 : n, i in 1 : m
            R[j, r[i]] += A[i, j] * z[i]
        end
    else
        l = size(zz, 1)
        for j in 1 : n, i in 1 : m
            roffset = (r[i] - 1) * l
            aij = A[i, j]
            for k in 1:l
                R[j, roffset + k] += aij * z[k, i]
            end
        end
    end
    R
end

Base.Ac_mul_B(A::DenseVecOrMat, B::ReMat) = Ac_mul_B!(Array(eltype(A), (size(A, 2), size(B, 2))), A, B)

(*){T}(D::Diagonal{T}, A::ScalarReMat{T}) = ScalarReMat(A.f, D * A.z, A.fnm, A.cnms)

function Base.A_mul_B!{T}(C::ScalarReMat{T}, A::Diagonal{T}, B::ScalarReMat{T})
    map!(*, C.z, A.diag, B.z)
    C
end

function Base.A_mul_B!{T}(A::Diagonal{T}, B::ScalarReMat{T})
    a, b = A.diag, B.z
    if length(a) ≠ length(b)
        throw(DimensionMismatch("A_mul_B!, A: diagonal $(size(A, 1)), B: ScalarReMat $(size(B))"))
    end
    for i in eachindex(a)
        b[i] *= a[i]
    end
    B
end

function Base.A_mul_B!{T}(A::Diagonal{T}, B::VectorReMat{T})
    a, b = A.diag, B.z
    if length(a) ≠ size(b, 2)
        throw(DimensionMismatch("A_mul_B!, A: diagonal $(size(A, 1)), B: ScalarReMat $(size(B))"))
    end
    k = size(b, 1)
    for j in eachindex(a)
        aj = a[j]
        for i in 1 : k
            b[i, j] *= aj
        end
    end
    B
end

(*){T}(D::Diagonal{T}, A::VectorReMat{T}) = VectorReMat(A.f, A.z * D, A.fnm, A.cnms)

function Base.A_mul_B!{T}(C::VectorReMat{T}, A::Diagonal{T}, B::VectorReMat{T})
    A_mul_B!(C.z, B.z, A)
    C
end


text:
# randn_in and rand_out
# =====================================

function randn_in(tm::TM) where {Tf<:Real, TM<:𝕎{Tf}}
    wx = randn(eltype_in(tm), size_in(tm)) ./ sqrt.(Ωpix(tm))
    wx
end

# This needs testing
function randn_out(tm::TM) where {Tf<:Complex, TM<:𝕎{Tf}}
    wk = randn(eltype_in(tm), size_in(tm)) ./ sqrt.(Ωfreq(tm))
    wk
end

# dot_in and dot_out
# =====================================

function dot_in(f::Xfield{FT}, g::Xfield{FT}) where {T<:Real, FT<:𝕎{T}}
    tm     = fieldtransform(f)
    fdata  = fielddata(MapField(f))
    gdata  = fielddata(MapField(g))
    Ω      = Ωpix(tm)
    return  sum_kbn(fdata .* gdata .* Ω)
end

function dot_in(f::Xfield{FT}, g::Xfield{FT}) where {T<:Complex, FT<:𝕎{T}}
    tm     = fieldtransform(f)
    fdata  = fielddata(MapField(f))
    gdata  = fielddata(MapField(g))
    Ω      = Ωpix(tm)
    return  sum_kbn(fdata .* conj.(gdata) .* Ω)
end

# TODO:
# function dot_out(f::Xfield{FT}, g::Xfield{FT}) where {T<:Real, FT<:𝕎{Complex{T}}}
#     tm     = fieldtransform(f)
#     fdata  = fielddata(FourierField(f))
#     gdata  = fielddata(FourierField(g))
#     Ω      = Ωfreq(tm)
#     return  sum_kbn(fdata .* conj.(gdata) .* Ω)
# end

# copied from https://github.com/JuliaMath/KahanSummation.jl
function sum_kbn(A)
    T = Base.@default_eltype(A)
    c = Base.reduce_empty(+, T)
    it = iterate(A)
    it === nothing && return c
    Ai, i = it
    s = Ai - c
    while (it = iterate(A, i)) !== nothing
        Ai, i = it
        t = s + Ai
        if abs(s) >= abs(Ai)
            c -= ((s-t) + Ai)
        else
            c -= ((Ai-t) + s)
        end
        s = t
    end
    s - c
end


text:
using Mimi
using Test

@testset "Mimi-MCS" begin

    @info("test_empirical.jl")
    include("test_empirical.jl")

    @info("test_defmcs.jl")
    include("test_defmcs.jl")

end


text:
using StanSample

#StanBase.set_cmdstan_home!(homedir() * "/Projects/StanSupport/cmdstan_stanc3")

ProjDir = @__DIR__
cd(ProjDir)

bernoulli_model = "
  functions{
    #include model_specific_funcs.stan
    #include shared_funcs.stan // a comment
    //#include shared_funcs.stan // a comment
  }
  data {
    int<lower=1> N;
    int<lower=0,upper=1> y[N];
  }
  parameters {
    real<lower=0,upper=1> theta;
  }
  model {
    model_specific_function();
    theta ~ beta(my_function(),1);
    y ~ bernoulli(theta);
  }
";

tmpdir = joinpath(ProjDir, "tmp")
stanmodel = SampleModel("bernoulli", bernoulli_model)

observeddata = Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1])

(sample_file, log_file) = stan_sample(stanmodel, data=observeddata)

if sample_file !== Nothing
  # Convert to an MCMCChains.Chains object
  chns = read_samples(stanmodel)

  # Describe the MCMCChains using MCMCChains statistics
  cdf = describe(chns)
  display(cdf)

  # Show cmdstan summary in DataFrame format
  sdf = read_summary(stanmodel)
  display(sdf)
  println()

  # Retrieve mean value of theta from the summary
  sdf[:theta, :mean]
end


text:
# Julia wrapper for header: /usr/include/scip/scip_bandit.h
# Automatically generated using Clang.jl wrap_c

function SCIPincludeBanditvtable(scip, banditvtable, name, banditfree, banditselect, banditupdate, banditreset)
    ccall((:SCIPincludeBanditvtable, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_BANDITVTABLE}}, Cstring, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}), scip, banditvtable, name, banditfree, banditselect, banditupdate, banditreset)
end

function SCIPfindBanditvtable(scip, name)
    ccall((:SCIPfindBanditvtable, libscip), Ptr{SCIP_BANDITVTABLE}, (Ptr{SCIP_}, Cstring), scip, name)
end

function SCIPfreeBandit(scip, bandit)
    ccall((:SCIPfreeBandit, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_BANDIT}}), scip, bandit)
end

function SCIPresetBandit(scip, bandit, priorities, seed)
    ccall((:SCIPresetBandit, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_BANDIT}, Ptr{Cdouble}, UInt32), scip, bandit, priorities, seed)
end


text:
ENV["GKS_ENCODING"]="utf-8"

using Test, Flux, NeuralNetDiffEq, Plots, Plots.PlotMeasures
using DiffEqDevTools

# Run a solve on vectors

linear = (u,p,t) -> [cos(2pi*t)]
tspan = (0.0f0, 1.0f0)
u0 = 0.0f0
prob = ODEProblem(linear, u0 ,tspan)
chain = Flux.Chain(Dense(1,5,σ),Dense(5,1))
opt = Flux.ADAM(0.1, (0.9, 0.95))
sol = solve(prob, NeuralNetDiffEq.NNODE(chain,opt), initial_condition = [u0], dt=1/20f0, abstol=1e-10,
           verbose = true, maxiters=200)

# Second-Order ODE - Dirichlet
linear = (u,p,t) -> cos(2pi*t)
tspan = (0.0f0, 1.0f0)
u0 = 0.0f0
u1 = 1.0f0
initial_condition = [u0, u1]
prob = ODEProblem(linear, u0, tspan)

chain = Flux.Chain(Dense(1,5,σ),Dense(5,1))
opt = Flux.ADAM(0.1, (0.9, 0.95))
sol = solve(prob, NeuralNetDiffEq.NNODE(chain,opt), initial_condition = [u0, u1], dt=1/20f0, verbose = true,
            abstol=1e-10, maxiters = 200)

analytical_sol = t -> ((4pi^2 * t) - cos(2pi*t) + 1) / (4 * pi^2)
plot(analytical_sol, label = "analytical", linestyle = :dot, linewidth = 8, legend=:bottomright, bottom_margin = 10mm, left_margin=10mm, top_margin=10mm, right_margin=10mm)
plot!(sol, title = "Second-Order Dirichlet : Cos(2𝝅t)", label = "phi", fillalpha = .5, linewidth = 4)
xlabel!("Time (t)")
display(ylabel!("f(x)", ylims = (0,1)))
savefig("plot.png")

#=
anat_sum = zeros(0)
append!(anat_sum, 0)
for i = 1:20
    k = i * .05
    append!(anat_sum, analytical_sol(k))
end

print("length of sol is: ", length(sol))
print("length of anat_sum is:", length(anat_sum))
#error = sol .- analytical_sol
#print("Error: ", error)
#@test error < 0.5

=#

# Run a solve on vectors
linear = (u,p,t) -> [cos(2pi*t)]
tspan = (0.0f0, 1.0f0)
u0 = [0.0f0]
prob = ODEProblem(linear, u0 ,tspan)
chain = Flux.Chain(Dense(1,5,σ),Dense(5,1))
opt = Flux.ADAM(0.1, (0.9, 0.95))
sol = solve(prob, NeuralNetDiffEq.NNODE(chain,opt), initial_condition = [u0], dt=1/20f0, abstol=1e-10,
            verbose = true, maxiters=200)

#Example 1
linear = (u,p,t) -> @. t^3 + 2*t + (t^2)*((1+3*(t^2))/(1+t+(t^3))) - u*(t + ((1+3*(t^2))/(1+t+t^3)))
linear_analytic = (u0,p,t) -> [exp(-(t^2)/2)/(1+t+t^3) + t^2]
prob = ODEProblem(ODEFunction(linear,analytic=linear_analytic),[1f0],(0.0f0,1.0f0))
chain = Flux.Chain(Dense(1,5,σ),Dense(5,1))
opt = Flux.ADAM(0.1, (0.9, 0.95))
sol  = solve(prob,NeuralNetDiffEq.NNODE(chain,opt),initial_condition = [prob.u0], verbose = true, dt=1/5f0)
err = sol.errors[:l2]
sol  = solve(prob,NeuralNetDiffEq.NNODE(chain,opt),initial_condition = [prob.u0], verbose = true, dt=1/20f0)
sol.errors[:l2]/err < 0.5

#=
dts = 1f0 ./ 2f0 .^ (6:-1:2)
sim = test_convergence(dts, prob, NeuralNetDiffEq.NNODE(chain, opt))
@test abs(sim.𝒪est[:l2]) < 0.1
@test minimum(sim.errors[:l2]) < 0.5
=#

#Example 2
linear = (u,p,t) -> -u/5 + exp(-t/5).*cos(t)
linear_analytic = (u0,p,t) ->  exp(-t/5)*(u0 + sin(t))
prob = ODEProblem(ODEFunction(linear,analytic=linear_analytic),0.0f0,(0.0f0,1.0f0))
chain = Flux.Chain(Dense(1,5,σ),Dense(5,1))
opt = Flux.ADAM(0.1, (0.9, 0.95))
sol  = solve(prob,NeuralNetDiffEq.NNODE(chain,opt), initial_condition = [prob.u0], verbose = true, dt=1/5f0)
err = sol.errors[:l2]
sol  = solve(prob,NeuralNetDiffEq.NNODE(chain,opt),initial_condition = [prob.u0], verbose = true, dt=1/20f0)
sol.errors[:l2]/err < 0.5

#=
dts = 1f0 ./ 2f0 .^ (6:-1:2)
sim = test_convergence(dts, prob, NeuralNetDiffEq.NNODE(chain, opt))
@test abs(sim.𝒪est[:l2]) < 0.5
@test minimum(sim.errors[:l2]) < 0.1
=#


text:
using Blink
using Test

# IMPORTANT: Window(...) cannot appear inside of a @testset for as-of-yet
# unknown reasons.
w = Window(Blink.@d(:show => false), async=false);
@testset "content! Tests" begin
    body!(w, "", async=false);
    @test (@js w document.querySelector("body").innerHTML) == ""

    # Test reloading body and a selector element.
    html = """
hi world

bye
"""
    body!(w, html, async=false);
    @test (@js w document.getElementById("a").innerHTML) == "hi world"
    @test (@js w document.getElementById("b").innerHTML) == "bye"
    content!(w, "div", "hello world", async=false);
    @test (@js w document.getElementById("a").innerHTML) == "hello world"
    # TODO(nhdaly): Is this right? Should content!(w,"div",...) change _all_ divs?
    @test (@js w document.getElementById("b").innerHTML) == "bye"
end

# Test `fade` parameter and scripts:
# Must create a new window to ensure javascript is reset.
w = Window(Blink.@d(:show => false), async=false);
fadeTestHtml = """
hi world
"""
@testset "Fade True" begin
    body!(w, fadeTestHtml; fade=true, async=false);
    @test (@js w testJS) == "test"
end

# Must create a new window to ensure javascript is reset.
w = Window(Blink.@d(:show => false), async=false);
@testset "Fade False" begin

    body!(w, fadeTestHtml; fade=false, async=false);
    @test (@js w testJS) == "test"
end

w = Window(Blink.@d(:show => false), async=false);
@testset "Sync/Async content reload tests" begin
    sleep_content(seconds) = """

      """

    @timed sleep(0.1);   # Throw-away statement to warm-up @sync and @async

    x, t = @timed body!(w, sleep_content(3); fade=true, async=false)
    #@test x == true  # TODO: What should it return?
    @test t >= 3.0 # seconds

    x, t = @timed body!(w, sleep_content(3); fade=false, async=false)
    @test t >= 3.0 # seconds

    x, t = @timed body!(w, sleep_content(3); fade=true, async=true);
    @test t < 3.0 # seconds
    sleep(3)  # (Wait until the end of the previous body! call.)

    x, t = @timed body!(w, sleep_content(3); fade=false, async=true);
    @test t < 3.0 # seconds
    sleep(3)  # (Wait until the end of the previous body! call.)

    @sync begin  # Throw-away block to warm-up @sync and @async
        @async sleep(0.1)
        @async sleep(0.1)
    end
    # Test using Julia's async mechanisms with synchronous `content!`.
    _, t = @timed @sync begin
      @async body!(w, sleep_content(4); async=false);
      sleep(4)
    end

    @test t >= 4.0
    @test t < 8
end


text:
@testset "non-mutating operator without keyward parameters" begin

    for Tv in [Float32, Float64, ComplexF32, ComplexF64]

        # size of the linear operator
        m = 33; n = 17;
        A = randn(Tv, m, n)

        # generate mutating forward operator
        function forward_generator(A)

            g = function (x)
                    y = A * x
                    return y
                end
            return g
        end
        f1 = forward_generator(A)

        # generate mutating adjoint operator
        function adjoint_generator(A)

            g = function (x)
                    y = A' * x
                    return y
                end
            return g
        end
        f2 = adjoint_generator(A)

        # create a linear operator
        L = LinearOperator(Tv, f1, f2, m, n)

        # apply forward operator to input x
        x = rand(Tv, n)
        y = Vector{Tv}(undef, m)
        y1= Vector{Tv}(undef, m)
        A_mul_b!(y, L, x)
        y1 = L * x
        y2 = A * x
        @test isapprox(y, y1)
        @test isapprox(y, y2)

        # apply adjoint operator to input x
        x = randn(Tv, m)
        y = Vector{Tv}(undef, n)
        Ac_mul_b!(y, L, x)
        y1 = A' * x
        @test isapprox(y, y1)

        # transpose opertor
        At_mul_b!(y, L, x)
        mul!(y1, transpose(A), x)
        @test isapprox(y, y1)

    end
end

@testset "non-mutating operator with keyword parameters" begin

    for Tv in [Float32, Float64, ComplexF32, ComplexF64]

        # size of the linear operator
        m = 33; n = 17;
        A = randn(Tv, m, n)

        # define forward operator which support mutating
        function forward_op(x; B=0)
            y = B * x
            return y
        end

        # define adjoint operator which support mutating
        function adjoint_op(x; B=0)
            y = B' * x
            return y
        end

        # create a linear operator
        L = LinearOperator(Tv, forward_op, adjoint_op, m, n; forward_params=(B=A,), adjoint_params=(B=A,))

        # apply forward operator to input x
        x = randn(Tv, n)
        y = Vector{Tv}(undef, m)
        y1= Vector{Tv}(undef, m)
        A_mul_b!(y, L, x)
        y1 = L * x
        y2 = A * x
        @test isapprox(y, y1)
        @test isapprox(y, y2)

        # apply adjoint operator to input x
        x = randn(Tv, m)
        y = Vector{Tv}(undef, n)
        Ac_mul_b!(y, L, x)
        y1 = A' * x
        @test isapprox(y, y1)

        # transpose opertor
        At_mul_b!(y, L, x)
        mul!(y1, transpose(A), x)
        @test isapprox(y, y1)

    end
end

@testset "mutating operator without keyward parameters" begin

    for Tv in [Float32, Float64, ComplexF32, ComplexF64]

        # size of the linear operator
        m = 127; n = 129;
        A = randn(Tv, m, n)

        # generate mutating forward operator
        function forward_generator(A)

            g = function (y, x)
                    mul!(y, A, x)
                    return nothing
                end
            return g
        end
        f1 = forward_generator(A)

        # generate mutating adjoint operator
        function adjoint_generator(A)

            g = function (y, x)
                    mul!(y, adjoint(A), x)
                    return nothing
                end
            return g
        end
        f2 = adjoint_generator(A)

        # create a linear operator
        L = LinearOperator(Tv, f1, f2, m, n)

        # apply forward operator to input x
        x = rand(Tv, n)
        y = Vector{Tv}(undef, m)
        y1= Vector{Tv}(undef, m)
        A_mul_b!(y, L, x)
        y1 = L * x
        y2 = A * x
        @test isapprox(y, y1)
        @test isapprox(y, y2)

        # apply adjoint operator to input x
        x = randn(Tv, m)
        y = Vector{Tv}(undef, n)
        Ac_mul_b!(y, L, x)
        y1 = A' * x
        @test isapprox(y, y1)

        # transpose opertor
        At_mul_b!(y, L, x)
        mul!(y1, transpose(A),x)
        @test isapprox(y, y1)

    end
end

@testset "mutating operator with keyward parameters" begin

    for Tv in [Float32, Float64, ComplexF32, ComplexF64]

        # size of the linear operator
        m = 100; n = 90;
        A = randn(Tv, m, n)

        # define forward operator which support mutating
        function forward_op!(y, x; B=0)
            mul!(y, B, x)
            return nothing
        end

        # define adjoint operator which support mutating
        function adjoint_op!(y, x; B=0)
            mul!(y, adjoint(B), x)
            return nothing
        end

        # create a linear operator
        L = LinearOperator(Tv, forward_op!, adjoint_op!, m, n; forward_params=(B=A,), adjoint_params=(B=A,))

        # apply forward operator to input x
        x = rand(Tv, n)
        y = Vector{Tv}(undef, m)
        y1= Vector{Tv}(undef, m)
        A_mul_b!(y, L, x)
        y1 = L * x
        y2 = A * x
        @test isapprox(y, y1)
        @test isapprox(y, y2)

        # apply adjoint operator to input x
        x = randn(Tv, m)
        y = Vector{Tv}(undef, n)
        Ac_mul_b!(y, L, x)
        y1 = A' * x
        @test isapprox(y, y1)

        # transpose opertor
        At_mul_b!(y, L, x)
        mul!(y1, transpose(A),x)
        @test isapprox(y, y1)

    end
end


text:
function heatbath_SU2!(U,NC,temps,β,Dim=4)

    temp1 = temps[1]
    temp2 = temps[2]
    V = temps[3]
    ITERATION_MAX = 10^5

    temps2 = Array{Matrix{ComplexF64},1}(undef,5)
    for i=1:5
        temps2[i] = zeros(ComplexF64,2,2)
    end

    mapfunc!(A,B) = SU2update_KP!(A,B,β,NC,temps2,ITERATION_MAX)

    for μ=1:Dim

        loops = loops_staple[(Dim,μ)]
        iseven = true

        evaluate_gaugelinks_evenodd!(V,loops,U,[temp1,temp2],iseven)
        map_U!(U[μ],mapfunc!,V,iseven)

        iseven = false
        evaluate_gaugelinks_evenodd!(V,loops,U,[temp1,temp2],iseven)
        map_U!(U[μ],mapfunc!,V,iseven)
    end

end

function heatbath_SU3!(U,NC,temps,β,Dim=4)
    temp1 = temps[1]
    temp2 = temps[2]
    V = temps[3]
    ITERATION_MAX = 10^5

    temps2 = Array{Matrix{ComplexF64},1}(undef,5)
    temps3 = Array{Matrix{ComplexF64},1}(undef,5)
    for i=1:5
        temps2[i] = zeros(ComplexF64,2,2)
        temps3[i] = zeros(ComplexF64,3,3)
    end

    mapfunc!(A,B) = SU3update_matrix!(A,B,β,NC,temps2,temps3,ITERATION_MAX)

    for μ=1:Dim

        loops = loops_staple[(Dim,μ)]
        iseven = true

        evaluate_gaugelinks_evenodd!(V,loops,U,[temp1,temp2],iseven)
        map_U!(U[μ],mapfunc!,V,iseven)

        iseven = false
        evaluate_gaugelinks_evenodd!(V,loops,U,[temp1,temp2],iseven)
        map_U!(U[μ],mapfunc!,V,iseven)
    end

end

function heatbath_SUN!(U,NC,temps,β,Dim = 4)
    #Dim = 4
    temp1 = temps[1]
    temp2 = temps[2]
    V = temps[3]
    ITERATION_MAX = 10^5

    temps2 = Array{Matrix{ComplexF64},1}(undef,5)
    temps3 = Array{Matrix{ComplexF64},1}(undef,5)
    for i=1:5
        temps2[i] = zeros(ComplexF64,2,2)
        temps3[i] = zeros(ComplexF64,NC,NC)
    end

    mapfunc!(A,B) = SUNupdate_matrix!(A,B,β,NC,temps2,temps3,ITERATION_MAX)

    for μ=1:Dim

        loops = loops_staple[(Dim,μ)]
        iseven = true

        evaluate_gaugelinks_evenodd!(V,loops,U,[temp1,temp2],iseven)
        map_U!(U[μ],mapfunc!,V,iseven)

        iseven = false
        evaluate_gaugelinks_evenodd!(V,loops,U,[temp1,temp2],iseven)
        map_U!(U[μ],mapfunc!,V,iseven)
    end

end

function heatbathtest_4D(NX,NY,NZ,NT,β,NC)
    Dim = 4
    Nwing = 1

    #=
    u1 = IdentityGauges(NC,Nwing,NX,NY,NZ,NT)
    U = Array{typeof(u1),1}(undef,Dim)
    U[1] = u1
    for μ=2:Dim
        U[μ] = IdentityGauges(NC,Nwing,NX,NY,NZ,NT)
    end
    =#

    U = Initialize_Gaugefields(NC,Nwing,NX,NY,NZ,NT,condition = "cold")

    h = Heatbath(U,β)

    #=
    gauge_action = GaugeAction(U)
    plaqloop = make_loops_fromname("plaquette",Dim=Dim)
    append!(plaqloop,plaqloop')
    βinp = β/2
    push!(gauge_action,βinp,plaqloop)
    =#

    temp1 = similar(U[1])
    temp2 = similar(U[1])
    temp3 = similar(U[1])

    comb = 6
    factor = 1/(comb*U[1].NV*U[1].NC)
    @time plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
    println("plaq_t = $plaq_t")
    poly = calculate_Polyakov_loop(U,temp1,temp2)
    println("polyakov loop = $(real(poly)) $(imag(poly))")
    #hnew = Heatbath_update(U,gauge_action)

    plaq_ave = 0.0
    numhb = 200
    for itrj = 1:numhb
        #heatbath!(U,hnew)
        heatbath!(U,h)
        #=
        if NC == 2
            heatbath_SU2!(U,NC,[temp1,temp2,temp3],β)
        elseif NC == 3
            heatbath_SU3!(U,NC,[temp1,temp2,temp3],β)
        else
            heatbath_SUN!(U,NC,[temp1,temp2,temp3],β)
        end
        =#
        plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
        plaq_ave += plaq_t

        if itrj % 40 == 0
            #@time plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
            println("$itrj plaq_t = $plaq_t average: $(plaq_ave/itrj)")
            #println("$itrj plaq_t = $plaq_t")
            poly = calculate_Polyakov_loop(U,temp1,temp2)
            println("$itrj polyakov loop = $(real(poly)) $(imag(poly))")
        end
        #=

        if itrj % 40 == 0
            @time plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
            println("$itrj plaq_t = $plaq_t")
            poly = calculate_Polyakov_loop(U,temp1,temp2)
            println("$itrj polyakov loop = $(real(poly)) $(imag(poly))")
        end
        =#
    end

    return plaq_ave/numhb

end

function heatbathtest_2D(NX,NT,β,NC)
    Dim = 2
    Nwing = 1

    #=
    u1 = RandomGauges(NC,Nwing,NX,NT)
    U = Array{typeof(u1),1}(undef,Dim)
    U[1] = u1
    for μ=2:Dim
        U[μ] = RandomGauges(NC,Nwing,NX,NT)
    end
    =#

    U = Initialize_Gaugefields(NC,Nwing,NX,NT,condition = "hot",randomnumber="Reproducible")

    temp1 = similar(U[1])
    temp2 = similar(U[1])
    temp3 = similar(U[1])

    #comb = 6
    if Dim == 4
        comb = 6 #4*3/2
    elseif Dim == 3
        comb = 3
    elseif Dim == 2
        comb = 1
    else
        error("dimension $Dim is not supported")
    end

    factor = 1/(comb*U[1].NV*U[1].NC)
    @time plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
    println("plaq_t = $plaq_t")
    poly = calculate_Polyakov_loop(U,temp1,temp2)
    println("polyakov loop = $(real(poly)) $(imag(poly))")

    #=
    gauge_action = GaugeAction(U)
    plaqloop = make_loops_fromname("plaquette",Dim=Dim)
    append!(plaqloop,plaqloop')
    βinp = β/2
    push!(gauge_action,βinp,plaqloop)
    hnew = Heatbath_update(U,gauge_action)

    =#

    numhb = 200
    plaq_ave = 0.0
    for itrj = 1:numhb
        #heatbath!(U,hnew)
        heatbath!(U,[temp1,temp2,temp3],β)
        #=
        if NC == 2
            heatbath_SU2!(U,NC,[temp1,temp2,temp3],β,Dim)
        elseif NC == 3
            heatbath_SU3!(U,NC,[temp1,temp2,temp3],β,Dim)
        else
            heatbath_SUN!(U,NC,[temp1,temp2,temp3],β,Dim)
        end
        =#
        plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
        plaq_ave += plaq_t

        if itrj % 40 == 0
            #@time plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
            println("$itrj plaq_t = $plaq_t average: $(plaq_ave/itrj)")
            #println("$itrj plaq_t = $plaq_t")
            poly = calculate_Polyakov_loop(U,temp1,temp2)
            println("$itrj polyakov loop = $(real(poly)) $(imag(poly))")
        end

        #=
        if itrj % 40 == 0
            @time plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
            println("$itrj plaq_t = $plaq_t")
            poly = calculate_Polyakov_loop(U,temp1,temp2)
            println("$itrj polyakov loop = $(real(poly)) $(imag(poly))")
        end
        =#
    end

    return plaq_ave/numhb

end

#eps = 1e-1

println("4D system")
@testset "4D" begin
    NX = 4
    NY = 4
    NZ = 4
    NT = 4
    Nwing = 1

    @testset "NC=2" begin
        β = 2.3
        NC = 2
        println("NC = $NC")
        #val =0.6414596466929057
        val = 0.6042874193905048
        #@time plaq_t = heatbathtest_2D(NX,NT,β,NC)
        @time plaq_ave = heatbathtest_4D(NX,NY,NZ,NT,β,NC)
        @test abs(plaq_ave-val)/abs(val) < eps

    end

    @testset "NC=3" begin
        β = 5.7
        NC = 3
        println("NC = $NC")
        #val = 0.5779454661484242
        val = 0.5616169101071591
        #@time plaq_t = heatbathtest_2D(NX,NT,β,NC)
        @time plaq_ave = heatbathtest_4D(NX,NY,NZ,NT,β,NC)
        @test abs(plaq_ave-val)/abs(val) < eps
    end

    @testset "NC=4" begin
        β = 5.7
        NC = 4
        println("NC = $NC")
        #val  =0.19127260002797497
        val = 0.1967198548214144
        #@time plaq_t = heatbathtest_2D(NX,NT,β,NC)
        @time plaq_ave = heatbathtest_4D(NX,NY,NZ,NT,β,NC)
        @test abs(plaq_ave-val)/abs(val) < eps
    end

end

println("2D system")
@testset "2D" begin
    NX = 4
    #NY = 4
    #NZ = 4
    NT = 4
    Nwing = 1

    @testset "NC=2" begin
        β = 2.3
        NC = 2
        println("NC = $NC")
        val = 0.47007878197368624
        #@time plaq_t = heatbathtest_2D(NX,NT,β,NC)
        @time plaq_ave = heatbathtest_2D(NX,NT,β,NC)
        @test abs(plaq_ave-val)/abs(val) < eps
    end

    @testset "NC=3" begin
        β = 5.7
        NC = 3
        println("NC = $NC")
        val = 0.40073421896125794
        #@time plaq_t = heatbathtest_2D(NX,NT,β,NC)
        @time plaq_ave = heatbathtest_2D(NX,NT,β,NC)
        @test abs(plaq_ave-val)/abs(val) < eps
    end

    @testset "NC=4" begin
        β = 5.7
        NC = 4
        println("NC = $NC")
        val = 0.17476796328668975
        #@time plaq_t = heatbathtest_2D(NX,NT,β,NC)
        @time plaq_ave = heatbathtest_2D(NX,NT,β,NC)
        @test abs(plaq_ave-val)/abs(val) < eps
        #@test abs(plaq_t-val)/abs(val) < eps
    end

end


text:
using NumericalMethodsforEngineers, Test

a = [2. -3. 1.; 1. -1. -2.; 3. 1. -1.]
b = [7., -2., 0.]

(lower, upper) = lufac(a)
println("NumericalMethodsforEngineers lufac(a) - lower:")
println(lower)
println()
println("NumericalMethodsforEngineers lufac(a) - upper:")
println(upper)
println()
@show lower * upper
println()

c1 = subfor(lower, b)
c2 = subbac(upper, c1)

subfor!(lower, b)
println("Updated RHS:")
println(b)
println()

@test ≈(c1,  b, atol=1e-10)
println()

c3 = subbac(upper, c1)

subbac!(upper, b)
println("Solution vector:")
println(b)
println()

c = a \ [7., -2., 0.]
c |> display
@test ≈(c2, c, atol=1e-10)
@test ≈(c3, c, atol=1e-10)

@test ≈(a * c, [7., -2., 0.], atol=1e-10)
println()

println("Julia lu(a):")
f = lu(a);
f |> display
println()
f.U |> display


text:
#################### NaV ###############################
mutable struct Prinz_NaV{S,T} <: IonChannel
    ḡNa::S
    mNa::T
    hNa::T
    #ENa::T
end

Prinz_NaV(x) = Prinz_NaV(x, 0.,0.)

m∞(::Prinz_NaV, V) =  1.0/(1.0+exp((V+25.5)/-5.29))
h∞(::Prinz_NaV, V) =  1.0/(1.0+exp((V+48.9)/5.18))
τm(::Prinz_NaV, V) =  1.32 - 1.26/(1+exp((V+120.0)/-25.0))
τh(::Prinz_NaV, V) =  (1.34/(1.0+exp((V+62.9)/-10.0)))*(1.5+1.0/(1.0+exp((V+34.9)/3.6)))
ionic_current(::Prinz_NaV, sys::ODESystem) = sys.INa
external_params(::Prinz_NaV) = (:ENa,)

function channel_dynamics(ch::Prinz_NaV, V, Ca, D, t)
    states = @variables mNa(t) hNa(t) INa(t)
    parameters = @parameters ḡNa ENa
    eqs = [ D(mNa) ~       (1/τm(ch, V))*(m∞(ch, V) - mNa),
            D(hNa) ~       (1/τh(ch, V))*(h∞(ch, V) - hNa),
            INa ~ ḡNa*mNa^3*hNa*(ENa - V)]
    current = [eqs[3]]
    defaultmap = [mNa => ch.mNa, hNa => ch.hNa, ḡNa => ch.ḡNa]
    return eqs, states, parameters, current, defaultmap
end

#################### Slow calcium current #############################
mutable struct Prinz_CaS{S,T} <: IonChannel
    ḡCaS::S
    mCaS::T
    hCaS::T
end

Prinz_CaS(x,y) = Prinz_CaS(x,y, 0.)
Prinz_CaS(x) = Prinz_CaS(x,20., 0.)

m∞(::Prinz_CaS, V) = 1.0/(1.0+exp((V+33.0)/-8.1))
h∞(::Prinz_CaS, V) = 1.0/(1.0+exp((V+60.0)/6.2))
τm(::Prinz_CaS, V) = 2.8 + 14.0/(exp((V+27.0)/10.0) + exp((V+70.0)/-13.0))
τh(::Prinz_CaS, V) = 120.0 + 300.0/(exp((V+55.0)/9.0) + exp((V+65.0)/-16.0))

ionic_current(::Prinz_CaS, sys::ODESystem) = sys.ICaS
calcium_current(::Prinz_CaS, sys::ODESystem) = sys.ICaS
ECa(::Prinz_CaS, Ca) = (500.0)*(8.6174e-5)*(283.15)*(log(max((3000.0/Ca), 0.001)))
external_params(::Prinz_CaS) = nothing

function channel_dynamics(ch::Prinz_CaS, V, Ca, D, t)
    states = @variables mCaS(t) hCaS(t) ICaS(t)
    parameters = @parameters ḡCaS
    eqs = [ D(mCaS) ~       (1/τm(ch, V))*(m∞(ch, V) - mCaS),
            D(hCaS) ~       (1/τh(ch, V))*(h∞(ch, V) - hCaS),
            ICaS    ~       ḡCaS*mCaS^3*hCaS*(ECa(ch, Ca) - V)
            ]
    current = [eqs[3]]
    defaultmap = [mCaS => ch.mCaS, hCaS => ch.hCaS, ḡCaS => ch.ḡCaS]
    return eqs, states, parameters, current, defaultmap
end

#################### Transient calcium current ######################

mutable struct Prinz_CaT{S,T} <: IonChannel
    ḡCaT::S
    mCaT::T
    hCaT::T
end
Prinz_CaT(x) = Prinz_CaT(x, 0., 0.)
m∞(::Prinz_CaT, V) = 1.0/(1.0 + exp((V+27.1)/-7.2))
h∞(::Prinz_CaT, V) = 1.0/(1.0 + exp((V+32.1)/5.5))
τm(::Prinz_CaT, V) = 43.4 - 42.6/(1.0 + exp((V+68.1)/-20.5))
τh(::Prinz_CaT, V) = 210.0 - 179.6/(1.0 + exp((V+55.0)/-16.9))

ionic_current(::Prinz_CaT, sys::ODESystem) = sys.ICaT
calcium_current(::Prinz_CaT, sys::ODESystem) = sys.ICaT
ECa(::Prinz_CaT, Ca) = (500.0)*(8.6174e-5)*(283.15)*(log(max((3000.0/Ca), 0.001)))
external_params(::Prinz_CaT) = nothing

function channel_dynamics(ch::Prinz_CaT, V, Ca, D, t)
    states = @variables mCaT(t) hCaT(t) ICaT(t)
    parameters = @parameters ḡCaT
    eqs = [ D(mCaT) ~       (1/τm(ch, V))*(m∞(ch, V) - mCaT),
            D(hCaT) ~       (1/τh(ch, V))*(h∞(ch, V) - hCaT),
            ICaT ~ ḡCaT*mCaT^3*hCaT*(ECa(ch, Ca) - V),
            ]
    current = [eqs[3]]
    defaultmap = [mCaT => ch.mCaT, hCaT => ch.hCaT, ḡCaT => ch.ḡCaT]
    return eqs, states, parameters, current, defaultmap
end

####################  #########################
"""
A-type potassium current
"""
mutable struct Prinz_Ka{S,T} <: IonChannel
    ḡKa::S
    mKa::T
    hKa::T
    #EK::T -80
end
Prinz_Ka(x) = Prinz_Ka(x,0.,0.)

m∞(::Prinz_Ka, V) = 1.0/(1.0+exp((V+27.2)/-8.7))
h∞(::Prinz_Ka, V) = 1.0/(1.0+exp((V+56.9)/4.9))
τm(::Prinz_Ka, V) = 23.2 - 20.8/(1.0+exp((V+32.9)/-15.2))
τh(::Prinz_Ka, V) = 77.2 - 58.4/(1.0+exp((V+38.9)/-26.5))

ionic_current(::Prinz_Ka, sys::ODESystem) = sys.IKa
external_params(::Prinz_Ka) = (:EK,)

function channel_dynamics(ch::Prinz_Ka, V, Ca, D, t)
    states = @variables mKa(t) hKa(t) IKa(t)
    parameters = @parameters ḡKa EK
    eqs = [ D(mKa) ~       (1/τm(ch, V))*(m∞(ch, V) - mKa),
            D(hKa) ~       (1/τh(ch, V))*(h∞(ch, V) - hKa),
            IKa ~ ḡKa*mKa^3*hKa*(EK - V)]
    current = [eqs[3]]
    defaultmap = [mKa => ch.mKa, hKa => ch.hKa, ḡKa => ch.ḡKa]
    return eqs, states, parameters, current, defaultmap
end

################### Calcium-activated potassium current ########

mutable struct Prinz_KCa{S,T} <: IonChannel
    ḡKCa::S
    mKCa::T
    #EK::T
end

Prinz_KCa(x) = Prinz_KCa(x,0.)
m∞(::Prinz_KCa, V, Ca) = (Ca/(Ca+3.0))/(1.0+exp((V+28.3)/-12.6));
τm(::Prinz_KCa, V) = 180.6 - 150.2/(1.0+exp((V+46.0)/-22.7))
ionic_current(::Prinz_KCa, sys::ODESystem) = sys.IKCa
external_params(::Prinz_KCa) = (:EK,)

function channel_dynamics(ch::Prinz_KCa, V, Ca, D, t)
    states = @variables mKCa(t) IKCa(t)
    parameters = @parameters ḡKCa EK
    eqs = [ D(mKCa) ~       (1/τm(ch, V))*(m∞(ch, V, Ca) - mKCa),
            IKCa ~ (ḡKCa*mKCa^4)*(EK - V)]
    current = [eqs[2]]
    defaultmap = [mKCa => ch.mKCa, ḡKCa => ch.ḡKCa]
    return eqs, states, parameters, current, defaultmap
end

"""
    Delayed rectifier potassium current
"""
mutable struct Prinz_Kdr{S,T} <: IonChannel
    ḡKdr::S
    mKdr::T
    #EK::T
end
Prinz_Kdr(x) = Prinz_Kdr(x,0.)

m∞(::Prinz_Kdr, V)=  1.0/(1.0+exp((V+12.3)/-11.8));
τm(::Prinz_Kdr, V)=  14.4 - 12.8/(1.0+exp((V+28.3)/-19.2))
ionic_current(::Prinz_Kdr, sys::ODESystem) = sys.IKdr
external_params(::Prinz_Kdr) = (:EK,)

function channel_dynamics(ch::Prinz_Kdr, V, Ca, D, t)
    states = @variables mKdr(t) IKdr(t)
    parameters = @parameters ḡKdr EK
    eqs = [ D(mKdr) ~       (1/τm(ch, V))*(m∞(ch, V) - mKdr),
            IKdr ~ (ḡKdr*mKdr^4)*(EK - V)]
    current = [eqs[2]]
    defaultmap = [mKdr => ch.mKdr, ḡKdr => ch.ḡKdr]
    return eqs, states, parameters, current, defaultmap
end

"""
H current
"""
mutable struct Prinz_H{S,T} <: IonChannel
    ḡH::S
    mH::T
    # EH::T -20.
end

Prinz_H(x) = Prinz_H(x, 0.)

m∞(::Prinz_H, V) = 1.0/(1.0+exp((V+75.0)/5.5))
τm(::Prinz_H, V) = (2/( exp((V+169.7)/(-11.6)) + exp((V- 26.7)/(14.3)) ))
ionic_current(::Prinz_H, sys::ODESystem) = sys.IH
external_params(::Prinz_H) = (:EH,)

function channel_dynamics(ch::Prinz_H, V, Ca, D, t)
    states = @variables mH(t) IH(t)
    parameters = @parameters ḡH EH
    eqs = [ D(mH) ~       (1/τm(ch, V))*(m∞(ch, V) - mH),
            IH ~ ḡH*mH*(EH - V)]
    current = [eqs[2]]
    defaultmap = [mH => ch.mH, ḡH => ch.ḡH]
    return eqs, states, parameters, current, defaultmap
end

"""
leak current
"""
mutable struct Prinz_Leak{S} <: IonChannel
    ḡLeak::S
    #ELeak::S -50
end

ionic_current(::Prinz_Leak, sys::ODESystem) = sys.ILeak
external_params(::Prinz_Leak) = (:ELeak,)

function channel_dynamics(ch::Prinz_Leak, V, Ca, D, t)
    states = @variables ILeak(t)
    parameters = @parameters ḡLeak ELeak
    eqs = [ILeak ~ ḡLeak*(ELeak - V)]
    current = [eqs[1]]
    defaultmap = [ḡLeak => ch.ḡLeak]
    return eqs, states, parameters, current, defaultmap
end


text:
using Base.Test

import LightTable: round3

pi

float(pi)

100pi
pi/10000

@test round3(1pi) == "3.142"
@test round3(100pi) == "314.159"
@test round3(pi/10000) == "0.000314"
@test round3(1/0.999) == "1.001"
@test round3(1/9.95) == "0.101"
@test round3(1/0.999) == "1.001"

# The rounder currently bails out on these
# round3(1.9999)
# round3(0.09999)


text:
# CUDA-specific operations on pointers with address spaces

## adrspace aliases

export AS

module AS

const Generic  = 0
const Global   = 1
const Shared   = 3
const Constant = 4
const Local    = 5

end

## ldg

# operand types supported by llvm.nvvm.ldg.global
# NOTE: CUDA 8.0 supports more caching modifiers, but those aren't supported by LLVM yet
const LDGTypes = Union{UInt8, UInt16, UInt32, UInt64,
                       Int8, Int16, Int32, Int64,
                       Float32, Float64}

# TODO: this functionality should throw <sm_32
@generated function pointerref_ldg(p::LLVMPtr{T,AS.Global}, i::Int,
                                   ::Val{align}) where {T<:LDGTypes,align}
    sizeof(T) == 0 && return T.instance
    JuliaContext() do ctx
        eltyp = convert(LLVMType, T, ctx)

        # TODO: ccall the intrinsic directly

        T_int = convert(LLVMType, Int, ctx)
        T_int32 = LLVM.Int32Type(ctx)
        T_ptr = convert(LLVMType, p, ctx)
        T_typed_ptr = LLVM.PointerType(eltyp, AS.Global)

        # create a function
        param_types = [T_ptr, T_int]
        llvm_f, _ = create_function(eltyp, param_types)

        # create the intrinsic
        intrinsic_name = let
            class = if isa(eltyp, LLVM.IntegerType)
                :i
            elseif isa(eltyp, LLVM.FloatingPointType)
                :f
            else
                error("Cannot handle $eltyp argument to unsafe_cached_load")
            end
            width = sizeof(T)*8
            typ = Symbol(class, width)
            "llvm.nvvm.ldg.global.$class.$typ.p1$typ"
        end
        mod = LLVM.parent(llvm_f)
        intrinsic_typ = LLVM.FunctionType(eltyp, [T_typed_ptr, T_int32])
        intrinsic = LLVM.Function(mod, intrinsic_name, intrinsic_typ)

        # generate IR
        Builder(ctx) do builder
            entry = BasicBlock(llvm_f, "entry", ctx)
            position!(builder, entry)

            typed_ptr = bitcast!(builder, parameters(llvm_f)[1], T_typed_ptr)
            typed_ptr = gep!(builder, typed_ptr, [parameters(llvm_f)[2]])
            ld = call!(builder, intrinsic,
                    [typed_ptr, ConstantInt(Int32(align), ctx)])

            metadata(ld)[LLVM.MD_tbaa] = tbaa_addrspace(AS.Global, ctx)

            ret!(builder, ld)
        end

        call_function(llvm_f, T, Tuple{LLVMPtr{T,AS.Global}, Int}, :((p, Int(i-one(i)))))
    end
end

# interface

export unsafe_cached_load

unsafe_cached_load(p::LLVMPtr{<:LDGTypes,1}, i::Integer=1, align::Val=Val(1)) =
    pointerref_ldg(p, Int(i), align)
# NOTE: fall back to normal unsafe_load for unsupported types. we could be smarter here,
#       e.g. destruct/load/reconstruct, but that's too complicated for what it's worth.
unsafe_cached_load(p::LLVMPtr, i::Integer=1, align::Val=Val(1)) =
    unsafe_load(p, i, align)


text:
# Dependency.
using BenchmarkTools

# Define a function to benchmark.
foo(x) = x^2

# Define a benchmark group and fill it with some tests.
regressions = BenchmarkGroup()
srand = 123 # Seed as in the tutorial.
regressions["floats"] = @benchmarkable foo.($(rand(100))) # Interpolate as in the tutorial.
regressions["ints"] = @benchmarkable foo.($(rand(1:1000, 100)))

# Get and compress results.
results = run(regressions)
results = median(results)

# Write results to file.
BenchmarkTools.save("benchmarks.json", results)

# Load results from file.
old = BenchmarkTools.load("benchmarks.json")[1]

# Make judgment.
@show judgment = judge(old, results)

# Inspect judgment element.
@show judgment["floats"]


text:
include("poly_factorization_project.jl")
println("For this repository, a new type was made for polynomials.")
println("The following are tests for the different functions that can be done with this new type")

x = x_poly()
@show p1 = 2x^3 + 4x^2 - 3x
@show p2 = 2x^4 - 4x^2 - 3x + 3

println("The first test is for polynomial mutliplication")
@show p1*p2
println("The following is a test for polynomial power")
@show p1^3

println("In the following 2 lines, it can be seen that the derivative function works on the polynomial class")
println("The first line is the derivative function on the polynomials multiplied togeher")
println("The second line is using the product rule")
d1 = derivative(p1*p2)
d2 = derivative(p1)*p2 + p1*derivative(p2);
println("Derivative of multiplied function: ", d1)
println("Derivative with product rule: ", d2)

prime = 17
p = mod((7x^3 + 2x^2 + 8x + 1)*(x^2+x+1),prime)
println("Will factor this polynomial (mod $prime): ", p)
factorization = factor(p,prime)
println("Here is the factorization: ", factorization)

pr = mod(expand_factorization(factorization),prime)
println("Reconstructing: ", pr)

println("The following is a test for polynomial addition")
x = x_poly()
@show p1 = 2x^2+3x +(-5) #Note we need +(-5)... try later without... how to fix?
@show p2 = -3x^2 - 4x +6
@show p1+p2

extended_euclid_alg(p1*p2,p2,101)


text:
using KitBase

@code_warntype PSpace2D()
@code_warntype VSpace2D()


text:
export group_by

macro group_by(data, ex)
    @capture(ex, args__)
    esc(
        :(
            group_by($data)
        )
    )
end

function group_by(args...; kwargs...)
    @info "group_by: The equivalent function in DataFrames.jl is `groupby`"
    groupby(args...; kwargs...)
end


text:
module RunTests
using Test
export runtests
function runtests()

    include("test_Tools.jl")
    include("test_DistSystem.jl")
    include("test_SimplePF.jl")
    include("test_Estimator.jl")
    include("test_HCEstimator.jl")

    @testset "example" begin
        include("../example/example.jl")
        @test true
    end

end
end

RunTests.runtests()
print("")


text:
function mountaincar_cmax_main()
    cmax_steps = []
    model = MountainCar(0.0)
    cmax_planner = MountainCarCMAXPlanner(model, 1000, true_params)
    generateHeuristic!(cmax_planner)
    for rock_c in range_of_values
        println("Rock_c is ", rock_c)
        mountaincar = MountainCar(rock_c)
        cmax_agent = MountainCarCMAXAgent(mountaincar, cmax_planner)
        n_steps = run(cmax_agent, max_steps = 1e4)
        push!(cmax_steps, n_steps)
        clearResiduals!(cmax_planner)
    end
    cmax_steps
end


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

using Markdown
using InteractiveUtils

# ╔═╡ b0723702-fe6d-11ea-3809-35a9826b419b
using Pkg, DrWatson

# ╔═╡ b0727334-fe6d-11ea-0e69-f563a75647df
begin
	#@quickactivate "StatisticalRethinkingStan"
	using StanSample
	using StatisticalRethinking
end

# ╔═╡ e626c09e-fe6c-11ea-1ca0-63c764e7a6fd
md"## Clip-06-07s.jl"

# ╔═╡ b072ebc0-fe6d-11ea-337b-a57fbb329ea5
# Snippet 6.1

begin
	N = 100
	df = DataFrame(
	  height = rand(Normal(10, 2), N),
	  leg_prop = rand(Uniform(0.4, 0.5), N),
	)
	df.leg_left = df.leg_prop .* df.height + rand(Normal(0, 0.02), N)
	df.leg_right = df.leg_prop .* df.height + rand(Normal(0, 0.02), N)
end;

# ╔═╡ d18a9308-fe6d-11ea-3c08-3f7840232a0a
Text(precis(df; io=String))

# ╔═╡ b0813202-fe6d-11ea-2b1a-a37c2943d8dc
md"### Snippet 6.2"

# ╔═╡ b081d1a8-fe6d-11ea-0e63-cbc35b0bbdfb
m6_2 = "
data {
  int <lower=1> N;
  vector[N] H;
  vector[N] LL;
}

parameters {
  real a;
  real bL;
  real <lower=0> sigma;
}

model {
  vector[N] mu;
  mu = a + bL * LL;
  a ~ normal(10, 100);
  bL ~ normal(2, 10);
  sigma ~ exponential(1);
  H ~ normal(mu, sigma);
}
";

# ╔═╡ b08c2a02-fe6d-11ea-05b2-5973a5b63804
begin
	m6_2s = SampleModel("m6.2s", m6_2)
	m6_2_data = Dict(:H => df.height, :LL => df.leg_left, :N => size(df, 1))
	rc6_2s = stan_sample(m6_2s, data=m6_2_data)
	success(rc6_2s) && (part6_2s = read_samples(m6_2s, :particles))
end

# ╔═╡ b09dac48-fe6d-11ea-3927-f5a982c7b715
if success(rc6_2s)
	chns6_2s = read_samples(m6_2s, :mcmcchains)
	CHNS(chns6_2s)
end

# ╔═╡ b0a8ce2a-fe6d-11ea-19e5-d7f5a9022680
success(rc6_2s) && plot(chns6_2s; seriestype=:traceplot)

# ╔═╡ bf05bf54-044f-11eb-1570-3b64f374127a
success(rc6_2s) && plot(chns6_2s; seriestype=:density)

# ╔═╡ b0aa516e-fe6d-11ea-33a6-e9a126093f43
md"## End of clip-06-07s.jl"

# ╔═╡ Cell order:
# ╟─e626c09e-fe6c-11ea-1ca0-63c764e7a6fd
# ╠═b0723702-fe6d-11ea-3809-35a9826b419b
# ╠═b0727334-fe6d-11ea-0e69-f563a75647df
# ╠═b072ebc0-fe6d-11ea-337b-a57fbb329ea5
# ╠═d18a9308-fe6d-11ea-3c08-3f7840232a0a
# ╟─b0813202-fe6d-11ea-2b1a-a37c2943d8dc
# ╠═b081d1a8-fe6d-11ea-0e63-cbc35b0bbdfb
# ╠═b08c2a02-fe6d-11ea-05b2-5973a5b63804
# ╠═b09dac48-fe6d-11ea-3927-f5a982c7b715
# ╠═b0a8ce2a-fe6d-11ea-19e5-d7f5a9022680
# ╠═bf05bf54-044f-11eb-1570-3b64f374127a
# ╟─b0aa516e-fe6d-11ea-33a6-e9a126093f43


text:
using DifferentialEquations
using Plots
using Roots
barrier=4
function condition(u,t,integrator) # Event when event_f(u,t) == 0
 u[3]
end

function affect!(integrator)
end
#cb1 = SavingCallback((u,t,integrator)->(u[4],u[2]), saved_values)
cb= ContinuousCallback(condition, affect!)

    Energy=.1;
    H=(2*Energy)^2
    h=sqrt(H)
    n_iter=3
    # Y_start=rand(1)*1e-6
    # Y_end=Y_start
    t_end=2.1e3
    # Y=range(Y_start,stop=Y_end,length=n_iter)
    # X=[rand(1)*1e-6]
    # P=zeros(n_iter)
    Final_H=zeros(n_iter,n_iter)
    H1(p,q,t)= ( (q[2]-q[1])^2+(p[2]-p[1])^2 )*( (q[2]+q[1])^2+(p[2]+p[1])^2 )
    H2(p,q,t)= (p[2]^4+2*p[2]^2*(q[1]^2-1)+(1+q[1]^2)^2 )*(q[2]^4+2*q[2]^2*(p[1]^2+1)+(p[1]^2-1)^2 )
    HAM(p,q,t)=H1(p,q,t)/H2(p,q,t)
    PofY_find(YY)=sqrt((h+YY.^2-h*YY.^2 )./(1+h-h*YY.^2))
    p=plot()
    DPI_VAL=1000

    Hsolve(XX,YY,PP)=((XX.^2+(PP-YY).^2)*(XX.^2+(PP+YY).^2)) ./((PP.^4+2*PP.^2*(XX.^2-1)+(1+XX.^2).^2)*(YY.^2-1).^2)
    A_start=1e-5
    A_end=.1
    t_end=2.1e3
    global
    B=range(A_start,stop=A_end,length=n_iter)

    for j=1:n_iter
    x=B[j]
      for k=1:n_iter
        y=B[k]
        q0,p0=[zeros(2) for i in 1:2]
        F(PP)=Hsolve(x,y,PP)-H
        pp=find_zero(F,.01)
        q0[1]=x
        p0[1]=pp
        q0[2]=0
        p0[2]=y
        Final_H[j,k]=HAM(p0,q0,0)
        prob= HamiltonianProblem{true}(HAM, p0, q0, (0., t_end))
        t,A=solve(prob, Vern9(),reltol=1e-8,abstol=1e-20,maxiters=1e15, callback=cb,save_end=false,save_everystep=false)
        plot!(A[2,:],A[4,:],seriestype=:scatter,marker=:circle,msw = 0,ms=1,seriescolor = :black,legend=false,dpi=DPI_VAL)
        plot!(A[2,:],-A[4,:],seriestype=:scatter,marker=:circle,msw = 0,ms=1,seriescolor = :black,legend=false)
        plot!(-A[2,:],A[4,:],dpi=DPI_VAL,seriestype=:scatter,marker=:circle,msw = 0,ms=1,seriescolor = :black,legend=false)
        plot!(-A[2,:],-A[4,:],dpi=DPI_VAL,seriestype=:scatter,marker=:circle,msw = 0,ms=1,seriescolor = :black,legend=false)
        plot!([A[2,end]],[A[4,end]],seriestype=:scatter,marker=:circle,msw = 0,ms= 1,seriescolor = :red,legend=false,dpi=DPI_VAL)
        plot!([≈],[-A[4,end]],dpi=DPI_VAL,seriestype=:scatter,marker=:circle,msw = 0,ms=1,seriescolor = :red,legend=false)
        plot!([A[2,end]],[A[4,end]],dpi=DPI_VAL,seriestype=:scatter,marker=:circle,msw = 0,ms=1,seriescolor = :red,legend=false)
        plot!([-A[2,end]],[-A[4,end]],dpi=DPI_VAL,seriestype=:scatter,marker=:circle,msw = 0,ms=1,seriescolor =:red,legend=false)
    # XX=A[2,end]
    # YY=A[4,end]
end
end
 # img = SVG("pss.svg", 14cm, 8cm)
 # draw(img, p)


text:
module Example

# Write your package code here.
using ForwardDiff

include("extra_file.jl")

export my_f
export derivative_of_my_f

end


text:
"""
    AbstractModelGenerator

Abstract type of problem generators. They are used to fill a CPModel which will be solved to train
a LearnedHeuristic. Some Generators are provided directly with the package but a user can also create
his own problem generator.
"""
abstract type AbstractModelGenerator end


text:
module DenseM
    include("../tools/glorotuniform.jl")
    using .GlorotUniform
    using LoopVectorization

    mutable struct Dense
        save_layer::Any
        load_layer::Any
        activate::Any
        initialize::Any
        update::Any

        input_shape::Tuple
        output_shape::Tuple
        layer_size::Int64
        activation_function::Module

        weights::Array{Float32}
        biases::Array{Float32}

        Vdw::Array{Float32}
        Sdw::Array{Float32}
        Vdb::Array{Float32}
        Sdb::Array{Float32}

        value::Array{Float32}
        output::Array{Float32}
        ∇biases::Array{Float32}
        δ::Array{Float32}

        function Dense(;input_shape::Tuple, layer_size::Int64, activation_function::Module, randomization::Bool=true)
            if randomization
                weights = Array{Float32}(GUMatrix(input_shape[1], layer_size, (layer_size, input_shape[1])))
                biases = Array{Float32}(GUMatrix(input_shape[1], layer_size, layer_size))
            else
                weights = zeros(Float32, (layer_size, input_shape[1]))
                biases = zeros(Float32, layer_size)
            end

            output_shape = (layer_size,)

            Vdw = zeros(Float32, size(weights))
            Sdw = zeros(Float32, size(weights))
            Vdb = zeros(Float32, size(biases))
            Sdb = zeros(Float32, size(biases))
            new(save_Dense, load_Dense, activate_Dense, init_Dense, update_Dense, input_shape, output_shape, layer_size, activation_function, weights, biases, Vdw, Sdw, Vdb, Sdb)
        end
    end

    function init_Dense(layer::Dense, mini_batch::Int64)
        layer.value = zeros(Float32, (layer.layer_size, mini_batch))
        layer.output = zeros(Float32, (layer.layer_size, mini_batch))
        layer.∇biases = zeros(Float32, (layer.layer_size, mini_batch))
        layer.δ = zeros(Float32, (layer.input_shape[1], mini_batch))
    end

    function activate_Dense(layer::Dense, input::Array{Float32})
        @avxt for x in axes(layer.weights, 1), y in axes(input, 2)
            c = 0.0f0
            for z in axes(layer.weights, 2)
                c += layer.weights[x,z]*input[z,y]
            end
            layer.value[x,y] = c+layer.biases[x]
        end
        # layer.value = layer.weights*input .+ layer.biases
        layer.activation_function.func!(layer.output, layer.value)
    end

    function update_Dense(layer::Dense, optimizer::String, input::Array{Float32}, δₗ₊₁::Array{Float32}, α::Float64, parameters::Tuple, direction::Int64=1)
        layer.activation_function.get_∇biases!(layer.∇biases, layer.value, δₗ₊₁)
        @avxt for x in axes(layer.weights, 2), y in axes(layer.∇biases, 2)
            c = 0.0f0
            for z in axes(layer.weights, 1)
                c += layer.weights[z,x]*layer.∇biases[z,y]
            end
            layer.δ[x,y] = c
        end
        # layer.δ = transpose(layer.weights)*∇biases

        if optimizer=="SGD"
            @avxt for x in axes(layer.∇biases, 1), y in axes(input, 1)
                layer.weights[x,y] -= α*layer.∇biases[x,1]*input[y,1]*direction
            end
            # layer.weights -= ∇biases*transpose(input).*α
            @avxt for i in 1:length(layer.biases)
                layer.biases[i] -= α*layer.∇biases[i,1]*direction
            end
            # layer.biases -= sum(∇biases, dims=2).*α

        elseif optimizer=="Minibatch_GD"
            batch_size = size(input, 2)
            @avxt for x in axes(layer.∇biases, 1), y in axes(input, 1)
                c = 0.0f0
                for b in axes(layer.∇biases, 2)
                    c += layer.∇biases[x,b]*input[y,b]
                end
                layer.weights[x,y] -= c*α*direction/batch_size
            end
            # layer.weights -= ∇biases*transpose(input).*α
            @avxt for i in 1:length(layer.biases)
                c = 0.0f0
                for b in axes(layer.∇biases, 2)
                    c += layer.∇biases[i,b]
                end
                layer.biases[i] -= c*α*direction/batch_size
            end
            # layer.biases -= sum(∇biases, dims=2).*α

        elseif optimizer=="Adam"
            t = parameters[1]
            β₁ = parameters[2]
            β₂ = parameters[3]
            ϵ = parameters[4]

            @avxt for x in axes(layer.∇biases, 1), y in axes(input, 1)
                layer.Vdw[x,y] = layer.Vdw[x,y]*β₁ + layer.∇biases[x,1]*input[y,1]*(1-β₁)
                layer.Sdw[x,y] = layer.Sdw[x,y]*β₂ + (layer.∇biases[x,1]*input[y,1])^2*(1-β₂)
                layer.weights[x,y] -= α*(layer.Vdw[x,y]/(1-β₁^t))/(sqrt(layer.Sdw[x,y]/(1-β₂^t))+ϵ)*direction
            end

            @avxt for i in 1:length(layer.biases)
                layer.Vdb[i] = layer.Vdb[i]*β₁ + layer.∇biases[i,1]*(1-β₁)
                layer.Sdb[i] = layer.Sdb[i]*β₂ + layer.∇biases[i,1]^2*(1-β₂)
                layer.biases[i] -= α*(layer.Vdb[i]/(1-β₁^t))/(sqrt(layer.Sdb[i]/(1-β₂^t))+ϵ)*direction
            end

        elseif optimizer=="AdaBelief"
            t = parameters[1]
            β₁ = parameters[2]
            β₂ = parameters[3]
            ϵ = parameters[4]

            @avxt for x in axes(layer.∇biases, 1), y in axes(input, 1)
                layer.Vdw[x,y] = layer.Vdw[x,y]*β₁ + layer.∇biases[x,1]*input[y,1]*(1-β₁)
                layer.Sdw[x,y] = layer.Sdw[x,y]*β₂ + (layer.∇biases[x,1]*input[y,1]-layer.Vdw[x,y])^2*(1-β₂) + ϵ
                layer.weights[x,y] -= α*(layer.Vdw[x,y]/(1-β₁^t))/(sqrt(layer.Sdw[x,y]/(1-β₂^t))+ϵ)*direction
            end

            @avxt for i in 1:length(layer.biases)
                layer.Vdb[i] = layer.Vdb[i]*β₁ + layer.∇biases[i,1]*(1-β₁)
                layer.Sdb[i] = layer.Sdb[i]*β₂ + (layer.∇biases[i,1]-layer.Vdb[i])^2*(1-β₂)
                layer.biases[i] -= α*(layer.Vdb[i]/(1-β₁^t))/(sqrt(layer.Sdb[i]/(1-β₂^t))+ϵ)*direction
            end
        end
    end

    function save_Dense(layer::Dense, file::Any, id::String)
        write(file, id, "Dense")
        write(file, id*"input_shape", collect(layer.input_shape))
        write(file, id*"layer_size", layer.layer_size)
        write(file, id*"activation_function", layer.activation_function.get_name())

        write(file, id*"weights", layer.weights)
        write(file, id*"biases", layer.biases)
    end

    function load_Dense(layer::Dense, file::Any, id::String)
        layer.weights = read(file, id*"weights")
        layer.biases = read(file, id*"biases")
    end

    function get_args(file::Any, id::String)
        input_shape = tuple(read(file, id*"input_shape")...)
        layer_size = read(file, id*"layer_size")
        return (input_shape=input_shape, layer_size=layer_size)
    end
end


text:
make_t1() = Template(data=attr(
    scatter=[attr(marker_color="red"), attr(marker_color="blue")]
))

make_t2() = Template(data=attr(
    scatter=[attr(marker_line_color="black")],
    bar=[attr(marker_color="cyan")]
))
make_t3() = merge(make_t1(), make_t2())

@testset "Templates" begin

    @testset "constructor with data::PlotlyAttribute" begin
        t1 = make_t1()
        @test t1 isa Template
        @test isempty(t1.layout)
        @test length(t1.data[:scatter]) == 2
    end

    @testset "merge templates" begin
        t1 = make_t1()
        t2 = make_t2()
        t3 = merge(t1, t2)

        # test that neither t1 nor t2 were modified
        @test length(t1.data[:scatter]) == 2
        @test length(t2.data[:scatter]) == 1
        @test !haskey(t1.data, :bar)

        # test that updates were applied
        @test length(t3.data[:scatter]) == 2
        @test t3.data[:scatter][1] == attr(marker=attr(color="red", line_color="black"))
        @test t3.data[:scatter][2] == attr(marker_color="blue", marker_line_color="black")
        @test t3.data[:bar][1] == attr(marker_color="cyan")
    end

    @testset "Can save Template to `templates` object" begin
        t1 = make_t1()
        templates.foobar = t1
        @test haskey(templates.templates, :foobar)

        # clean up
        pop!(templates.templates, :foobar)
    end

    @testset "Multiple templates via `+` syntax" begin
        t1 = make_t1()
        t2 = make_t2()
        templates.t1 = t1
        templates.t2 = t2

        @test templates["t1+t2"] == merge(t1, t2)
        # clean up
        pop!(templates.templates, :t1)
        pop!(templates.templates, :t2)
    end

    @testset "loading known templates from file" begin
        @test templates["ggplot2"] isa Template
        @test templates["plotly_dark"] isa Template

    end

end


text:
function spindown_convergence_plot(dtlist, l2err, maxerr, sdc_order)
    orderstring = latexstring("\\mathcal{O}(\\Delta t^$sdc_order)")
    h = loglog(dtlist, maxerr, ".-", label=LaTeXString("\$||\\tilde u^{($sdc_order)} - u||_\\infty\$"))
    color = h[1][:get_color]()
    loglog(dtlist, l2err, "*-", label=LaTeXString("\$||\\tilde u^{($sdc_order)} - u||_2\$"), color=color)
    matchidx = 4
    C = ( maxerr[matchidx] + l2err[matchidx] )/2
    @show dtlist_ext = vcat(dtlist[1]*2, dtlist, dtlist[end]/2)
    loglog(dtlist_ext, C*(dtlist_ext/dtlist[matchidx]).^sdc_order, "--", label=orderstring, color=color)
    #axis("equal")
    grid("on")
    #xlim(left=3e-5, right=1e-2)
    axis("tight")
    ylim([1e-13, 1e-2])
    xlabel(L"Time step $\Delta t$")
    ylabel("Relative error")
    legend(bbox_to_anchor=(1.0,1.03), ncol=2)
    tight_layout(0.1)
end

close("all")
figure(figsize=(9,3))
spindown_convergence_plot(dtlist1, l2err1, maxerr1, 1)
spindown_convergence_plot(dtlist2, l2err2, maxerr2, 2)
spindown_convergence_plot(dtlist3, l2err3, maxerr3, 3)
spindown_convergence_plot(dtlist4, l2err4, maxerr4, 4)

savefig("../docs/paper/figs/convergence_time.pdf")


text:
using MPI, Test

include("../testhelpers.jl")

@testset "Driver" begin
    runmpi(joinpath(@__DIR__, "cr_unit_tests.jl"), ntasks = 1)
    runmpi(joinpath(@__DIR__, "driver_test.jl"), ntasks = 1)

    #TODO: these take too long on Azure CI -- need to go in weekly
    # runmpi(`$(joinpath(@__DIR__, mms3.jl)) --checkpoint-walltime 60 --checkpoint-keep-all`; ntasks=1)
    # runmpi(`$(joinpath(@__DIR__, mms3.jl)) --restart-from-num 3`; ntasks=1)
end


text:
# SEIR Model

using DifferentialEquations

beta = 0.215; gamma = 0.07; N = 7000;
alpha = 0.25;

"""SEIR Model
Extension of the basic SIR model

"""
function SEIR!(du,u,p,t)
    du[1] = -beta* u[1]/N*u[3]
    du[2] = beta*u[1]/N*u[3] - alpha*u[2]
    du[3] = alpha* u[2] - gamma*u[3]
    du[4] = gamma*u[3]
end

   u0 = [N-15.0;15.0;0.0;0.0]
   tspan = (0.0,100.0)
   prob = ODEProblem(SEIR!,u0,tspan)
   sol = solve(prob)

   plot(sol)


text:
using Test
import Omega.Lang: AppExpr, LetExpr
import Omega.Lang: unparse

function test1()
  op = """
  μ : uniform 0 1
  x : normal μ 1
  rand x + μ
  """
  oe = LetExpr(:μ, AppExpr(:uniform, [0, 1]),
              LetExpr(:x, AppExpr(:uniform, [0, 1]),
              AppExpr(:rand, )
              ))

  @test unparse(oe) == op
end

function test2()
  expr = quote
  function f(w)
    x = sqrt(w) + 1
    y = sqrt(x) + 3
    if

  end
end


text:
include("tsne_imports.jl")

get_simdat_path(path_to_output) = joinpath(splitdir(path_to_output)[1], "simdat.bson")
results_dir = "/scratchdata/alui2/cytof/results/repfam/test-sim-6-8-4"

for zind in 1:3
  for pmiss in (0.0, 0.6)
    # Get simulation name
    simname = "pmiss$(pmiss)-phi0-zind$(zind)"

    println("Processing: ", simname)

    # Get path to simuldate data.
    path = joinpath(results_dir, simname, "simdat.bson")

    # Fit TSNE jointly on all samples.
    tsne, sample_ind, Y, M, true_labels = let
      use_complete_data = (pmiss == 0.0)
      compute_combined_tsne(path, use_complete_data=use_complete_data,
                            seed=1, verbose=2)
    end

    # Create data frame
    df = let
      num_markers = size(Y, 2)
      Ynames = [Symbol("Y$(j)") for j in 1:num_markers]
      Mnames = [Symbol("M$(j)") for j in 1:num_markers]
      colnames = [:tsne1, :tsne2, :sample_ind, :true_labels,
                  Ynames..., Mnames...]
      DataFrame([tsne sample_ind true_labels Y M], colnames)
    end

    # Write output to file.
    csvdir = "viz/csv"
    mkpath(csvdir)
    CSV.write("$(csvdir)/tsne-$(simname).csv", df)
  end
end


text:
using Test, Sockets, LanguageServer, CSTParser, SymbolServer, SymbolServer.Pkg, StaticLint, JSON
using LanguageServer: Document, get_text, get_offset, get_line_offsets, get_position_at, get_open_in_editor, set_open_in_editor, is_workspace_file, applytextdocumentchanges
const LS = LanguageServer
const Range = LanguageServer.Range

@testset "LanguageServer" begin

@testset "document" begin
include("test_document.jl")
end
@testset "communication" begin
include("test_communication.jl")
end
@testset "hover" begin
include("test_hover.jl")
end
@testset "edit" begin
include("test_edit.jl")
end
@testset "actions" begin
include("test_actions.jl")
end
@testset "paths" begin
include("test_paths.jl")
end

end


text:
function acceleration(time::Float64, pos::Vector{Float64})
    moon_pos = moon_position(time)

    distance_from_earth = pos
    distance_to_moon = pos - moon_pos
    mag_e = norm(distance_from_earth)
    mag_m = norm(distance_to_moon)
    return -G * (ME * distance_from_earth / mag_e^3 + MM * distance_to_moon / mag_m^3)
end

function update(me::Command_Module, time::Float64, h::Float64)
    acceleration0 = acceleration(time, me.position)
    velocityE = me.velocity + h * acceleration0
    positionE = me.position + h * me.velocity
    velocityH = me.velocity + h * 0.5 * (acceleration0 + acceleration(time + h, positionE))
    positionH = me.position + h * 0.5 * (me.velocity + velocityH)

    me.velocity = velocityH
    me.position = positionH

    me.positionH = positionH
    me.velocityH = velocityH
    me.positionE = positionE
    me.velocityE= velocityE

    me
end

max_stars_repo_name: rodrigolece/grupo-fisica-compu

max_issues_repo_name: rodrigolece/grupo-fisica-compu

max_forks_repo_name: rodrigolece/grupo-fisica-compu


text:
abstract Animal
speak(a::Animal) = "$(a.name) says $(sound(a))"

type Horse <: Animal
  name
end
sound(h::Horse) = "neigh"

type Cow <: Animal
  name
end
sound(c::Cow) = "moooo"

type Sheep <: Animal
  name
end
sound(s::Sheep) = "baaaa"

s = Horse("CJ")
@assert speak(s) == "CJ says neigh"
c = Cow("Bessie")
@assert speak(c) == "Bessie says moooo"
@assert speak(Sheep("Little Lamb")) == "Little Lamb says baaaa"


text:
"""
    PlanResult{S <: MAPFState, A <: MAPFAction, C <: Number}

Stores the result of an individual agent's search.

Attributes:
    - `states::Vector{Tuple{S,C}}` The list of (state, cost-to-go) on the path
    - `actions::Vector{Tuple{A,C}}` The list of (actions, cost) on the path.
    - `cost::C` The cumulative cost of the path
    - `fmin::C` The minimum f-value expanded during search
"""
@with_kw struct PlanResult{S <: MAPFState, A <: MAPFAction, C <: Number}
    states::Vector{Tuple{S,C}}  = Vector{Tuple{S,C}}(undef, 0)
    actions::Vector{Tuple{A,C}} = Vector{Tuple{A,C}}(undef, 0)
    cost::C                     = zero(C)
    fmin::C                     = zero(C)
end

Base.isempty(pln::PlanResult) = isempty(pln.states)

"""
    get_mapf_state_from_idx(env::E, idx::Int64) where {E <: MAPFEnvironment}

My A* search implementation returns a list of graph indices; this function maps an index
to the corresponding state for that environment.
"""
function get_mapf_state_from_idx end

"""
    get_mapf_action(env::E, source::Int64, target::Int64) where {E <: MAPFEnvironment}

Returns the MAPFAction while going from source index to target index.
"""
function get_mapf_action end

"""
    get_plan_result_from_astar(env::E, a_star_dists::Dict, a_star_parent_indices::Dict,
                               start_idx::Int64, goal_idx::Int64,
                               best_fvalue::D) where {E <: MAPFEnvironment, D <: Number}

Converts the result of an A*-search (from my Graphs.jl implementation) to a corresponding PlanResult instance.
"""
function get_plan_result_from_astar(env::E, a_star_dists::Dict, a_star_parent_indices::Dict,
                                    start_idx::Int64, goal_idx::Int64,
                                    best_fvalue::D) where {E <: MAPFEnvironment, D <: Number}

    # First ensure the goal is reachable
    if ~(haskey(a_star_dists, goal_idx))
        return nothing
    end

    # Set both cost and fmin
    cost = a_star_dists[goal_idx]
    if best_fvalue == zero(D)
        fmin = a_star_dists[goal_idx]
    else
        fmin = best_fvalue
    end

    # Insert last elements to states and actions arrays
    goal_state = get_mapf_state_from_idx(env, goal_idx)
    states = [(goal_state, cost)]

    curr_idx = a_star_parent_indices[goal_idx]

    goal_action = get_mapf_action(env, curr_idx, goal_idx)
    action_cost = a_star_dists[goal_idx] - a_star_dists[curr_idx]
    actions = [(goal_action, action_cost)]

    # Now walk back to start
    while curr_idx != start_idx
        prev_idx = a_star_parent_indices[curr_idx]

        # Update states
        pushfirst!(states, (get_mapf_state_from_idx(env, curr_idx), a_star_dists[curr_idx]))

        # Update actions
        pushfirst!(actions, (get_mapf_action(env, prev_idx, curr_idx),
                        a_star_dists[curr_idx] - a_star_dists[prev_idx]))

        curr_idx = prev_idx
    end

    pushfirst!(states, (get_mapf_state_from_idx(env, start_idx), 0))

    return PlanResult(states=states, actions=actions, cost=cost, fmin=fmin)
end


text:
using LinearAlgebra
function test1(Nwing)
    NX = 4
    NY = 4
    NZ = 4
    NT = 4
    #Nwing = 1
    Dim = 4
    NC = 3

    u1 = IdentityGauges(NC,Nwing,NX,NY,NZ,NT)
    U = Array{typeof(u1),1}(undef,Dim)
    U[1] = u1
    for μ=2:Dim
        U[μ] = IdentityGauges(NC,Nwing,NX,NY,NZ,NT)
    end

    gauge_action = GaugeAction(U)
    plaqloop = make_loops_fromname("plaquette")
    append!(plaqloop,plaqloop')
    β = 1
    push!(gauge_action,β,plaqloop)

    show(gauge_action)

    Uout = evaluate_GaugeAction_untraced(gauge_action,U)
    println(tr(Uout))

end

@testset "initialization" begin
    println("with Nwing")
    Nwing = 1
    test1(Nwing)
    println("without Nwing")
    Nwing = 0
    test1(Nwing)
end


text:
"""
Computes random samples of functions from aligned data using Gaussian model

    gauss_model(warp_data; n=1, sort_samples=false)
    :param warp_data: fdawarp type from srsf_align of aligned data
    :param n: number of samples
    :param sort_samples: sort samples

    Returns warp_data containing
    :return fs: random aligned functions
    :return gams: random warping functions
    :return ft: random functions
"""
function gauss_model(warp_data::fdawarp; n=1, sort_samples=false)
    fn = warp_data.fn
    timet = warp_data.time
    qn = warp_data.qn
    gam = warp_data.gam
    # Parameters
    eps1 = eps(Float64);
    binsize = mean(diff(timet));
    M = length(timet);

    # compute mean and covariance in q-domain
    mq_new = mean(qn,dims=2);
    mididx = round(Integer, M/2);
    m_new = sign.(fn[mididx,:]) .* sqrt.(abs.(fn[mididx,:]));
    mqn = vec([mq_new; mean(m_new)]);
    qn2 = vcat(qn, m_new');
    C = Statistics.covm(qn2, mean(qn2,dims=2), 2);

    tmp = cholesky(C,check=false)
    if (!issuccess(tmp))
        for i = 1:M+1
            C[i,i] = C[i,i] + 1e-7
        end
    end
    di = MvNormal(mqn, C);
    q_s = rand(di, n);

    # compute the correspondence to the orignal function domain
    fs = zeros(M, n);
    for k in 1:n
        fs[:, k] = cumtrapzmid(timet, q_s[1:M,k].*abs.(q_s[1:M,k]),
                               sign(q_s[M+1,k])*(q_s[M+1,k]^2), mididx);
    end
    fbar = mean(fn,dims=2)
    fsbar = mean(fs,dims=2)
    err = repeat(fbar-fsbar,1,n)
    fs = fs + err

    # random warping generation
    rgam = random_gamma(gam, n);
    gams = zeros(M, n);
    for k in 1:n
        gams[:,k] = invert_gamma(rgam[:, k]);
    end

    if sort_samples
        mx = maximum(fs,1);
        seq1 = sortperm(vec(mx));

        # compute the psi-function
        psi = zeros(size(rgam));
        for ii in 1:size(rgam, 2)
            fy = gradient(rgam[:, ii], binsize);
            psi[:,ii] = fy ./ sqrt(abs(fy) + eps1);
        end
        ip = zeros(n);
        len = zeros(n);

        for i in 1:n
            tmp = ones(1,M);
            psi_tmp = (psi[:,i]/M);
            ip[i] = (tmp * psi_tmp)[1];
            len[i] = acos(ip[i]);
        end

        seq2 = sortperm(len);

        # combine x-variability and y-variability
        ft = zeros(M, n);
        time0 = [0:(M-1)]./(M-1);
        for k in 1:n
            ft[:,k] = approx(time0, fs[:,seq1[k]], gams[:,seq2[k]]);
            tmp = isnan(ft[:, k]);
            while any(tmp)
                rgam2 = random_gamma(gam, 1);
                gam_tmp = invert_gamma(rgam2);
                ft[:, k] = approx(time0, fs[:,seq1[k]], gam_tmp);
            end
        end
    else
        # combine x-variability and y-variability
        ft = zeros(M, n);
        time0 = collect(0:(M-1))./(M-1);
        for k in 1:n
            ft[:,k] = approx(time0, fs[:,k], gams[:,k]);
            tmp = isnan.(ft[:, k]);
            while any(tmp)
                rgam2 = random_gamma(gam, 1);
                gam_tmp = invert_gamma(rgam2);
                ft[:, k] = approx(time0, fs[:,k], gam_tmp);
            end
        end
    end

    out = fdawarp(warp_data.f, timet, fn, qn, warp_data.q0, warp_data.fmean,
                  warp_data.mqn, gam, warp_data.orig_var, warp_data.amp_var,
                  warp_data.phase_var, warp_data.cost, warp_data.lambda,
                  warp_data.method, warp_data.omethod, warp_data.gamI, true,
                  fs, rgam, ft, q_s)

    return out

end


text:
function initws(res_x::Int, res_y::Int)
    clearws()

    ds_width_meters, ds_height_meters, ds_width_px, ds_height_px = inqdspsize()
    ds_px_width_meters = ds_width_meters / ds_width_px
    ds_px_height_meters = ds_height_meters / ds_height_px
    ws_width_meters = res_x * ds_px_width_meters
    ws_height_meters = res_y * ds_px_height_meters

    setwsviewport(0., ws_width_meters, 0., ws_height_meters)
    setwswindow(0., 1., 0., 1.)
    setviewport(0., 1., 0., 1.)
    setwindow(0., 1., 0., 1.)

    updatews()
end

function render(gr_img::GRImage)
    clearws()
    drawimage(0., 1., 0., 1., gr_img.res_x, gr_img.res_y, gr_img.abgr, 0)
    updatews()
end