text:
using Random

@testset "DataLoader" begin
    X = reshape([1:10;], (2, 5))
    Y = [1:5;]

    d = DataLoader(X, batchsize=2)
    @inferred first(d)
    batches = collect(d)
    @test eltype(batches) == eltype(d) == typeof(X)
    @test length(batches) == 3
    @test batches[1] == X[:,1:2]
    @test batches[2] == X[:,3:4]
    @test batches[3] == X[:,5:5]

    d = DataLoader(X, batchsize=2, partial=false)
    @inferred first(d)
    batches = collect(d)
    @test eltype(batches) == eltype(d) == typeof(X)
    @test length(batches) == 2
    @test batches[1] == X[:,1:2]
    @test batches[2] == X[:,3:4]

    d = DataLoader((X,), batchsize=2, partial=false)
    @inferred first(d)
    batches = collect(d)
    @test eltype(batches) == eltype(d) == Tuple{typeof(X)}
    @test length(batches) == 2
    @test batches[1] == (X[:,1:2],)
    @test batches[2] == (X[:,3:4],)

    d = DataLoader((X, Y), batchsize=2)
    @inferred first(d)
    batches = collect(d)
    @test eltype(batches) == eltype(d) == Tuple{typeof(X), typeof(Y)}
    @test length(batches) == 3
    @test length(batches[1]) == 2
    @test length(batches[2]) == 2
    @test length(batches[3]) == 2
    @test batches[1][1] == X[:,1:2]
    @test batches[1][2] == Y[1:2]
    @test batches[2][1] == X[:,3:4]
    @test batches[2][2] == Y[3:4]
    @test batches[3][1] == X[:,5:5]
    @test batches[3][2] == Y[5:5]

    # test with NamedTuple
    d = DataLoader((x=X, y=Y), batchsize=2)
    @inferred first(d)
    batches = collect(d)
    @test eltype(batches) == eltype(d) == NamedTuple{(:x, :y), Tuple{typeof(X), typeof(Y)}}
    @test length(batches) == 3
    @test length(batches[1]) == 2
    @test length(batches[2]) == 2
    @test length(batches[3]) == 2
    @test batches[1][1] == batches[1].x == X[:,1:2]
    @test batches[1][2] == batches[1].y == Y[1:2]
    @test batches[2][1] == batches[2].x == X[:,3:4]
    @test batches[2][2] == batches[2].y == Y[3:4]
    @test batches[3][1] == batches[3].x == X[:,5:5]
    @test batches[3][2] == batches[3].y == Y[5:5]

    # test interaction with `train!`
    θ = ones(2)
    X = zeros(2, 10)
    loss(x) = sum((x .- θ).^2)
    d  = DataLoader(X)
    Flux.train!(loss, [θ], ncycle(d, 10), Descent(0.1))
    @test norm(θ) < 1e-4

    # test interaction with `train!`
    θ = zeros(2)
    X = ones(2, 10)
    Y = fill(2, 10)
    loss(x, y) = sum((y - x'*θ).^2)
    d  = DataLoader((X, Y))
    Flux.train!(loss, [θ], ncycle(d, 10), Descent(0.1))
    @test norm(θ .- 1) < 1e-10

    # specify the rng
    d = map(identity, DataLoader(X, batchsize=2; shuffle=true, rng=Random.seed!(Random.default_rng(), 5)))
end


text:
export modnewton

"""
modnewton(f,df,H,x)

Newton method with Hessian shift for solving min_x f(x)

"""
function modnewton(f::Function,df::Function,H::Function,x::Vector;maxIter=20,atol=1e-8,out::Int=0,
	storeInterm::Bool=false,beta=1e-1,
		lineSearch::Function=(f,df,fk,dfk,xk,pk)->armijo(f,fk,dfk,xk,pk,maxIter=10))

    his = zeros(maxIter,4)
    X = (storeInterm) ? zeros(length(x),maxIter) : []
    i = 1; flag = -1; LL = []
    while i<=maxIter
        fk = f(x)
        dfk = df(x)
        his[i,1:2] = [norm(fk) norm(dfk)]
        if storeInterm; X[:,i] = x; end;

        if(norm(dfk)<atol)
            flag = 0
            his = his[1:i,:]
            break
        end

        # get search direction
        d2f = H(x)
        tau = (minimum(diag(d2f)) > 0) ? 0 :  -minimum(diag(d2f))+beta
        nfac = 0;
        for k=1:20
            try
                nfac +=1
                LL    = cholfact(d2f + tau*speye(length(x)))
                break;
            catch err
                if isa(err,Base.LinAlg.PosDefException) && k<20
                    tau = max(2*tau,beta)
                elseif isa(err,Base.LinAlg.PosDefException)
                    flag = -2
                    his[i,4] = nfac
                    his = his[1:i,:]
                    break; break;
                else
                    throw(err)
                end
            end
        end
        his[i,4] = nfac
        pk    = - (LL\dfk)

        # line search
		ak,his[i,3] = lineSearch(f,df,fk,dfk,x,pk)
        if his[i,3]==-1
            flag = -3
            his  = his[1:i,:]
            break
        end
        if out>0
            @printf "iter=%04d\t|f|=%1.2e\t|df|=%1.2e\tLS=%d\n" i his[i,1] his[i,2] his[i,3]
        end
        # update x and H
        x  += ak*pk

        i+=1
    end
    i = min(maxIter,i)

    if out>=0
        if flag==-1
            println(@sprintf("modnewton iterated maxIter (=%d) times but reached only atol of %1.2e instead of tol=%1.2e",i,his[i,2],atol))
        elseif flag==-2
            println(@sprintf("modnewton stopped because the Hessian at iteration %d was not positive definite.",i))
        elseif flag==-3
            println(@sprintf("modnewton stopped of a line search fail at iteration %d.",i))
		elseif out>1
            println(@sprintf("modnewton achieved desired atol of %1.2e at iteration %d.",atol,i))
        end
    end
    if storeInterm; X = X[:,1:i]; end
    return x,flag,his,X
end

ext: jl


text:
function get_first_visible_neighbours(layout, ix, jy)
    on_grid(x, y) = 1 <= x <= nr_rows && 1 <= y <= nr_cols

    nr_rows, nr_cols = size(layout)
    rays = [
        map(x -> layout[x, jy], (ix-1):-1:1),
        map(((x, y),) -> layout[x, y], filter(((x, y),) -> on_grid(x, y), map(dx -> (ix - dx, jy + dx), 1:nr_rows))),
        map(y -> layout[ix, y], (jy+1:nr_cols)),
        map(((x, y),) -> layout[x, y], filter(((x, y),) -> on_grid(x, y), map(dx -> (ix + dx, jy + dx), 1:nr_rows))),
        map(x -> layout[x, jy], (ix+1:nr_rows)),
        map(((x, y),) -> layout[x, y], filter(((x, y),) -> on_grid(x, y), map(dx -> (ix + dx, jy - dx), 1:nr_rows))),
        map(y -> layout[ix, y], (jy-1):-1:1),
        map(((x, y),) -> layout[x, y], filter(((x, y),) -> on_grid(x, y), map(dx -> (ix - dx, jy - dx), 1:nr_rows)))
    ]

    function find(itr)
        for item in itr
            if item == 'L' || item == '#'
                return item
            end
        end
        '.'
    end

    map(find, rays)
end

function single_step(layout)
    updated = Array{Char, 2}(undef, size(layout))
    nr_changed = 0
    for jy in 1:size(layout, 2)
        for ix in 1:size(layout, 1)
            item = layout[ix, jy]
            if item == '.'
                updated[ix, jy] = '.'
            else
                neighbours = get_first_visible_neighbours(layout, ix, jy)
                if item == 'L' && all(neighbour != '#' for neighbour in neighbours)
                    updated[ix, jy] = '#'
                    nr_changed += 1
                elseif item == '#' && count(n -> n == '#', neighbours) >= 5
                    updated[ix, jy] = 'L'
                    nr_changed += 1
                else
                    updated[ix, jy] = item
                end
            end
        end
    end
    updated, nr_changed
end

function seating_system(seat_layout)
    while true
        seat_layout, nr_changed_seats = single_step(seat_layout)
        if nr_changed_seats == 0
            return count(x -> x == '#', seat_layout)
        end
    end
end

function main()
    lines = readlines("../resources/input_11.txt")
    seat_layout = permutedims(hcat(map(collect, lines)...))
    result = seating_system(seat_layout)
    println(result)
end

main()


text:
function integrate_icp(ξ0, p::BezierCurve, ω, θ)
    exp_ω_θ = exp(ω * θ)
    return ξ0 * exp_ω_θ - ω * exp_ω_θ * exponential_integral(p, -ω, θ)
end

const COP_TRAJ_DEGREE = 4

struct ICPTrajectoryGenerator{T, N, O<:MOI.AbstractOptimizer, L}
    ω::T
    model::Model{T, O}
    cop_pieces::Vector{BezierCurve{COP_TRAJ_DEGREE + 1, SVector{2, Variable}}}
    icp_knots::Vector{SVector{2, Variable}}
    Δts::Vector{T}
    cop_polyhedra::Vector{MHRep{N, 2, T, L}}
    preferred_cops::Vector{SVector{2, T}}
    initial_cop::Base.RefValue{SVector{2, T}}
    initial_icp::Base.RefValue{SVector{2, T}}
    final_icp::Base.RefValue{SVector{2, T}}
    num_active_segments::Base.RefValue{Int}

    function ICPTrajectoryGenerator{T, N}(optimizer::O, num_segments::Integer, ω::Number) where {T, N, O<:MOI.AbstractOptimizer}
        # TODO:
        # * find the right objective function

        # Variables: (2 * COP_TRAJ_DEGREE + 4) * n + 2:
        # * CoP control points: 2 * (COP_TRAJ_DEGREE + 1) * n
        # * ICP knots: 2 * (n + 1)
        # Equality constraints: 6 * n + 4:
        # * initial CoP: 2
        # * CoP C⁰ continuity: 2 * (n - 1)
        # * CoP C¹ continuity: 2 * (n - 1)
        # * ICP dynamics: 2 * n
        # * initial ICP: 2
        # * final ICP: 2
        # * final CoP == final ICP: 2

        n = num_segments
        L = N * 2
        model = Model(optimizer)

        # Parameters
        Δts = Parameter(model, val=zeros(T, n))
        cop_polyhedra = Parameter(model, val=[zero(MHRep{N, 2, T, L}) for i = 1 : n])
        preferred_cops = Parameter(model, val=[zero(SVector{2, T}) for i = 1 : n])
        initial_cop = Parameter(model, val=zero(SVector{2, T}))
        initial_icp = Parameter(model, val=zero(SVector{2, T}))
        final_icp = Parameter(model, val=zero(SVector{2, T}))

        # CoP trajectory
        cop_pieces = [BezierCurve(ntuple(i -> SVector(Variable(model), Variable(model)), Val(COP_TRAJ_DEGREE + 1))) for i = 1 : n]

        # Initial CoP
        @constraint model first(first(cop_pieces).coeffs) == initial_cop

        # CoP trajectory continuity
        for i = 1 : n - 1
            # C⁰
            pᵢ = cop_pieces[i]
            pᵢ₊₁ = cop_pieces[i + 1]
            @constraint model first(pᵢ₊₁.coeffs) == last(pᵢ.coeffs)

            # C¹. Note: clearing denominators.
            pᵢ′ = derivative(pᵢ)
            pᵢ₊₁′ = derivative(pᵢ₊₁)
            # @constraint model first(pᵢ₊₁′.coeffs) * Δts[i] == last(pᵢ′.coeffs) * Δts[i + 1] #  # TODO: allocates. Fix missing optimization for SVector * scalar.
            @constraint model Vector(first(pᵢ₊₁′.coeffs)) * Δts[i] == Vector(last(pᵢ′.coeffs)) * Δts[i + 1]
        end

        # CoP bounds
        for i = 1 : n
            Aᵢ = @expression cop_polyhedra[i].A
            bᵢ = @expression cop_polyhedra[i].b
            pᵢ = cop_pieces[i]
            for point in pᵢ.coeffs
                @constraint model Aᵢ * point <= bᵢ
            end
        end

        # ICP dynamics
        icp_knots = [SVector(Variable(model), Variable(model)) for i = 1 : n + 1]
        for i = 1 : n
            # Use an indicator variable to make icp knots equal if the Δt is zero
            z = @expression Δts[i] != 0
            next_icp = icp_knots[i] * exp(ω) - ω * exp(ω) * exponential_integral(cop_pieces[i] , -ω)
            @constraint model icp_knots[i + 1] == z * Vector(next_icp) + (1 - z) * Vector(icp_knots[i]) # TODO: don't use Vector
        end
        @constraint model first(icp_knots) == initial_icp
        @constraint model last(icp_knots) == final_icp
        @constraint model last(last(cop_pieces).coeffs) == final_icp # implies that final ICP velocity is zero

        # Objective
        objective = Parametron.LazyExpression(identity, zero(QuadraticFunction{Float64}))
        for i = 1 : n
            pᵢ = cop_pieces[i]
            z = @expression Δts[i] != 0 # inactive segments shouldn't incur cost
            for point in pᵢ.coeffs
                e = @expression point - preferred_cops[i]
                objective = @expression objective + z * (e ⋅ e)
            end
        end
        @objective model Minimize objective

        num_active_segments = Ref(0)

        new{T, N, O, L}(
            ω,
            model,
            cop_pieces, icp_knots,
            Δts.val[], cop_polyhedra.val[], preferred_cops.val[],
            initial_cop.val, initial_icp.val, final_icp.val,
            num_active_segments
        )
    end
end

function Base.empty!(generator::ICPTrajectoryGenerator)
    for hrep in generator.cop_polyhedra
        hrep.A .= 0
        hrep.b .= 0
    end
    generator.Δts .= 0
    generator.num_active_segments[] = 0
    generator
end

# TODO: name?
function initialize!(generator::ICPTrajectoryGenerator, initial_cop::AbstractVector, initial_icp::AbstractVector, final_icp::AbstractVector)
    generator.initial_cop[] = initial_cop
    generator.initial_icp[] = initial_icp
    generator.final_icp[] = final_icp
    generator
end

function push_segment!(generator::ICPTrajectoryGenerator,
        Δt::Number, cop_polyhedron::ConvexHull, preferred_cop::SVector = centroid(cop_polyhedron))
    i = generator.num_active_segments[] += 1
    @boundscheck i <= length(generator.Δts) || error()
    @inbounds begin
        generator.Δts[i] = Δt
        hrep = generator.cop_polyhedra[i]
        hrep!(hrep.A, hrep.b, cop_polyhedron)
        generator.preferred_cops[i] = preferred_cop
    end
    generator
end

function solve!(generator::ICPTrajectoryGenerator)
    Parametron.solve!(generator.model)
    checkstatus(generator.model)
end

function checkstatus(model::Parametron.Model)
    ok = terminationstatus(model) == MOI.OPTIMAL && primalstatus(model) == MOI.FEASIBLE_POINT
    if !ok
        okish = terminationstatus(model) == MOI.ALMOST_OPTIMAL && primalstatus(model) == MOI.UNKNOWN_RESULT_STATUS
        if !okish
            """
            Solve failed.
            * termination status: $(terminationstatus(model))
            * primal status: $(primalstatus(model))
            * dual status: $(dualstatus(model))
            """ |> error
        end
    end
end

function initial_icp(generator::ICPTrajectoryGenerator, segment_number::Integer)
    Parametron.value.(generator.model, generator.icp_knots[segment_number])
end

function cop_piece(generator::ICPTrajectoryGenerator, segment_number::Integer)
    BezierCurve(map(x -> Parametron.value.(generator.model, x), generator.cop_pieces[segment_number].coeffs))
end

final_time(generator::ICPTrajectoryGenerator) = sum(generator.Δts)

function find_segment(generator::ICPTrajectoryGenerator{T}, t::Number) where T
    i = 1
    t0 = zero(T)
    for Δt in generator.Δts
        if t <= t0 + Δt || i >= generator.num_active_segments[]
            break
        else
            i += 1
            t0 += Δt
        end
    end
    i, t0
end

num_segments(generator::ICPTrajectoryGenerator) = length(generator.Δts)

function (generator::ICPTrajectoryGenerator{T})(t::Number) where T
    t < 0 && throw(ArgumentError("Cannot evaluate at a negative time."))
    if t < sum(generator.Δts)
        # Find segment index and segment start time
        i, t0 = find_segment(generator, t)

        # Segment phase
        θ = (t - t0) / generator.Δts[i]

        # CoP
        pᵢ = cop_piece(generator, i)
        p = pᵢ(θ)

        # ICP
        ω = generator.ω
        ξ0 = initial_icp(generator, i)
        ξ = integrate_icp(ξ0, pᵢ, ω, θ)

        # ICP derivative. Use (12) in Capturability-based Analysis Part 1.
        ξd = ω * (ξ - p)
    else
        # Hold at the final ICP
        ξ = generator.final_icp[]
        ξd = zero(ξ)
    end
    ξ, ξd
end

function cop(generator::ICPTrajectoryGenerator, t::Number)
    ξ, ξd = generator(t)
    ω = generator.ω
    return ξ - ξd / ω
end


text:
using CLArrays, CLArrays.Shorthands
using GPUArrays.TestSuite, Base.Test

for dev in CLArrays.devices()
    # we only test gpu devices for now - cpu drivers use different
    # inconsistent struct alignment, which has low priority to be fixed right now
    is_gpu(dev) || continue
    CLArrays.init(dev)
    @testset "Device: $dev" begin

        TestSuite.run_tests(CLArray)

        @testset "muladd & abs" begin
            a = rand(Float32, 32) - 0.5f0
            A = CLArray(a)
            x = abs.(A)
            @test Array(x) == abs.(a)
            y = muladd.(A, 2f0, x)
            @test Array(y) == muladd(a, 2f0, abs.(a))
            u0 = CLArray(ones(Float32, 32, 32)); u1 = zeros(u0);
            @test muladd.(2, Array(u0), Array(u1)) == Array(muladd.(2, u0, u1))
            ###########
            # issue #20
            against_base(a-> abs.(a), CLArray{Float32}, (10, 10))

            #### bools in kernel:
        end
        @testset "bools" begin
            A, B = rand(Bool, 10), rand(Bool, 10)
            Ag, Bg = CLArray(A), CLArray(B)
            res = A .& B
            resg = Ag .& Bg
            @test res == Array(resg)
            # this version needs to have a fix in GPUArrays, since it uses T.(array)
            # in copy to convert to array type, but that actually convert Array{Bool} to BitArray
            # against_base((a, b)-> a .& b, CLArray{Bool}, (10,), (10,))
        end

        @testset "Shorthand Test" begin
            GPUArrays.allowslow(true)
            @test collect(cl([1,2])) == [1,2]
            @test collect(cl([1 2;3 4])) == [1 2;3 4]
            @test cl([1,2,3]) == CLArray([1,2,3])
        end

        @testset "Compilation from String" begin
            copy_source = """
                __kernel void copy(
                    __global float *dest,
                    __global float *source
                ){
                int gid = get_global_id(0);
                dest[gid] = source[gid];
                }
            """
            source = rand(CLArray{Float32}, 1023, 11)
            dest = zeros(CLArray{Float32}, size(source))
            gpu_call(:copy => copy_source, dest, (dest, source))
            Array(dest) == Array(source)
        end
    end
end

# Indexing with
# Issue
#   ([CartesianIndex(2,2), CartesianIndex(2,1)],) # Array{CartesianIndex} # FAIL

# #The above is equal to:
# Typ = CLArray
# GPUArrays.allowslow(false)
# TestSuite.run_gpuinterface(Typ)
# TestSuite.run_base(Typ)
# TestSuite.run_blas(Typ)
# TestSuite.run_broadcasting(CLArray)
# TestSuite.run_construction(Typ)
# TestSuite.run_fft(Typ)
# TestSuite.run_linalg(Typ)
# TestSuite.run_mapreduce(Typ)
# TestSuite.run_indexing(Typ)
#
# function test_sizes6(state, out)
#     x1 = (1, 2, 3)
#     out[1] = sizeof(x1)
#
#     x2 = (1f0, 2f0, 3f0)
#     out[2] = sizeof(x2)
#
#     x3 = ((1f0, 2f0), 1.0)
#     out[3] = sizeof(x3)
#
#     x4 = ((1f0, 2f0, 3f0), 1.0)
#     out[4] = sizeof(x4)
#
#     x5 = ((1f0, 2f0, 3f0, 4f0), 1.0)
#     out[5] = sizeof(x5)
#
#     x6 = (1.0, (1f0, 2f0))
#     out[6] = sizeof(x6)
#
#     x7 = (1.0, (1f0, 2f0, 3f0))
#     out[7] = sizeof(x7)
#
#     x8 = (1.0, (1f0, 2f0, 3f0, 4f0))
#     out[8] = sizeof(x8)
#
#     x9 = ((1f0, 2f0), UInt32(1), 1.0)
#     out[9] = sizeof(x9)
#
#     return
# end
#
# function test_sizes6(state, out)
#     x1 = 1
#     out[1] = sizeof(x1)
#
#     x2 = (1, 2)
#     out[2] = sizeof(x2)
#
#     x3 = (1, 2, 3)
#     out[3] = sizeof(x3)
#
#     x4 = 1f0
#     out[4] = sizeof(x4)
#
#     x5 = (1f0, 2f0)
#     out[5] = sizeof(x5)
#
#     x6 = (1f0, 2f0, 3f0)
#     out[6] = sizeof(x6)
#
#     x7 = (1f0, 2f0, 3f0, 4f0, 5f0)
#     out[7] = sizeof(x7)
#
#     x8 =  (1f0, 2f0, 3f0, 4f0, 5f0, 6f0)
#     out[8] = sizeof(x8)
#
#     x9 = (1, 2, 3, 4, 5, 6)
#     out[9] = sizeof(x9)
#
#     x10 = UInt8(1)
#     out[10] = sizeof(x10)
#
#     x11 = (x10, x10)
#     out[11] = sizeof(x11)
#
#     x12 = (x10, x10, x10)
#     out[12] = sizeof(x12)
#
#
#     x13 = CLArrays.DeviceArray{Complex{Float64},3,CLArrays.HostPtr{Complex{Float64}}}(
#         CLArrays.HostPtr{Complex{Float64}}(),
#         (Cuint(2), Cuint(2), Cuint(2))
#     )
#     out[13] = sizeof(x13)
#     #
#     # x14 = (x13, x13)
#     # out[14] = sizeof(x14)
#     #
#     # x15 = (x13, x13, x13)
#     # out[15] = sizeof(x15)
#     return
# end


text:
#
# The whole submodule thing is weird, so ...
#

# TODO sub2ind for assignments
# TODO broadcast fallback for empty arrays
# TODO broadcast_reduce tag-team

include("factors_main.jl") # first because aux needs it
include("errors.jl")
include("auxiliary.jl")
include("factors_dims.jl")
include("factors_access.jl")
include("factors_dataframes.jl")
include("factors_io.jl")


text:
function load(::Day{5}, path)
    # Shifting indices by 1 for Julia's 1-based indexing
    return [
        [(Tuple(parse.(Int, split(t, ","))) .+ 1) for t in split(l, " -> ")] for
        l in readlines(path)
    ]
end

# Part 1
getx(t::NTuple{2,<:Integer}) = first(t)
gety(t::NTuple{2,<:Integer}) = last(t)
xmax(data) = maximum(maximum(getx.(line) for line in data))
ymax(data) = maximum(maximum(gety.(line) for line in data))
not_diagonal(line) = ==(getx.(line)...) || ==(gety.(line)...)
_range(a, b) = a:sign(b - a):b

function get_indices(line)
    line = sort(line)
    if not_diagonal(line)
        return CartesianIndex(line[1]):CartesianIndex(line[2])
    end
    return CartesianIndex.(_range(getx.(line)...), _range(gety.(line)...))
end

function solve1(::Day{5}, data) # count count overlaps
    counts = zeros(Int, xmax(data), ymax(data))
    @inbounds @simd for line in data
        if not_diagonal(line)
            counts[get_indices(line)] .+= 1
        end
    end
    return sum(counts .> 1)
end

# Part 2 – we just need to remove the `not_diagonal` check
function solve2(::Day{5}, data)
    counts = zeros(Int, xmax(data), ymax(data))
    @inbounds @simd for line in data
        counts[get_indices(line)] .+= 1
    end
    return sum(counts .> 1)
end

testresult1(::Day{5}) = 5
testresult2(::Day{5}) = 12


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

## FEM benchmark from https://github.com/JuliaLang/julia/issues/9668

# assemble the finite-difference laplacian
function fdlaplacian(N)
    # create a 1D laplacian and a sparse identity
    fdl1 = sparse(SymTridiagonal(fill(-2., N), fill(1., N-1)))
    # laplace operator on the full grid
    return kron(sparse(1.0I, N, N), fdl1) + kron(fdl1, sparse(1.0I, N, N))
end

# get the list of boundary dof-indices
function get_free(N)
    L = zeros(Int, N, N)
    L[2:N-1, 2:N-1] = 1
    return findall(L)
end

# timing of assembly, slice and solve
function run_fem(N)
    Ifree = get_free(N)
    # assembly
    A = fdlaplacian(N)
    # boundary condition
    B = A[Ifree, Ifree]
    # solver
    lufact(B)
end

function fem_perf()
    # run tests once to compile
    run_fem(10)

    # run
    @timeit run_fem(256) "finite_elements" ""

    # runs the tests
    #NN = 2.^(3:8)
    #TT = zeros(3, length(NN))
    #@printf("(All times are seconds)\n")
    #@printf("     N   |   assembly |   slice   |   lufact  | slice / N^4 \n")
    #@printf("---------|------------|-----------|-----------|-------------\n")
    #for n = 1:length(NN)
    #    t1, t2, t3 = timings(NN[n])
    #    @printf(" %4.1e |  %4.2e  | %4.2e  |  %4.2e |   %4.2e \n",
    #            NN[n], t1, t2, t3, t2 / NN[n]^4)
    #end
end

max_stars_repo_name: blegat/julia

max_issues_repo_name: blegat/julia

max_forks_repo_name: blegat/julia


text:
# Gating variables
abstract type AbstractGatingVariable end

output(x::AbstractGatingVariable) = getfield(x, :output)
timeconstant(x::AbstractGatingVariable) = getfield(x, :tau)
steadystate(x::AbstractGatingVariable) = getfield(x, :steadystate)
forward_rate(x::AbstractGatingVariable) = getfield(x, :alpha)
reverse_rate(x::AbstractGatingVariable) = getfield(x, :beta)
hasexponent(x::AbstractGatingVariable) = hasfield(typeof(x), :p) ? getfield(x, :p) !== one(typeof(x.p)) : false
exponent(x::AbstractGatingVariable) = hasexponent(x) ? getfield(x, :p) : nothing

@enum GateVarType NullType Activation Inactivation
@enum GateVarForm AlphaBeta SteadyStateTau

struct GatingVariable <: AbstractGatingVariable
    output::Num
    alpha::Num
    beta::Num
    steadystate::Num
    tau::Num
    p::Real
    type::GateVarType # for future use
end

const Gate = GatingVariable

function GatingVariable(T::GateVarForm, x::Num, y::Num, p::Real = 1; name = Base.gensym("GateVar"))
    if T == AlphaBeta
        alpha, beta = x, y
        ss = alpha/(alpha + beta)
        tau = inv(alpha + beta)
    elseif T == SteadyStateTau
        ss, tau = x, y
        alpha = ss/tau
        beta = inv(tau) - alpha
    end

    out = only(@variables $name(t) = ss)
    GatingVariable(out, alpha, beta, ss, tau, p, NullType)
end

macro gate(ex::Expr, p::Expr = :(p=1))
    name = Base.gensym("GateVar")
    _make_gate_variable(name, ex, p)
end

macro gate(name::Symbol, ex::Expr, p::Expr = :(p=1))
    _make_gate_variable(name, ex, p)
end

function _make_gate_variable(name::Symbol, ex::Expr, p::Expr = :(p=1))
    ex = MacroTools.striplines(ex)
    length(ex.args) !== 2 && throw("Invalid number of input equations.")
    p.args[1] != :p && throw("Please use `p` to define a gate exponent.")
    gate_pow = p.args[2]

    rate_vars = map(eq -> eq.args[1], ex.args)

    #= Find paramters in rate equations =#
    param_list = Symbol[]
    foreach(rate->MacroTools.postwalk(x -> extract_symbols(x, param_list), rate.args[2]), ex.args)
    filter!(x -> x != :Vₘ, param_list)

    if issetequal([:α, :β], rate_vars)  # Equation order doesn't matter
        rate_type = AlphaBeta
        gate_expr = :(GatingVariable($AlphaBeta, α, β, $gate_pow; name=$(QuoteNode(name))))
    elseif issetequal([:ss, :τ], rate_vars)  # Equation order doesn't matter
        rate_type = SteadyStateTau
        gate_expr = :(GatingVariable($SteadyStateTau, ss, τ, $p.args[2]; name=$(QuoteNode(name))))
    else
        throw("invalid keywords")
    end
    pushfirst!(ex.args, :(Vₘ = MembranePotential()))
    foreach(param -> pushfirst!(ex.args, :(@parameters $param)), param_list)  # generate parameters
    push!(ex.args, gate_expr)
    ex
end

# Likely replace this with macro
# function GatingVariable(T::GateVarForm; p = one(Int64), kwargs...)
#     syms = keys(kwargs)
#     length(syms) !== 2 && throw("Invalid number of input equations.")
#     if T == SteadyStateTau
#         if issetequal([:m∞, :τₘ], syms)
#             return Gate(SteadyStateTau, kwargs[:m∞], kwargs[:τₘ], p; name = :m)
#         elseif issetequal([:h∞, :τₕ], syms)
#             return Gate(SteadyStateTau, kwargs[:h∞], kwargs[:τₕ], p; name = :h)
#         elseif isesetequal([:n∞, :τₙ], syms)
#             return Gate(SteadyStateTau, kwargs[:n∞], kwargs[:τₙ], p; name = :n)
#         else
#             throw("invalid keywords")
#         end
#     elseif T == AlphaBeta
#         if issetequal([:αₘ, :βₘ], syms)
#             return Gate(AlphaBeta, kwargs[:αₘ], kwargs[:βₘ], p; name = :m)
#         elseif issetequal([:αₕ, :βₕ], syms)
#             return Gate(AlphaBeta, kwargs[:αₕ], kwargs[:βₕ], p; name = :h)
#         elseif issetequal([:αₙ, :βₙ], syms)
#             return Gate(AlphaBeta, kwargs[:αₙ], kwargs[:βₙ], p; name = :n)
#         else
#             throw("invalid keywords")
#         end
#     end
# end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg
using Base.BinaryPlatforms
const YGGDRASIL_DIR = "../.."
include(joinpath(YGGDRASIL_DIR, "platforms", "mpi.jl"))

name = "MPItrampoline"
version = v"4.0.2"

mpich_version_str = "4.0.2"
mpiwrapper_version = v"2.8.1"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/eschnett/MPItrampoline/archive/refs/tags/v$(version).tar.gz",
                  "89abda0526dba9e52a3b6334d1ac86709c12567ff114acd610471e66c6190b89"),
    ArchiveSource("https://www.mpich.org/static/downloads/$(mpich_version_str)/mpich-$(mpich_version_str).tar.gz",
                  "5a42f1a889d4a2d996c26e48cbf9c595cbf4316c6814f7c181e3320d21dedd42"),
    ArchiveSource("https://github.com/eschnett/MPIwrapper/archive/refs/tags/v$(mpiwrapper_version).tar.gz",
                  "e6fc1c08ad778675e5b58b91b4658b12e3f985c6d4c5c2c3e9ed35986146780e"),
]

# Bash recipe for building across all platforms
script = raw"""
################################################################################
# MPItrampoline
################################################################################

cd $WORKSPACE/srcdir/MPItrampoline-*
mkdir build
cd build
cmake \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_FIND_ROOT_PATH=$prefix \
    -DCMAKE_INSTALL_PREFIX=$prefix \
    -DBUILD_SHARED_LIBS=ON \
    -DMPITRAMPOLINE_DEFAULT_LIB="@MPITRAMPOLINE_DIR@/lib/libmpiwrapper.so" \
    -DMPITRAMPOLINE_DEFAULT_MPIEXEC="@MPITRAMPOLINE_DIR@/bin/mpiwrapperexec" \
    ..
cmake --build . --config RelWithDebInfo --parallel $nproc
cmake --build . --config RelWithDebInfo --parallel $nproc --target install

################################################################################
# Install MPICH
################################################################################

cd ${WORKSPACE}/srcdir/mpich*

EXTRA_FLAGS=()
if [[ "${target}" != i686-linux-gnu ]] || [[ "${target}" != x86_64-linux-* ]]; then
    # Define some obscure undocumented variables needed for cross compilation of
    # the Fortran bindings.  See for example
    # * https://stackoverflow.com/q/56759636/2442087
    # * https://github.com/pmodels/mpich/blob/d10400d7a8238dc3c8464184238202ecacfb53c7/doc/installguide/cfile
    export CROSS_F77_SIZEOF_INTEGER=4
    export CROSS_F77_SIZEOF_REAL=4
    export CROSS_F77_SIZEOF_DOUBLE_PRECISION=8
    export CROSS_F77_FALSE_VALUE=0
    export CROSS_F77_TRUE_VALUE=1

    if [[ ${nbits} == 32 ]]; then
        export CROSS_F90_ADDRESS_KIND=4
    else
        export CROSS_F90_ADDRESS_KIND=8
    fi
    export CROSS_F90_OFFSET_KIND=8
    export CROSS_F90_INTEGER_KIND=4
    export CROSS_F90_INTEGER_MODEL=9
    export CROSS_F90_REAL_MODEL=6,37
    export CROSS_F90_DOUBLE_MODEL=15,307
    export CROSS_F90_ALL_INTEGER_MODELS=2,1,4,2,9,4,18,8,
    export CROSS_F90_INTEGER_MODEL_MAP={2,1,1},{4,2,2},{9,4,4},{18,8,8},

    if [[ "${target}" == i686-linux-musl ]]; then
        # Our `i686-linux-musl` platform is a bit rotten: it can run C programs,
        # but not C++ or Fortran.  `configure` runs a C program to determine
        # whether it's cross-compiling or not, but when it comes to running
        # Fortran programs, it fails.  In addition, `configure` ignores the
        # above exported variables if it believes it's doing a native build.
        # Small hack: edit `configure` script to force `cross_compiling` to be
        # always "yes".
        sed -i 's/cross_compiling=no/cross_compiling=yes/g' configure
        EXTRA_FLAGS+=(ac_cv_sizeof_bool="1")
    fi
fi

# Building with an external hwloc leads to problems loading the
# resulting libraries and executable via MPIwrapper, because this
# happens outside of Julia's control.

# Produce static libraries since these are easier to handle
# downstream. (MPIwrapper doesn't need to find its dependencies.) Also
# make sure that these static libraries can be linked into shared
# libraries.

# Disable OpenCL since it requires the option `-framework OpenCL` on
# macOS, and cmake does not handle this option well. cmake assumes
# that this is one option and one file name and splits them up.

export CFLAGS='-fPIC -DPIC'
export CXXFLAGS='-fPIC -DPIC'
export FFLAGS='-fPIC -DPIC'
export FCFLAGS='-fPIC -DPIC'

if [[ "${target}" == aarch64-apple-* ]]; then
    export FFLAGS="$FFLAGS -fallow-argument-mismatch"
fi

if [[ "${target}" == *-apple-* ]]; then
    # MPICH uses the link options `-flat_namespace` on Darwin. This
    # conflicts with MPItrampoline, which requires the option
    # `-twolevel_namespace`.
    EXTRA_FLAGS+=(--enable-two-level-namespace)
fi

if [[ "${target}" == aarch64-apple-* ]]; then
    EXTRA_FLAGS+=(
        FFLAGS=-fallow-argument-mismatch
        FCFLAGS=-fallow-argument-mismatch
    )
fi

./configure \
    --build=${MACHTYPE} \
    --host=${target} \
    --disable-dependency-tracking \
    --docdir=/tmp \
    --enable-shared=no \
    --enable-static=yes \
    --enable-threads=multiple \
    --enable-opencl=no \
    --with-device=ch3 \
    --prefix=${prefix}/lib/mpich \
    "${EXTRA_FLAGS[@]}"

# Remove empty `-l` flags from libtool
# (Why are they there? They should not be.)
# Run the command several times to handle multiple (overlapping) occurrences.
sed -i 's/"-l /"/g;s/ -l / /g;s/-l"/"/g' libtool
sed -i 's/"-l /"/g;s/ -l / /g;s/-l"/"/g' libtool
sed -i 's/"-l /"/g;s/ -l / /g;s/-l"/"/g' libtool

make -j${nproc}
make -j${nproc} install

# Delete duplicate file
if ar t $prefix/lib/mpich/lib/libpmpi.a | grep -q setbotf.o; then
    ar d $prefix/lib/mpich/lib/libmpifort.a setbotf.o
fi

################################################################################
# Install MPIwrapper
################################################################################

cd $WORKSPACE/srcdir/MPIwrapper-*
mkdir build
cd build
# Yes, this is tedious. No, without being this explicit, cmake will
# not properly auto-detect the MPI libraries on Darwin.
if [[ "${target}" == *-apple-* ]]; then
    ext='a'
    cmake \
        -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
        -DCMAKE_FIND_ROOT_PATH="${prefix}/lib/mpich;${prefix}" \
        -DCMAKE_INSTALL_PREFIX=${prefix} \
        -DBUILD_SHARED_LIBS=ON \
        -DMPI_C_COMPILER=cc \
        -DMPI_CXX_COMPILER=c++ \
        -DMPI_Fortran_COMPILER=gfortran \
        -DMPI_C_LIB_NAMES='mpi;pmpi' \
        -DMPI_CXX_LIB_NAMES='mpicxx;mpi;pmpi' \
        -DMPI_Fortran_LIB_NAMES='mpifort;mpi;pmpi' \
        -DMPI_pmpi_LIBRARY=${prefix}/lib/mpich/lib/libpmpi.${ext} \
        -DMPI_mpi_LIBRARY=${prefix}/lib/mpich/lib/libmpi.${ext} \
        -DMPI_mpicxx_LIBRARY=${prefix}/lib/mpich/lib/libmpicxx.${ext} \
        -DMPI_mpifort_LIBRARY=${prefix}/lib/mpich/lib/libmpifort.${ext} \
        -DMPIEXEC_EXECUTABLE=${prefix}/lib/mpich/bin/mpiexec \
        ..
else
    cmake \
        -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
        -DCMAKE_FIND_ROOT_PATH="${prefix}/lib/mpich;${prefix}" \
        -DMPIEXEC_EXECUTABLE=${prefix}/lib/mpich/bin/mpiexec \
        -DBUILD_SHARED_LIBS=ON \
        -DCMAKE_INSTALL_PREFIX=${prefix} \
        ..
fi

cmake --build . --config RelWithDebInfo --parallel $nproc
cmake --build . --config RelWithDebInfo --parallel $nproc --target install

################################################################################
# Install licenses
################################################################################

install_license $WORKSPACE/srcdir/MPItrampoline-*/LICENSE.md $WORKSPACE/srcdir/mpich*/COPYRIGHT
"""

augment_platform_block = """
    using Base.BinaryPlatforms
    $(MPI.augment)
    augment_platform!(platform::Platform) = augment_mpi!(platform)
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; experimental=true)

# MPItrampoline requires `RTLD_DEEPBIND` for `dlopen`, and thus does
# not support musl or BSD.
# FreeBSD: https://reviews.freebsd.org/D24841
platforms = filter(p -> !(Sys.iswindows(p) || libc(p) == "musl"), platforms)
platforms = filter(!Sys.isfreebsd, platforms)
platforms = expand_gfortran_versions(platforms)
# libgfortran3 does not support `!GCC$ ATTRIBUTES NO_ARG_CHECK`. (We
# could in principle build without Fortran support there.)
platforms = filter(p -> libgfortran_version(p) ≠ v"3", platforms)

# Add `mpi+mpitrampoline` platform tag
foreach(p -> (p["mpi"] = "MPItrampoline"), platforms)

# The products that we will ensure are always built
products = [
    ExecutableProduct("mpicc", :mpicc),
    ExecutableProduct("mpicxx", :mpicxx),
    ExecutableProduct("mpifc", :mpifc),
    ExecutableProduct("mpifort", :mpifort),
    ExecutableProduct("mpiexec", :mpiexec),

    # We need to call this library `:libmpi` in Julia so that Julia's
    # `MPI.jl` will find it
    LibraryProduct("libmpitrampoline", :libmpi),

    # MPICH
    ExecutableProduct("mpiexec", :mpich_mpiexec, "lib/mpich/bin"),

    # MPIwrapper
    ExecutableProduct("mpiwrapperexec", :mpiwrapperexec),
    # `libmpiwrapper` is a plugin, not a library, and thus has the
    # suffix `.so` on macOS. This confuses BinaryBuilder.
    # LibraryProduct("libmpiwrapper", :libmpiwrapper; dont_dlopen=true),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"), v"0.5.2"),
    Dependency(PackageSpec(name="MPIPreferences", uuid="3da0fdf6-3ccc-4f1b-acd9-58baa6c99267"); compat="0.1"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               augment_platform_block, julia_compat="1.6")

max_stars_repo_name: staticfloat/BinaryBuilderRepos

max_issues_repo_name: staticfloat/BinaryBuilderRepos

max_forks_repo_name: staticfloat/BinaryBuilderRepos


text:
import Dates
using Dates

include("../util/util.jl")

include("../../conf.jl")

include("base.jl")

include("qGroups.jl")
include("qStacks.jl")
include("qStyles.jl")
include("qAttributes.jl")

include("../providers/data.jl")

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

println("Loading Plots...")

@time begin
    using Plots
    using PlotThemes
    using Plots.PlotMeasures
    # Don't try to open windows for the plots
    ENV["GKSwstype"] = "nul"
    gr()
    # init it now
    theme(PLOT_THEME)
    plot([0], [0])
    logger(clibraries(), "clibraries()")
end

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

function addTicks!(datesI::Array{DateTime,1}, ctx::PlotContext)
    # Add x-tick-labels
    local datesTicks = Array{DateTime,1}()
    local datesLabels = Array{String,1}()
    local dates = sort(datesI) # copy and sort
    local lastDay = 0
    local tickDistance = (dates[end] - dates[1]) ÷ NUMBER_OF_TICKS
    local lastTick = 0

    # add interval and day based ticks
    for date in dates
        local day = Dates.format(date, "dd")
        if day != lastDay
            push!(datesLabels, Dates.format(date, "yy/mm/dd")) # HH:MM"))
            push!(datesTicks, date)
            lastDay = day
            lastTick = date
        elseif NUMBER_OF_TICKS > 0 && date - lastTick >= tickDistance
            push!(datesLabels, Dates.format(date, "HH:MM"))
            push!(datesTicks, date)
            lastTick = date
        end
    end
    # add hour based ticks
    if NUMBER_OF_TICKS < 0
        local cpDatesTicks = copy(datesTicks)
        for i in 1:length(cpDatesTicks)
            local date = cpDatesTicks[i]
            local now = date + Dates.Hour(-NUMBER_OF_TICKS)
            while now < date + Dates.Hour(24) && now < dates[end] && now < ctx.maxDate && (
                i == length(cpDatesTicks) ||
                now < cpDatesTicks[i+1] + Dates.Hour(round(-NUMBER_OF_TICKS / 2))
            )
                push!(datesLabels, Dates.format(now, "HH:MM"))
                push!(datesTicks, now)
                now += Dates.Hour(-NUMBER_OF_TICKS)
            end
        end
    end
    plot!(ctx.plot, xticks = (datesTicks, datesLabels))
end

function addBackground!(dates::Array{DateTime,1}, ctx::PlotContext)
    local lastDay = 0
    for date in dates
        local day = Dates.format(date, "dd")
        if day != lastDay
            local range = [
                Dates.DateTime(Dates.format(date, "yyyy/mm/dd 00:00"), "yyyy/mm/dd HH:MM"),
                min(
                    Dates.DateTime(
                        Dates.format(date, "yyyy/mm/dd 06:00"),
                        "yyyy/mm/dd HH:MM",
                    ),
                    ctx.maxDate,
                ),
            ]
            vspan!(ctx.plot, range, color = :blue, alpha = 0.05, labels = "")
            lastDay = day
        end
    end
end

function addPoints!(points::Dict{Tuple{String,String},Array{DateTime,1}}, ctx::PlotContext)
    # plot points

    for (label, ps) in points
        vline!(ctx.plot, ps, label = label[2], color = nextColor!(ctx.palette))

        # TODO: Scatter is more beautiful. But where on the yaxis do we place the points?
        # _on_ the xaxis would look nice. But that is not possible (?)

        #=scatter!(
            ctx.plot,
            ps,
            [.0 for _ in 1:length(ps)],
            label = label[2],
            color = nextColor!(ctx.palette),
        )=#
    end
end

global globalDataStore = DataStore()

function analyze(widget::QPlot, ctx::PlotContext)
    local dates = Set{DateTime}()
    local srcCtx = Dict{String,DataAttributeContext}()
    local mini = typemax(Float64)
    local maxi = typemin(Float64)
    local titles = Array{String,1}()

    for attr in widget.data
        local _dates, _srcCtx, _mini, _maxi, _titles = analyze(attr, ctx)

        for c in _dates
            push!(dates, c)
        end
        for (k, v) in _srcCtx
            srcCtx[k] = v
        end
        mini = min(mini, _mini)
        maxi = max(maxi, _maxi)
        titles = [titles..., _titles...]
    end

    return sort(collect(dates)), srcCtx, mini, maxi, titles
end

function renderWidget(widget::QPlot, today::DateTime, saveTo::String)
    logger("", "Plotting $(widget.title)", true)

    local from = today - Dates.Day(widget.days)
    local to = today

    local size = widget.size
    if size === nothing
        size = DEFAULT_SIZE
    end

    # Get the DataAttributeContext
    local _, srcCtx, _, __, titles = analyze(
        widget,
        PlotContext(today = today, store = globalDataStore, maxDate = to, minDate = from),
    )

    # aggregate the context
    local aggregatedCtx = Dict{Tuple{String,String},Array{String,1}}()
    for (k, v) in srcCtx
        for w in v.ctx
            if !((w[1], w[2]) in keys(aggregatedCtx))
                aggregatedCtx[(w[1], w[2])] = []
            end
            push!(aggregatedCtx[(w[1], w[2])], w[3])
        end
    end

    # multiply the context
    local dataAttrCtxs = Array{DataAttributeContext,1}()
    for (k, v) in aggregatedCtx
        local nDataAttrCtxs = Array{DataAttributeContext,1}()
        for w in v
            if length(dataAttrCtxs) == 0
                push!(nDataAttrCtxs, DataAttributeContext(Set([(k..., w)])))
            else
                for c in dataAttrCtxs
                    local cc = copy(c)
                    push!(cc.ctx, (k..., w))
                    push!(nDataAttrCtxs, cc)
                end
            end
        end
        dataAttrCtxs = nDataAttrCtxs
    end

    local results = []

    if length(dataAttrCtxs) == 0
        push!(dataAttrCtxs, DataAttributeContext(Set()))
    end

    # Iterate all elements of the product
    for dataAttrCtx in dataAttrCtxs
        local ctxString = ""
        local description = ""
        # Generate Description from Specialization Context!
        for v in dataAttrCtx.ctx
            ctxString *= "$(v[1])($(v[2])) = $(v[3]) "
            local desc = ""
            if (v[1], v[2]) in keys(globalDataStore.descriptions)
                desc = globalDataStore.descriptions[(v[1], v[2])]
            end
            description *= "$(v[1])($(v[2])) = $(v[3])$desc
\n"
        end
        description = strip(description)
        ctxString = strip(ctxString)
        local title = "$(widget.title) $ctxString"

        logger(title, "Plotting with context", true)

        local ctx = PlotContext(
            palette = loadPalette(PLOT_THEME),
            today = today,
            store = globalDataStore,
            maxDate = to,
            minDate = from,
        )

        # init the plot
        ctx.plot = plot(
            left_margin = 5mm,
            title = title,
            size = size,
            legend = :topleft,
            yformatter = Y_FORMATTER,
        )

        ctx.styles["source context"] = dataAttrCtx
        local dates, _, _, _, titles = analyze(widget, ctx)

        if DRAW_NIGHT_BACKGROUND && Dates.days(to - from) <= 10
            addBackground!(dates, ctx)
            # vspan overwrites the default format. Restore it.
            plot!(ctx.plot, yformatter = Y_FORMATTER)
        end

        for el in widget.data
            ctx.styles = defaultStyles()
            ctx.styles["source context"] = dataAttrCtx
            renderWidget!(el, ctx)
        end

        local points = getPoints(globalDataStore, from, to)

        addPoints!(points, ctx)
        addTicks!(dates, ctx)
        plot!(ctx.plot, yformatter = Y_FORMATTER)

        # save to plot
        mkpath(saveTo)
        local path = joinpath(saveTo, "$(sanitizeFile(title)).png")
        savefig(ctx.plot, path)

        titles = map(titles) do x
            local str = ""
            if x in keys(globalDataStore.descriptions)
                str = globalDataStore.descriptions[x]
            end
            (x..., str)
        end

        @assert isfile("$path") "Error saving the Plot"

        # Apply strong png compression to make e-mail smaller
        if USE_PNGQUANT
            exe(`pngquant --quality=60-80 --force --output $path $path`)
        end

        push!(results, (path, titles, title, description))
    end

    return results
end

max_stars_repo_name: Agapanthus/server-diary

max_issues_repo_name: Agapanthus/server-diary

max_forks_repo_name: Agapanthus/server-diary


text:
# ## Line style

#md # ![FILE_NAME.png](images/FILE_NAME.png)

Plots.supported_styles()

#-

## ref https://jp.mathworks.com/matlabcentral/fileexchange/35230-matlab-plot-gallery-standard-line-styles?focused=6792901&tab=example

using SpecialFunctions

style = Plots.supported_styles()[2:end]
style = reshape(style, 1, length(style))
x = 0:0.2:10;
y0 = besselj.(0,x);
y1 = besselj.(1,x);
y2 = besselj.(2,x);
y3 = besselj.(3,x);
y4 = besselj.(4,x);

plot(x, [y0 y1 y2 y3 y4], line=(3, style))


text:
struct WhiteRBM{V,H,W,Av,Ah}
    visible::V
    hidden::H
    w::W
    affine_v::Av
    affine_h::Ah
    """
        WhiteRBM(rbm, affine_v, affine_h)

    Creates a whitened RBM with the given transforms. Note that the resulting model is
    not equivalent to the original `rbm`. To obtain an equivalent model (where energies
    differ by a constant), see [`whiten`](@ref).
    """
    function WhiteRBM(rbm::RBM{V,H,W}, affine_v::Affine, affine_h::Affine) where {V,H,W}
        @assert length(visible(rbm)) == length(affine_v.u)
        @assert length(hidden(rbm)) == length(affine_h.u)
        return new{V, H, W, typeof(affine_v), typeof(affine_h)}(
            visible(rbm), hidden(rbm), weights(rbm), affine_v, affine_h
        )
    end
end

function WhiteRBM(visible, hidden, w, affine_v, affine_h)
    rbm = RBM(visible, hidden, w)
    return WhiteRBM(rbm, affine_v, affine_h)
end

const AffineRBM{Av,Ah,V,H,W} = WhiteRBM{V,H,W,Av,Ah}

RBMs.visible(white_rbm::WhiteRBM) = white_rbm.visible
RBMs.hidden(white_rbm::WhiteRBM) = white_rbm.hidden
RBMs.weights(white_rbm::WhiteRBM) = white_rbm.w

"""
    WhiteRBM(rbm)

Creates a WhiteRBM with identity transforms.
"""
function WhiteRBM(rbm::RBM)
    T = eltype(weights(rbm))
    N = length(visible(rbm))
    M = length(hidden(rbm))
    affine_v = Affine(Diagonal(ones(T, N)), zeros(T, N))
    affine_h = Affine(Diagonal(ones(T, M)), zeros(T, M))
    return WhiteRBM(rbm, affine_v, affine_h)
end

"""
    RBM(white_rbm::WhiteRBM)

Returns an (unwhitened) `RBM` which neglects the affine transforms of `white_rbm`.
The resulting model is *not* equivalent to the original `white_rbm`.
To construct an equivalent model, use the function
`blacken(white_rbm)` instead (see [`blacken`](@ref)).
"""
function RBMs.RBM(white_rbm::WhiteRBM)
    return RBM(visible(white_rbm), hidden(white_rbm), weights(white_rbm))
end

function whiten_v(white_rbm::WhiteRBM, v::AbstractArray)
    return reshape(white_rbm.affine_v * flatten(visible(white_rbm), v), size(v))
end

function whiten_h(white_rbm::WhiteRBM, h::AbstractArray)
    return reshape(white_rbm.affine_h * flatten(hidden(white_rbm), h), size(h))
end

function RBMs.energy(white_rbm::WhiteRBM, v::AbstractArray, h::AbstractArray)
    Ev = RBMs.energy(visible(white_rbm), v)
    Eh = RBMs.energy(hidden(white_rbm), h)
    Ew = RBMs.interaction_energy(white_rbm, v, h)
    return Ev .+ Eh .+ Ew
end

function RBMs.interaction_energy(white_rbm::WhiteRBM, v::AbstractArray, h::AbstractArray)
    white_v = whiten_v(white_rbm, v)
    white_h = whiten_h(white_rbm, h)
    return RBMs.interaction_energy(RBM(white_rbm), white_v, white_h)
end

function RBMs.free_energy(white_rbm::WhiteRBM, v::AbstractArray)
    inputs = RBMs.inputs_v_to_h(white_rbm, v)
    E_vis = RBMs.energy(visible(white_rbm), v)
    F_hid = RBMs.free_energy(hidden(white_rbm), inputs)
    b = reshape(white_rbm.affine_h.u, size(hidden(white_rbm)))
    ΔE = RBMs.energy(RBMs.Binary(b), inputs)
    return E_vis - ΔE + F_hid
end

function RBMs.inputs_v_to_h(white_rbm::WhiteRBM, v::AbstractArray)
    white_v = whiten_v(white_rbm, v)
    inputs = inputs_v_to_h(RBM(white_rbm), white_v)
    I_flat = flatten(hidden(white_rbm), inputs)
    return reshape(white_rbm.affine_h.A' * I_flat, size(inputs))
end

function RBMs.inputs_h_to_v(white_rbm::WhiteRBM, h::AbstractArray)
    white_h = whiten_h(white_rbm, h)
    inputs = inputs_h_to_v(RBM(white_rbm), white_h)
    I_flat = flatten(visible(white_rbm), inputs)
    return reshape(white_rbm.affine_v.A' * I_flat, size(inputs))
end

function RBMs.sample_h_from_v(white_rbm::WhiteRBM, v::AbstractArray)
    inputs = inputs_v_to_h(white_rbm, v)
    return RBMs.sample_from_inputs(hidden(white_rbm), inputs)
end

function RBMs.sample_v_from_h(white_rbm::WhiteRBM, h::AbstractArray)
    inputs = inputs_h_to_v(white_rbm, h)
    return RBMs.sample_from_inputs(visible(white_rbm), inputs)
end

function RBMs.sample_v_from_v(white_rbm::WhiteRBM, v::AbstractArray; steps::Int = 1)
    @assert size(visible(white_rbm)) == size(v)[1:ndims(visible(white_rbm))]
    for _ in 1:steps
        v = oftype(v, RBMs.sample_v_from_v_once(white_rbm, v))
    end
    return v
end

function RBMs.sample_h_from_h(white_rbm::WhiteRBM, h::AbstractArray; steps::Int = 1)
    @assert size(hidden(white_rbm)) == size(h)[1:ndims(hidden(white_rbm))]
    for _ in 1:steps
        h = oftype(h, RBMs.sample_h_from_h_once(white_rbm, h))
    end
    return h
end

function RBMs.sample_v_from_v_once(white_rbm::WhiteRBM, v::AbstractArray)
    h = RBMs.sample_h_from_v(white_rbm, v)
    v = RBMs.sample_v_from_h(white_rbm, h)
    return v
end

function RBMs.sample_h_from_h_once(white_rbm::WhiteRBM, h::AbstractArray)
    v = RBMs.sample_v_from_h(white_rbm, h)
    h = RBMs.sample_h_from_v(white_rbm, v)
    return h
end

function RBMs.mean_h_from_v(white_rbm::WhiteRBM, v::AbstractArray)
    inputs = RBMs.inputs_v_to_h(white_rbm, v)
    return RBMs.mean_from_inputs(hidden(white_rbm), inputs)
end

function RBMs.mean_v_from_h(white_rbm::WhiteRBM, h::AbstractArray)
    inputs = RBMs.inputs_h_to_v(white_rbm, h)
    return RBMs.mean_from_inputs(visible(white_rbm), inputs)
end

function RBMs.mode_v_from_h(white_rbm::WhiteRBM, h::AbstractArray)
    inputs = RBMs.inputs_h_to_v(white_rbm, h)
    return RBMs.mode_from_inputs(visible(white_rbm), inputs)
end

function RBMs.mode_h_from_v(white_rbm::WhiteRBM, v::AbstractArray)
    inputs = RBMs.inputs_v_to_h(white_rbm, v)
    return RBMs.mode_from_inputs(hidden(white_rbm), inputs)
end

function RBMs.reconstruction_error(white_rbm::WhiteRBM, v::AbstractArray; steps::Int = 1)
    @assert size(visible(white_rbm)) == size(v)[1:ndims(visible(white_rbm))]
    v1 = RBMs.sample_v_from_v(white_rbm, v; steps)
    ϵ = Statistics.mean(abs.(v .- v1); dims = 1:ndims(visible(white_rbm)))
    if ndims(v) == ndims(visible(white_rbm))
        return only(ϵ)
    else
        return reshape(ϵ, size(v)[end])
    end
end

function RBMs.mirror(white_rbm::WhiteRBM)
    function p(i)
        if i ≤ ndims(visible(white_rbm))
            return i + ndims(hidden(white_rbm))
        else
            return i - ndims(visible(white_rbm))
        end
    end
    perm = ntuple(p, ndims(weights(white_rbm)))
    w = permutedims(weights(white_rbm), perm)
    rbm = RBM(hidden(white_rbm), visible(white_rbm), w)
    return WhiteRBM(rbm, white_rbm.affine_h, white_rbm.affine_v)
end

function RBMs.∂free_energy(
    white_rbm::WhiteRBM, v::AbstractArray; wts = nothing,
    stats = RBMs.suffstats(visible(white_rbm), v; wts)
)
    inputs = RBMs.inputs_v_to_h(white_rbm, v)
    h = RBMs.mean_from_inputs(hidden(white_rbm), inputs)
    ∂v = RBMs.∂energy(visible(white_rbm), stats)
    ∂h = RBMs.∂free_energy(hidden(white_rbm), inputs; wts)
    ∂w = RBMs.∂interaction_energy(white_rbm, v, h; wts)
    return (visible = ∂v, hidden = ∂h, w = ∂w)
end

function RBMs.∂interaction_energy(
    white_rbm::WhiteRBM, v::AbstractArray, h::AbstractArray; wts = nothing
)
    white_v = whiten_v(white_rbm, v)
    white_h = whiten_h(white_rbm, h)
    ∂w = RBMs.∂interaction_energy(RBM(white_rbm), white_v, white_h; wts)
    return ∂w
end

function RBMs.log_pseudolikelihood(rbm::WhiteRBM, v::AbstractArray)
    return RBMs.log_pseudolikelihood(blacken(rbm), v)
end


text:
#=
# 5. Viscous flow about a moving body
In this notebook we will demonstrate the simulation of a moving body. It is straightforward
to set up a moving body. The main caveat is that the simulation is slower,
because the integrator must update the operators continuously throughout the simulation.

We will demonstrate this on an oscillating flat plate.
=#

using ViscousFlow
#-
using Plots

#=
### Problem specification and discretization
=#
Re = 200; # Reynolds number
U = 1.0; # Free stream velocity
U∞ = (U,0.0);
#-
xlim = (-1.0,3.0)
ylim = (-1.0,1.0);
Δx, Δt = setstepsizes(Re,gridRe=4.0)

#=
### Set up body
Set up the plate and place it at the origin
=#
body = Plate(1.0,1.0Δx)
T = RigidTransform((0.,0.),0.)
T(body)

#=
### Set the body motion
Now we specify the body motion. We will use oscillatory pitch-heave kinematics for this:
=#
a = 0.25 # location of pitch axis, a = 0.5 is leading edge
ϕp = -π/2  # phase lag of pitch
ϕh = 0.0  # phase lag of heave
A = 0.25  # amplitude/chord
fstar = 1/π # fc/U
α₀ = 0 # mean angle of attack
Δα = 10π/180 # amplitude of pitching
U₀ = 0.0 # translational motion (set to zero in place of free stream)
K = π*fstar # reduced frequency, K = πfc/U

oscil1 = RigidBodyTools.PitchHeave(U₀,a,K,ϕp,α₀,Δα,A,ϕh)
motion = RigidBodyMotion(oscil1)

# We can inspect the kinematics in this `motion` by plotting them:
plot(motion)

#=
### Construct the system structure
Here, we supply the motion as an another argument.
=#
sys = NavierStokes(Re,Δx,xlim,ylim,Δt,body,motion,freestream = U∞)
#-
u0 = newstate(sys)
tspan = (0.0,10.0)
integrator = init(u0,tspan,sys)

#=
### Solve
This takes longer than it does for stationary bodies. Here, we only run it
for a little while just to demonstrate it.
=#
step!(integrator,0.1)

# ### Examine the solution
plot(vorticity(integrator),sys,clim=(-10,10),levels=range(-10,10,length=30),color = :RdBu)

# and the forces
sol = integrator.sol
fx, fy = force(sol,sys,1);
#-
plot(
plot(sol.t,2*fx,xlim=(0,Inf),ylim=(-3,3),xlabel="Convective time",ylabel="\$C_D\$",legend=:false),
plot(sol.t,2*fy,xlim=(0,Inf),ylim=(-6,6),xlabel="Convective time",ylabel="\$C_L\$",legend=:false),
    size=(800,350)
)


text:
"""
TagType is a parameterized type which allows
making Index tags into Julia types. One use case
is overloading functions such as `op` which
generates physics operators for indices
with certain tags such as "S=1/2".

To make a TagType, you can use the string
macro notation: `TagType"MyTag"`
"""
struct TagType{T}
end

macro TagType_str(s)
  TagType{Tag(s)}
end

function _call_op(s::Index,
                  opname::AbstractString;
                  kwargs...)
  use_tag = 0
  nfound = 0
  for n=1:length(tags(s))
    TType = TagType{tags(s)[n]}
    if hasmethod(op,Tuple{TType,Index,AbstractString})
      use_tag = n
      nfound += 1
    end
  end
  if nfound == 0
    error("Overload of \"op\" function not found for Index tags $(tags(s))")
  elseif nfound > 1
    error("Multiple tags from $(tags(s)) overload the function \"op\"")
  end
  TType = TagType{tags(s)[use_tag]}
  return op(TType(),s,opname;kwargs...)
end

function op(s::Index,
            opname::AbstractString;
            kwargs...)::ITensor

  opname = strip(opname)

  if opname == "Id"
    Op = emptyITensor(dag(s),s')
    for n=1:dim(s)
      Op[dag(s)(n),s'(n)] = 1.0
    end
    return Op
  end

  # Interpret operator names joined by *
  # as acting sequentially on the same site
  starpos = findfirst("*",opname)
  if !isnothing(starpos)
    op1 = opname[1:starpos.start-1]
    op2 = opname[starpos.start+1:end]
    return product(op(s,op1;kwargs...),op(s,op2;kwargs...))
  end

  return _call_op(s,opname;kwargs...)
end

function op(s::Vector{<:Index},
            opname::AbstractString,
            n::Int;
            kwargs...)::ITensor
  return op(s[n],opname;kwargs...)
end

state(s::Index,n::Integer) = s[n]

function state(s::Index,
               str::String)::IndexVal
  use_tag = 0
  nfound = 0
  for n=1:length(tags(s))
    TType = TagType{tags(s)[n]}
    if hasmethod(state,Tuple{TType,AbstractString})
      use_tag = n
      nfound += 1
    end
  end
  if nfound == 0
    error("Overload of \"state\" function not found for Index tags $(tags(s))")
  elseif nfound > 1
    error("Multiple tags from $(tags(s)) overload the function \"state\"")
  end
  TType = TagType{tags(s)[use_tag]}
  sn = state(TType(),str)
  return s[sn]
end

function state(sset::Vector{<:Index},
               j::Integer,
               st)::IndexVal
  return state(sset[j],st)
end

function siteinds(d::Integer,
                  N::Integer; kwargs...)
  return [Index(d,"Site,n=$n") for n=1:N]
end

function siteinds(str::String,
                  N::Integer; kwargs...)
  TType = TagType{Tag(str)}
  if !hasmethod(siteinds,Tuple{TType,Int})
    error("Overload of \"siteinds\" function not found for tag type \"$str\"")
  end
  return siteinds(TType(),N; kwargs...)
end

function has_fermion_string(s::Index,
                            opname::AbstractString;
                            kwargs...)::Bool
  opname = strip(opname)
  use_tag = 0
  nfound = 0
  for n=1:length(tags(s))
    TType = TagType{tags(s)[n]}
    if hasmethod(has_fermion_string,Tuple{TType,Index,AbstractString})
      use_tag = n
      nfound += 1
    end
  end
  if nfound == 0
    return false
  elseif nfound > 1
    error("Multiple tags from $(tags(s)) overload the function \"has_fermion_string\"")
  end
  TType = TagType{tags(s)[use_tag]}
  return has_fermion_string(TType(),s,opname;kwargs...)
end


text:
module diffevo

using BlackBoxOptimizationBenchmarking
using Plots
using Statistics
include("plot_utils.jl")

const CR = 0.941343846536435
const F = 0.5184306753627251
const NP = 100
const DIM = 10
const MAX_GEN = 800
const NEXPERIMENTS = 10
const EXPERIMENT_NAME = "diffevo_time";
const COMPARE_TO = "nothing";

const GENS_PER_FRAME = 100

mutable struct Agent
    pos::Vector{Float64}
    obj::Float64
    fit::Float64
    Agent() = begin
        pos = [randu(-4.0, 4.0) for _ in 1:DIM]
        return new(pos, Inf, -Inf)
    end
end

function randu(lo::Float64, hi::Float64)
    return (hi - lo)*rand() + lo
end

function copy!(dest::Agent, src::Agent)
    dest.pos .= src.pos
    dest.obj = src.obj
    dest.fit = src.fit
end

const Population = Vector{Agent}

function find_best(pop::Population)::Agent
    if length(pop) == 0
        @error("Empty population!")
    end
    winner = pop[1]
    for x in pop
        if x.fit > winner.fit
            winner = x
        end
    end
    return winner
end

function eval!(x::Agent, fun::BBOBFunction)
    val = fun(x.pos)
    x.fit = -val
    x.obj = val - fun.f_opt + 1e-15
end

function breed!(parents::Population, children::Population, fun::BBOBFunction)
    Threads.@threads for i in 1:NP
        x = parents[i]
        y = children[i]
        a = rand(parents)
        b = rand(parents)
        c = rand(parents)
        R = rand(1:DIM)
        for i in 1:DIM
            if R == i || rand() < CR
                y.pos[i] = a.pos[i] + F*(b.pos[i] - c.pos[i])
            else
                y.pos[i] = x.pos[i]
            end
        end
        eval!(y, fun)
    end
end

function exchange!(parents::Population, children::Population)
    for i in 1:NP
        if parents[i].fit <= children[i].fit
            temp = children[i]
            children[i] = parents[i]
            parents[i] = temp
        end
    end
end

function eva(fun::BBOBFunction)
    log_best = Vector{Float64}()
    #random initial population
    parents = [Agent() for _ in 1:NP]
    children = [Agent() for _ in 1:NP]
    best = parents[1]

    #eval them
    for x in parents
        eval!(x, fun)
    end
    #then iterate:
    for gen in 1:MAX_GEN
        breed!(parents, children, fun)
        exchange!(parents, children)
        best = find_best(parents)
        #save results
        if (gen % GENS_PER_FRAME == 0)
            println(gen, '\t', best.obj)
        end
        push!(log_best, best.obj)
    end
    return (best, log_best)
end

function main()
    bbob_funs = BlackBoxOptimizationBenchmarking.list_functions()
    for n in [2,6,8,10]
        log = zeros(NEXPERIMENTS, MAX_GEN)
        time_avg = 0.0
        for experiment in 1:NEXPERIMENTS
            println("experiment #", experiment)
            stats = @timed (result, log_best) = eva(bbob_funs[n])
            for gen in 1:MAX_GEN
                log[experiment, gen] = log_best[gen]
            end
            time_avg += stats.time/NEXPERIMENTS
            #uncomment this if you want to save the best solution
            #saveind(result, "julia-partition-results/"*EXPERIMENT_NAME*"-"*string(k)*".best")
        end
        t = LinRange(0.0, time_avg, MAX_GEN)
        new_path = "julia-cont_opt-results/"*EXPERIMENT_NAME*string(n)*".log"
        old_path = "julia-cont_opt-results/"*COMPARE_TO*string(n)*".log"
        fig_path = "julia-cont_opt-results/"*EXPERIMENT_NAME*string(n)*".pdf"
        plot_utils.save_log(t, log, new_path)
        if ispath(old_path)
            p = plot_utils.plot_log(old_path, new_path)
            savefig(p, fig_path)
        else
            p = plot_utils.plot_log(new_path)
            savefig(p, fig_path)
        end
    end
end

end

max_stars_repo_name: OndrejKincl/evaTeaching-python

max_issues_repo_name: OndrejKincl/evaTeaching-python

max_forks_repo_name: OndrejKincl/evaTeaching-python


text:
ByteArr = Vector{UInt8} # byte array

"""
    BinDescr

Describe (a chunk) of binary data, e.g. numeric array or image in the LVServer format used
for exchange between Julia and LabVIEW. Gets converted from/to JSON for transport.
"""
mutable struct BinDescr
    nofbytes::Int
    numtype::String
    start::Int
    arrdims::Vector{Int}
    kwarg_name::Union{String,Symbol}
    category::String
end

BinDescr() = BinDescr(0, "", 1, [], "", "numarrays")
JSON3.StructTypes.StructType(::Type{BinDescr}) = JSON3.StructTypes.Struct()


text:
abstract type Loop end

const LoopBody = Vector{<:Recurrence}

struct EmptyLoop <: Loop end

struct SingleLoop <: Loop
    body::LoopBody
    lc::Sym
    vars::Array{Sym,1}
end

struct MultiLoop <: Loop
    branches::Array{SingleLoop}
    vars::Array{Sym,1}
end

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

function Base.show(io::IO, body::LoopBody)
    print(io, "[")
    join(io, body, ", ")
    println(io, "]")
end

function Base.show(io::IO, loop::SingleLoop)
    if get(io, :compact, false)
        Base.show(io, loop.body)
    else
        println(io, "$(length(loop.body))-element $(typeof(loop)):")
        for l in loop.body
            print(io, "  $(l)\n")
        end
    end
end

function Base.show(io::IO, loop::MultiLoop)
    println(io, "$(length(loop.branches))-element $(typeof(loop)):")
    for l in loop.branches
        print(io, " ")
        summary(io, l)
        # print(io, "\n")
    end
end

function closed_forms(loop::MultiLoop)
    [ClosedFormSystem(rec_solve(l.body), l.lc, l.vars) for l in loop.branches]
end

function closed_forms(loop::SingleLoop)
    ClosedFormSystem(rec_solve(loop.body), loop.lc, loop.vars)
end


text:
function WeightingOp(in,adj;w=1)

	return in.*w

end

function WeightingOp(m::AbstractString,d::AbstractString,adj;w="NULL")

	if (adj==true)
		d1,h1,e1 = SeisRead(d)
		d2,h2,e2 = SeisRead(w)
		SeisWrite(m,d1[:,:].*d2[:,:],h1,e1)
	else
		d1,h1,e1 = SeisRead(m)
		d2,h2,e2 = SeisRead(w)
		SeisWrite(d,d1[:,:].*d2[:,:],h1,e1)
	end

end


text:
function average_loss(m::Model, loss::Loss, inputs::Vector, weights::Vector)
    l, n = 0.0, 0
    n_inputs = length(inputs)
    for i in 1:n_inputs
        @inbounds x2, t2, x1, t1, y = inputs[i]
        @inbounds w = weights[i]
        p = predict(m, x2, t2, x1, t1)
        l += w * loss_value(loss, p, y)
        n += w
    end
    return l / n
end


text:
using Test
using ImageInTerminal, TestImages, ImageCore, ImageTransformations
using Random

if isinteractive()
    @info ("In interactive use, one should respond \"n\" when the program"
           * " offers to create or replace files associated with some tests.")
else
    @info ("Nine tests should correctly report failure in the transcript"
           * " (but not the test summary).")
end
# check for ambiguities
refambs = detect_ambiguities(ImageInTerminal, Base, Core)

using ReferenceTests
ambs = detect_ambiguities(ReferenceTests, ImageInTerminal, Base, Core)

@testset "ReferenceTests" begin

@test Set(setdiff(ambs, refambs)) == Set{Tuple{Method,Method}}()

# load/create some example images
lena = testimage("lena_color_256")
camera = testimage("cameraman")
cameras = similar(camera, size(camera)..., 2)
copyto!(view(cameras,:,:,1), camera)
copyto!(view(cameras,:,:,2), camera)
square = Gray{N0f8}[0.1 0.2 0.3; 0.4 0.5 0.6; 0.7 0.6 0.9]
rgb_rect = rand(RGB{N0f8}, 2, 3)

@testset "io2str" begin
    @test_throws LoadError eval(@macroexpand @io2str(::IO))
    @test_throws ArgumentError @io2str(2)
    @test_throws ArgumentError @io2str(string(2))
    @test @io2str(print(::IO, "foo")) == "foo"
    @test @io2str(println(::IO, "foo")) == "foo\n"
    @test @io2str(show(::IO, "foo")) == "\"foo\""
    A = ones(30,30)
    @test @io2str(show(IOContext(::IO, :limit => true, :displaysize => (5,5)), A)) == "[1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0; … ; 1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0]"
end

@testset "withcolor" begin
    @test_throws ArgumentError @withcolor throw(ArgumentError("foo"))
    @test @withcolor Base.have_color == true
    @test @withcolor @io2str(printstyled(IOContext(::IO, :color => true), "test", color=:green)) == "\e[32mtest\e[39m"
end

@testset "string as txt" begin
    foo = "foo"
    @test_reference "references/string1.txt" foo * "bar"
    @test_reference "references/string1.txt" [foo * "bar"]
    A = ones(30,30)
    @test_reference "references/string2.txt" @io2str show(IOContext(::IO, :limit=>true, :displaysize=>(5,5)), A)
    @test_reference "references/string3.txt" 1337
    @test_reference "references/string3.txt" 1338 by=(ref, x)->isapprox(ref, x; atol=10)
    @test_reference "references/string4.txt" @io2str show(::IO, MIME"text/plain"(), Int64.(collect(1:5)))

    @test_reference "references/string5.txt" """
        This is a
        multiline string that does not end with a new line."""

    @test_reference "references/string6.txt" """
        This on the other hand is a
        multiline string that does indeed end with a new line.
    """

    @test_throws ErrorException @test_reference "references/string1.txt" "intentionally wrong to check that this message prints"
    @test_throws ErrorException @test_reference "references/wrong.txt" "intentional error to check that this message prints"
    @test_throws ErrorException @test_reference "references/string5.txt" """
        This is an incorrect
        multiline string that does not end with a new line."""
end

@testset "string as unknown file type" begin
    @test_reference "references/string1.nottxt" "This is not a .txt file, but it should be treated as such.\n"
end

@testset "images as txt using ImageInTerminal" begin
    #@test_throws MethodError @test_reference "references/fail.txt" rand(2,2)
    @test_reference "references/camera.txt" camera size=(5,10)
    @test_reference "references/lena.txt" lena
end

@testset "plain ansi string" begin
    @test_reference(
        "references/ansii.txt",
        @io2str(printstyled(IOContext(::IO, :color=>true), "this should be blue", color=:blue)),
        render = ReferenceTests.BeforeAfterFull()
    )
    @test_throws ErrorException @test_reference(
        "references/ansii.txt",
        @io2str(printstyled(IOContext(::IO, :color=>true), "this should be red", color=:red)),
        render = ReferenceTests.BeforeAfterFull()
    )
end

@testset "string as SHA" begin
    @test_reference "references/number1.sha256" 1337
    foo = "foo"
    @test_reference "references/string1.sha256" foo * "bar"
    A = ones(30,30)
    @test_reference "references/string2.sha256" @io2str show(IOContext(::IO, :limit=>true, :displaysize=>(5,5)), A)
end

@testset "images as SHA" begin
    @test_reference "references/camera.sha256" camera
    @test_reference "references/lena.sha256" convert(Matrix{RGB{Float64}}, lena)
end

@testset "images as PNG" begin
    @test_reference "references/camera.png" imresize(camera, (64,64))
    @test_reference "references/camera.png" imresize(camera, (64,64)) by=psnr_equality(25)
    @test_throws Exception @test_reference "references/wrongfilename.png" imresize(camera, (64,64))
    @test_throws ErrorException @test_reference "references/camera.png" imresize(lena, (64,64))
    @test_throws Exception @test_reference "references/camera.png" camera # unequal size
end

using DataFrames, CSVFiles
@testset "DataFrame as CSV" begin
    @test_reference "references/dataframe.csv" DataFrame(v1=[1,2,3], v2=["a","b","c"])
    @test_throws ErrorException @test_reference "references/wrongfilename.csv" DataFrame(v1=[1,2,3], v2=["a","b","c"])
    @test_throws ErrorException @test_reference "references/dataframe.csv" DataFrame(v1=[1,2,3], v2=["c","b","c"])
end

end


text:
#KMER ORDER/SEQUENCE INTEGER CODING UTILITIES
#higher order DNA alphabet
struct CompoundAlphabet
  symbols::Dict{Mer{DNAAlphabet{2}},Integer}
  integers::Dict{Integer,Mer{DNAAlphabet{2}}}
  #build a CompoundAlphabet for DNA of some order_no
    function CompoundAlphabet(alphabet::Tuple, order_no::Integer)
        symbols = Dict{Mer{DNAAlphabet{2}},Integer}()
        integers = Dict{Integer,Mer{DNAAlphabet{2}}}()

        tuples = Array{Tuple}
        if order_no > 0
            tuples = collect(Iterators.product([alphabet for order in 0:order_no]...))
        else #0th order compound product of an alphabet is that alphabet
            tuples = collect(Iterators.product(alphabet))
        end

        @inbounds for index in eachindex(tuples)
            tuple_seq = Mer{DNAAlphabet{2}}(collect(tuples[index]))
            integers[index] = tuple_seq
            symbols[tuple_seq] = index
        end

        new(symbols,integers)
    end
end

struct N_Order_ntSequence
    alphabet::CompoundAlphabet
    seq_lengths::Vector{Integer}
    order_kmers::Vector{Vector{Mer{DNAAlphabet{2}}}}
  end

function code_job_obs(job_ids::Vector{Chain_ID}, obs_sets::Dict{String,Vector{LongSequence{DNAAlphabet{2}}}})
    code_jobs=Vector{Tuple{String,Integer}}()
    for id in job_ids #assemble a vector of observation encoding jobs
        code_job=(id.obs_id, id.order)
        !(code_job in code_jobs) && push!(code_jobs, code_job)
    end

    code_dict = Dict{Tuple{String,Integer}, AbstractArray}()
    @showprogress 1 "Encoding observations..." for (obs_id, order) in code_jobs #build the appropriate sample sets once
        order_seqs = get_order_n_seqs(obs_sets[obs_id],order) #get the kmer sequences at the appropriate order
        coded_seqs = code_seqs(order_seqs, sorted=true) #numerically code the sequences in trainable format
        code_dict[(obs_id, order)] = coded_seqs
    end
    return code_dict
end

#from a vector of LongSequences, get
function get_order_n_seqs(seqs::Vector{LongSequence{DNAAlphabet{2}}}, order_no::Integer, base_tuple::Tuple=ACGT)
    kmer_vecs = Vector{Vector{Mer{DNAAlphabet{2}}}}()
    length_vec = Vector{Integer}()
    window = order_no + 1

    for seq in seqs
        kmer_vec = Vector{Mer{DNAAlphabet{2}}}()
        @inbounds for (i, kmer) in collect(each(Mer{DNAAlphabet{2},window},seq))
            push!(kmer_vec, kmer)
        end

        push!(kmer_vecs, kmer_vec)
        push!(length_vec, length(seq))
    end

    return nordseqs = N_Order_ntSequence(CompoundAlphabet(base_tuple, order_no), length_vec, kmer_vecs)
end

#convert tuple kmers to symbol codes
function code_seqs(input::N_Order_ntSequence, offsets::Vector=[0 for i in 1:length(input.order_kmers)]; sorted::Bool=false)
    symbol_no=length(input.alphabet.symbols)
    if symbol_no <= typemax(UInt8)
        integer_type = UInt8
    elseif typemax(UInt8) < symbol_no < typemax(UInt16)
        integer_type = UInt16
    else
        integer_type = UInt32
    end

    alphabet = input.alphabet
    output = zeros(integer_type,  length(input.order_kmers), (maximum([length(seq) for seq in input.order_kmers])+1)) #leave 1 missing value after the longest sequence forindexing sequence length in CLHMM messages
    sorted && (sort_idxs = sortperm(input.seq_lengths,rev=true))
    sorted ? input_seqs = input.order_kmers[sort_idxs] : input_seqs = input.order_kmers

    for (i, seq) in enumerate(input_seqs)
        for t in 1:length(seq)
            curr_kmer = seq[t]
            curr_code = alphabet.symbols[curr_kmer]
            output[i,t+offsets[i]]=curr_code
        end
    end
    return output
end


text:
"""
    MutateAudioGenerator

An audio generator that wraps another.  It mutates its audio with the given
function.  This is useful for applying an envelope or distortion to a signal.

# Example
```julia
gen = MutateGen(
    SineAudioGenerator(0, 0.5, 400),
    (times, samples) -> (samples .* (1 .- exp.(times.*-80))
                                 .* exp.(times.*-3))
)
```
"""
mutable struct MutateAudioGenerator <: AbstractAudioGenerator
    gen::AbstractAudioGenerator
    env::Function  # Applies the envelope to the given (times, samples)
    _t::Float64
    MutateAudioGenerator(gen, env) = new(gen, env, 0.)
end

function reset!(gen::MutateAudioGenerator)
    reset!(gen.gen)
    gen._t = 0.
    nothing
end
function has_samples(gen::MutateAudioGenerator, sample_rate::Float64)
    has_samples(gen.gen, sample_rate)
end
function next_samples!(gen::MutateAudioGenerator, max_n::Int,
                       sample_rate::Float64)
    times = StepRangeLen(gen._t, 1/sample_rate, max_n)
    gen._t += max_n / sample_rate
    samples = next_samples!(gen.gen, max_n, sample_rate)
    gen.env(times, samples)
end


text:
"""
    StackedTimeSolver

A module that is part of StateSpaceEcon package.
Contains methods for solving the dynamic system of equations
for the model and running simulations.
"""
module StackedTimeSolver

using SparseArrays
using LinearAlgebra
using ForwardDiff
using DiffResults
using Printf

using ModelBaseEcon
using TimeSeriesEcon

using ..Plans

include("stackedtime/abstract.jl")
include("stackedtime/misc.jl")
include("stackedtime/solverdata.jl")
include("stackedtime/simulate.jl")

end # module

using .StackedTimeSolver
export FCType, fcgiven, fclevel, fcslope, fcrate, fcnatural
export simulate
export seriesoverlay, dictoverlay
export dict2array, array2dict, array2data


text:
"""
    mutable struct Deterministic <: AbstractDeterministic
        name::String
        data::Union{
            SortedDict{Dates.DateTime, Vector{CONSTANT}},
            SortedDict{Dates.DateTime, Vector{POLYNOMIAL}},
            SortedDict{Dates.DateTime, Vector{PWL}},
        }
        resolution::Dates.Period
        scaling_factor_multiplier::Union{Nothing, Function}
        internal::InfrastructureSystemsInternal
    end

A deterministic forecast for a particular data field in a Component.

# Arguments
- `name::String`: user-defined name
- `data::Union{SortedDict{Dates.DateTime, Vector{CONSTANT}}, SortedDict{Dates.DateTime, Vector{POLYNOMIAL}}, SortedDict{Dates.DateTime, Vector{PWL}}}`: timestamp - scalingfactor
- `resolution::Dates.Period`: forecast resolution
- `scaling_factor_multiplier::Union{Nothing, Function}`: Applicable when the time series
  data are scaling factors. Called on the associated component to convert the values.
- `internal::InfrastructureSystemsInternal`
"""
mutable struct Deterministic <: AbstractDeterministic
    "user-defined name"
    name::String
    "timestamp - scalingfactor"
    data::Union{
        SortedDict{Dates.DateTime, Vector{CONSTANT}},
        SortedDict{Dates.DateTime, Vector{POLYNOMIAL}},
        SortedDict{Dates.DateTime, Vector{PWL}},
    }
    "forecast resolution"
    resolution::Dates.Period
    "Applicable when the time series data are scaling factors. Called on the associated component to convert the values."
    scaling_factor_multiplier::Union{Nothing, Function}
    internal::InfrastructureSystemsInternal
end

function Deterministic(;
    name,
    data,
    resolution,
    scaling_factor_multiplier = nothing,
    normalization_factor = 1.0,
    internal = InfrastructureSystemsInternal(),
)
    data = handle_normalization_factor(convert_data(data), normalization_factor)
    return Deterministic(name, data, resolution, scaling_factor_multiplier, internal)
end

function Deterministic(
    name::AbstractString,
    data::AbstractDict,
    resolution::Dates.Period;
    normalization_factor::NormalizationFactor = 1.0,
    scaling_factor_multiplier::Union{Nothing, Function} = nothing,
)
    return Deterministic(
        name = name,
        data = data,
        resolution = resolution,
        scaling_factor_multiplier = scaling_factor_multiplier,
        internal = InfrastructureSystemsInternal(),
    )
end

"""
Construct Deterministic from a Dict of TimeArrays.

# Arguments
- `name::AbstractString`: user-defined name
- `input_data::AbstractDict{Dates.DateTime, TimeSeries.TimeArray}`: time series data.
- `normalization_factor::NormalizationFactor = 1.0`: optional normalization factor to apply
  to each data entry
- `scaling_factor_multiplier::Union{Nothing, Function} = nothing`: If the data are scaling
  factors then this function will be called on the component and applied to the data when
  [`get_time_series_array`](@ref) is called.
- `timestamp = :timestamp`: If the values are DataFrames is passed then this must be the
  column name that contains timestamps.
"""
function Deterministic(
    name::AbstractString,
    input_data::AbstractDict{Dates.DateTime, <:TimeSeries.TimeArray};
    normalization_factor::NormalizationFactor = 1.0,
    scaling_factor_multiplier::Union{Nothing, Function} = nothing,
)
    data_type = eltype(TimeSeries.values(first(values(input_data))))
    data = SortedDict{Dates.DateTime, Vector{data_type}}()
    resolution =
        TimeSeries.timestamp(first(values(input_data)))[2] -
        TimeSeries.timestamp(first(values(input_data)))[1]
    for (k, v) in input_data
        if length(size(v)) > 1
            throw(ArgumentError("TimeArray with timestamp $k has more than one column)"))
        end
        data[k] = TimeSeries.values(v)
    end

    return Deterministic(
        name = name,
        data = data,
        resolution = resolution,
        normalization_factor = normalization_factor,
        scaling_factor_multiplier = scaling_factor_multiplier,
    )
end

"""
Construct Deterministic from a CSV file. The first column must be a timestamp in
DateTime format and the columns the values in the forecast window.

# Arguments
- `name::AbstractString`: user-defined name
- `filename::AbstractString`: name of CSV file containing data
- `component::InfrastructureSystemsComponent`: component associated with the data
- `normalization_factor::NormalizationFactor = 1.0`: optional normalization factor to apply
  to each data entry
- `scaling_factor_multiplier::Union{Nothing, Function} = nothing`: If the data are scaling
  factors then this function will be called on the component and applied to the data when
  [`get_time_series_array`](@ref) is called.
"""
function Deterministic(
    name::AbstractString,
    filename::AbstractString,
    component::InfrastructureSystemsComponent,
    resolution::Dates.Period;
    normalization_factor::NormalizationFactor = 1.0,
    scaling_factor_multiplier::Union{Nothing, Function} = nothing,
)
    component_name = get_name(component)
    raw_data = read_time_series(Deterministic, filename, component_name)
    return Deterministic(
        name,
        raw_data,
        resolution;
        normalization_factor = normalization_factor,
        scaling_factor_multiplier = scaling_factor_multiplier,
    )
end

"""
Construct Deterministic from RawTimeSeries.
"""
function Deterministic(
    name::AbstractString,
    series_data::RawTimeSeries,
    resolution::Dates.Period;
    normalization_factor::NormalizationFactor = 1.0,
    scaling_factor_multiplier::Union{Nothing, Function} = nothing,
)
    return Deterministic(
        name = name,
        data = series_data.data,
        resolution = resolution,
        normalization_factor = normalization_factor,
        scaling_factor_multiplier = scaling_factor_multiplier,
    )
end

function Deterministic(ts_metadata::DeterministicMetadata, data::SortedDict)
    return Deterministic(
        name = get_name(ts_metadata),
        resolution = get_resolution(ts_metadata),
        data = data,
        scaling_factor_multiplier = get_scaling_factor_multiplier(ts_metadata),
        internal = InfrastructureSystemsInternal(get_time_series_uuid(ts_metadata)),
    )
end

function Deterministic(info::TimeSeriesParsedInfo)
    return Deterministic(
        info.name,
        info.data,
        info.resolution;
        normalization_factor = info.normalization_factor,
        scaling_factor_multiplier = info.scaling_factor_multiplier,
    )
end

"""
Construct a new Deterministic from an existing instance and a subset of data.
"""
function Deterministic(forecast::Deterministic, data)
    vals = Dict{Symbol, Any}()
    for (fname, ftype) in zip(fieldnames(Deterministic), fieldtypes(Deterministic))
        if ftype <: SortedDict
            val = data
        elseif ftype <: InfrastructureSystemsInternal
            # Need to create a new UUID.
            val = InfrastructureSystemsInternal()
        else
            val = getfield(forecast, fname)
        end

        vals[fname] = val
    end

    return Deterministic(; vals...)
end

convert_data(data::AbstractDict{Dates.DateTime, Vector{T}}) where {T} =
    SortedDict{Dates.DateTime, Vector{CONSTANT}}(data...)
convert_data(data::AbstractDict{Dates.DateTime, Vector{T}}) where {T <: Tuple} =
    SortedDict{Dates.DateTime, Vector{POLYNOMIAL}}(data...)
convert_data(data::AbstractDict{Dates.DateTime, Vector{Vector{T}}}) where {T <: Tuple} =
    SortedDict{Dates.DateTime, Vector{PWL}}(data...)
convert_data(
    data::Union{
        SortedDict{Dates.DateTime, Vector{CONSTANT}},
        SortedDict{Dates.DateTime, Vector{POLYNOMIAL}},
        SortedDict{Dates.DateTime, Vector{PWL}},
    },
) = data

# Workaround for a bug/limitation in SortedDict. If a user tries to construct
# SortedDict(i => ones(2) for i in 1:2)
# it won't discern the types and will return SortedDict{Any,Any,Base.Order.ForwardOrdering}
# https://github.com/JuliaCollections/DataStructures.jl/issues/239
# This will only work for the most common use case of Vector{CONSTANT}.
# For other types the user will need to create SortedDict with explicit key-value types.
convert_data(data::AbstractDict) = SortedDict{Dates.DateTime, Vector{CONSTANT}}(data...)

function get_array_for_hdf(forecast::Deterministic)
    return transform_array_for_hdf(forecast.data)
end

"""Get [`Deterministic`](@ref) `name`."""
get_name(value::Deterministic) = value.name
"""Get [`Deterministic`](@ref) `data`."""
get_data(value::Deterministic) = value.data
"""Get [`Deterministic`](@ref) `resolution`."""
get_resolution(value::Deterministic) = value.resolution
"""Get [`Deterministic`](@ref) `scaling_factor_multiplier`."""
get_scaling_factor_multiplier(value::Deterministic) = value.scaling_factor_multiplier
"""Get [`Deterministic`](@ref) `internal`."""
get_internal(value::Deterministic) = value.internal

"""Set [`Deterministic`](@ref) `name`."""
set_name!(value::Deterministic, val) = value.name = val
"""Set [`Deterministic`](@ref) `data`."""
set_data!(value::Deterministic, val) = value.data = val
"""Set [`Deterministic`](@ref) `resolution`."""
set_resolution!(value::Deterministic, val) = value.resolution = val
"""Set [`Deterministic`](@ref) `scaling_factor_multiplier`."""
set_scaling_factor_multiplier!(value::Deterministic, val) =
    value.scaling_factor_multiplier = val
"""Set [`Deterministic`](@ref) `internal`."""
set_internal!(value::Deterministic, val) = value.internal = val

eltype_data(forecast::Deterministic) = eltype_data_common(forecast)
get_count(forecast::Deterministic) = get_count_common(forecast)
get_horizon(forecast::Deterministic) = get_horizon_common(forecast)
get_initial_times(forecast::Deterministic) = get_initial_times_common(forecast)
get_initial_timestamp(forecast::Deterministic) = get_initial_timestamp_common(forecast)
get_interval(forecast::Deterministic) = get_interval_common(forecast)
iterate_windows(forecast::Deterministic) = iterate_windows_common(forecast)
get_window(f::Deterministic, initial_time; len = nothing) =
    get_window_common(f, initial_time; len = len)

function make_time_array(forecast::Deterministic)
    # Artificial limitation to reduce scope.
    @assert_op get_count(forecast) == 1
    timestamps = range(
        get_initial_timestamp(forecast);
        step = get_resolution(forecast),
        length = get_horizon(forecast),
    )
    data = first(values(get_data(forecast)))
    return TimeSeries.TimeArray(timestamps, data)
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export daq_graduatedbar

"""
    daq_graduatedbar(;kwargs...)

A GraduatedBar component.
A graduated bar component that displays
a value within some range as a
percentage.
Keyword arguments:
- `id` (String; optional): The ID used to identify this compnent in Dash callbacks
- `className` (String; optional): Class to apply to the root component element.
- `color` (optional): Color configuration for the graduated bar's
progress blocks.. color has the following type: String | lists containing elements 'default', 'gradient', 'ranges'.
Those elements have the following types:
  - `default` (String; optional): Fallback color to use when color.ranges
has gaps.
  - `gradient` (Bool; optional): Display ranges as a gradient between given colors.
Requires color.ranges to be contiguous along
the entirety of the graduated bar's range of values.
  - `ranges` (optional): Define multiple color ranges on the graduated bar's track.
The key determines the color of the range and
the value is the start,end of the range itself.. ranges has the following type: lists containing elements 'color'.
Those elements have the following types:
  - `color` (Array of Reals; optional)
- `label` (optional): Description to be displayed alongside the control. To control styling, pass an object with label and style properties.. label has the following type: String | lists containing elements 'style', 'label'.
Those elements have the following types:
  - `style` (Dict; optional)
  - `label` (String; optional)
- `labelPosition` (a value equal to: 'top', 'bottom'; optional): Where the component label is positioned.
- `max` (Real; optional): The maximum value of the graduated bar
- `min` (Real; optional): The minimum value of the graduated bar
- `showCurrentValue` (Bool; optional): If true, the current percentage
 of the bar will be displayed
- `size` (Real; optional): The size (length) of the graduated bar in pixels
- `step` (Real; optional): Value by which progress blocks appear
- `style` (Dict; optional): Style to apply to the root component element.
- `theme` (Dict; optional): Theme configuration to be set by a ThemeProvider
- `value` (Real; optional): The value of the graduated bar
- `vertical` (Bool; optional): If true, will display bar vertically instead of horizontally
"""
function daq_graduatedbar(; kwargs...)
        available_props = Symbol[:id, :className, :color, :label, :labelPosition, :max, :min, :showCurrentValue, :size, :step, :style, :theme, :value, :vertical]
        wild_props = Symbol[]
        return Component("daq_graduatedbar", "GraduatedBar", "dash_daq", available_props, wild_props; kwargs...)
end


text:
# ---
# title: 30. Substring with Concatenation of All Words
# id: problem30
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Hash Table, Two Pointers, String
# link: <https://leetcode.com/problems/substring-with-concatenation-of-all-words/description/>
# hidden: true
# ---
#
# You are given a string `s` and an array of strings `words` of **the same
# length**. Return all starting indices of substring(s) in `s` that is a
# concatenation of each word in `words` **exactly once** , **in any order** ,
# and **without any intervening characters**.
#
# You can return the answer in **any order**.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "barfoothefoobarman", words = ["foo","bar"]
#     Output: [0,9]
#     Explanation: Substrings starting at index 0 and 9 are "barfoo" and "foobar" respectively.
#     The output order does not matter, returning [9,0] is fine too.
#
#
# **Example 2:**
#
#
#
#     Input: s = "wordgoodgoodgoodbestword", words = ["word","good","best","word"]
#     Output: []
#
#
# **Example 3:**
#
#
#
#     Input: s = "barfoofoobarthefoobarman", words = ["bar","foo","the"]
#     Output: [6,9,12]
#
#
#
#
# **Constraints:**
#
#   * `1 <= s.length <= 104`
#   * `s` consists of lower-case English letters.
#   * `1 <= words.length <= 5000`
#   * `1 <= words[i].length <= 30`
#   * `words[i]` consists of lower-case English letters.
#
#
## @lc code=start
using LeetCode

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


text:
# Autogenerated wrapper script for Libbpf_jll for i686-linux-musl
export libbpf

using Elfutils_jll
JLLWrappers.@generate_wrapper_header("Libbpf")
JLLWrappers.@declare_library_product(libbpf, "libbpf.so.0")
function __init__()
    JLLWrappers.@generate_init_header(Elfutils_jll)
    JLLWrappers.@init_library_product(
        libbpf,
        "lib/libbpf.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
struct LexicalException <: Exception
  msg::String
end

@enum TokenClass begin
  identifier
  times
  plus
  minus
  divide
  modulo
  identity_op
  equals
  not_equal
  less_than
  greater_than
  less_than_or_equal
  greater_than_or_equal
  open_paren
  close_paren
  open_sqr_bracket
  close_sqr_bracket
  rng
  semicolon
  colon
  comma
  dot
  doublequote
  int_literal
  string_literal
  real_literal
  assign
  kw_program
  kw_begin
  kw_end
  kw_procedure
  kw_function
  kw_return
  kw_var
  kw_while
  kw_for
  kw_if
  kw_then
  kw_else
  kw_in
  kw_do
  kw_read
  kw_writeln
  kw_array
  kw_of
  kw_size
  kw_and
  kw_or
  kw_not
  kw_int
  kw_real
  kw_string
  kw_boolean
  kw_nothing
  kw_assert
  kw_true
  kw_false
  eoi
  undefined
end

whitespace = [' ', '\t', '\n']
end_of_input_symbol = '$'
delimiters = union(whitespace, end_of_input_symbol)
ident_or_kw_initial = union('a':'z', 'A':'Z')
ident_or_kw_body = union(ident_or_kw_initial, '0':'9', '_')
keywords = Dict([
  "program" => kw_program,
  "begin" => kw_begin,
  "end" => kw_end,
  "procedure" => kw_procedure,
  "function" => kw_function,
  "return" => kw_return,
  "var" => kw_var,
  "while" => kw_while,
  "for" => kw_for,
  "if" => kw_if,
  "then" => kw_then,
  "else" => kw_else,
  "in" => kw_in,
  "do" => kw_do,
  "read" => kw_read,
  "writeln" => kw_writeln,
  "array" => kw_array,
  "of" => kw_of,
  "size" => kw_size,
  "int" => kw_int,
  "real" => kw_real,
  "string" => kw_string,
  "boolean" => kw_boolean,
  "assert" => kw_assert,
  "true" => kw_true,
  "false" => kw_false,
  "and" => kw_and,
  "or" => kw_or,
  "not" => kw_not
  ])
symbol_initials = ['*', '+', '-', '/', '(', ')', '[', ']', '.', ',', ';', ':', '!', '=', '<', '>', '&', '%']
digits = '0':'9'
unary_ops = Dict(
#   log_not => '!'
)
binary_ops = Dict(
  times => '*',
  plus => '+',
  minus => '-',
  divide => '/',
  equals => '=',
  less_than => '<',
  modulo => '%'
)
syntactic_symbols = Dict(
  open_paren => '(',
  close_paren => ')',
  open_sqr_bracket => '[',
  close_sqr_bracket => ']',
  colon => ':',
  semicolon => ';',
  comma => ',',
  dot => '.',
  doublequote => '"'
)
symbol_to_character = union(unary_ops, binary_ops, syntactic_symbols)
character_to_symbol = Dict([(sym => op) for (op, sym) ∈ symbol_to_character])

mutable struct Token
  class::TokenClass
  lexeme::String
  line::Int
end

Token(class::TokenClass, lexeme::String) = Token(class, lexeme, 0)

function getfirst(predicate, iterable)
  for item in iterable
    predicate(item) && return item
  end
  return nothing
end

"""
The main function of the scanner. The "next" index points at the
character to scan next. The outer while loop skips whitespace, keeps
track of the line number and the comment nesting level, calls the
getToken function, stores the resulting tokens and increments the index.
"""
function scan_input(input::AbstractString, next = 1)
  input *= end_of_input_symbol
  tokens = Vector{Token}()
  lineNumber = 1
  while next <= length(input)
    commentNesting = 0
    while (input[next] in whitespace
      || commentNesting > 0
      || (input[next] == '{' && input[next+1] == '*'))
      # || (input[next] == '/' && input[next+1] == '*'))
      if input[next] == '\n'
        lineNumber += 1
        next += 1
      elseif input[next] == '{' && input[next+1] == '*'
        commentNesting += 1
        next += 2
      elseif input[next] == '*' && input[next+1] == '}'
        commentNesting -= 1
        next += 2
      else
        next += 1
      end
    end

    if next <= length(input)
      token, next = getToken(input, next, lineNumber)
      token.line = lineNumber
      push!(tokens, token)
    end
  end
  return tokens
end

"""
Looks at the first character of a potential lexeme and
calls the corresponding function to produce a token.
"""
function getToken(input, next, lineNumber)
  c = input[next]
  if c ∈ ident_or_kw_initial
    return getIdentOrKw(input, next, lineNumber)
  end
  if c ∈ symbol_initials
    return getOperator(input, next, lineNumber)
  end
  if c ∈ digits
    return getNumber(input, next, lineNumber)
  end
  if c == '"'
    return getString(input, next, lineNumber)
  end
  if c == '$'
    return Token(eoi, string(end_of_input_symbol)), next+1
  end
  throw(LexicalException(
    "Characted $(c) on or near line $(lineNumber) is not part of a legal token."))
end

function getIdentOrKw(input, next, lineNumber)
  initial = next
  while input[next] ∈ ident_or_kw_body next += 1 end
  str = lowercase(input[initial:next-1])
  token = str ∈ keys(keywords) ? Token(keywords[str], str) : Token(identifier, str)
  return token, next
end

function getOperator(input, next, lineNumber)
  # if input[next] == '.'
  #   input[next+1] != '.' && throw(LexicalException(
  #     "A single dot (on line $(lineNumber)) is not a valid token. Did you mean '..'?"))
  #   return Token(rng, ".."), next+2
  # end
  if input[next] == ':'
    input[next+1] == '=' && return Token(assign, ":="), next+2
    return Token(colon, ":"), next+1
  end
  if input[next] == '<'
    input[next+1] == '=' && return Token(less_than_or_equal, "<="), next+2
    input[next+1] == '>' && return Token(not_equal, "<>"), next+2
    return Token(less_than, "<"), next+1
  end
  if input[next] == '>'
    input[next+1] == '=' && return Token(greater_than_or_equal, ">="), next+2
    return Token(greater_than, ">"), next+1
  end
  tokenClass = character_to_symbol[input[next]]
  return Token(tokenClass, string(input[next])), next+1
end

function getNumber(input, next, lineNumber)
  initial = next
  integer_block, next = get_number_block(input, next)
  if input[next] == '.'
    next += 1
    decimal_block, next = get_number_block(input, next)
    if input[next] ∈ ['e','E']
      next += 1
      sign = '+'
      if input[next] ∈ ['+','-']
        sign = input[next]
        next += 1
      end
      exp_block, next = get_number_block(input, next)
      return Token(real_literal, integer_block * '.' * decimal_block * 'e' * sign * exp_block), next
    end
    return Token(real_literal, integer_block * '.' * decimal_block), next
  end
  return Token(int_literal, integer_block), next
end

function get_number_block(input, next)
  initial = next
  while input[next] ∈ digits next += 1 end
  return input[initial:next-1], next
end

# function getNumber(input, next, lineNumber)
#   initial = next
#   while input[next] ∈ digits next += 1 end
#   return Token(int_literal, input[initial:next-1]), next
# end

function getString(input, next, lineNumber)
  initial = next
  next += 1
  str = Vector{Char}()
  while (input[next] != '"' ||
    (input[next] == '"' && input[next-1] == '\\'))
    if input[next] == '\\'
      next += 1
      if input[next] == 'n'
        push!(str, '\n')
      end
      if input[next] == 't'
        push!(str, '\t')
      end
      if input[next] == '"'
        push!(str, '"')
      end
      if input[next] == '\\'
        push!(str, '\\')
      end
      next += 1
    else
      push!(str, input[next])
      next += 1
    end
    next >= length(input) && throw(LexicalException("Reached the end of the program while
      scanning a string literal. Did you forget the closing quote?"))
  end
  return Token(string_literal, join(str)), next+1
end

max_stars_repo_name: juhamyllari/minipascal-julia

max_issues_repo_name: juhamyllari/minipascal-julia

max_forks_repo_name: juhamyllari/minipascal-julia


text:
function equilibrium(
    q::Quadrature,
    ρ::T,
    u::VT,
    temperature::T,
) where {T <: Real, VT <: AbstractVector{T}}
    f = zeros(length(q.weights))
    hermite_based_equilibrium!(q, ρ, u, temperature, f)
    return f
end

function equilibrium!(
    q::Quadrature,
    ρ::T,
    u::VT,
    temperature::T,
    f,
) where {T <: Real, VT <: AbstractVector{T}}
    u_squared = 0.0
    for d in 1:dimension(q)
        u_squared += u[d] .^ 2
    end

    for idx in 1:length(q.weights)
        u_dot_xi = q.abscissae[1, idx] .* u[1] .+ q.abscissae[2, idx] .* u[2]

        f[idx] = _equilibrium(
            q,
            ρ,
            q.weights[idx],
            u_dot_xi,
            u_squared,
            temperature,
            q.abscissae[1, idx]^2 + q.abscissae[2, idx]^2,
        )
    end

    return
end

# Truncated upto order 2
# function _equilibrium(q::Quadrature, ρ::T, weight::T, u_dot_xi::T, u_squared::T, temperature::T, xi_squared::T) where { T <: Real }
function _equilibrium(
    q::Quadrature,
    ρ,
    weight,
    u_dot_xi,
    u_squared,
    temperature,
    xi_squared,
)
    cs = q.speed_of_sound_squared
    # a = 0.0
    a_H_0 = 1.0
    a_H_1 = cs * u_dot_xi

    # H_2_temperature = cs * ( temperature .- 1) .* (xi_squared * cs - dimension(q))
    H_2_temperature = 0.0 #zero(T)

    a_H_2 = cs^2 * (u_dot_xi .* u_dot_xi) .+ H_2_temperature .+ -cs * u_squared
    # return ρ .* weight .+ weight * (
    #     a_H_1 .+
    #     (1 / 2) * a_H_2
    # )
    return ρ .* weight .* (a_H_0 .+ a_H_1 .+ (1 / 2) * a_H_2)
end


text:
using LyceumMuJoCo, MuJoCo

m = jlModel("humanoid.xml")