text:
include("../utils.jl")

struct TestVectorSpaceType <: VectorSpaceType end

@testset "Tangent bundle" begin
    M = Sphere(2)

    @testset "Nice access to vector bundle components" begin
        TB = TangentBundle(M)
        p = ProductRepr([1.0, 0.0, 0.0], [0.0, 2.0, 4.0])
        @test p[TB, :point] === p.parts[1]
        @test p[TB, :vector] === p.parts[2]
        p[TB, :vector] = [0.0, 3.0, 1.0]
        @test p.parts[2] == [0.0, 3.0, 1.0]
        p[TB, :point] = [0.0, 1.0, 0.0]
        @test p.parts[1] == [0.0, 1.0, 0.0]
        @test_throws DomainError p[TB, :error]
        @test_throws DomainError p[TB, :error] = [1, 2, 3]
    end

    types = [Vector{Float64}]
    TEST_FLOAT32 && push!(types, Vector{Float32})
    TEST_STATIC_SIZED && push!(types, MVector{3,Float64})

    for T in types
        p = convert(T, [1.0, 0.0, 0.0])
        TB = TangentBundle(M)
        TpM = TangentSpaceAtPoint(M, p)
        @test sprint(show, TB) == "TangentBundle(Sphere(2, ℝ))"
        @test base_manifold(TB) == M
        @test manifold_dimension(TB) == 2 * manifold_dimension(M)
        @test representation_size(TB) == (6,)
        CTB = CotangentBundle(M)
        @test sprint(show, CTB) == "CotangentBundle(Sphere(2, ℝ))"
        @test sprint(show, VectorBundle(TestVectorSpaceType(), M)) ==
              "VectorBundle(TestVectorSpaceType(), Sphere(2, ℝ))"
        @testset "Type $T" begin
            pts_tb = [
                ProductRepr(convert(T, [1.0, 0.0, 0.0]), convert(T, [0.0, -1.0, -1.0])),
                ProductRepr(convert(T, [0.0, 1.0, 0.0]), convert(T, [2.0, 0.0, 1.0])),
                ProductRepr(convert(T, [1.0, 0.0, 0.0]), convert(T, [0.0, 2.0, -1.0])),
            ]
            @inferred ProductRepr(
                convert(T, [1.0, 0.0, 0.0]),
                convert(T, [0.0, -1.0, -1.0]),
            )
            for pt in pts_tb
                @test bundle_projection(TB, pt) ≈ pt.parts[1]
            end
            diag_basis = DiagonalizingOrthonormalBasis(log(TB, pts_tb[1], pts_tb[2]))
            basis_types = (
                DefaultOrthonormalBasis(),
                get_basis(TB, pts_tb[1], DefaultOrthonormalBasis()),
                diag_basis,
                get_basis(TB, pts_tb[1], diag_basis),
            )
            test_manifold(
                TB,
                pts_tb,
                test_injectivity_radius=false,
                test_reverse_diff=isa(T, Vector),
                test_forward_diff=isa(T, Vector),
                test_tangent_vector_broadcasting=false,
                test_vee_hat=true,
                test_project_tangent=true,
                test_project_point=true,
                test_default_vector_transport=true,
                vector_transport_methods=[ParallelTransport()],
                basis_types_vecs=basis_types,
                projection_atol_multiplier=4,
                test_inplace=true,
            )

            # tangent space at point
            pts_TpM = map(
                p -> convert(T, p),
                [[0.0, 0.0, 1.0], [0.0, 2.0, 0.0], [0.0, -1.0, 1.0]],
            )
            basis_types = (
                DefaultOrthonormalBasis(),
                get_basis(TpM, pts_TpM[1], DefaultOrthonormalBasis()),
            )
            @test get_basis(TpM, pts_TpM[1], basis_types[2]) === basis_types[2]
            test_manifold(
                TpM,
                pts_TpM,
                test_injectivity_radius=true,
                test_reverse_diff=isa(T, Vector),
                test_forward_diff=isa(T, Vector),
                test_tangent_vector_broadcasting=true,
                test_vee_hat=false,
                test_project_tangent=true,
                test_project_point=true,
                test_default_vector_transport=true,
                basis_types_vecs=basis_types,
                projection_atol_multiplier=4,
                test_inplace=true,
            )
        end
    end

    @test TangentBundle{ℝ,Sphere{2,ℝ}} ==
          VectorBundle{ℝ,Manifolds.TangentSpaceType,Sphere{2,ℝ}}
    @test CotangentBundle{ℝ,Sphere{2,ℝ}} ==
          VectorBundle{ℝ,Manifolds.CotangentSpaceType,Sphere{2,ℝ}}

    @test base_manifold(TangentBundle(M)) == M
    @testset "spaces at point" begin
        p = [1.0, 0.0, 0.0]
        t_p = TangentSpaceAtPoint(M, p)
        t_p2 = TangentSpace(M, p)
        @test t_p == t_p2
        ct_p = CotangentSpaceAtPoint(M, p)
        t_ps = sprint(show, "text/plain", t_p)
        sp = sprint(show, "text/plain", p)
        sp = replace(sp, '\n' => "\n ")
        t_ps_test = "Tangent space to the manifold $(M) at point:\n $(sp)"
        @test t_ps == t_ps_test
        @test base_manifold(t_p) == M
        @test base_manifold(ct_p) == M
        @test t_p.fiber.manifold == M
        @test ct_p.fiber.manifold == M
        @test t_p.fiber.fiber == TangentSpace
        @test ct_p.fiber.fiber == CotangentSpace
        @test t_p.point == p
        @test ct_p.point == p
        @test injectivity_radius(t_p) == Inf
        @test representation_size(t_p) == representation_size(M)
        v = [0.0, 0.0, 1.0]
        @test embed(t_p, v) == v
        @test embed(t_p, v, v) == v
        # generic vector space at
        fiber = VectorBundleFibers(TestVectorSpaceType(), M)
        v_p = VectorSpaceAtPoint(fiber, p)
        v_ps = sprint(show, "text/plain", v_p)
        fiber_s = sprint(show, "text/plain", fiber)
        v_ps_test = "$(typeof(v_p))\nFiber:\n $(fiber_s)\nBase point:\n $(sp)"
        @test v_ps == v_ps_test
    end

    @testset "tensor product" begin
        TT = Manifolds.TensorProductType(TangentSpace, TangentSpace)
        @test sprint(show, TT) == "TensorProductType(TangentSpace, TangentSpace)"
        @test vector_space_dimension(VectorBundleFibers(TT, Sphere(2))) == 4
        @test vector_space_dimension(VectorBundleFibers(TT, Sphere(3))) == 9
        @test base_manifold(VectorBundleFibers(TT, Sphere(2))) == M
        @test sprint(show, VectorBundleFibers(TT, Sphere(2))) ==
              "VectorBundleFibers(TensorProductType(TangentSpace, TangentSpace), Sphere(2, ℝ))"
    end

    @testset "Error messages" begin
        vbf = VectorBundleFibers(TestVectorSpaceType(), Euclidean(3))
        @test_throws ErrorException inner(vbf, [1, 2, 3], [1, 2, 3], [1, 2, 3])
        @test_throws ErrorException Manifolds.project!(vbf, [1, 2, 3], [1, 2, 3], [1, 2, 3])
        @test_throws ErrorException zero_vector!(vbf, [1, 2, 3], [1, 2, 3])
        @test_throws MethodError vector_space_dimension(vbf)
        a = fill(0.0, 6)
        @test_throws ErrorException get_coordinates!(
            TangentBundle(M),
            a,
            ProductRepr([1.0, 0.0, 0.0], [0.0, 0.0, 0.0]),
            ProductRepr([1.0, 0.0, 0.0], [0.0, 0.0, 1.0]),
            CachedBasis(DefaultOrthonormalBasis(), []),
        )
    end

    @testset "log and exp on tangent bundle for power and product manifolds" begin
        M = PowerManifold(Circle(ℝ), 2)
        N = TangentBundle(M)
        p1 = ProductRepr([0.0, 0.0], [0.0, 0.0])
        p2 = ProductRepr([-1.047, -1.047], [0.0, 0.0])
        X1 = log(N, p1, p2)
        @test isapprox(N, p2, exp(N, p1, X1))
        @test is_vector(N, p2, vector_transport_to(N, p1, X1, p2))

        M2 = ProductManifold(Circle(ℝ), Euclidean(2))
        N2 = TangentBundle(M2)
        p1_2 = ProductRepr(ProductRepr([0.0], [0.0, 0.0]), ProductRepr([0.0], [0.0, 0.0]))
        p2_2 = ProductRepr(
            ProductRepr([-1.047], [1.0, 0.0]),
            ProductRepr([-1.047], [0.0, 1.0]),
        )
        @test isapprox(N2, p2_2, exp(N2, p1_2, log(N2, p1_2, p2_2)))

        ppt = PowerVectorTransport(ParallelTransport())
        tbvt = Manifolds.VectorBundleVectorTransport(ppt, ppt)
        @test TangentBundle(M, tbvt).vector_transport === tbvt
        @test CotangentBundle(M, tbvt).vector_transport === tbvt
        @test VectorBundle(TangentSpace, M, tbvt).vector_transport === tbvt
    end
end


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

struct ElementContainer{NDIMS, RealT<:Real, uEltype<:Real, NDIMSP1, NDIMSP2, NDIMSP3}
  # Physical coordinates at each node
  node_coordinates      ::Array{RealT, NDIMSP2}   # [orientation, node_i, node_j, node_k, element]
  # ID of neighbor element in negative direction in orientation
  left_neighbors        ::Array{Int, 2}           # [orientation, elements]
  # Jacobian matrix of the transformation
  # [jacobian_i, jacobian_j, node_i, node_j, node_k, element] where jacobian_i is the first index of the Jacobian matrix,...
  jacobian_matrix       ::Array{RealT, NDIMSP3}
  # Contravariant vectors, scaled by J, in Kopriva's blue book called Ja^i_n (i index, n dimension)
  contravariant_vectors ::Array{RealT, NDIMSP3}   # [dimension, index, node_i, node_j, node_k, element]
  # 1/J where J is the Jacobian determinant (determinant of Jacobian matrix)
  inverse_jacobian      ::Array{RealT, NDIMSP1}   # [node_i, node_j, node_k, element]
  # Buffer for calculated surface flux
  surface_flux_values   ::Array{uEltype, NDIMSP2} # [variable, i, j, direction, element]
end

# Create element container and initialize element data
function init_elements(mesh::StructuredMesh{NDIMS, RealT},
                       equations::AbstractEquations,
                       basis, ::Type{uEltype}) where {NDIMS, RealT<:Real, uEltype<:Real}

  nelements = prod(size(mesh))
  node_coordinates      = Array{RealT, NDIMS+2}(undef, NDIMS, ntuple(_ -> nnodes(basis), NDIMS)..., nelements)
  left_neighbors        = Array{Int, 2}(undef, NDIMS, nelements)
  jacobian_matrix       = Array{RealT, NDIMS+3}(undef, NDIMS, NDIMS, ntuple(_ -> nnodes(basis), NDIMS)..., nelements)
  contravariant_vectors = similar(jacobian_matrix)
  inverse_jacobian      = Array{RealT, NDIMS+1}(undef, ntuple(_ -> nnodes(basis), NDIMS)..., nelements)
  surface_flux_values   = Array{uEltype, NDIMS+2}(undef, nvariables(equations),
                                                  ntuple(_ -> nnodes(basis), NDIMS-1)..., NDIMS*2, nelements)

  elements = ElementContainer{NDIMS, RealT, uEltype, NDIMS+1, NDIMS+2, NDIMS+3}(
      node_coordinates, left_neighbors, jacobian_matrix, contravariant_vectors,
      inverse_jacobian, surface_flux_values)

  init_elements!(elements, mesh, basis)
  return elements
end

@inline nelements(elements::ElementContainer) = size(elements.left_neighbors, 2)
@inline Base.ndims(::ElementContainer{NDIMS}) where NDIMS = NDIMS

Base.eltype(::ElementContainer{NDIMS, RealT, uEltype}) where {NDIMS, RealT, uEltype} = uEltype

include("containers_1d.jl")
include("containers_2d.jl")
include("containers_3d.jl")

end # @muladd

ext: jl


text:
julia> 3
3
julia> ans * ans, ans - 1
(9, 2)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
function load_bed(bed_file::AbstractString, assembly)
    io = open(bed_file)
    bed_file = ""
    # TODO: work around for readall missing in master
    if isdefined(Base, :readstring)
        bed_file = readstring(io)
    else
        bed_file = readall(io)
    end
    lines = split(bed_file, '\n')
    panel = []
    start_in_panel = 1
    panel_size = 0
    for line in lines
        line = rstrip(line, '\n')
        cols = split(line)
        if length(cols)<3
            continue
        end
        chr = ASCIIString(cols[1])
        if !haskey(assembly, chr)
            warn("$chr is not in assembly")
            continue
        end
        chrlen = length(assembly[chr])
        from = parse(Int64, ASCIIString(cols[2]))
        to = parse(Int64, ASCIIString(cols[3]))

        # extend from/to by 50*2 bp
        from += sign(from - to) * 50
        to += sign(to - from) * 50

        # clamp from 1 ~ chromo length
        from = max(1, min(chrlen, from))
        to = max(1, min(chrlen, to))

        len = abs(to - from) + 1
        contig_name = "unknown"
        if length(cols)>=4
            contig_name = ASCIIString(cols[4])
        end
        panel_size = start_in_panel+len-1
        rec = Dict(
            "chr"=>chr,
            "name"=>contig_name,
            "from"=>from, "to"=>to,
            "start_in_panel"=>start_in_panel,
            "end_in_panel"=> panel_size
            )
        push!(panel, rec)
        start_in_panel += len
    end
    return panel, panel_size
end


text:
using FastTransforms, Compat
using Compat.Test
if VERSION ≥ v"0.7-"
    using FFTW
end

@testset "BigFloat FFT and DCT" begin

    c = collect(range(-big(1.0),stop=1,length=16))
    @test norm(fft(c) - fft(Float64.(c))) < 3Float64(norm(c))*eps()
    @test norm(ifft(c) - ifft(Float64.(c))) < 3Float64(norm(c))*eps()

    c = collect(range(-big(1.0),stop=1.0,length=201))
    @test norm(ifft(fft(c))-c) < 200norm(c)eps(BigFloat)

    p = plan_dct(c)
    @test norm(FastTransforms.generic_dct(c) - p*c) == 0

    pi = plan_idct!(c)
    @test norm(pi*dct(c) - c) < 1000norm(c)*eps(BigFloat)

    @test norm(dct(c)-dct(map(Float64,c)),Inf) < 10eps()

    cc = cis.(c)
    @test norm(dct(cc)-dct(map(Complex{Float64},cc)),Inf) < 10eps()

    c = big.(rand(100)) + im*big.(rand(100))
    @test norm(dct(c)-dct(map(ComplexF64,c)),Inf) < 10eps()
    @test norm(idct(c)-idct(map(ComplexF64,c)),Inf) < 10eps()
    @test norm(idct(dct(c))-c,Inf) < 1000eps(BigFloat)
    @test norm(dct(idct(c))-c,Inf) < 1000eps(BigFloat)

    # Make sure we don't accidentally hijack any FFTW plans
    for T in (Float32, Float64)
        @test plan_fft(rand(BigFloat,10)) isa FastTransforms.DummyPlan
        @test plan_fft(rand(BigFloat,10), 1:1) isa FastTransforms.DummyPlan
        @test plan_fft(rand(Complex{BigFloat},10)) isa FastTransforms.DummyPlan
        @test plan_fft(rand(Complex{BigFloat},10), 1:1) isa FastTransforms.DummyPlan
        @test plan_fft!(rand(Complex{BigFloat},10)) isa FastTransforms.DummyPlan
        @test plan_fft!(rand(Complex{BigFloat},10), 1:1) isa FastTransforms.DummyPlan
        @test !( plan_fft(rand(T,10)) isa FastTransforms.DummyPlan )
        @test !( plan_fft(rand(T,10), 1:1) isa FastTransforms.DummyPlan )
        @test !( plan_fft(rand(Complex{T},10)) isa FastTransforms.DummyPlan )
        @test !( plan_fft(rand(Complex{T},10), 1:1) isa FastTransforms.DummyPlan )
        @test !( plan_fft!(rand(Complex{T},10)) isa FastTransforms.DummyPlan )
        @test !( plan_fft!(rand(Complex{T},10), 1:1) isa FastTransforms.DummyPlan )

        @test plan_ifft(rand(T,10)) isa FFTW.ScaledPlan
        @test plan_ifft(rand(T,10), 1:1) isa FFTW.ScaledPlan
        @test plan_ifft(rand(Complex{T},10)) isa FFTW.ScaledPlan
        @test plan_ifft(rand(Complex{T},10), 1:1) isa FFTW.ScaledPlan
        @test plan_ifft!(rand(Complex{T},10)) isa FFTW.ScaledPlan
        @test plan_ifft!(rand(Complex{T},10), 1:1) isa FFTW.ScaledPlan

        @test plan_bfft(rand(BigFloat,10)) isa FastTransforms.DummyPlan
        @test plan_bfft(rand(BigFloat,10), 1:1) isa FastTransforms.DummyPlan
        @test plan_bfft(rand(Complex{BigFloat},10)) isa FastTransforms.DummyPlan
        @test plan_bfft(rand(Complex{BigFloat},10), 1:1) isa FastTransforms.DummyPlan
        @test plan_bfft!(rand(Complex{BigFloat},10)) isa FastTransforms.DummyPlan
        @test plan_bfft!(rand(Complex{BigFloat},10), 1:1) isa FastTransforms.DummyPlan
        @test !( plan_bfft(rand(T,10)) isa FastTransforms.DummyPlan )
        @test !( plan_bfft(rand(T,10), 1:1) isa FastTransforms.DummyPlan )
        @test !( plan_bfft(rand(Complex{T},10)) isa FastTransforms.DummyPlan )
        @test !( plan_bfft(rand(Complex{T},10), 1:1) isa FastTransforms.DummyPlan )
        @test !( plan_bfft!(rand(Complex{T},10)) isa FastTransforms.DummyPlan )
        @test !( plan_bfft!(rand(Complex{T},10), 1:1) isa FastTransforms.DummyPlan )

        @test plan_dct(rand(BigFloat,10)) isa FastTransforms.DummyPlan
        @test plan_dct(rand(BigFloat,10), 1:1) isa FastTransforms.DummyPlan
        @test plan_dct(rand(Complex{BigFloat},10)) isa FastTransforms.DummyPlan
        @test plan_dct(rand(Complex{BigFloat},10), 1:1) isa FastTransforms.DummyPlan
        @test plan_dct!(rand(Complex{BigFloat},10)) isa FastTransforms.DummyPlan
        @test plan_dct!(rand(Complex{BigFloat},10), 1:1) isa FastTransforms.DummyPlan
        @test !( plan_dct(rand(T,10)) isa FastTransforms.DummyPlan )
        @test !( plan_dct(rand(T,10), 1:1) isa FastTransforms.DummyPlan )
        @test !( plan_dct(rand(Complex{T},10)) isa FastTransforms.DummyPlan )
        @test !( plan_dct(rand(Complex{T},10), 1:1) isa FastTransforms.DummyPlan )
        @test !( plan_dct!(rand(Complex{T},10)) isa FastTransforms.DummyPlan )
        @test !( plan_dct!(rand(Complex{T},10), 1:1) isa FastTransforms.DummyPlan )

        @test plan_idct(rand(BigFloat,10)) isa FastTransforms.DummyPlan
        @test plan_idct(rand(BigFloat,10), 1:1) isa FastTransforms.DummyPlan
        @test plan_idct(rand(Complex{BigFloat},10)) isa FastTransforms.DummyPlan
        @test plan_idct(rand(Complex{BigFloat},10), 1:1) isa FastTransforms.DummyPlan
        @test plan_idct!(rand(Complex{BigFloat},10)) isa FastTransforms.DummyPlan
        @test plan_idct!(rand(Complex{BigFloat},10), 1:1) isa FastTransforms.DummyPlan
        @test !( plan_idct(rand(T,10)) isa FastTransforms.DummyPlan )
        @test !( plan_idct(rand(T,10), 1:1) isa FastTransforms.DummyPlan )
        @test !( plan_idct(rand(Complex{T},10)) isa FastTransforms.DummyPlan )
        @test !( plan_idct(rand(Complex{T},10), 1:1) isa FastTransforms.DummyPlan )
        @test !( plan_idct!(rand(Complex{T},10)) isa FastTransforms.DummyPlan )
        @test !( plan_idct!(rand(Complex{T},10), 1:1) isa FastTransforms.DummyPlan )

        @test plan_rfft(rand(BigFloat,10)) isa FastTransforms.DummyPlan
        @test plan_rfft(rand(BigFloat,10), 1:1) isa FastTransforms.DummyPlan
        @test plan_brfft(rand(Complex{BigFloat},10), 19) isa FastTransforms.DummyPlan
        @test plan_brfft(rand(Complex{BigFloat},10), 19, 1:1) isa FastTransforms.DummyPlan
        @test !( plan_rfft(rand(T,10)) isa FastTransforms.DummyPlan )
        @test !( plan_rfft(rand(T,10), 1:1) isa FastTransforms.DummyPlan )
        @test !( plan_brfft(rand(Complex{T},10), 19) isa FastTransforms.DummyPlan )
        @test !( plan_brfft(rand(Complex{T},10), 19, 1:1) isa FastTransforms.DummyPlan )

    end

end


text:
using ExtensibleEffects
using DataTypesBasic
using TypeClasses
using Distributed

# Generic implementations
# -----------------------

# standard pure - fallback to TypeClass.pure
ExtensibleEffects.eff_pure(T, value) = TypeClasses.pure(T, value)

# standard eff_flatmap - fallback to map, flip_types, and flatten
function ExtensibleEffects.eff_flatmap(continuation, a)
  a_of_eff_of_a = map(continuation, a)
  eff_of_a_of_a = flip_types(a_of_eff_of_a)
  eff_of_a = map(flatten, eff_of_a_of_a)
  eff_of_a
end

# NoEffect
# --------
ExtensibleEffects.eff_applies(handler::Type{<:NoEffect}, effectful::NoEffect) = true
ExtensibleEffects.eff_pure(handler::Type{<:NoEffect}, value) = value
ExtensibleEffects.eff_flatmap(continuation, effectful::NoEffect) = continuation(effectful.value)

# Identity
# --------
# we choose to Identity{T} instead of plain T to be in accordance with behaviour of syntax_flatmap
ExtensibleEffects.eff_applies(handler::Type{<:Identity}, effectful::Identity) = true
ExtensibleEffects.eff_pure(handler::Type{<:Identity}, value) = Identity(value)
# Extra handling of Const so that order of executing Const or Identity handler does not matter
# This is especially important for ExtensibleEffects.autorun, as here it might be "random" whether we first see
# an Identity or a Const
ExtensibleEffects.eff_pure(handler::Type{<:Identity}, value::Const) = value
ExtensibleEffects.eff_flatmap(continuation, effectful::Identity) = continuation(effectful.value)

# Const
# -----
ExtensibleEffects.eff_applies(handler::Type{<:Const}, effectful::Const) = true
# usually Const does not have a pure, however within Eff, it is totally fine,
# as continuations on Const never get evaluated anyways,
# (and eff_pure is only called at the very end, when literal values are reached)
ExtensibleEffects.eff_pure(handler::Type{<:Const}, value) = value
ExtensibleEffects.eff_flatmap(continuation, effectful::Const) = effectful

# Either
# ------
# with this you can use Option/Try/Either as explicit handlers within `@runhandlers` calls
ExtensibleEffects.eff_applies(handler::Type{<:Either}, effectful::Either) = true
# eff_flatmap follows completely from Const and Identity
ExtensibleEffects.eff_pure(handler::Type{<:Either}, value) = ExtensibleEffects.eff_pure(Identity, value)  # Const would never reach this

# Iterable
# --------
ExtensibleEffects.eff_applies(handler::Type{<:Iterable}, effectful::Iterable) = true
# everything else follows from the generic implementations of eff_autohandler, eff_pure and eff_flatmap

# AbstractVector
# --------------
ExtensibleEffects.eff_applies(handler::Type{T}, effectful::T) where {T<:AbstractArray} = true
# eff_flatmap, eff_pure follow the generic implementation

# Future/Task
# -----------
ExtensibleEffects.eff_applies(handler::Type{<:Task}, effectful::Task) = true
ExtensibleEffects.eff_applies(handler::Type{<:Future}, effectful::Future) = true
# we directly interprete Task and Future
# finally surround interpreted expression by `@async`/`@spawnnat :any` to get back the Task/Future context
ExtensibleEffects.eff_pure(handler::Type{<:Union{Task, Future}}, value) = value
function ExtensibleEffects.eff_flatmap(continuation, effectful::Union{Task, Future})
  continuation(fetch(effectful))
end

# Writer
# ------
ExtensibleEffects.eff_applies(handler::Type{<:Writer}, effectful::Writer) = true
function ExtensibleEffects.eff_flatmap(continuation, a::Writer)
  eff_of_writer = continuation(a.value)
  map(eff_of_writer) do b
    Writer(a.acc ⊕ b.acc, b.value)
  end
end

"""
    WriteHandler(pure_accumulator=neutral)

Handler for generic Writers. The default accumulator works with Option values.
"""
struct WriterHandler{Acc}
  pure_acc::Acc
end
WriterHandler() = WriterHandler(neutral)  # same default pure-accumulator which is also used in TypeClasses
ExtensibleEffects.eff_applies(handler::WriterHandler, effectful::Writer) = true
ExtensibleEffects.eff_pure(handler::WriterHandler, value) = Writer(handler.pure_acc, value)
# eff_flatmap is the same

# Callable
# --------

"""
    CallableHandler(args...; kwargs...)

Handler for functions, providing the arguments and keyword arguments for calling the functions.
"""
struct CallableHandler{Args, Kwargs}
  args::Args
  kwargs::Kwargs
  CallableHandler(args...; kwargs...) = new{typeof(args), typeof(kwargs)}(args, kwargs)
end
ExtensibleEffects.eff_applies(handler::CallableHandler, effectful::Callable) = true
# we interpret callable by adding an extra Functor from the top outside, so that internally we can interpret each call
# by just getting args and kwargs from the context
ExtensibleEffects.eff_pure(handler::CallableHandler, value) = value
function ExtensibleEffects.eff_flatmap(handler::CallableHandler, continuation, effectful::Callable)
  continuation(effectful(handler.args...; handler.kwargs...))
end

"""
    @runcallable(eff)

translates to

    Callable(function(args...; kwargs...)
      @insert_into_runhandlers CallableHandler(args...; kwargs...) eff
    end)

Thanks to `@insert_into_runhandlers` this outer runner can compose well with other outer runners.
"""
macro runcallable(expr)
  esc(:(ExtensibleEffects.TypeClasses.Callable(function(args...; kwargs...)
    ExtensibleEffects.@insert_into_runhandlers(ExtensibleEffects.CallableHandler(args...; kwargs...), $expr)
  end)))
end

# ContextManager
# --------------

"""
    ContextManagerHandler(continuation)

Handler for `DataTypesBasic.ContextManager`.

The naive handler implementation for contextmanager would immediately run the continuation within the contextmanager.
However this does not work, as handling one effect does not mean that all "inner" effects are already handled.
Hence, such a handler would actually initialize and finalize the contextmanager, without its value being
processed already. When the other "inner" effects are run later on, they would find an already destroyed
contextmanager session.
We need to make sure, that the contextmanager is really the last Effect run. Therefore we create a custom
handler.
"""
struct ContextManagerHandler{F}
  cont::F
  # extra performance support for using Types as continuations.
  ContextManagerHandler(cont) = new{Core.Typeof(cont)}(cont)
end

ExtensibleEffects.eff_applies(handler::ContextManagerHandler, effectful::ContextManager) = true
ExtensibleEffects.eff_pure(handler::ContextManagerHandler, value) = handler.cont(value)
function ExtensibleEffects.eff_flatmap(::ContextManagerHandler, continuation, effectful::ContextManager)
  result = effectful(continuation)
  @assert(result.effectful isa NoEffect,
    "ContextManager should be run after all other effects,"*
    " however found result `$(result)` of type $(typeof(result))")
  result
end

"""
    @runcontextmanager(eff)

translates to

    ContextManager(function(cont)
      @insert_into_runhandlers ContextManagerHandler(cont) eff
    end)

Thanks to `@insert_into_runhandlers` this outer runner can compose well with other outer runners.
"""
macro runcontextmanager(expr)
  esc(:(ExtensibleEffects.DataTypesBasic.ContextManager(function(cont)
    ExtensibleEffects.@insert_into_runhandlers(ExtensibleEffects.ContextManagerHandler(cont), $expr)
  end)))
end

"""
    @runcontextmanager_(eff)

like `@runcontextmanager(eff)`, but immediately runs the final ContextManager
"""
macro runcontextmanager_(expr)
  esc(:(
    ExtensibleEffects.@insert_into_runhandlers(ExtensibleEffects.ContextManagerHandler(Base.identity), $expr)
  ))
end

# Combine ContextManager with Vector
# ----------------------------------

"""
It happens that the plain contextmanager handler nests all found ContextManager, i.e. the first ContextManager will
only be finalized if all other ContextManagers run, even if it is completely independent of them.
For example consider
```julia
@runcontextmanager_ @syntax_eff begin
  a = [1,2,3]
  b = mycontextmanager(a)
end
```
then the execution is
```
start mycontextmanager(1)
start mycontextmanager(2)
start mycontextmanager(3)
finish mycontextmanager(3)
finish mycontextmanager(2)
finish mycontextmanager(1)
```
That is still semantically fine, however you may use ContextManager in order to save memory and really would like to
have resources released as soon as possible. So you wish to have the following execution order
```
start mycontextmanager(1)
finish mycontextmanager(1)
start mycontextmanager(2)
finish mycontextmanager(2)
finish mycontextmanager(3)
start mycontextmanager(3)
```
This is not easily possible by a mere ContextManager handler, because of several reasons. To give an intution, note
two things:
1. if you run a handler in ExtensibleEffects, it is automatically applied everywhere where possible
2. if you run a continuation, it is run until it finds the next non-run Effect, or returns the final value.

So if we just do `@runcontextmanager_` as above, it will try to evaluate the given continuation on the first
contextmanager, where the continuation itself also already run the contextmanager. Hence there is no non-run Effect left
in the continuation, and everything is run until the one final result. This results in the execution order seen above
for `@runcontextmanager_`.

To circumvent this we only have to combine the ContextManagerHandler with the Vector handler, as it is the Vector
handler, which merges all different branches into the one single continuation which makes the problems above.
If we just run both handlers at once, the continuation which is seen by the ContextManager is the continuation of the
current branch, and hence is only run until the final value of the current branch.
This corrects the execution order successfully.

The only downside is that you cannot use autorun, but you have to default to using `@syntax_eff_noautorun` and then
run your handlers manually.

IMPORTANT: As this runs the ContextManager it is crucial that all other handlers are run beforehand.
"""

"""
  ContextManagerCombinedHandler(otherhandler)

We can combine ContextManager with any other Handler. This is possible because ContextManager, within eff_flatmap,
does not constrain the returned eff of the continuation.
"""
struct ContextManagerCombinedHandler{OtherHandler, Func}
  other_handler::OtherHandler
  contextmanager_handler::ContextManagerHandler{Func}
end
function ContextManagerCombinedHandler(other_handler, func::Union{Type, Function} = identity)
  ContextManagerCombinedHandler{Core.Typeof(other_handler), Core.Typeof(func)}(other_handler, ContextManagerHandler(func))
end

function ExtensibleEffects.eff_applies(handler::ContextManagerCombinedHandler, effectful)
  eff_applies(handler.other_handler, effectful) || eff_applies(handler.contextmanager_handler, effectful)
end

function ExtensibleEffects.eff_pure(handler::ContextManagerCombinedHandler, value)
  value′ = ExtensibleEffects.eff_pure(handler.contextmanager_handler, value)
  ExtensibleEffects.eff_pure(handler.other_handler, value′)
end
function ExtensibleEffects.eff_flatmap(handler::ContextManagerCombinedHandler, continuation, effectful)
  if eff_applies(handler.other_handler, effectful)
    ExtensibleEffects.eff_flatmap(handler.other_handler, continuation, effectful)
  elseif eff_applies(handler.contextmanager_handler, effectful)
    ExtensibleEffects.eff_flatmap(handler.contextmanager_handler, continuation, effectful)
  else
    error("ContextManagerCombinedHandler should only be eff_flatmap on values which can either be handled "*
    "by ContextManagerHandler or by other_handler = `$(handler.other_handler)`. However got effectful `$effectful`")
  end
end

# State
# -----

"""
    StateHandler(state)

Handler for running State. Gives the initial state.
"""
struct StateHandler{T}
  state::T
end
ExtensibleEffects.eff_applies(handler::StateHandler, effectful::State) = true
ExtensibleEffects.eff_pure(handler::StateHandler, value) = (value, handler.state)

# The updating of the state cannot be described by plain `eff_flatmap`.
# We need to define our own runhandler instead. It is a bit more complex, but still straightforward and compact.
function ExtensibleEffects.runhandler(handler::StateHandler, eff::Eff)
  eff_applies(handler, eff.effectful) || return runhandler_not_applies(handler, eff)

  value, nextstate = eff.effectful(handler.state)
  nexthandler = StateHandler(nextstate)
  if isempty(eff.cont)
    _eff_pure(nexthandler, value)
  else
    runhandler(nexthandler, eff.cont(value))
  end
end

"""
    @runstate eff

Note that unlike Callable, a State has to ensure that it is always the first outer Eff being run,
as it returns the inner state as an additional argument.

If you would nest it with runcallable, e.g. like `@runstate @runcallable eff` it wouldn't work,
as now the appended state is within the Callable and not directly within the State.
"""
macro runstate(expr)
  # Note, that we have to use `runhandlers` explicitly, such that other outer handlers using `@insert_into_runhandlers`
  # can interact well with this outer handler.
  esc(quote
    let eff = $expr
      isa(eff, ExtensibleEffects.Eff) || error("""
        `@runstate` only works for `Eff` type, got `$(typeof(eff))`.
        Try to use `@runstate` as you first outer handler, which is directly applied to the `Eff`.
        """)
      ExtensibleEffects.State() do state
        ExtensibleEffects.runhandlers(ExtensibleEffects.StateHandler(state), eff)
      end
    end
  end)
end

max_stars_repo_name: schlichtanders/ExtensibleEffects

max_issues_repo_name: schlichtanders/ExtensibleEffects

max_forks_repo_name: schlichtanders/ExtensibleEffects


text:
using fpzip_jll

# lib = "/usr/local/lib/libfpzip.1.3.0.dylib"

const FPZIP_TYPE_FLOAT = 0
const FPZIP_TYPE_DOUBLE = 1

# array meta data and stream handle
mutable struct FPZ
    type::Cint
    prec::Cint
    nx::Cint
    ny::Cint
    nz::Cint
    nw::Cint
    nf::Cint
end

FPZ(meta::Ptr) = unsafe_load(Ptr{FPZ}(meta))

const FPZIP_MEDIUM_PRECISION = 16
const FPZIP_HIGH_PRECISION = 24

depth = 768
spectrum = Float32(100.0) * randn(Float32, depth)

function fpzip_compress(src::Array{Float32,1}, precision::Integer)

    ndims = length(size(src))
    ndims in [1] || throw(DimensionMismatch("FPZIP compression only for a 1D array (for the time being)."))

    local status, outbytes

    bufsize = 1024 + length(src) * sizeof(eltype(src))
    dest = Vector{UInt8}(undef, bufsize)
    # dest = zeros(UInt8, bufsize)

    # compress to memory
    fpz = ccall((:fpzip_write_to_buffer, libfpzip), Ptr{FPZ}, (Ptr{Cvoid}, Csize_t), dest, bufsize)

    # convert a pointer into Julia struct to access nx,ny,nz,nf etc.
    meta = FPZ(fpz)

    println("an empty fpz: ", meta)

    # fill in the metadata
    meta.type = FPZIP_TYPE_FLOAT
    meta.prec = precision
    meta.nx = length(src)
    meta.ny = 1
    meta.nz = 1
    meta.nf = 1

    println("after filling: ", meta)

    # sync the Julia structure back to C
    unsafe_store!(Ptr{FPZ}(fpz), meta)

    meta2 = FPZ(fpz)
    println("meta copied to fpz; reconstructed from fpz: ", meta2)

    println("got here#1")

    # write header
    status = ccall((:fpzip_write_header, libfpzip), Cint, (Ptr{Cvoid},), fpz)

    println("got here#2")

    if status == 0
        throw(error("Writing FPZIP header failed."))
    else
        # compress the data
        outbytes = ccall((:fpzip_write, libfpzip), Csize_t, (Ptr{Cvoid}, Ptr{Cvoid}), fpz, src)
        println("outbytes = $outbytes")
    end

    # release the metadata structure
    ccall((:fpzip_write_close, libfpzip), Cvoid, (Ptr{Cvoid},), fpz)

    if outbytes > 0
        return dest[1:outbytes]
    else
        return Nothing
    end
end

fpzip_compress(spectrum, FPZIP_MEDIUM_PRECISION)


text:
#-------------------------------------------------------------------------------
# defining the radiative domain and associated functions

export AtmosphericDomain

"""
Structure defining the temperature and pressure ranges over which absorption cross-sections are generated when constructing gas objects. `AtmosphericDomain` objects store the temperature and pressure coordinates of cross-section interpolation grids. More points lead to higher accuracy interpolation. Generally, about 12 temperature points and 24 pressure points results in maximum error of ~1 % and much smaller average error.

| Field | Type | Description |
| ----- | :--- | :---------- |
| `T` | `Vector{Float64}` | temperature coordinates of grid [K] |
| `Tmin` | `Float64` | lowest temperature value |
| `Tmax` | `Float64` | highest temperature value |
| `nT` | `Int64` | number of temperature coordinates |
| `P` | `Vector{Float64}` | pressure coordinates of grid [Pa] |
| `Pmin` | `Float64` | lowest pressure value |
| `Pmax` | `Float64` | highest pressure value |
| `nP` | `Int64` | number of pressure coordinates |

# Constructor

    AtmosphericDomain(Trange, nT, Prange, nP)

Creates a domain with the given temperature/pressure ranges and numbers of points. `Trange` and `Prange` should be tuples of two values. `nT` and `nP` indicate the number of points to use.
"""
struct AtmosphericDomain
    #temperature samples for interpolator [K]
    T::Vector{Float64}
    Tmin::Float64
    Tmax::Float64
    nT::Int64
    #pressure samples for interpolator [atm]
    P::Vector{Float64}
    Pmin::Float64
    Pmax::Float64
    nP::Int64
end

function Base.show(io::IO, Ω::AtmosphericDomain)
    print(io, "AtmosphericDomain:\n")
    print(io, "  $(Ω.nT) temperature nodes ∈ [$(Ω.Tmin),$(Ω.Tmax)] K\n")
    print(io, "  $(Ω.nP) pressure nodes ∈ [$(Ω.Pmin),$(Ω.Pmax)] Pa")
end

function AtmosphericDomain(Trange::NTuple{2,Real}, nT::Int,
                           Prange::NTuple{2,Real}, nP::Int)
    #check for negatives
    @assert all(Trange .> 0) "temperature range must be positive"
    @assert all(Prange .> 0) "pressure range must be positive"
    #check the Qref/Q range
    @assert all(Trange .>= TMIN) "minimum temperature with Qref/Q accuracy is $TMIN K"
    @assert all(Trange .<= TMAX) "maximum temperature with Qref/Q accuracy is $TMAX K"
    #order
    @assert Trange[1] < Trange[2] "Trange[1] ($(Trange[1])) can't be greater than Trange[2] ($(Trange[2]))"
    @assert Prange[1] < Prange[2] "Prange[1] ($(Prange[1])) can't be greater than Prange[2] ($(Prange[2]))"
    #generate grid points
    T = chebygrid(Trange[1], Trange[2], nT)
    P = exp.(chebygrid(log(Prange[1]), log(Prange[2]), nP))
    #assemble!
    AtmosphericDomain(T, Trange[1], Trange[2], nT, P, Prange[1], Prange[2], nP)
end

#-------------------------------------------------------------------------------
#wrapper type for BichebyshevInterpolators used for cross-sections

export OpacityTable

struct OpacityTable{M,N}
    Φ::BichebyshevInterpolator{M,N,Float64}
end

# T: temperature grid coordinates [K]
# P: pressure grid coordinates [Pa]
# σ: cross-section grid [cm^2/molecule]
function OpacityTable(T, P, σ)
    tiny = floatmin(eltype(T))
    ζ = all(σ .<= tiny) #whether all values are effectively zero
    lnP = log.(P)
    lnσ = ζ ? fill(log(tiny), size(σ)) : log.(σ)
    Φ = BichebyshevInterpolator(T, lnP, lnσ)
    OpacityTable(Φ)
end

#gets cross-section out of interpolators, un-logged [cm^2/molecule]
(Π::OpacityTable)(T, P) = exp(Π.Φ(T, log(P)))

#-------------------------------------------------------------------------------
#function for building gas opacity tables

function checkν(ν::AbstractVector)::Nothing
    #check wavenumbers for problems
    @assert all(diff(ν) .> 0) "wavenumbers must be unique and in ascending order"
    @assert all(ν .>= 0) "wavenumbers must be positive"
    nothing
end

function bake(sl::SpectralLines,
              fC::F,
              shape!::G,
              Δνcut::Real,
              ν::Vector{Float64},
              Ω::AtmosphericDomain,
              progress::Bool=true
              ) where {F, G<:Function}
    checkν(ν)
    #number of wavenumbers
    nν = length(ν)
    #create a single block of cross-sections
    σ = zeros(nν, Ω.nT, Ω.nP)
    #initialize progress meter if desired
    if progress
        prg = Progress(Ω.nT*Ω.nP, 0.1, "evaluating $(sl.formula) cross-sections ")
    end
    #fill σ by evaluating in batches of wavenumbers (slow part)
    @threads for i ∈ eachindex(Ω.T)
        for j ∈ eachindex(Ω.P)
            #get a view into the big σ array
            σᵢⱼ = @view σ[:,i,j]
            #get temperature, pressure, concentration
            T = Ω.T[i]
            P = Ω.P[j]
            C = fC(T, P)
            #make sure concentration isn't wacky
            @assert 0 <= C <= 1 "gas molar concentrations must be in [0,1], not $C (encountered @ $T K, $P Pa)"
            #evaluate line shapes (slow part)
            shape!(σᵢⱼ, ν, sl, T, P, C*P, Δνcut)
            #update progress meter
            progress && next!(prg)
        end
    end
    #check for weirdness
    z = zeros(Bool, nν)
    for i ∈ eachindex(ν)
        σᵥ = @view σ[i,:,:]
        if (minimum(σᵥ) == 0) & (maximum(σᵥ) > 0)
            z[i] = true
        end
    end
    if any(z)
        @info "Zero cross-section values are mixed with non-zero values for the following wavenumbers for $(sl.name):\n\n$(ν[z])\n\n Likely, absorption is extremely weak in these regions, causing underflow. Absorption is being set to zero for all temperatures and pressures at those wavenumbers to avoid non-smooth and inaccurate interpolation tables."
        σ[z,:,:] .= 0.0
    end
    #split the block and create opacity tables for each ν
    [OpacityTable(Ω.T, Ω.P, σ[i,:,:]) for i ∈ eachindex(ν)]
end

#-------------------------------------------------------------------------------
#for testing opacity table errors

export opacityerror

function opacityerror(Π::OpacityTable,
                      Ω::AtmosphericDomain,
                      sl::SpectralLines,
                      ν::Real,
                      C::F, #C(T,P)
                      shape::G=voigt,
                      N::Int=50) where {F,G}
    #create T and P grids from the domain
    T = LinRange(Ω.Tmin, Ω.Tmax, N)
    P = 10 .^ LinRange(log10(Ω.Pmin), log10(Ω.Pmax), N)
    #compute exact and approximate cross-sections
    σop = zeros(N,N)
    σex = zeros(N,N)
    @threads for i = 1:N
        for j = 1:N
            σop[i,j] = Π(T[i], P[j])
            σex[i,j] = shape(ν, sl, T[i], P[j], C(T[i],P[j])*P[j])
        end
    end
    #compute error and relative error
    aerr = σop .- σex
    rerr = aerr./σex
    return T, P, aerr, rerr
end

#-------------------------------------------------------------------------------
#defining Gas struct and access to cross-sections

abstract type AbstractGas end

export Gas
export rawσ, concentration, reconcentrate

"""
Gas type for real-gas, radiatively active, atmospheric constituents. Must be constructed from a `.par` file or a [`SpectralLines`](@ref) object.

# Constructors

    Gas(sl::SpectralLines, fC, ν, Ω, shape!=voigt!, Δνcut=25; progress=true)

* `sl`: a [`SpectralLines`](@ref) object
* `fC`: molar concentration of the constituent [mole/mole] as a function of temperature and pressure `fC(T,P)`
* `ν`: vector of wavenumber samples [cm``^{-1}``]
* `Ω`: [`AtmosphericDomain`](@ref)
* `shape!`: line shape to use, must be the in-place version ([`voigt!`](@ref), [`lorentz!`](@ref), etc.)
* `Δνcut`: profile truncation distance [cm``^{-1}``]
* `progress`: whether to display the progress meter

    Gas(par::String, fC, ν, Ω, shape!=voigt!, Δνcut=25; progress=true)

Same arguments as the first constructor, but reads a `par` file directly into the gas object. Keyword arguments are passed through to [`readpar`](@ref).
"""
struct Gas{T,F} <: AbstractGas
    name::String
    formula::String
    μ::Float64 #mean molar mass
    ν::Vector{Float64}
    Ω::AtmosphericDomain
    Π::Vector{T} #cross-section interpolators
    fC::F #molar concentration [mole/mole] function fC(T,P)
end

function Base.show(io::IO, g::Gas)
    print(io, "$(g.name) ($(g.formula))\n")
    print(io, "  μ = $(round(g.μ, sigdigits=8)) kg/mole\n")
    print(io, "  ν: $(minimum(g.ν)) - $(maximum(g.ν)) cm^-1 ($(length(g.ν)) samples)\n")
    Tmin, Tmax = g.Ω.Tmin, g.Ω.Tmax
    Pmin, Pmax = g.Ω.Pmin, g.Ω.Pmax
    print(io, "  T: $Tmin - $Tmax K ($(g.Ω.nT) samples)\n")
    print(io, "  P: $Pmin - $Pmax Pa ($(g.Ω.nP) samples)\n")
end

function Gas(sl::SpectralLines,
             fC::F,
             ν::AbstractVector{<:Real},
             Ω::AtmosphericDomain,
             shape!::Function=voigt!,
             Δνcut::Real=25;
             progress::Bool=true) where {F}
    @assert length(ν) > 0
    μ = sum(sl.A .* sl.μ)/sum(sl.A) # mean molar mass [kg/mole]
    ν = collect(Float64, ν)
    Π = bake(sl, fC, shape!, Δνcut, ν, Ω, progress)
    T = eltype(Π)
    Gas{T,F}(sl.name, sl.formula, μ, ν, Ω, Π, fC)
end

function Gas(filename::String,
             fC::F,
             ν::AbstractVector{<:Real},
             Ω::AtmosphericDomain,
             shape!::Function=voigt!,
             Δνcut=25;
             kwargs...) where {F}
    sl = SpectralLines(filename; kwargs...)
    Gas(sl, fC, ν, Ω, shape!, Δνcut; kwargs...)
end

"""
    rawσ(g::Gas, i, T, P)

Retrieve the absorption cross-section at wavenumber index `i`, temperature `T` [K], and pressure `P` [Pa], **without** the concentration scaling.
"""
rawσ(g::Gas, i::Int, T, P) = g.Π[i](T, P)

"""
    rawσ(g::Gas, T, P)

Retrieve the absorption cross-section for all wavenumbers at temperature `T` [K], and pressure `P` [Pa], **without** the concentration scaling.
"""
rawσ(g::Gas, T, P) = [Π(T, P) for Π ∈ g.Π]

"""
    concentration(g::Gas, T, P)

Furnishes the molar concentration [mole/mole] of a [`Gas`](@ref) object at a particular temperature and pressure. Identical to `g.C(T,P)`.
"""
concentration(g::Gas, T, P) = g.fC(T,P)

function concentration(g::Gas, T::AbstractVector, P::AbstractVector)
    @assert length(T) == length(P)
    [concentration(g, T[i], P[i]) for i ∈ eachindex(P)]
end

#single concentration-scaled cross-section
(g::Gas)(i::Int, T, P) = concentration(g, T, P)*rawσ(g, i, T, P)

#for full vectors of concentration-scaled cross-sections with whatever gas
(g::Gas)(T, P) = concentration(g, T, P)*rawσ(g, T, P)

"""
    reconcentrate(g::Gas, fC)

Create a copy of a [`Gas`](@ref) object with a new molar concentration function, `fC(T,P)` [mole/mole].

!!! warning

    The self-broadening component of the line shape is not recomputed when using the `reconcentrate` function. This component is generally small when partial pressure is low, but may be appreciable if the concentration changes significantly.
"""
function reconcentrate(g::Gas, fC)::Gas
    #check for invalid concentrations or unexpected input
    for P ∈ g.Ω.P, T ∈ g.Ω.T
        C = try
            fC(T,P)
        catch
            error("concentration function must have the form fC(T,P), where T is temperature [K] and P is pressure [Pa], but something went wrong when the function was evaluated")
        end
        @assert 0 <= C <= 1.0 "gas molar concentrations must be in [0,1], not $C, which was encountered at T=$T P=$P"
    end
    #construct a new gas WITHOUT COPYING
    Gas(g.name, g.formula, g.μ, g.ν, g.Ω, g.Π, fC)
end

"""
    reconcentrate(g::Gas, C)

Create a copy of a [`Gas`](@ref) object with a new, **uniform** molar concentration [mole/mole].

!!! warning

    The self-broadening component of the line shape is not recomputed when using the `reconcentrate` function. This component is generally small when partial pressure is low, but may be appreciable if the concentration changes significantly.
"""
function reconcentrate(g::Gas, C::Real)::Gas where {U,V}
    #check for invalid concentrations
    @assert 0 <= C <= 1.0 "gas molar concentrations must be in [0,1], not $C"
    #construct a new gas WITHOUT COPYING
    Gas(g.name, g.formula, g.μ, g.ν, g.Ω, g.Π, (T,P)->float(C))
end

function Base.getindex(g::Gas, x)
    Gas(
        g.name[:],
        g.formula[:],
        g.μ,
        g.ν[x],
        g.Ω,
        g.Π[x],
        g.fC
    )
end

Base.firstindex(g::Gas) = 1

Base.lastindex(g::Gas) = length(g.ν)

#---------------------------------------

export GrayGas

struct GrayGas{T} <: AbstractGas
    name::String
    formula::String
    μ::Float64
    ν::Vector{Float64}
    σ::T
end

function Base.show(io::IO, g::GrayGas)
    print(io, "Gray Gas\n")
    print(io, "  ν: $(minimum(g.ν)) - $(maximum(g.ν)) cm^-1 ($(length(g.ν)) samples)\n")
    print(io, "  σ = $(g.σ) cm^2/molecule")
end

function GrayGas(σ::Real, ν::AbstractVector{<:Real})
    GrayGas("Gray", "Gray", NaN, collect(Float64, ν), σ)
end

(g::GrayGas)(x...) = g.σ

#---------------------------------------

export SemiGrayGas

struct SemiGrayGas{T} <: AbstractGas
    name::String
    formula::String
    μ::Float64
    ν::Vector{Float64}
    νcut::Float64
    σ::T
end

function Base.show(io::IO, g::SemiGrayGas)
    println(io, "Semi-Gray Gas")
    println(io, "  ν: $(minimum(g.ν)) - $(maximum(g.ν)) cm^-1 ($(length(g.ν)) samples)")
    println(io, "  σ = $(g.σ) cm^2/molecule")
    print(io, "  absorption cutoff = $(g.νcut) cm^-1")
end

function SemiGrayGas(σ::Real, ν::AbstractVector{<:Real}, νcut::Real)
    SemiGrayGas("SemiGray", "SemiGray", NaN, collect(Float64, ν), Float64(νcut), σ)
end

(g::SemiGrayGas)(i::Int, _...) = (g.ν[i] <= g.νcut) ? g.σ : zero(g.σ)


text:
module BERT

export
    BertPreTraining,
    BertClassification,
    BertConfig,
    load_from_torch_base,
    load_from_torch_pretraining,
    load_from_torch_classification,
    BertAdam,
    bert_tokenize

using Knet, SpecialFunctions, LinearAlgebra
include("model.jl")
include("optimizer.jl")
include("preprocess.jl")

end # module

max_issues_repo_name: OsmanMutlu/BERT-in-Knet

max_forks_repo_name: OsmanMutlu/BERT-in-Knet


text:
function return_cache(::typeof(*),a::AbstractArray{<:Number},b::AbstractArray{<:Number})
  c = a*b
  CachedArray(c)
end

@inline function evaluate!(cache,::typeof(*),a::AbstractMatrix{<:Number},b::AbstractVector{<:Number})
  m = size(a,1)
  setsize!(cache,(m,))
  c = cache.array
  mul!(c,a,b)
  c
end

@inline function evaluate!(cache,::typeof(*),a::AbstractMatrix{<:Number},b::AbstractMatrix{<:Number})
  m = size(a,1)
  n = size(b,2)
  setsize!(cache,(m,n))
  c = cache.array
  mul!(c,a,b)
  c
end

struct MulAddMap{T} <: Map
  α::T
  β::T
end

function return_cache(k::MulAddMap,a,b,c)
  d = a*b+c
  CachedArray(d)
end

@inline function evaluate!(cache,k::MulAddMap,a,b,c)
  setsize!(cache,size(c))
  d = cache.array
  copyto!(d,c)
  mul!(d,a,b,k.α,k.β)
  d
end

struct AddEntriesMap{F} <: Map
  combine::F
end

function evaluate!(cache,k::AddEntriesMap,A,v,i,j)
  add_entries!(k.combine,A,v,i,j)
end

function evaluate!(cache,k::AddEntriesMap,A,v,i)
  add_entries!(k.combine,A,v,i)
end

struct TouchEntriesMap <: Map end

function evaluate!(cache,k::TouchEntriesMap,A,v,i,j)
  add_entries!(+,A,nothing,i,j)
end

function evaluate!(cache,k::TouchEntriesMap,A,v,i)
  add_entries!(+,A,nothing,i)
end


text:
module GalerkinSparseGrids

# The prerequisite packages as of May 2018 are Cubature.jl and ODE.jl

using Cubature
using ODE

# The following script files are used

include("Schemes.jl")					# CartesianIndex Manipulation Methods
include("Additional_Tools.jl")			# Helper functions not specific to any part of the package
# include("Hat_Methods.jl")				# Using non-galerkin elementary 'hat' basis functions
include("DG_Basis.jl")					# Gram-Schmidt procedure for DG basis functions in 1-D
include("1D_DG_Functions.jl") 			# Explicitly building the 1-D Galerkin Basis
include("DG_Methods.jl")				# Multidimensional hierarchical & sparse coefficients
include("DG_vMethods.jl")				# Going between a dictionary & a vector of coeffs
include("Derivative_Matrix_Elements.jl")# 1-D symbolic piecewise derivative
include("Derivative_Precompute.jl")		# Precomputing derivative matrix for coeff vect
include("1D_Derivative.jl")				# Constructing ideal 1D derivative matrix using boundary terms
include("Multidim_Derivative.jl")		# Multidimensional DG Derivatives in full & sparse bases
include("PDEs.jl")						# Solving the n-D wave equation with periodic boundary
include("1D_Nodal_Basis.jl")			# Construct 1D nodal basis for multiplication
include("Multidim_Nodal_Basis.jl")		# Construct Multidimensional nodal basis for multiplication

include("Error_Measure.jl")				# Monte Carlo Methods to measure error
include("Tensor_Construct.jl")			# Quickly calculates coeffs of simple tensors of functions
include("Traveling_Wave_Example.jl")	# Construct and evolve traveling waves

export
# standard_coeffs,
# standard_reconstruct,
# coeffs_hat,
# reconstruct_hat,

coeffs_DG,
reconstruct_DG,

get_size,
V2D,
D2V,
V2Dref,
D2Vref,
vcoeffs_DG,

D_matrix,
grad_matrix,
laplacian_matrix,

wave_evolve_1D,
wave_evolve,

nodal2points_1D,
points2nodal_1D,
nodal2heir_1D,
hier2points_1D,
points2hier_1D,

mcerr,
mcerr2,

tensor_construct,

cos_coeffs,
sin_coeffs,
traveling_wave_solver

end # module


text:
@assert ARGS[1] in ["chimera", "square"]
@assert ARGS[2] in ["forwarddiff", "nilang"]

const AD = ARGS[2]

using CUDAnative, CuArrays
device!(parse(Int, get(ARGS, 3, "0")))
using CuYao
using DelimitedFiles, TropicalTensors
using BenchmarkTools
using ForwardDiff

if ARGS[1] == "chimera"
    const LATTICE = ChimeraLattice
else
    const LATTICE = SquareLattice
end

suite = BenchmarkGroup()

function bfunc(::Type{T}, L::Int) where T
    lt = LATTICE(L, L)
    if T <: Integer
        Js = rand(T[-1, 1], length(sgbonds(lt)))
    else
        Js = randn(T, length(sgbonds(lt)))
    end
    hs = zeros(T, length(lt))
    sg = Spinglass(lt, Js, hs)
    if AD == "forwarddiff"
        CuArrays.@sync ForwardDiff.gradient(x->solve(lt,
            convert.(eltype(x), Js),
            x; usecuda=true).n, hs)
    else
        opt_config(sg)
    end
end

if LATTICE === SquareLattice
    const RANGE = 4:2:25
    const TYPE = Float32
else
    const RANGE = 2:6
    const TYPE = Int16
end

for L = RANGE
    suite[L] = @benchmarkable bfunc($TYPE, $L)
end

println("loaded")

tune!(suite)
res = run(suite)

function analyze_res(res)
    times = zeros(length(res))
    for (k,L) = enumerate(RANGE)
        times[k] = minimum(res[L].times)
    end
    return times
end

times = analyze_res(res)
fname = joinpath(@__DIR__, "data", "bench_$(LATTICE)_$AD.dat")
println("Writing benchmark results to file: $fname.")
writedlm(fname, times)


text:
export iauPvup
"""
Update a pv-vector, discarding the velocity component.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  vector/matrix support function.

Given:
   dt       double            time interval
   pv       double[2][3]      pv-vector

Returned:
   p        double[3]         p-vector

Notes:

   1. "Update" means "refer the position component of the vector to a
   new date dt time units from the existing date".

   2. The time units of dt must match those of the velocity.

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauPvup(double dt, double pv[2][3], double p[3])
function iauPvup(dt::Real, pv::AbstractMatrix{<:Real})
   # Preallocate return values
   p = zeros(Float64, 3)

   ccall((:iauPvup, libsofa_c), Cvoid,
         (Cdouble, Ptr{Cdouble}, Ref{Cdouble}),
          convert(Float64, dt),
          convert(Array{Float64, 2}, pv'),
          p)

   return SVector{3}(p)
end


text:
import Base: consume, foldl, @deprecate, @deprecate_binding
export lift, consume, foldl, keepwhen, keepif, dropif, dropwhen

@deprecate_binding Input Signal

@deprecate lift(f, s::Signal...; kwargs...) map(f,s...; kwargs...)
@deprecate consume(f::Union{Function, DataType}, s::Signal...;kwargs...) map(f, s...;kwargs...)
@deprecate foldl(f, x, s::Signal...;kwargs...) foldp(f, x, s...;kwargs...)
@deprecate keepwhen filterwhen
@deprecate keepif filter
@deprecate dropif(f, default, signal) filter(x -> !f(x), default, signal)
@deprecate dropwhen(predicate, x, signal) filterwhen(map(!, predicate), x, signal)


text:
using Test
using Liblet

# Constructor

test_automaton_from_grammar() = begin
    A = Automaton(Grammar("""
        S -> a A
        S -> a B
        A -> b B
        A -> b C
        B -> c A
        B -> c C
        C -> a
    """))
    s = "Automaton(N={A, B, C, S, ◇}, T={a, b, c}, transitions=(S-a->A, S-a->B, A-b->B, A-b->C, B-c->A, B-c->C, C-a->◇), F={◇}, q0=S)"
    s == string(A)
end

test_automaton_from_grammar_fail3() = Automaton(Grammar("S -> a b c"))

test_automaton_from_grammar_fail2a() = Automaton(Grammar("S -> a b"))

test_automaton_from_grammar_fail2b() = Automaton(Grammar("S -> B B"))

test_automaton_from_ε_grammar() = begin
    A = Automaton(Grammar("""
    S -> A
    S -> a B
    A -> a A
    A -> ε
    B -> b B
    B -> b
    """))
    s = "Automaton(N={A, B, S, ◇}, T={a, b}, transitions=(S-ε->A, S-a->B, A-a->A, A-ε->◇, B-b->B, B-b->◇), F={◇}, q0=S)"
    s == string(A)
end

test_automaton_from_string() = begin
    A = Automaton("""
        A, 0, B
        A, 1, F
        B, 0, G
        B, 1, C
        C, 1, C
        D, 0, C
        D, 1, G
        E, 0, H
        E, 1, F
        F, 0, C
        F, 1, G
        G, 0, G
        H, 0, G
        H, 1, C
    """, Set(["C"]))
    s = "Automaton(N={A, B, C, D, E, F, G, H}, T={0, 1}, transitions=(A-0->B, A-1->F, B-0->G, B-1->C, C-1->C, D-0->C, D-1->G, E-0->H, E-1->F, F-0->C, F-1->G, G-0->G, H-0->G, H-1->C), F={C}, q0=A)"
    s == string(A)
end

test_automaton_overlapTN() = Automaton(Set(["A", "B"]), Set(["B", "C"]), [], "A",  Set())

test_automaton_q0N() = Automaton(Set(["A", "B"]), Set(["b", "c"]), [], "X", Set())

test_automaton_FN() = Automaton(Set(["A", "B"]), Set(["b", "c"]), [], "A", Set(["B", "C"]))

# δ

test_automaton_δ() = begin
    states = Automaton(Grammar("""
        S -> a A
        S -> a B
        A -> b B
        A -> b C
        B -> c A
        B -> c C
        C -> a
    """)) |> x -> δ(x, "S", "a")
    states == Set(["A","B"])
end

# Tests

function runautomatontests()
    @testset "automaton_tests" begin
        #@test test_automaton_from_grammar()
        @test test_automaton_δ()
        @test_throws ArgumentError test_automaton_from_grammar_fail3()
        @test_throws ArgumentError test_automaton_from_grammar_fail2a()
        @test_throws ArgumentError test_automaton_from_grammar_fail2b()
        #@test test_automaton_from_ε_grammar()
        #@test test_automaton_from_string()
        @test_throws ArgumentError test_automaton_overlapTN()
        @test_throws ArgumentError test_automaton_q0N()
        @test_throws ArgumentError test_automaton_FN()
    end
end


text:
# each leaf node is a partition.

import Random
import PyPlot
import Statistics

import AbstractTrees

using LinearAlgebra
#import Interpolations
#using Revise

import Colors

#include("../src/RKHSRegularization.jl")
import RKHSRegularization # https://github.com/RoyCCWang/RKHSRegularization

#PyPlot.matplotlib["rcParams"][:update](["font.size" => 22, "font.family" => "serif"])

include("./helpers/visualization.jl")

PyPlot.close("all")

PyPlot.matplotlib["rcParams"][:update](["font.size" => 22, "font.family" => "serif"])

Random.seed!(25)

fig_num = 1

### plotting.

D = 2
N = 500
X = collect( randn(D) for n = 1:N )
levels = 2 # 2^(levels-1) leaf nodes. Must be larger than 1.

root, X_parts, _ = RKHSRegularization.setuppartition(X, levels)

# # print using AbstractTrees.
# AbstractTrees.printnode(io::IO, node::RKHSRegularization.BinaryNode) = print(io, node.data)
# AbstractTrees.print_tree(root)

### visualize tree.
# traverse all leaf nodes.
# X_parts = Vector{Vector{Vector{Float64}}}(undef, 0)
# RKHSRegularization.buildXpart!(X_parts, root)

centroid = Statistics.mean( Statistics.mean(X_parts[i]) for i = 1:length(X_parts) )
max_dist = maximum( maximum( norm(X_parts[i][j]-centroid) for j = 1:length(X_parts[i])) for i = 1:length(X_parts) ) * 1.1

y_set = Vector{Vector{Float64}}(undef, 0)
t_set = Vector{Vector{Float64}}(undef, 0)
min_t = -5.0
max_t = 5.0
max_N_t = 5000
RKHSRegularization.getpartitionlines!(y_set, t_set, root, levels, min_t, max_t, max_N_t, centroid, max_dist)

fig_num, _ = visualize2Dpartition(X_parts, y_set, t_set, fig_num, "levels = $(levels)")

levels = 3 # 2^(levels-1) leaf nodes. Must be larger than 1.

root, X_parts, _ = RKHSRegularization.setuppartition(X, levels)

# # print using AbstractTrees.
# AbstractTrees.printnode(io::IO, node::RKHSRegularization.BinaryNode) = print(io, node.data)
# AbstractTrees.print_tree(root)

### visualize tree.
centroid = Statistics.mean( Statistics.mean(X_parts[i]) for i = 1:length(X_parts) )
max_dist = maximum( maximum( norm(X_parts[i][j]-centroid) for j = 1:length(X_parts[i])) for i = 1:length(X_parts) ) * 1.1

y_set = Vector{Vector{Float64}}(undef, 0)
t_set = Vector{Vector{Float64}}(undef, 0)
min_t = -5.0
max_t = 5.0
max_N_t = 5000
RKHSRegularization.getpartitionlines!(y_set, t_set, root, levels, min_t, max_t, max_N_t, centroid, max_dist)

fig_num, _ = visualize2Dpartition(X_parts, y_set, t_set, fig_num, "levels = $(levels)")

levels = 4 # 2^(levels-1) leaf nodes. Must be larger than 1.

root, X_parts, _ = RKHSRegularization.setuppartition(X, levels)

# # print using AbstractTrees.
# AbstractTrees.printnode(io::IO, node::RKHSRegularization.BinaryNode) = print(io, node.data)
# AbstractTrees.print_tree(root)

### visualize tree.
centroid = Statistics.mean( Statistics.mean(X_parts[i]) for i = 1:length(X_parts) )
max_dist = maximum( maximum( norm(X_parts[i][j]-centroid) for j = 1:length(X_parts[i])) for i = 1:length(X_parts) ) * 1.1

y_set = Vector{Vector{Float64}}(undef, 0)
t_set = Vector{Vector{Float64}}(undef, 0)
min_t = -5.0
max_t = 5.0
max_N_t = 5000
RKHSRegularization.getpartitionlines!(y_set, t_set, root, levels, min_t, max_t, max_N_t, centroid, max_dist)

fig_num, ax4 = visualize2Dpartition(X_parts, y_set, t_set, fig_num, "levels = $(levels)")

levels = 5 # 2^(levels-1) leaf nodes. Must be larger than 1.

root, X_parts, _ = RKHSRegularization.setuppartition(X, levels)

# # print using AbstractTrees.
# AbstractTrees.printnode(io::IO, node::RKHSRegularization.BinaryNode) = print(io, node.data)
# AbstractTrees.print_tree(root)

### visualize tree.
centroid = Statistics.mean( Statistics.mean(X_parts[i]) for i = 1:length(X_parts) )
max_dist = maximum( maximum( norm(X_parts[i][j]-centroid) for j = 1:length(X_parts[i])) for i = 1:length(X_parts) ) * 1.1

y_set = Vector{Vector{Float64}}(undef, 0)
t_set = Vector{Vector{Float64}}(undef, 0)
min_t = -5.0
max_t = 5.0
max_N_t = 5000
RKHSRegularization.getpartitionlines!(y_set, t_set, root, levels, min_t, max_t, max_N_t, centroid, max_dist)

fig_num, ax = visualize2Dpartition(X_parts, y_set, t_set, fig_num, "levels = $(levels)")

#### test partition search given points. Write as unit stress test later.

x = [0.42; 2.05] # leaf node 11
x = [3.99; 3.35] # leaf node 16
x = [0.1; 0.16] # leaf node 9
x_region_ind = RKHSRegularization.findpartition(x, root, levels)

#

#@assert 1==2

####### point, find all hyperplanes intersection with sphere.

levels = 4 # 2^(levels-1) leaf nodes. Must be larger than 1.

root, X_parts, _ = RKHSRegularization.setuppartition(X, levels)

# # print using AbstractTrees.
# AbstractTrees.printnode(io::IO, node::RKHSRegularization.BinaryNode) = print(io, node.data)
# AbstractTrees.print_tree(root)

### visualize tree.
centroid = Statistics.mean( Statistics.mean(X_parts[i]) for i = 1:length(X_parts) )
max_dist = maximum( maximum( norm(X_parts[i][j]-centroid) for j = 1:length(X_parts[i])) for i = 1:length(X_parts) ) * 1.1

y_set = Vector{Vector{Float64}}(undef, 0)
t_set = Vector{Vector{Float64}}(undef, 0)
min_t = -5.0
max_t = 5.0
max_N_t = 5000
RKHSRegularization.getpartitionlines!(y_set, t_set, root, levels, min_t, max_t, max_N_t, centroid, max_dist)

fig_num, ax = visualize2Dpartition(X_parts, y_set, t_set, fig_num, "levels = $(levels)")

### work on a training set that is different than the partition.
ε = 0.3
X_set, X_set_inds, region_list_set,
problematic_inds = RKHSRegularization.organizetrainingsets(root, levels, X, ε)

fig_num = visualizesingleregions(X_set, y_set, t_set, fig_num)

N2 = collect( length(X_set[n]) for n = 1:length(X_set) )
N1 = collect( length(X_parts[n]) for n = 1:length(X_parts) )
println("N: X_parts: ", N1)
println("N: X_set:   ", N2)
println()

#p = randn(D) # test point.
p = x

# find the region for p.
p_region_ind = RKHSRegularization.findpartition(p, root, levels)
println("p_region_ind = ", p_region_ind)
println()

# get all hyperplanes.
hps = RKHSRegularization.fetchhyperplanes(root)
@assert length(hps) == length(X_parts) - 1 # sanity check.

radius = ε #0.3
δ = 1e-5
region_inds, ts, zs, hps_keep_flags = RKHSRegularization.findneighbourpartitions(p, radius, root, levels, hps, p_region_ind; δ = δ)

# debug.
hps_kept = hps[hps_keep_flags]
hp = hps_kept[1]
u = hp.v
c = hp.c

z_kept = zs[hps_keep_flags]
t_kept = ts[hps_keep_flags]
#z = z_kept[1]

dists = collect( norm(z_kept[i]-p) for i = 1:length(z_kept) )
@assert norm( dists - abs.(t_kept) ) < 1e-10 # sanity check.

# visualize on fresh plot.
fig_num, ax = visualize2Dpartition(X_parts, y_set, t_set, fig_num, "levels = $(levels)")

PyPlot.scatter([p[1];], [p[2];], marker = "x", s = 600.0, label = "p")

for i = 1:length(z_kept)
    PyPlot.scatter([z_kept[i][1];], [z_kept[i][2];], marker = "d", s = 600.0, label = "z = $(i)")
end

PyPlot.axis([p[1]-radius; p[1]+radius; p[2]-radius; p[2]+radius])
PyPlot.axis("scaled")

PyPlot.legend()

max_stars_repo_name: RoyCCWang/RKHSRegularization

max_issues_repo_name: RoyCCWang/RKHSRegularization

max_forks_repo_name: RoyCCWang/RKHSRegularization


text:
# Copyright 2018 IBM Corp.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

module IntegrationTest
using Test

using Catlab.WiringDiagrams
import Catlab.WiringDiagrams: add_box!
using SemanticFlowGraphs

# Convenience methods to read raw and semantic flow graphs.
read_raw_graph(name::String) =
  read_raw_graphml(joinpath(@__DIR__, "$name.raw.graphml"))
read_semantic_graph(name::String) =
  read_semantic_graphml(joinpath(@__DIR__, "$name.semantic.graphml"))

# Convenience method for building semantic flow graphs.
add_box!(d::WiringDiagram, f::Monocl.Hom) = add_box!(d, Box(f))

# Load all concepts at the outset.
const db = OntologyDB()
load_concepts(db)

@testset "PyFlowGraphs" begin
  include("PyFlowGraphs.jl")
end

@testset "RFlowGraphs" begin
  include("RFlowGraphs.jl")
end

end

max_stars_repo_name: epatters/flowgraphs

max_issues_repo_name: epatters/flowgraphs

max_forks_repo_name: epatters/flowgraphs


text:
using StaticArrays
using Base.Threads: @spawn

abstract type AbstractCoordinationStrategy end

struct VarEl <: AbstractCoordinationStrategy
end

Base.@kwdef struct MaxPlus <:AbstractCoordinationStrategy
    message_iters::Int64 = 10
    message_norm::Bool = true
    use_agent_utils::Bool = false
    node_exploration::Bool = true
    edge_exploration::Bool = true
end

"""
Factored Value Monte Carlo Tree Search solver datastructure

Fields:
    n_iterations::Int64
        Number of iterations during each action() call.
        default: 100

    max_time::Float64
        Maximum CPU time to spend computing an action.
        default::Inf

    depth::Int64
        Number of iterations during each action() call.
        default: 100

    exploration_constant::Float64:
        Specifies how much the solver should explore. In the UCB equation, Q + c*sqrt(log(t/N)), c is the exploration constant.
        The exploration terms for FV-MCTS-Var-El and FV-MCTS-Max-Plus are different but the role of c is the same.
        default: 1.0

    rng::AbstractRNG:
        Random number generator

    estimate_value::Any (rollout policy)
        Function, object, or number used to estimate the value at the leaf nodes.
        If this is a function `f`, `f(mdp, s, depth)` will be called to estimate the value.
        If this is an object `o`, `estimate_value(o, mdp, s, depth)` will be called.
        If this is a number, the value will be set to that number
        default: RolloutEstimator(RandomSolver(rng))

    init_Q::Any
        Function, object, or number used to set the initial Q(s,a) value at a new node.
        If this is a function `f`, `f(mdp, s, a)` will be called to set the value.
        If this is an object `o`, `init_Q(o, mdp, s, a)` will be called.
        If this is a number, Q will be set to that number
        default: 0.0

    init_N::Any
        Function, object, or number used to set the initial N(s,a) value at a new node.
        If this is a function `f`, `f(mdp, s, a)` will be called to set the value.
        If this is an object `o`, `init_N(o, mdp, s, a)` will be called.
        If this is a number, N will be set to that number
        default: 0

    reuse_tree::Bool
        If this is true, the tree information is re-used for calculating the next plan.
        Of course, clear_tree! can always be called to override this.
        default: false

    coordination_strategy::AbstractCoordinationStrategy
        The specific strategy with which to compute the best joint action from the current MCTS statistics.
        default: VarEl()
"""
Base.@kwdef mutable struct FVMCTSSolver <: AbstractMCTSSolver
    n_iterations::Int64 = 100
    max_time::Float64 = Inf
    depth::Int64 = 10
    exploration_constant::Float64 = 1.0
    rng::AbstractRNG = Random.GLOBAL_RNG
    estimate_value::Any = RolloutEstimator(RandomSolver(rng))
    init_Q::Any = 0.0
    init_N::Any = 0
    reuse_tree::Bool = false
    coordination_strategy::AbstractCoordinationStrategy = VarEl()
end

mutable struct FVMCTSTree{S,A,CS<:CoordinationStatistics}

    # To map the multi-agent state vector to the ID of the node in the tree
    state_map::Dict{S,Int64}

    # The next two vectors have one for each node ID in the tree
    total_n::Vector{Int}                # The number of times the node has been tried
    s_labels::Vector{S} # The state corresponding to the node ID

    # List of all individual actions of each agent for coordination purposes.
    all_agent_actions::Vector{A}

    coordination_stats::CS
    lock::ReentrantLock
end

function FVMCTSTree(all_agent_actions::Vector{A},
                    coordination_stats::CS,
                    init_state::S,
                    lock::ReentrantLock,
                    sz::Int64=10000) where {S, A, CS <: CoordinationStatistics}

    return FVMCTSTree{S,A,CS}(Dict{S,Int64}(),
                                 sizehint!(Int[], sz),
                                 sizehint!(S[], sz),
                                 all_agent_actions,
                                 coordination_stats,
                                 lock
                                 )
end # function

Base.isempty(t::FVMCTSTree) = isempty(t.state_map)
state_nodes(t::FVMCTSTree) = (FVStateNode(t, id) for id in 1:length(t.total_n))

struct FVStateNode{S}
    tree::FVMCTSTree{S}
    id::Int64
end

# Accessors for state nodes
@inline state(n::FVStateNode) = n.tree.s_labels[n.id]
@inline total_n(n::FVStateNode) = n.tree.total_n[n.id]

## No need for `children` or ActionNode just yet

mutable struct FVMCTSPlanner{S, A, SE, CS <: CoordinationStatistics, RNG <: AbstractRNG} <: AbstractMCTSPlanner{JointMDP{S,A}}
    solver::FVMCTSSolver
    mdp::JointMDP{S,A}
    tree::FVMCTSTree{S,A,CS}
    solved_estimate::SE
    rng::RNG
end

"""
Called internally in solve() to create the FVMCTSPlanner where Var-El is the specific action coordination strategy.
Creates VarElStatistics internally with the CG components and the minimum degree ordering heuristic.
"""
function varel_joint_mcts_planner(solver::FVMCTSSolver,
                                  mdp::JointMDP{S,A},
                                  init_state::S,
                                  ) where {S,A}

    # Get coordination graph components from maximal cliques
    #adjmat = coord_graph_adj_mat(mdp)
    #@assert size(adjmat)[1] == n_agents(mdp) "Adjacency Matrix does not match number of agents!"

    #adjmatgraph = SimpleGraph(adjmat)
    adjmatgraph = coordination_graph(mdp)

    coord_graph_components = maximal_cliques(adjmatgraph)
    min_degree_ordering = sortperm(degree(adjmatgraph))

    # Initialize full agent actions
    all_agent_actions = Vector{(actiontype(mdp))}(undef, n_agents(mdp))
    for i = 1:n_agents(mdp)
        all_agent_actions[i] = agent_actions(mdp, i)
    end

    ve_stats = VarElStatistics{S}(coord_graph_components, min_degree_ordering,
                                                   Dict{typeof(init_state),Vector{Vector{Int64}}}(),
                                                   Dict{typeof(init_state),Vector{Vector{Int64}}}(),
                                                   )

    # Create tree from the current state
    tree = FVMCTSTree(all_agent_actions, ve_stats,
                         init_state, ReentrantLock(), solver.n_iterations)
    se = convert_estimator(solver.estimate_value, solver, mdp)

    return FVMCTSPlanner(solver, mdp, tree, se, solver.rng)
end # end function

"""
Called internally in solve() to create the FVMCTSPlanner where Max-Plus is the specific action coordination strategy.
Creates MaxPlusStatistics and assumes the various MP flags are sent down from the CoordinationStrategy object given to the solver.
"""
function maxplus_joint_mcts_planner(solver::FVMCTSSolver,
                                    mdp::JointMDP{S,A},
                                    init_state::S,
                                    message_iters::Int64,
                                    message_norm::Bool,
                                    use_agent_utils::Bool,
                                    node_exploration::Bool,
                                    edge_exploration::Bool,
                                    ) where {S,A}

    @assert (node_exploration || edge_exploration) "At least one of nodes or edges should explore!"

#=     adjmat = coord_graph_adj_mat(mdp)
    @assert size(adjmat)[1] == n_agents(mdp) "Adjacency Mat does not match number of agents!" =#

    #adjmatgraph = SimpleGraph(adjmat)
    adjmatgraph = coordination_graph(mdp)
    @assert size(adjacency_matrix(adjmatgraph))[1] == n_agents(mdp)

    # Initialize full agent actions
    # TODO(jkg): this is incorrect? Or we need to override actiontype to refer to agent actions?
    all_agent_actions = Vector{(actiontype(mdp))}(undef, n_agents(mdp))
    for i = 1:n_agents(mdp)
        all_agent_actions[i] = agent_actions(mdp, i)
    end

    mp_stats = MaxPlusStatistics{S}(adjmatgraph,
                                                    message_iters,
                                                    message_norm,
                                                    use_agent_utils,
                                                    node_exploration,
                                                    edge_exploration,
                                                    Dict{S,PerStateMPStats}())

    # Create tree from the current state
    tree = FVMCTSTree(all_agent_actions, mp_stats,
                      init_state, ReentrantLock(), solver.n_iterations)
    se = convert_estimator(solver.estimate_value, solver, mdp)

    return FVMCTSPlanner(solver, mdp, tree, se, solver.rng)
end

# Reset tree.
function clear_tree!(planner::FVMCTSPlanner)

    # Clear out state map dict entirely
    empty!(planner.tree.state_map)

    # Empty state vectors with state hints
    sz = min(planner.solver.n_iterations, 100_000)

    empty!(planner.tree.s_labels)
    sizehint!(planner.tree.s_labels, planner.solver.n_iterations)

    # Don't touch all_agent_actions and coord graph component
    # Just clear comp stats dict
    clear_statistics!(planner.tree.coordination_stats)
end

MCTS.init_Q(n::Number, mdp::JointMDP, s, c, a) = convert(Float64, n)
MCTS.init_N(n::Number, mdp::JointMDP, s, c, a) = convert(Int, n)

# No computation is done in solve; the solver is just given the mdp model that it will work with
# and in case of MaxPlus, the various flags for the MaxPlus behavior
function POMDPs.solve(solver::FVMCTSSolver, mdp::JointMDP)
    if typeof(solver.coordination_strategy) == VarEl
        return varel_joint_mcts_planner(solver, mdp, rand(solver.rng, initialstate(mdp)))
    elseif typeof(solver.coordination_strategy) == MaxPlus
        return maxplus_joint_mcts_planner(solver, mdp, rand(solver.rng, initialstate(mdp)),
                                          solver.coordination_strategy.message_iters,
                                          solver.coordination_strategy.message_norm,
                                          solver.coordination_strategy.use_agent_utils,
                                          solver.coordination_strategy.node_exploration,
                                          solver.coordination_strategy.edge_exploration)
    else
        throw(error("Not Implemented"))
    end
end

# IMP: Overriding action for FVMCTSPlanner here
# NOTE: Hardcoding no tree reuse for now
function POMDPs.action(planner::FVMCTSPlanner, s)
    clear_tree!(planner) # Always call this at the top
    plan!(planner, s)
    action =  coordinate_action(planner.mdp, planner.tree, s)
    return action
end

function POMDPModelTools.action_info(planner::FVMCTSPlanner, s)
    clear_tree!(planner) # Always call this at the top
    plan!(planner, s)
    action = coordinate_action(planner.mdp, planner.tree, s)
    return action, nothing
end

function plan!(planner::FVMCTSPlanner, s)
    planner.tree = build_tree(planner, s)
end

# build_tree can be called on the assumption that no reuse AND tree is reinitialized
function build_tree(planner::FVMCTSPlanner, s::S) where S

    n_iterations = planner.solver.n_iterations
    depth = planner.solver.depth

    root = insert_node!(planner.tree, planner, s)

    # Simulate can be multi-threaded
    @sync for n = 1:n_iterations
        @spawn simulate(planner, root, depth)
    end
    return planner.tree
end

function simulate(planner::FVMCTSPlanner, node::FVStateNode, depth::Int64)

    mdp = planner.mdp
    rng = planner.rng
    s = state(node)
    tree = node.tree

    # once depth is zero return
    if isterminal(planner.mdp, s)
        return 0.0
    elseif depth == 0
        return estimate_value(planner.solved_estimate, planner.mdp, s, depth)
    end

    # Choose best UCB action (NOT an action node as in vanilla MCTS)
    ucb_action = coordinate_action(mdp, planner.tree, s, planner.solver.exploration_constant, node.id)

    # Monte Carlo Transition
    sp, r = @gen(:sp, :r)(mdp, s, ucb_action, rng)

    spid = lock(tree.lock) do
        get(tree.state_map, sp, 0) # may be non-zero even with no tree reuse
    end
    if spid == 0
        spn = insert_node!(tree, planner, sp)
        spid = spn.id

        q = r .+ discount(mdp) * estimate_value(planner.solved_estimate, planner.mdp, sp, depth - 1)
    else
        q = r .+ discount(mdp) * simulate(planner, FVStateNode(tree, spid) , depth - 1)
    end

    # NOTE: Not bothering with tree visualization right now
    # Augment N(s)
    lock(tree.lock) do
        tree.total_n[node.id] += 1
    end

    # Update component statistics! (non-trivial)
    # This is related but distinct from initialization
    update_statistics!(mdp, tree, s, ucb_action, q)

    return q
end

@POMDP_require simulate(planner::FVMCTSPlanner, s, depth::Int64) begin
    mdp = planner.mdp
    P = typeof(mdp)
    @assert P <: JointMDP
    SV = statetype(P)
    @req iterate(::SV)
    #@assert typeof(SV) <: AbstractVector
    AV = actiontype(P)
    @assert typeof(AV) <: AbstractVector
    @req discount(::P)
    @req isterminal(::P, ::SV)
    @subreq insert_node!(planner.tree, planner, s)
    @subreq estimate_value(planner.solved_estimate, mdp, s, depth)
    @req gen(::P, ::SV, ::AV, ::typeof(planner.rng)) # XXX this is not exactly right - it could be satisfied with transition

    ## Requirements from MMDP Model
    @req agent_actions(::P, ::Int64)
    @req agent_actions(::P, ::Int64, ::eltype(SV))
    @req n_agents(::P)
    @req coordination_graph(::P)

    # TODO: Should we also have this requirement for SV?
    @req isequal(::S, ::S)
    @req hash(::S)
end

function insert_node!(tree::FVMCTSTree{S,A,CS}, planner::FVMCTSPlanner,
                      s::S) where {S,A,CS <: CoordinationStatistics}

    lock(tree.lock) do
        push!(tree.s_labels, s)
        tree.state_map[s] = length(tree.s_labels)
        push!(tree.total_n, 1)

        # NOTE: Could actually make actions state-dependent if need be
        init_statistics!(tree, planner, s)
    end

    # length(tree.s_labels) is just an alias for the number of state nodes
    ls = lock(tree.lock) do
        length(tree.s_labels)
    end
    return FVStateNode(tree, ls)
end

@POMDP_require insert_node!(tree::FVMCTSTree, planner::FVMCTSPlanner, s) begin

    P = typeof(planner.mdp)
    AV = actiontype(P)
    A = eltype(AV)
    SV = typeof(s)
    #S = eltype(SV)

    # TODO: Review IQ and IN
    IQ = typeof(planner.solver.init_Q)
    if !(IQ <: Number) && !(IQ <: Function)
        @req init_Q(::IQ, ::P, ::SV, ::Vector{Int64}, ::AbstractVector{A})
    end

    IN = typeof(planner.solver.init_N)
    if !(IN <: Number) && !(IN <: Function)
        @req init_N(::IQ, ::P, ::SV, ::Vector{Int64}, ::AbstractVector{A})
    end

    @req isequal(::S, ::S)
    @req hash(::S)
end


text:
function extract_nodedbcs!(node_dbcs::Dict{String, Vector{Tuple{TI,TF}}}, file) where {TI, TF}
    pattern_zero = r"([^,]+)\s*,\s*(\d)"
    pattern_range = r"([^,]+)\s*,\s*(\d)\s*,\s*(\d)"
    pattern_other = r"([^,]+)\s*,\s*(\d)\s*,\s*(\d)\s*,\s*(\-?\d+\.\d*)"
    line = readline(file)
    m = match(stopping_pattern, line)
    while m isa Nothing
        m = match(pattern_other, line)
        if m != nothing
            nodesetname = m[1]
            dof1 = parse(TI, m[2])
            dof2 = parse(TI, m[3])
            val = parse(TF, m[4])
            if haskey(node_dbcs, nodesetname)
                for dof in dof1:dof2
                    push!(node_dbcs[nodesetname], (dof, val))
                end
            else
                node_dbcs[nodesetname] = [(dof1, val)]
                for dof in dof1+1:dof2
                    push!(node_dbcs[nodesetname], (dof, val))
                end
            end
            line = readline(file)
            m = match(stopping_pattern, line)
            continue
        end
        m = match(pattern_range, line)
        if m != nothing
            nodesetname = m[1]
            dof1 = parse(TI, m[2])
            dof2 = parse(TI, m[3])
            if haskey(node_dbcs, nodesetname)
                for dof in dof1:dof2
                    push!(node_dbcs[nodesetname], (dof, zero(TF)))
                end
            else
                node_dbcs[nodesetname] = [(dof1, zero(TF))]
                for dof in dof1+1:dof2
                    push!(node_dbcs[nodesetname], (dof, zero(TF)))
                end
            end
            line = readline(file)
            m = match(stopping_pattern, line)
            continue
        end
        m = match(pattern_zero, line)
        if m != nothing
            nodesetname = m[1]
            dof = parse(TI, m[2])
            if haskey(node_dbcs, nodesetname)
                push!(node_dbcs[nodesetname], (dof, zero(TF)))
            else
                node_dbcs[nodesetname] = [(dof, zero(TF))]
            end
            line = readline(file)
            m = match(stopping_pattern, line)
            continue
        end
        line = readline(file)
        m = match(stopping_pattern, line)
    end
    return line
end


text:
using Boltzmann
using MNIST

X, y = traindata()
X = X[:, 1:1000]                     # take only 1000 observations for speed
X = X / (maximum(X) - (minimum(X)))  # normalize to [0..1]

layers = [("vis", BaseRBM(Units(784, :gaussian), Units(256))),
          ("hid1", BaseRBM(Units(256), Units(100))),
          ("hid2", BaseRBM(Units(100), Units(100)))]
dbn = DBN(layers)
fit(dbn, X)
transform(dbn, X)


text:
#========================================================================================#
#	Laboratory 11
#
# Deterministic CHAOS and the course project.
#
# Author: Niall Palfreyman, 23/03/2022
#========================================================================================#
[
	Activity(
		"""
		This is the last laboratory in this course, and it presents the assessed project for
		the course. In this laboratory, you will implement a program to demonstrate chaos in a
		gravitational system containing three bodies of equal mass in two dimensions. You will use
		Julia to analyse the execution of a complex simulation program, and adapt this program to
		use Runge-Kutta integration to (a) simulate 3-body motion and (b) evaluate its chaoticity.
		You will present your methods and results in ONE deliverable: A complete Julia module with
		associated demonstration code that presents a pedagogically pleasing description of the
		n-body problem, its solution using Runge-Kutta-2 integration and an analysis of the
		dynamics of your chosen system. This deliverable will also contain a title, information on
		authors and references, and diagrams.
		""",
		"",
		x -> true
	),
	Activity(
		"""
		In your project, you will investigate DETERMINISTIC CHAOS. Chaos is important for biology
		because it is the source of the spontaneity that we observe in living organisms. To see
		how chaos works, first define a Julia version of the Mathematica function nestlist():

			function nestlist( f::Function, x0, n::Integer)
				(n ≤ 0) ? [x0] : begin
					list = Vector{typeof(x0)}(undef, n+1)
					list[1] = x0
					for i in 1:n
						list[i+1] = f(list[i])
					end
					list
				end
			end

		Experiment with nestlist(). For example, what is the result of the following invocation?

			nestlist(x->2x,3,5)
		""",
		"",
		x -> x == [3,6,12,24,48,96]
	),
	Activity(
		"""
		nestlist() applies the function f repeatedly to the initial value x0, creating a list of
		the values that it generates in this way. You can see how this might be useful if we want
		to create a simulation of timesteps that repeat themselves over time. We'll do that now ...

		Imagine a population of wasps living on an island with renewable but limited resources.
		They can fill the island's carrying capacity (population x = 1), they can die out (x = 0),
		or else their population can have any value in the range 0 ≤ x ≤ 1. Also, they all die each
		winter, so their population takes on a new value each year: [x0,x1,x2,...,xn].

		For specific birthrate r, the wasp population's growth is governed by this logistic equation:

			x[i+1] = Λ(r)(x[i]), where
			Λ(r) is the function: (x -> r x (1-x))

		Implement this model now and experiment with it, then tell me the population two years
		after an initial population of x0 = 0.3, if r = 0.5
		""",
		"Λ(r) = (x -> r .* x .* (1 .- x))",
		x -> (0.04698 < x < 0.04699)
	),
	Activity(
		"""
		Use nestlist() to create a list of 5 simulation steps of the wasp population, starting from
		an initial value of 0.3 and using specific growth rate r = 0.5. You will see that the
		population x is tending towards a particular limiting value - what is that limit value?
		""",
		"nestlist(Λ(1),0.3,5)",
		x -> x==0
	),
	Activity(
		"""
		You should now be able to do the following: Use nestlist() to create a list of 20
		simulation steps of the wasp population with x0=0.3, r=0.9, then plot them in a graph.
		What is the limit point of this sequence?
		""",
		"lines(nestlist(Λ(0.9),0.3,20))",
		x -> x==0
	),
	Activity(
		"""
		OK, now we have the equipment, we can investigate the onset of chaos. We will slowly
		increase the value of the specific birthrate r to discover how this affects the
		developmental motion of the wasp population. Use the initial value x0=0.3 for all of the
		folloiwng excercises until I tell you otherwise.

		Just to make our language clear: a LIMIT POINT of the wasp population's motion is any
		value of the population that stays the same from one generation to the next: Λ(r)(x) == x.

		What was the limit point of the motion Λ(0.9)?
		""",
		"",
		x -> x==0
	),
	Activity(
		"""
		What is the approximate value of the limit point of the motion Λ(1.1)?
		""",
		"",
		x -> 0.07<x<0.1
	),
	Activity(
		"""
		What is the EXACT value of the limit point of the motion Λ(1.5)?
		""",
		"You can work it out for yourself by solving the equation Λ(r)(x) == x, or: r x (1-x) = x",
		x -> x == 1//3
	),
	Activity(
		"""
		Use the technique from the hint in the previous activity to calculate the limit point of
		the motion Λ(2), then check this value using your simulator:
		""",
		"",
		x -> x==0.5
	),
	Activity(
		"""
		Now investigate the motion Λ(2.1). In which direction does the population change from x[6]
		to x[7]: positive (+) or negative (-)?
		""",
		"",
		x -> (x == -)
	),
	Activity(
		"""
		This is interesting 00! Up to now, the development of the wasps has been monotone: the
		population has only either shrunk or grown. But now it grows above the value, then drops
		back down again. Let's investigate this further: find the limit point of the motion
		Λ(2.5) both by calculating and by simulating:
		""",
		"",
		x -> x==0.6
	),
	Activity(
		"""
		Notice that now we have an oscillating motion that dies away as x comes to rest at the
		limit point. Try out various values of r between 2.5 and 2.95. Does the oscillation
		still die away?
		""",
		"",
		x -> occursin("y",lowercase(x))
	),
	Activity(
		"""
		Now find the limiting motion of the population for Λ(3). Does the oscillation die away?
		""",
		"",
		x -> occursin("n",lowercase(x))
	),
	Activity(
		"""
		It now no longer makes sense to speak of a limit VALUE, since the motion oscillates
		forever around the value 2/3. Instead, we speak of a limit CYCLE: the motion Λ(3) displays
		a limit cycle around the value 2/3. "cycle" means the value oscillates forever, and
		"limit" means that it will converge to this cycle, no matter where we start the motion.
		Check this for yourself by changing the value of x0. Does the motion always converge on
		a cycle around 2/3?
		""",
		"",
		x -> occursin("y",lowercase(x))
	),
	Activity(
		"""
		How many oscillatory periods are contained within one complete limit cycle of the motion
		Λ(3)? That is, how many times does the population wobble up and down within ONE cycle? If
		you have difficulty understanding exactly what I am asking, try looking at the hint for
		this activity.
		""",
		"We call this limit cycle a \"period-1 limit-cycle\"!",
		x -> x==1
	),
	Activity(
		"""
		Let's investigate further. Look at the limit-cycle of the motion Λ(3.2) (You may want
		to extend the simulation's duration to 50). How many periods are in this limit cycle?
		""",
		"The number of periods should not (yet) have changed",
		x -> x==1
	),
	Activity(
		"""
		How many periods are in the limit-cycle of the motion Λ(3.45)? That is, how many complete
		oscillations does the population make before the motion repeats itself?
		""",
		"Count very carefully: This should be a period-2 limit-cycle!",
		x -> x==2
	),
	Activity(
		"""
		How many periods are in the limit-cycle of the motion Λ(3.55)? By now, things will move
		quite quickly and the wasp population will need a while to stabilise towards the limit-
		cycle. You may like to save time by using a line like this:

			lines(nestlist(Λ(3.55),0.5,5000)[end-40:end])
		""",
		"Look very carefully!",
		x -> x==4
	),
	Activity(
		"""
		How many periods are in the limit-cycle of the motion Λ(3.567)?
		""",
		"Count ve-ery carefully - remembering to inspect the smaller oscillations as well!",
		x -> x==8
	),
	Activity(
		"""
		How many periods are in the limit-cycle of the motion Λ(3.5695)? You will find this
		one difficult to count, but it is actually a period-16 limit-cycle. If you can't
		distinguish the oscillations, don't worry - just go on to the next activity.
		""",
		"",
		x -> x==16
	),
	Activity(
		"""
		We have seen that as r increases, period-doubling occurs, so the motion takes longer
		and longer before it repeats. Now check out the motion Λ(3.7). How long does it take
		before this motion repeats itself? Or in other words: What is the period-length of
		this limit-cycle?
		""",
		"Find out how to write \"Infinity\" in Julia :)",
		x -> x==Inf
	),
	Activity(
		"""
		We have come a long way, and what we have learned is that some naturally occurring systems
		can enter a type of motion in which they require an infinite amount of time to repeat, and
		so we call them CHAOTIC. Notice that this phenomenon has nothing to do with randomness. The
		population of Vespula Island is perfectly deterministic - we can always predict what will
		happen in the next step, but ONLY by actually simulating it! We cannot calculate in advance
		what will happen after 1000 iterations without actually performing ALL of the 1000 steps
		that lead up to it.

		Actually, the situation is even worse than this. Use the following command to inspect the
		behaviour of Λ(3.7) for the different starting conditions x0 ∈ [0.5,0.50001,4.99999]

			lines(nestlist(Λ(3.7),0.5,5000)[end-40:end])

		Are these three graphs at all similar to each other?
		""",
		"",
		x -> occursin("n",lowercase(x))
	),
	Activity(
		"""
		So you see, chaotic motion meanns that even tinsey-tiny changes in the initial conditions
		of a chaotic system lead to completely different behaviour. So even if we measured the
		current wasp population, we could never be sure that we had done it sufficiently accurately
		to be ABSOLUTELY sure that we are accurately predicting the development of the system!

		So. What has all this to do with your course project? Henri Poincaré was the first to
		notice deterministic chaos. In 1908, he showed that the gravitational motion of three
		orbiting bodies cannot be solved exactly. Weather, dripping taps and driven pendula have a
		similar problem. We cannot predict the story of their future motion without simulation,
		which as you know is never precise! As an introduction to chaos in the three-body problem,
		please view the following video-clip now, and then proceed to the next activity:

		https://www.youtube.com/watch?v=LwkvO3t1b30&t=113s
		""",
		"",
		x -> true
	),
	Activity(
		"""
		The demo() function of the module NBodies creates an instance of the type NBody to simulate
		the motion of TWO bodies in ONE spatial dimension. Run it now: do the blue and red bodies
		orbit each other? Investigate this further by increasing the resolution of the simulation
		and of the figure. Try focussing your graphic on the time-snapshots around where the two
		bodies meet, and magnify the image. Do they orbit now? What exactly is happening here?
		This is the origin of chaotic motion!

		Check that the vector x returned by the method NBodies.simulate() has three levels: the
		elements of x span the individual time snapshots of the simulation; each snapshot contains
		a list of the bodies participating in the motion; and each body contains a list of that
		body's coordinates. What do the contents of p represent?
		""",
		"",
		x -> occursin("time",lowercase(x))
	),
	Activity(
		"""
		Software development is always driven by the needs of its CLIENT program - in our case,
		this is the function NBodies.demo(). As you work with demo(), think about what kinds of
		behaviour you as a user demand of the type NBody, and include that functionality in your
		NBodies module design. Your use-cases and newly developed design will form an important
		part of your deliverables for this project.

		Notice that the two bodies (blue and red) do not orbit around each other. Why not? It may
		help your thinking to consider the following question: Is energy conserved in this system?
		""",
		"Think about the bodies' positions and momenta during the simulation",
		x -> occursin("no",lowercase(x))
	),
	Activity(
		"""
		Be aware that I have not given much attention to the naming of variables - you will need to
		comment, reformat and document the module much better than I have! Doing this will help you
		to understand the matrix programming techniques I have used, but they are not easy! You
		will need to thoroughly analyse the code and use Help and internet search to look up
		EVERYTHING you do not yet understand. Here is how to do this ...

		First remember: Matrices work very naturally in Julia, so if a=[1 2;3 5] and b=[2 3;4 5],
		then a*b and a.*b will deliver very different results. First test this idea at the Julia
		prompt, and then tell me the value of a*b-a.*b :
		""",
		"",
		x -> x == [8 7;14 9]
	),
	Activity(
		"""
		The secret of matrix programming is this: Whenever you think you need a for-loop to
		manipulate some values, build them together instead into a matrix that does the job for you.
		For example define and test an anonymous version of the factorial function that requires no
		recursion or iteration, but uses the prod() function.
		""",
		"",
		x -> occursin("->prod(1:",replace(x," "=>""))
	),
	Activity(
		"""
		Analyse the internal helper function relpos(), which calculates the relative positions
		x[i]-x[j] between all pairs of particles. Instead of looping over the particles, I compute
		the relative positions more efficiently in a matrix operation. Find where relpos() is
		called, and notice that its argument is a Vector - for example [1,2].

		At the Julia prompt, define relpos() and then call relpos([1,2]). What does the result look
		like? What does relpos() do? At the Julia prompt, analyse the calls

			locations = [1,2]
			locPerBody = repeat(locations,1,length(locations))
			permutedims(locPerBody)

		Discuss your findings with others: How, exactly, does relPos() work?

		Now tell me the result of repeat([1,2],2,3)?
		""",
		"",
		x -> x==[1 1 1;2 2 2;1 1 1;2 2 2]
	),
	Activity(
		"""
		Set a breakpoint and use the Debugger to investigate how simulate() uses relpos() to
		calculate the gravitational forces acting between the different bodies in the system.
		Which variable is used to calculate Newton's inverse-square law that falls with
		increased distance between the sources?
		""",
		"",
		x -> x=="invSq"
	),
	Activity(
		"""
		OK, now we can start on your actual project. First, extend my simulator to 2 bodies in
		2 dimensions and get the two bodies to orbit each other. This is your first major
		success - well done!

		Remember to keep your solution very neat, so that others can understand what you have
		done by reading your code.
		""",
		"Just move on when you have completed this activity",
		x -> true
	),
	Activity(
		"""
		Notice that your two bodies spiral outwards on their orbits because of my inaccurate
		simulator. Improve my integration method to achieve (practically) closed orbits.

		What is the name of the integration technique you are using?
		""",
		"",
		x -> occursin("runge-kutta",lowercase(x))
	),
	Activity(
		"""
		Now introduce a third body into your system (use the same mass for all three bodies). Try
		to create a situation in which the three bodies circle around each other without any of
		them getting thrown out of the system to infinity. Ensure that your system is conserving
		energy; that is, the three bodies taken together should neither gain nor lose energy!
		""",
		"Just move on when you have completed this activity",
		x -> true
	),
	Activity(
		"""
		You will need to find a way of testing your simulation for chaotic motion. There are
		several ways to do this - some are easy, some more difficult; some are approximate, some
		are very accurate. Here is a paper containing an accurate method of detecting chaos:

		https://www.nature.com/articles/s42003-019-0715-9.pdf

		Alternatively, in chaotic dynamics, small differences in initial conditions grow
		exponentially over time in the short term. The easiest way to test for exponential growth
		is by plotting a logarithm graph over time to see if it is linear. Decide how you want to
		measure chaoticity in your system.
		""",
		"",
		x -> true
	),
	Activity(
		"""
		Your project deliverable is a Julia module that must fulfil the following constraints:
			- the demonstration program MUST execute correctly within its own folder;
			- simulates the motion of three or more gravitational bodies;
			- executes an aesthetically pleasing chaotic trajectory;
			- demonstrates that the motion conserves energy;
			- demonstrates that the motion is indeed chaotic;
			- displays explanatory text so that anyone can understand it WITHOUT studying the code!
			- the code must be neat, well documented and easily readable and comprehensible.
		""",
		"Good luck!!!",
		x -> true
	),
]