text:
using Dates
using DelimitedFiles
using LinearAlgebra
using MultivariateStats
using NearestNeighbors # https://github.com/KristofferC/NearestNeighbors.jl
using Printf
using Statistics
using StatsBase

mutable struct PointCloud

    x::Vector{Float64}
    y::Vector{Float64}
    z::Vector{Float64}

    nx::Vector{Float64}
    ny::Vector{Float64}
    nz::Vector{Float64}
    planarity::Vector{Float64}

    no_points::Int64
    sel::Vector{Int64}

    function PointCloud(x, y, z)
        no_points = length(x)
        new(x,
            y,
            z,
            fill(NaN, no_points),
            fill(NaN, no_points),
            fill(NaN, no_points),
            fill(NaN, no_points),
            no_points,
            collect(1:no_points))
    end

end

function select_in_range!(pc::PointCloud, X::Array, max_range::Number)

    size(X)[2] == 3 || error(""""X" must have 3 columns""")
    max_range > 0 || error(""""max_range" must be > 0""")

    kdtree = KDTree(X')

    query_points = [pc.x[pc.sel]'; pc.y[pc.sel]'; pc.z[pc.sel]']

    _, distances = nn(kdtree, query_points)

    keep = [d <= max_range for d in distances]

    pc.sel = pc.sel[keep]

end

function select_n_points!(pc::PointCloud, n)

    no_selected_points = length(pc.sel)

    if no_selected_points > n
        idx = round.(Int, range(1, no_selected_points, length=n))
        pc.sel = pc.sel[idx]
    end

end

function estimate_normals!(pc::PointCloud, neighbors)

    kdtree = KDTree([pc.x'; pc.y'; pc.z'])
    query_points = [pc.x[pc.sel]'; pc.y[pc.sel]'; pc.z[pc.sel]']
    idxNN_all_qp, = knn(kdtree, query_points, neighbors, false)

    for (i, idxNN) in enumerate(idxNN_all_qp)
        selected_points = [pc.x[idxNN]'; pc.y[idxNN]'; pc.z[idxNN]']
        # P = fit(PCA, selected_points, pratio=1.0)
        # pc.nx[pc.sel[i]] = projection(P)[1,3]
        # pc.ny[pc.sel[i]] = projection(P)[2,3]
        # pc.nz[pc.sel[i]] = projection(P)[3,3]
        C = cov(selected_points, dims=2)
        F = eigen(C) # eigenvalues are in ascending order
        pc.nx[pc.sel[i]] = F.vectors[1,1]
        pc.ny[pc.sel[i]] = F.vectors[2,1]
        pc.nz[pc.sel[i]] = F.vectors[3,1]
        pc.planarity[pc.sel[i]] = (F.values[2]-F.values[1])/F.values[3];
    end

end

function transform!(pc, H)

    XInH = euler_coord_to_homogeneous_coord([pc.x pc.y pc.z])
    XOutH = transpose(H*XInH')
    XOut = homogeneous_coord_to_euler_coord(XOutH)

    pc.x = XOut[:,1]
    pc.y = XOut[:,2]
    pc.z = XOut[:,3]

    return pc

end

function matching!(pcmov::PointCloud, pcfix)

    kdtree = KDTree([pcmov.x'; pcmov.y'; pcmov.z'])
    query_points = [pcfix.x[pcfix.sel]'
                    pcfix.y[pcfix.sel]'
                    pcfix.z[pcfix.sel]']
    idxNN, = knn(kdtree, query_points, 1)
    pcmov.sel = vcat(idxNN...)

    dx = pcmov.x[pcmov.sel] - pcfix.x[pcfix.sel]
    dy = pcmov.y[pcmov.sel] - pcfix.y[pcfix.sel]
    dz = pcmov.z[pcmov.sel] - pcfix.z[pcfix.sel]

    nx = pcfix.nx[pcfix.sel]
    ny = pcfix.ny[pcfix.sel]
    nz = pcfix.nz[pcfix.sel]

    distances = [dx[i]*nx[i] + dy[i]*ny[i] + dz[i]*nz[i] for i in 1:length(pcmov.sel)]

    return distances

end

function reject!(pcmov::PointCloud, pcfix::PointCloud, min_planarity, distances)

    planarity = pcfix.planarity[pcfix.sel]

    med = median(distances)
    sigmad = mad(distances, normalize=true)

    keep_distance = [abs(d-med) <= 3*sigmad for d in distances]
    keep_planarity = [p > min_planarity for p in planarity]

    keep = keep_distance .& keep_planarity

    pcmov.sel = pcmov.sel[keep]
    pcfix.sel = pcfix.sel[keep]
    deleteat!(distances, .!keep)

    return nothing

end

function estimate_rigid_body_transformation(x_fix, y_fix, z_fix, nx_fix, ny_fix, nz_fix,
                                            x_mov, y_mov, z_mov)

    A = hcat(-z_mov.*ny_fix + y_mov.*nz_fix,
              z_mov.*nx_fix - x_mov.*nz_fix,
             -y_mov.*nx_fix + x_mov.*ny_fix,
             nx_fix,
             ny_fix,
             nz_fix)

    l = nx_fix.*(x_fix-x_mov) + ny_fix.*(y_fix-y_mov) + nz_fix.*(z_fix-z_mov)

    x = A\l

    residuals = A*x-l

    R = euler_angles_to_linearized_rotation_matrix(x[1], x[2], x[3])

    t = x[4:6]

    H = create_homogeneous_transformation_matrix(R, t)

    return H, residuals

end

function euler_angles_to_linearized_rotation_matrix(α1, α2, α3)

    dR = [  1 -α3  α2
           α3   1 -α1
          -α2  α1   1]

end

function create_homogeneous_transformation_matrix(R, t)

    H = [R          t
         zeros(1,3) 1]

end

function euler_coord_to_homogeneous_coord(XE)

    no_points = size(XE, 1)
    XH = [XE ones(no_points,1)]

end

function homogeneous_coord_to_euler_coord(XH)

    XE = XH[:,1:3]./XH[:,4]

end

function check_convergence_criteria(distances_new, distances_old, min_change)

    change(new, old) = abs((new-old)/old*100)

    change_of_mean = change(mean(distances_new), mean(distances_old))
    change_of_std = change(std(distances_new), std(distances_old))

    return change_of_mean < min_change && change_of_std < min_change ? true : false

end

function simpleicp(X_fix::Array, X_mov::Array;
                   correspondences::Integer=1000,
                   neighbors::Integer=10,
                   min_planarity::Number=0.3,
                   max_overlap_distance::Number=Inf,
                   min_change::Number=3,
                   max_iterations::Integer=100)

    size(X_fix)[2] == 3 || error(""""X_fix" must have 3 columns""")
    size(X_mov)[2] == 3 || error(""""X_mov" must have 3 columns""")
    correspondences >= 10 || error(""""correspondences" must be >= 10""")
    min_planarity >= 0 && min_planarity < 1 || error(""""min_planarity" must be >= 0 and < 1""")
    neighbors >= 2 || error(""""neighbors" must be >= 2""")
    min_change > 0 || error(""""min_change" must be > 0""")
    max_iterations > 0 || error(""""max_iterations" must be > 0""")

    dt = @elapsed begin
        @info "Create point cloud objects ..."
        pcfix = PointCloud(X_fix[:,1], X_fix[:,2], X_fix[:,3])
        pcmov = PointCloud(X_mov[:,1], X_mov[:,2], X_mov[:,3])

        if isfinite(max_overlap_distance)
            @info "Consider partial overlap of point clouds ..."
            select_in_range!(pcfix, X_mov, max_overlap_distance)
            if length(pcfix.sel) == 0
                error(@sprintf("Point clouds do not overlap within max_overlap_distance = %.5f! ",
                max_overlap_distance) * "Consider increasing the value of max_overlap_distance.")
            end
        end

        @info "Select points for correspondences in fixed point cloud ..."
        select_n_points!(pcfix, correspondences)
        sel_orig = pcfix.sel

        @info "Estimate normals of selected points ..."
        estimate_normals!(pcfix, neighbors)

        H = Matrix{Float64}(I,4,4)
        residual_distances = Any[]

        @info "Start iterations ..."
        for i in 1:max_iterations

            initial_distances = matching!(pcmov, pcfix)

            reject!(pcmov, pcfix, min_planarity, initial_distances)

            dH, residuals = estimate_rigid_body_transformation(
                pcfix.x[pcfix.sel], pcfix.y[pcfix.sel], pcfix.z[pcfix.sel],
                pcfix.nx[pcfix.sel], pcfix.ny[pcfix.sel], pcfix.nz[pcfix.sel],
                pcmov.x[pcmov.sel], pcmov.y[pcmov.sel], pcmov.z[pcmov.sel])

            push!(residual_distances, residuals)

            transform!(pcmov, dH)

            H = dH*H
            pcfix.sel = sel_orig

            if i > 1
                if check_convergence_criteria(residual_distances[i], residual_distances[i-1],
                                              min_change)
                    @info "Convergence criteria fulfilled -> stop iteration!"
                    break
                end
            end

            if i == 1
                @info @sprintf(" %9s | %15s | %15s | %15s", "Iteration", "correspondences",
                               "mean(residuals)", "std(residuals)")
                @info @sprintf(" %9s | %15d | %15.4f | %15.4f", "orig:0", length(initial_distances),
                               mean(initial_distances), std(initial_distances))
            end
            @info @sprintf(" %9d | %15d | %15.4f | %15.4f", i, length(residual_distances[i]),
                           mean(residual_distances[i]), std(residual_distances[i]))

        end

    end
    @info "Estimated transformation matrix H:\n" *
    @sprintf("[%12.6f %12.6f %12.6f %12.6f]\n", H[1,1], H[1,2], H[1,3], H[1,4]) *
    @sprintf("[%12.6f %12.6f %12.6f %12.6f]\n", H[2,1], H[2,2], H[2,3], H[2,4]) *
    @sprintf("[%12.6f %12.6f %12.6f %12.6f]\n", H[3,1], H[3,2], H[3,3], H[3,4]) *
    @sprintf("[%12.6f %12.6f %12.6f %12.6f]\n", H[4,1], H[4,2], H[4,3], H[4,4])
    @info "Finished in " * @sprintf("%.3f", dt) * " seconds!"

    X_mov_transformed = [pcmov.x pcmov.y pcmov.z]

    return H, X_mov_transformed

end

max_stars_repo_name: mitjap/simpleICP

max_issues_repo_name: mitjap/simpleICP

max_forks_repo_name: mitjap/simpleICP


text:
"""
For integers \$a\$, \$b\$, and \$n\$ efficiently compute
\$a^b \\bmod n\$.
"""
@fastmath function power_mod(a, b, n)
  r = 1
  while b > 0
    if isodd(b)
      r = r * a % n
    end
    b = div(b, 2)
    a = a * a % n
  end
  return r
end

"""
For integer \$d\$, \$k\$, and \$j\$ compute
\\[ ( 16^{d-k} \\bmod 8k + j ) / (8k + j). \\]
"""
ratio(d, k, j) = power_mod(16, d - k, 8k + j) / (8k + j)

"""
Fractional part of a real number \$x\$.
"""
frac_part(x) = x - floor(x)

"""
For integer \$d\$ and \$j\$ computes the fractional part of
\\[ \\sum_{k=0}^d \\frac{16^{d-k} \\bmod 8k+j}{8k+j} + \\sum_{k=d+1}^\\infty \\frac{1}{16^{k-d} (8k+j)}. \\]
"""
function sigma(d, j)
  r = 0.
  @fastmath @simd for k in 0:d
    r += ratio(d, k, j)
    if r > 1.0
      r = frac_part(r)
    end
  end

  kk  = 8. * (d + 1) + j
  u   = 1 / 16
  rat = u / kk

  @fastmath while rat > eps(Float64)
    r   += rat
    kk  += 8.
    u   /= 16.
    rat  = u / kk
  end

  return frac_part(r)
end

"""
Compute \$d+1\$th hexadecimal digit of \$\\pi\$.
"""
function digpi(d)
  r = 4*sigma(d, 1) - 2*sigma(d, 4) - sigma(d, 5) - sigma(d, 6)
  return floor(Int, 16*frac_part(r))
end

ext: jl


text:
using Test
using JuMP
using Polyhedra
@testset "PolyhedraToLPBridge{$T}" for T in [Float64, Int, BigInt, Rational{BigInt}]
    h = HalfSpace(ones(T, 2), one(T)) ∩ HalfSpace(-ones(T, 2), -one(T))
    set = Polyhedra.PolyhedraOptSet(h)
    S = typeof(set)
    Fs = [MOI.VectorOfVariables, MOI.VectorAffineFunction{T}, MOI.VectorQuadraticFunction{T}]
    mock = MOI.Utilities.MockOptimizer(MOI.Utilities.Model{T}())
    for F in Fs
        @test !MOI.supports_constraint(mock, F, S)
    end

    @testset "LazyBridgeOptimizer" begin
        bridged = MOI.Bridges.LazyBridgeOptimizer(mock)
        for F in Fs
            @test !MOI.supports_constraint(bridged, F, S)
        end
        MOI.Bridges.add_bridge(bridged, Polyhedra.PolyhedraToLPBridge{T})
        for F in Fs
            @test MOI.supports_constraint(bridged, F, S)
        end
        x = MOI.add_constrained_variables(bridged, set)
     end

    @testset "LazyBridgeOptimizer" begin
        bridged = MOI.Bridges.full_bridge_optimizer(mock, T)
        for F in Fs
            @test !MOI.supports_constraint(bridged, F, S)
        end
        MOI.Bridges.add_bridge(bridged, Polyhedra.PolyhedraToLPBridge{T})
        for F in Fs
            @test MOI.supports_constraint(bridged, F, S)
        end
        x = MOI.add_constrained_variables(bridged, set)
     end
end


text:
struct PostProcessor{T, R <: BEMResult{T}}
    bem     ::R
    elements::TriangleVector{T}
end

@inline function PostProcessor(bem::R) where {T, R <: BEMResult{T}}
    PostProcessor(bem, TriangleVector(bem.model.elements))
end

function φ⃰Ω(Ξ::Vector{Vector{T}}, proc::PostProcessor{T}) where T
    Ξ = PositionVector(Ξ)
    Ṽ = LaplacePotentialMatrix{SingleLayer}(Ξ, proc.elements)
    W = LaplacePotentialMatrix{DoubleLayer}(Ξ, proc.elements)

    (NESSie.potprefactor(T) / T(4π)) .* (Ṽ * proc.bem.q .- W * proc.bem.u)
end

function φΩ(Ξ::Vector{Vector{T}}, proc::PostProcessor{T}) where T
    εΩ = proc.bem.model.params.εΩ
    φm  = φmol(Ξ, proc.bem.model.charges)

    φ⃰Ω(Ξ, proc) .+ (NESSie.potprefactor(T) / εΩ .* φm)
end

@inline φΩ(
    Ξ  ::Vector{Vector{T}},
    bem::R
) where {T, R <: BEMResult{T}} = φΩ(Ξ, PostProcessor(bem))

function rfenergy(proc::PostProcessor{T}) where T
    qposs = [c.pos for c in proc.bem.model.charges]
    qvals = [c.val for c in proc.bem.model.charges]

    φ⃰Ω(qposs, proc) ⋅ qvals * T(48.242647748524156)
end

@inline rfenergy(
    bem::R
) where {T, R <: BEMResult{T}} = rfenergy(PostProcessor(bem))

function φΣ(
    Ξ::Vector{Vector{T}},
    proc::PostProcessor{T, LocalBEMResult{T, Triangle{T}}}
) where T
    εΩ = proc.bem.model.params.εΩ
    εΣ = proc.bem.model.params.εΣ
    Ξ  = PositionVector(Ξ)
    Ṽ  = LaplacePotentialMatrix{SingleLayer}(Ξ, proc.elements)
    W  = LaplacePotentialMatrix{DoubleLayer}(Ξ, proc.elements)

    (NESSie.potprefactor(T) / T(4π)) .*
        (W * (proc.bem.u .+ proc.bem.umol) .-
        ((εΩ / εΣ) .* (Ṽ * (proc.bem.q .+ proc.bem.qmol))))
end

function φΣ(
    Ξ::Vector{Vector{T}},
    proc::PostProcessor{T, NonlocalBEMResult{T, Triangle{T}}}
) where T
    εΩ   = proc.bem.model.params.εΩ
    εΣ   = proc.bem.model.params.εΣ
    ε∞   = proc.bem.model.params.ε∞
    yuk  = yukawa(proc.bem.model.params)
    Ξ    = PositionVector(Ξ)
    Ṽ    = LaplacePotentialMatrix{SingleLayer}(Ξ, proc.elements)
    W    = LaplacePotentialMatrix{DoubleLayer}(Ξ, proc.elements)
    ṼʸmṼ = ReYukawaPotentialMatrix{SingleLayer}(Ξ, proc.elements, yuk)
    WʸmW = ReYukawaPotentialMatrix{DoubleLayer}(Ξ, proc.elements, yuk)

    (NESSie.potprefactor(T) / T(4π)) .* (
        (Ṽ * ((-εΩ/ε∞) .* (proc.bem.q .+ proc.bem.qmol)))
        .+ (W * (proc.bem.u .+ proc.bem.umol))
        .+ (ṼʸmṼ * ((εΩ * (1/εΣ - 1/ε∞)) .* (proc.bem.q .+ proc.bem.qmol)))
        .+ (WʸmW * (proc.bem.u .- ((ε∞/εΣ) * proc.bem.w) .+ ((1 - εΩ/εΣ) * proc.bem.umol)))
    )
end

@inline φΣ(
    Ξ  ::Vector{Vector{T}},
    bem::R
) where {T, R <: BEMResult{T}} = φΣ(Ξ, PostProcessor(bem))


text:
@testset "Height functions for elliptic curves" begin
  E = EllipticCurve([1, -1, 1, -2063758701246626370773726978, 32838647793306133075103747085833809114881])

  P = E([-30987785091199, 258909576181697016447])
  R = ArbField(100)
  @testset "Naive height" begin

    @test contains( R(3106461421344758392571240556046391//100000000000000000000000000000000) , @inferred naive_height(P))

  end

  @testset "Local heights" begin
    @test contains( R(-1888808896037477386833023078582084//1000000000000000000000000000000000) , @inferred local_height(P, 17))
    @test contains( R(-1098612288668109691395245236922525//1000000000000000000000000000000000) , @inferred local_height(P, 3))
  end

  @testset "Real height" begin
     @test contains( R(31812867857980432756826379798033//1000000000000000000000000000000) , @inferred real_height(P))
  end

  @testset "Canonical height" begin
    @test  contains( R(25860317067546190743868840740735//1000000000000000000000000000000) , @inferred canonical_height(P))
  end

end


text:
expected_coskew_tensor = reshape(expected_coskew, 3, 3, 3)
expected_cokurt_tensor = reshape(expected_cokurt, 3, 3, 3, 3)

@test tolerance(_cov(data, bias=1), cov(data, corrected=true))
@test tolerance(_cov(data), cov(data, corrected=false))

ltm = _cov(data, dense=false)
@test tolerance(ltm + ltm' - diagm(diag(ltm)), _cov(data, dense=true))

@test tolerance(coskew(data, flatten=true), expected_coskew)
@test tolerance(coskew(data), expected_coskew_tensor)
@test tolerance(coskew(data, flatten=true)[1,4], coskew(data, flatten=true, dense=false)[1,4])
@test tolerance(coskew(data, flatten=true, dense=false)[1,4], coskew(data, dense=false)[2,1,1])
@test tolerance(coskew(data, dense=false)[2,1,1], coskew(data)[2,1,1])

@test tolerance(cokurt(data, flatten=true), expected_cokurt)
@test tolerance(cokurt(data), expected_cokurt_tensor)
@test tolerance(cokurt(data, flatten=true)[1,10], cokurt(data, flatten=true, dense=false)[1,10])
@test tolerance(cokurt(data, flatten=true, dense=false)[1,10], cokurt(data, dense=false)[2,1,1,1])


text:
# using LinearAlgebra, Plots, JSON
# include("../src/SFFM.jl")
include("METools.jl")

include("stochasticApproximations.jl")
p = plot!()

# # define SFM
# T = [-2.0 2.0; 1.0 -1.0]#[-2.0 2.0 0; 1.0 -2.0 1; 1 1 -2]
# C = [1.0; -2.0]#; -1]
# fn(x) = [ones(size(x)) ones(size(x))]# ones(size(x))]
# model = SFFM.Model(;T=T,C=C,r=(r=fn,R=fn),Bounds=[0 10;-Inf Inf])
# N₋ = sum(C.<=0)
# N₊ = sum(C.>=0)
# NPhases = length(C)
#
# t = 4.0
# τ = SFFM.FixedTime(T=t)
# NSim = 100_000
# sims = SFFM.SimSFM(model=model,StoppingTime=τ,InitCondition=(φ=2*ones(Int,NSim),X=zeros(NSim)))

function MakeGlobalApprox(;NCells = 3,up, down,T,C,bkwd=false,D=[],NCounters=2)
    αup,Qup = up
    αdown,Qdown = down
    N₋ = sum(C.<=0)
    N₊ = sum(C.>=0)
    NPhases = length(C)
    NBases = length(αup)
    qup = -sum(Qup,dims=2)
    qdown = -sum(Qdown,dims=2)
    Q = zeros(NCells*NBases*NPhases*NCounters,NCells*NBases*NPhases*NCounters)
    for n in 1:NCells
        for i in 1:NPhases
            for c in 1:NCounters
                tempmat = zeros(NBases,NBases)
                for nn in 1:(c-1)
                    phase = mod(1+i+nn,2)+1
                    tempmat += T[phase,phase].*I(NBases)
                end
                idx = (1:NBases) .+ (n-1)*(NBases*NPhases*NCounters) .+ (i-1)*NBases .+ (c-1)*NBases*NPhases
                if C[i]>0
                    Q[idx,idx] = (Qup+tempmat)*abs(C[i])
                    if n<NCells
                        qup = -sum((Qup+tempmat),dims=2)
                        Q[idx,(1:NBases) .+ (n)*NBases*NPhases*NCounters .+ NBases*(i-1)] = qup*αup*abs(C[i])
                    end
                elseif C[i]<0
                    Q[idx,idx] = (Qdown+tempmat)*abs(C[i])
                    if n>1
                        qdown = -sum((Qdown+tempmat),dims=2)
                        Q[idx,(n-2)*NBases*NPhases*NCounters .+ (1:NBases) .+ NBases*(i-1)] = qdown*αdown*abs(C[i])
                    end
                end
            end
        end
    end
    T₋₋ = T[C.<=0,C.<=0]
    T₊₋ = T[C.>=0,:].*((C.<0)')
    T₋₊ = T[C.<=0,:].*((C.>0)')
    T₊₊ = T[C.>=0,C.>=0]
    negidx = kron(
        ones(NCells*NCounters),
        kron((C.<0)[:],ones(NBases))
        )
    inLower = -sum(Q,dims=2).*negidx
    outLower = [kron(1,kron(T₋₊,αup)) zeros(N₋,N₊+(NCells-1)*NPhases*NBases*NCounters+(NCounters-1)*NPhases*NBases)]
    inUpper = -sum(Q,dims=2).*(1 .- negidx)
    outUpper = [zeros(N₊,N₋+(NCells-1)*NPhases*NBases*NCounters) kron(1,kron(T₊₋,αdown)) zeros(N₊,(NCounters-1)*NPhases*NBases)]
    # display(inLower)
    # display(outLower)
    # display(inUpper)
    # display(outUpper)
    # display(NCells)
    # display(NPhases)
    # display(NBases)
    # display(N₊)
    # display(N₋)
    # display(kron(I(NCells),kron(T,I(NBases)))+Q)

    # idx = [1; [3:2:NBases 2:2:NBases]'[:]]
    Tdiag = diagm(diag(T))
    Toff = T-diagm(diag(T))
    T₊₋ = T.*((C.>0)*(C.<0)')
    T₋₊ = T.*((C.<0)*(C.>0)')
    Tchange = T₊₋ + T₋₊
    Tnochange = T-Tchange
    # πME = -αup*Qup^-1
    # μ = sum(πME)
    # πME = πME./μ
    I2 = I(NBases)#[:,idx]#diagm(πME[:])#repeat(πME,length(πME),1)#
    D = zeros(size(Q))
    for c in 1:NCells, i in 1:length(C), j in 1:length(C), n in 1:(NCounters)
        if n < NCounters
            idxi = (1:NBases) .+ (i-1)*NBases .+ (n-1)*NBases*NPhases .+ (c-1)*NBases*NPhases*NCounters
            idxj = (1:NBases) .+ (j-1)*NBases .+ (n)*NBases*NPhases .+ (c-1)*NBases*NPhases*NCounters
        else
            idxi = (1:NBases) .+ (i-1)*NBases .+ (n-1)*NBases*NPhases .+ (c-1)*NBases*NPhases*NCounters
            idxj = (1:NBases) .+ (j-1)*NBases .+ (n-1)*NBases*NPhases .+ (c-1)*NBases*NPhases*NCounters
        end
        # idxi = (n-1)*NBases .+ ((i-1)*NBases*NCounters) .+ (1:NBases) .+ ((c-1)*NBases*NPhases)
        # idxj = (n-1)*NBases .+ ((j-1)*NBases*NCounters) .+ (1:NBases) .+ ((c-1)*NBases*NPhases)
        if i!=j
            tempmat = zeros(NBases,NBases)
            for nn in 1:(n)
                phase = mod(1+i+nn,2)+1
                tempmat += T[phase,phase].*I(NBases)
            end
            if C[i]>0 && C[j]<0
                # display(Qup+T[i,i]*I(NBases)*plusI)
                # display(Qup)
                if n < NCounters
                    πtemp = -αup*(abs(C[i])*Qup+tempmat)^-1
                    πtemp = πtemp./sum(πtemp)
                    D[idxi,idxj] = T[i,j].*repeat(πtemp,NBases,1)
                else
                    πtemp = αup#*(-Qup^-1)
                    # πtemp = -αup*(abs(C[i])*Qup+tempmat)^-1
                    πtemp = πtemp./sum(πtemp)
                    D[idxi,idxj] = T[i,j].*repeat(πtemp,NBases,1)
                end
            elseif C[i]<0 && C[j]>0
                if n < NCounters
                    πtemp = -αdown*(abs(C[i])*Qdown+tempmat)^-1
                    πtemp = πtemp./sum(πtemp)
                    D[idxi,idxj] = T[i,j].*repeat(πtemp,NBases,1)
                else
                    πtemp = αdown#*(-Qdown)^-1
                    # πtemp = -αdown*(abs(C[i])*Qdown+tempmat)^-1
                    πtemp = πtemp./sum(πtemp)
                    D[idxi,idxj] = T[i,j].*repeat(πtemp,NBases,1)
                end
            else
                D[idxi,idxj] = T[i,j].*I(NBases)
            end
        else
                if n < NCounters
                    D[idxi,idxj.-NBases*NPhases] = T[i,i].*I(NBases)
                else
                    D[idxi,idxj] = T[i,i].*I(NBases)
                end
        end
    end
    B = [
        T₋₋ outLower;
        inLower D+Q inUpper;
        outUpper T₊₊;
    ]
    return Q, B
end

for NCounters in 2:2
    globalerrME = []
    globalerrDG = []
    for NBases in 1:2:21
        # NBases = 2
        Nodes = collect(0:Δ:10)
        NCells = length(Nodes)-1
        Erlang = MakeErlang(NBases,mean=Δ)#MakeME(CMEParams[NBases], mean = Δ)#
        display(Erlang.Q)
        Q, B = MakeGlobalApprox(
            NCells = NCells,
            up = Erlang,
            down = Erlang,
            T = T,
            C = C,
            NCounters = NCounters,
        )
        display(B)
        display(Q)

        DGMesh = SFFM.MakeMesh(model=model,NBases=1,Nodes=collect(Nodes[1]:Δ/NBases:Nodes[end]),Basis="lagrange")
        All = SFFM.MakeAll(model=model,mesh=DGMesh)

        initDist = zeros(1,size(All.B.B,1))
        initDist[1] = 1

        temp = initDist*exp(Matrix(All.B.B)*t)#SFFM.EulerDG(D=All.B.B,y=t,x0=initDist)#
        DGdist_t = SFFM.Coeffs2Dist(model=model,mesh=DGMesh,Coeffs=temp,type="probability")

        initDist = zeros(1,size(B,1))
        initDist[1] = 1
        dist_t = initDist*exp(B*t)#SFFM.EulerDG(D=B,y=t,x0=initDist)#
        pm_t = dist_t[[1:N₋;(end-N₊+1):end]]
        dist_t = sum(reshape(dist_t[N₋+1:end-N₊],NBases,NCounters*NPhases*NCells),dims=1)
        dist_t = sum(reshape(dist_t,NPhases,NCounters,NCells),dims=2)
        mesh = SFFM.MakeMesh(model=model,NBases=NBases,Nodes=Nodes,Basis="lagrange")
        simDist = SFFM.Sims2Dist(model=model,mesh=mesh,sims=sims,type="probability")

        localerrME = sum(abs.(pm_t-simDist.pm))
        localerrDG = sum(abs.(DGdist_t.pm-simDist.pm))

        # q = plot(layout = (2,1))
        for n in 1:NCells
            x = simDist.x[n]
            for i in 1:NPhases
                if n==1
                    label1 = "DG"
                    label2 = "ME"
                else
                    label1 = false
                    label2 = false
                end
                yvalsME = [sum(dist_t[i,1,n])]
                yvalsDG = [sum(DGdist_t.distribution[:,(1:NBases).+NBases*(n-1),i])]
                # scatter!([x],yvalsDG,label=label1,subplot=i,color=:blue,markershape=:rtriangle)
                # scatter!([x],yvalsME,label=label2,subplot=i,color=:black,markershape=:ltriangle)
                localerrME += abs(sum(yvalsME-simDist.distribution[:,n,i]))
                localerrDG += abs(sum(yvalsDG-simDist.distribution[:,n,i]))
            end
        end
        # display(q)
        push!(globalerrME,localerrME)
        push!(globalerrDG,localerrDG)
        # for i in 1:NPhases
        #     scatter!(simDist.x,simDist.distribution[:,:,i][:],subplot=i,label="sim")
        # end
        # display(plot!())
    end
    plot!(p,1:2:21,log.(globalerrME),label=false,linestyle=:dash,colour=1+NCounters)
    plot!(p,1:2:21,log.(globalerrDG),label=false,linestyle=:solid,colour=3)
    display(p)
end

max_stars_repo_name: angus-lewis/SFFM

max_issues_repo_name: angus-lewis/SFFM

max_forks_repo_name: angus-lewis/SFFM


text:
mutable struct Snake
    body::Vector
    dir::Vector{<:Real}
    len::Int
    Snake() = begin
        new([ [0, 0] ], [1, 0], 1)
    end
end

function update(snake::Snake)
    head = copy(last(snake.body))
    head += snake.dir
    popfirst!(snake.body)
    push!(snake.body, head)
end

function show(snake::Snake)
    for pos in snake.body
        sethue("white")
        rect(point(pos), 1, 1, :fill)
    end
end

function grow(snake::Snake)
    head = last(snake.body)
    snake.len += 1
    push!(snake.body, head)
end

function eat(snake::Snake, food::Vector)
    if food == last(snake.body)
        grow(snake)
        return true
    else
        return false
    end
end

function gameover(this::Snake)
    x, y = last(this.body)

    if x > w-1 || x < 0 || y > h-1 || y < 0
        return true
    end

    for part in this.body[1:end-1]
        if part == [x, y]
            return true
        end
    end

    return false
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 6a8d7cf4-08d2-11ec-2bfe-77d314082c71
using StatsModels

# ╔═╡ b7bc91a0-c24b-4c18-a37d-9d82ea06174b
using StatsModels: hasresponse, hasintercept, omitsintercept, drop_term

# ╔═╡ 5ecccbc3-37db-4326-b363-c15d2c1c55fb
using Random

# ╔═╡ 45ae7bfb-c4cf-425c-9541-a4410f6bbf02
using Tables

# ╔═╡ a68d1c2d-cedd-4fd7-923b-69bf70c91b50
using DataFrames

# ╔═╡ ea2b9175-ebd3-4f33-b2c4-d947651281c7
using GLM

# ╔═╡ 3b93d4d1-9d32-4e6e-a517-87ad8815f071
using StableRNGs

# ╔═╡ b6c4bdd5-378c-4338-8690-d155abe5eda5
Random.seed!(1234)

# ╔═╡ 43502b67-c04e-4d4b-91dd-1b9fec574097
rng = StableRNG(1)

# ╔═╡ 57233fca-dd70-4113-88e7-3699f2987a10
f1 = @formula(y ~ 1 + a + b + c + b&c) # can only allow ints as intercept and also can be of any form such as @formula(y ~ 5 + a - b * c / b|c)

# ╔═╡ 4504aa01-5bee-4d96-afeb-098e6ce706cf
hasresponse(f1)

# ╔═╡ 565d4afd-1dbf-4668-bd4d-c2112e74d9b2
StatsModels.implicit_intercept(f1)

# ╔═╡ 81046045-8818-4dd5-826b-93162cf41010
StatsModels.drop_intercept(f1)

# ╔═╡ 7d24f2de-f734-4c7b-bbeb-a39353c639c1
omitsintercept(f1)

# ╔═╡ d525b311-cae2-4edd-8585-3051de11761f
f1.rhs == 0 #checks rhs and lhs of the formulae variable

# ╔═╡ 9ebcd948-32c8-499d-bdff-ab15632ec118
terms(f1)

# ╔═╡ 639a83b6-4b11-44ec-9d39-07de19f59900
typeof(f1)

# ╔═╡ 7b99d6b6-b046-4b8a-adb6-51365daa3522
d1 = drop_term(f1, term(:1))

# ╔═╡ 5763afa6-41a8-47dc-b8b0-4a8370bd4340
hasintercept(d1)

# ╔═╡ 5014e8fd-15a7-4202-a550-95c518dda936
a, b = term.((:a, :b))

# ╔═╡ 3be813d9-35e1-4111-ab90-f9b075bad927
@formula(0 ~ a * b).rhs == (a, b, a&b) # This is Star property of Formula Rewrite.

# ╔═╡ 59a1d516-1349-40cf-a63b-ba397c51429c
df = DataFrame(y = rand(rng, 9), a = 1:9, b = rand(rng, 9), c = repeat(["d","e","f"], 3))

# ╔═╡ a0faf1c1-bbb6-4ddb-ae52-9e618b0541e3
f = apply_schema(f1, schema(f1, df))

# ╔═╡ c774b110-cac9-433b-aca8-6e2ae4de7509
resp, pred = modelcols(f, df)

# ╔═╡ 41bc742c-8347-4ea8-9d20-73b53c65cec2
resp

# ╔═╡ 7b823b2d-8152-432b-bed0-8e956b370750
pred

# ╔═╡ de56af73-fb18-4a65-b0bd-10eae75d2861
coefnames(f)

# ╔═╡ d4c50010-b5ff-4e2a-8f7e-7e3638a08d7c
modelmatrix(@formula(y ~ 1 + a + log(1+a)), df)

# ╔═╡ 9de377d7-65bb-4acc-a381-e3ac5641bcc1
lm(@formula(log(y) ~ 1 + a + b), df)

# ╔═╡ 4a115e3e-0393-44f2-965d-5b28ce45844d
# Creating formulae programmatically

# ╔═╡ 77facf95-32cd-46f0-92e7-36d366bbf6e9
Term(:y) ~ ConstantTerm(1) + Term(:a) + Term(:b) + Term(:a) & Term(:b)

# ╔═╡ e1559ad2-81ac-42c9-9773-aa1198bfcd05
ts = term.((5, :a, "b"))

# ╔═╡ 8c2764f2-d3e6-4ccc-8564-5fd70ec79ea6
f2 = term(:y) ~ foldl(+, ts)

# ╔═╡ 6d1f3e8b-f6ad-4b70-9161-3946b4dae736
# Fitting a model from formulae

# ╔═╡ 7ccea56a-9287-44d8-97ea-ca3346273f81
data = DataFrame(a = rand(rng, 100), b = repeat(["w", "x", "y", "z"], 25))

# ╔═╡ 4c78b70c-e582-416d-9601-6d30d82d7819
X = StatsModels.modelmatrix(@formula(y ~ 1 + a*b).rhs, data)

# ╔═╡ 219f7dcc-880e-4dd3-b565-515f7ab76570
data.y = X*(1:8) .+ randn(rng, 100)*0.1 # be careful with the operator sequence

# ╔═╡ 200dfb5a-5cf8-4107-bb88-a928c4c50cf7
mod = fit(LinearModel, @formula(y ~ 1 + a*b * a+b), data)

# ╔═╡ 54f3f851-7858-42c2-b9ed-0cb821983a58
formula(mod)

# ╔═╡ 0b413bf4-4d4f-4d5a-8b86-eb0411bb8a20
y = rand(200)

# ╔═╡ 4cda786f-80f1-415d-931f-96399c737dc4
x1 = rand(200)

# ╔═╡ b932797a-8f5f-437e-8374-a839b31987ea
x2 = rand(200)

# ╔═╡ 9b6b93d1-43c3-44ea-8c5a-f2b0f48afae6
 d = DataFrame(x1=x1, x2=x2, y=y)

# ╔═╡ b341db66-b885-4eeb-a43b-5bc9329c1021
mod1 = fit(LinearModel, @formula(y ~ 1 + x1), d)

# ╔═╡ 9ceea998-0072-4d5a-8eb4-0efa2ddb6b5a
mod2 = fit(LinearModel, @formula(y ~ 1 + x1 * x2), d)

# ╔═╡ cd39d09b-1851-4d22-91b4-5bc08815ed44
StatsModels.isnested(mod1, mod2) # not yet defined

# ╔═╡ 8354c1ba-ef08-4932-bf40-846ddfac262e
StatsModels.lrtest(mod1, mod2)

# ╔═╡ 0c37727f-614d-47dc-b820-68adef96e97a
#  Julia's @formula is meant to be extendable as feasible through Julian's usual techniques of multiple dispatch.

# ╔═╡ e9096f5b-147e-4727-89e2-6ec58d6abcd5
# Modelling categorical data

# ╔═╡ 40a59d50-4507-4e84-a1d1-9df3feeb9d71
#=
StatsModels uses a number of contrast coding techniques to transform category data into a numerical representation appropriate for modelling.
One of the most common contrast-coding systems is mapping a categorical vector with k levels onto a model matrix with k-1 linearly independent rows.
=#

# ╔═╡ 1fdedb94-6afa-4b1a-836c-c7fbdbf1b762
#=
Some of the most contrast coding techniques implemented in StatsModels are:

1. DummyCoding: Code each non-base level as a 0-1 indicator column

2. EffectsCoding: Code each non-base level as 1, and base as -1

3. HelmertCoding: Code each non-base level as the difference from the mean of the lower levels

4. HypothesisCoding: Manually specify contrasts via a hypothesis matrix, which gives the weighting for the average response for each level

5. SeqDiffCoding:  Code for differences between sequential levels of the variable

6. ContrastsCoding: Manually specify contrasts matrix, which is directly copied into the model matrix

The default coding system is DummyCoding. To change the contrast coding system in the model frame, set `contrasts` argument to the desired coding system.

=#

# ╔═╡ 85c7b430-7d36-4246-829b-4548b131892c
f_def = @formula(y ~ 1 + log(w) + x)

# ╔═╡ 3864f3ed-43a6-4df5-a258-1e5ab2d4f6c2
df_def = DataFrame(y = rand(rng, 8), w = rand(rng, 8), x = repeat(["a", "b"], 4))

# ╔═╡ cac06c6b-dd5f-44cb-b3c7-c6e0c6d2002a
schema(f_def, df_def)

# ╔═╡ be6b2a96-7202-4d24-bd2b-95b3342305ff
f_def1 = apply_schema(f_def, schema(f_def, df_def))

# ╔═╡ 62aa8fd4-ec01-4c7b-94e0-931099fb2438
#=
On the basis of the given contrasts and the ModelFrame data, a new schema is computed, and then applied to the FormulaTerm of the ModelFrame.

Only the ModelFrame is mutated as AbstractTerms are immutable and any changes to Terms will produce a copy.
=#

# ╔═╡ 00bda272-74c8-4f14-a2e3-df87d6e3b873
mf = ModelFrame(f_def, df_def, contrasts = Dict(:x => SeqDiffCoding())) # when constructing ModelFrame

# ╔═╡ c6ce6ae3-5f3a-4608-95f0-e0256ffc2584
setcontrasts!(mf, Dict(:x => EffectsCoding())) # when trying to change the already created ModelFrame

# ╔═╡ 7593f656-c72b-47d2-bda0-8203aab852d6
#=

There's a distinction between HypothesisCoding and StatsModels.
ContrastsCoding and other contrast coding systems. These allow you to set a contrast matrix manually.
Contrast matrices M are used for variables that have more than one level, and transfer those levels onto the columns of a model matrix.
This is the full-rank indicator matrix (X), which is equal to 1 in the case of x[i], and zero in all other cases, based on the fact that levels(x)[j] == levels(x)[i].
Therefore, Y = X * M is the model matrix column created by the contrasts matrix, M.

In the context of categorical data levels, a ContrastsMatrix object represents the instantiation of a contrast coding scheme.

=#

# ╔═╡ 30fcd878-7e21-458f-af39-40249dea9dd7
#=
There is a non-zero mean in the columns, and they are orthogonal to each other but collinear with an intercept column (and lower-order columns).
For example, in regression models, the intercept is the mean of the dependant variable at a certain base level when using dummy coding.

It is similar to One-Hot encoding.
=#

StatsModels.ContrastsMatrix(DummyCoding(), ["a", "c", "d", "f"]).matrix

# ╔═╡ 3ef9ab8d-b3dd-4de7-a9c5-f1f3affeca13
#=
Model matrix columns centered on the mean are generated when all levels are equally common.
The resulting columns are not orthogonal when there are more than two levels.The intercept is the grand mean of the regression model with an effects-coded argument.

Similar to DummyCoding, but with a base level of -1 instead of 0.
=#
StatsModels.ContrastsMatrix(EffectsCoding(), ["a", "c", "d", "f"]).matrix

# ╔═╡ ba812640-89f4-4c4c-b978-77f6d082e608
#=
According to Helmert coding, when all levels are equally common, orthogonal columns are generated that are mean-centered.

In Helmert coding, each level is coded as the difference from the average of the levels below it.
=#
StatsModels.ContrastsMatrix(HelmertCoding(), ["a", "c", "d", "f"]).matrix

# ╔═╡ b444f628-95ee-46dd-938a-aaeddbae964c
#=
In order to test "sequential difference" hypotheses, code each level such that it can be compared to the one below.
nth predictor examines a null hypothesis that the difference between levels n and n+1 is zero in this case.

=#
seqdiff = StatsModels.ContrastsMatrix(SeqDiffCoding(), ["a", "c", "d", "f"]).matrix

# ╔═╡ 655f9232-c85e-4094-b757-9be7901b4c61
StatsModels.hypothesis_matrix(seqdiff) # this shows a better picture of the corresponding hypothesis matrix

# ╔═╡ 246223cb-df6e-4dcc-b590-f6ec5e08636b
#=
In Full Dummy coding, each level including the base level, has one indication (1 or 0) column. This is exactly One-hot encoding. This is not exported.
=#

StatsModels.ContrastsMatrix(StatsModels.FullDummyCoding(), ["a", "b", "c"]).matrix

# ╔═╡ d4d4e674-4591-41f3-b544-4c6c4bf87fbd
# Temporal Terms

# ╔═╡ b4c15ce5-4f87-4f0a-af83-b45cd7ac9a89
df1 = DataFrame(y=-10:10, x=-20:2:20) # or use ``y, x = term(:y), term(:x)``

# ╔═╡ d419deee-45b4-4910-a6a7-32b9670d41e4
fun = @formula(y ~ x + lag(x, 2) + lead(x, 2))

# ╔═╡ 1546bb8c-e34b-4a8d-89fa-979cb93b5bf9
func = apply_schema(fun, schema(fun, df1))

# ╔═╡ c10e0eac-a057-4a85-8c19-137ee8e26f63
modelmatrix(func, df1)

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
GLM = "38e38edf-8417-5370-95a0-9cbb8c7f171a"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
StatsModels = "3eaba693-59b7-5ba5-a881-562e759f1c8d"
Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"

[compat]
DataFrames = "~1.2.2"
GLM = "~1.5.1"
StableRNGs = "~1.0.0"
StatsModels = "~0.6.24"
Tables = "~1.5.0"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "bdc0937269321858ab2a4f288486cb258b9a0af7"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.3.0"

[[Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "727e463cfebd0c7b999bbf3e9e7e16f254b94193"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.34.0"

[[CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[Crayons]]
git-tree-sha1 = "3f71217b538d7aaee0b69ab47d9b7724ca8afa0d"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.0.4"

[[DataAPI]]
git-tree-sha1 = "ee400abb2298bd13bfc3df1c412ed228061a2385"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.7.0"

[[DataFrames]]
deps = ["Compat", "DataAPI", "Future", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrettyTables", "Printf", "REPL", "Reexport", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"]
git-tree-sha1 = "d785f42445b63fc86caa08bb9a9351008be9b765"
uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
version = "1.2.2"

[[DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "7d9d316f04214f7efdbb6398d545446e246eff02"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.10"

[[DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[Distributions]]
deps = ["FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SparseArrays", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns"]
git-tree-sha1 = "f389cb8974e02d7eaa6ae2ccedbbfb43174cd8e8"
uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
version = "0.25.14"

[[DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "a32185f5428d3986f47c2ab78b1f216d5e6cc96f"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.5"

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[FillArrays]]
deps = ["LinearAlgebra", "Random", "SparseArrays", "Statistics"]
git-tree-sha1 = "7c365bdef6380b29cfc5caaf99688cd7489f9b87"
uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
version = "0.12.2"

[[Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

[[Future]]
deps = ["Random"]
uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"

[[GLM]]
deps = ["Distributions", "LinearAlgebra", "Printf", "Reexport", "SparseArrays", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns", "StatsModels"]
git-tree-sha1 = "f564ce4af5e79bb88ff1f4488e64363487674278"
uuid = "38e38edf-8417-5370-95a0-9cbb8c7f171a"
version = "1.5.1"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[InvertedIndices]]
deps = ["Test"]
git-tree-sha1 = "15732c475062348b0165684ffe28e85ea8396afc"
uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
version = "1.0.0"

[[IrrationalConstants]]
git-tree-sha1 = "f76424439413893a832026ca355fe273e93bce94"
uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
version = "0.1.0"

[[IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.3.0"

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[LogExpFunctions]]
deps = ["DocStringExtensions", "IrrationalConstants", "LinearAlgebra"]
git-tree-sha1 = "3d682c07e6dd250ed082f883dc88aee7996bf2cc"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.3.0"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "2ca267b08821e86c5ef4376cffed98a46c2cb205"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.1"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[OpenSpecFun_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
version = "0.5.5+0"

[[OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[PDMats]]
deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"]
git-tree-sha1 = "4dd403333bcf0909341cfe57ec115152f937d7d8"
uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
version = "0.11.1"

[[Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[PooledArrays]]
deps = ["DataAPI", "Future"]
git-tree-sha1 = "cde4ce9d6f33219465b55162811d8de8139c0414"
uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720"
version = "1.2.1"

[[Preferences]]
deps = ["TOML"]
git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.2"

[[PrettyTables]]
deps = ["Crayons", "Formatting", "Markdown", "Reexport", "Tables"]
git-tree-sha1 = "0d1245a357cc61c8cd61934c07447aa569ff22e6"
uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
version = "1.1.0"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[QuadGK]]
deps = ["DataStructures", "LinearAlgebra"]
git-tree-sha1 = "12fbe86da16df6679be7521dfb39fbc861e1dc7b"
uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
version = "2.4.1"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[Rmath]]
deps = ["Random", "Rmath_jll"]
git-tree-sha1 = "bf3188feca147ce108c76ad82c2792c57abe7b1f"
uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa"
version = "0.7.0"

[[Rmath_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "68db32dff12bb6127bac73c209881191bf0efbb7"
uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f"
version = "0.3.0+0"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[ShiftedArrays]]
git-tree-sha1 = "22395afdcf37d6709a5a0766cc4a5ca52cb85ea0"
uuid = "1277b4bf-5013-50f5-be3d-901d8477a67a"
version = "1.0.0"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

[[SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[SpecialFunctions]]
deps = ["ChainRulesCore", "LogExpFunctions", "OpenSpecFun_jll"]
git-tree-sha1 = "a322a9493e49c5f3a10b50df3aedaf1cdb3244b7"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "1.6.1"

[[StableRNGs]]
deps = ["Random", "Test"]
git-tree-sha1 = "3be7d49667040add7ee151fefaf1f8c04c8c8276"
uuid = "860ef19b-820b-49d6-a774-d7a799459cd3"
version = "1.0.0"

[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[StatsAPI]]
git-tree-sha1 = "1958272568dc176a1d881acb797beb909c785510"
uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
version = "1.0.0"

[[StatsBase]]
deps = ["DataAPI", "DataStructures", "LinearAlgebra", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"]
git-tree-sha1 = "8cbbc098554648c84f79a463c9ff0fd277144b6c"
uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
version = "0.33.10"

[[StatsFuns]]
deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "Reexport", "Rmath", "SpecialFunctions"]
git-tree-sha1 = "46d7ccc7104860c38b11966dd1f72ff042f382e4"
uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
version = "0.9.10"

[[StatsModels]]
deps = ["DataAPI", "DataStructures", "LinearAlgebra", "Printf", "ShiftedArrays", "SparseArrays", "StatsBase", "StatsFuns", "Tables"]
git-tree-sha1 = "a209a68f72601f8aa0d3a7c4e50ba3f67e32e6f8"
uuid = "3eaba693-59b7-5ba5-a881-562e759f1c8d"
version = "0.6.24"

[[SuiteSparse]]
deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"]
git-tree-sha1 = "d0c690d37c73aeb5ca063056283fde5585a41710"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.5.0"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╠═6a8d7cf4-08d2-11ec-2bfe-77d314082c71
# ╠═b7bc91a0-c24b-4c18-a37d-9d82ea06174b
# ╠═5ecccbc3-37db-4326-b363-c15d2c1c55fb
# ╠═45ae7bfb-c4cf-425c-9541-a4410f6bbf02
# ╠═a68d1c2d-cedd-4fd7-923b-69bf70c91b50
# ╠═ea2b9175-ebd3-4f33-b2c4-d947651281c7
# ╠═3b93d4d1-9d32-4e6e-a517-87ad8815f071
# ╠═b6c4bdd5-378c-4338-8690-d155abe5eda5
# ╠═43502b67-c04e-4d4b-91dd-1b9fec574097
# ╠═57233fca-dd70-4113-88e7-3699f2987a10
# ╠═4504aa01-5bee-4d96-afeb-098e6ce706cf
# ╠═565d4afd-1dbf-4668-bd4d-c2112e74d9b2
# ╠═81046045-8818-4dd5-826b-93162cf41010
# ╠═7d24f2de-f734-4c7b-bbeb-a39353c639c1
# ╠═d525b311-cae2-4edd-8585-3051de11761f
# ╠═9ebcd948-32c8-499d-bdff-ab15632ec118
# ╠═639a83b6-4b11-44ec-9d39-07de19f59900
# ╠═7b99d6b6-b046-4b8a-adb6-51365daa3522
# ╠═5763afa6-41a8-47dc-b8b0-4a8370bd4340
# ╠═5014e8fd-15a7-4202-a550-95c518dda936
# ╠═3be813d9-35e1-4111-ab90-f9b075bad927
# ╠═59a1d516-1349-40cf-a63b-ba397c51429c
# ╠═a0faf1c1-bbb6-4ddb-ae52-9e618b0541e3
# ╠═c774b110-cac9-433b-aca8-6e2ae4de7509
# ╠═41bc742c-8347-4ea8-9d20-73b53c65cec2
# ╠═7b823b2d-8152-432b-bed0-8e956b370750
# ╠═de56af73-fb18-4a65-b0bd-10eae75d2861
# ╠═d4c50010-b5ff-4e2a-8f7e-7e3638a08d7c
# ╠═9de377d7-65bb-4acc-a381-e3ac5641bcc1
# ╠═4a115e3e-0393-44f2-965d-5b28ce45844d
# ╠═77facf95-32cd-46f0-92e7-36d366bbf6e9
# ╠═e1559ad2-81ac-42c9-9773-aa1198bfcd05
# ╠═8c2764f2-d3e6-4ccc-8564-5fd70ec79ea6
# ╠═6d1f3e8b-f6ad-4b70-9161-3946b4dae736
# ╠═7ccea56a-9287-44d8-97ea-ca3346273f81
# ╠═4c78b70c-e582-416d-9601-6d30d82d7819
# ╠═219f7dcc-880e-4dd3-b565-515f7ab76570
# ╠═200dfb5a-5cf8-4107-bb88-a928c4c50cf7
# ╠═54f3f851-7858-42c2-b9ed-0cb821983a58
# ╠═0b413bf4-4d4f-4d5a-8b86-eb0411bb8a20
# ╠═4cda786f-80f1-415d-931f-96399c737dc4
# ╠═b932797a-8f5f-437e-8374-a839b31987ea
# ╠═9b6b93d1-43c3-44ea-8c5a-f2b0f48afae6
# ╠═b341db66-b885-4eeb-a43b-5bc9329c1021
# ╠═9ceea998-0072-4d5a-8eb4-0efa2ddb6b5a
# ╠═cd39d09b-1851-4d22-91b4-5bc08815ed44
# ╠═8354c1ba-ef08-4932-bf40-846ddfac262e
# ╠═0c37727f-614d-47dc-b820-68adef96e97a
# ╠═e9096f5b-147e-4727-89e2-6ec58d6abcd5
# ╠═40a59d50-4507-4e84-a1d1-9df3feeb9d71
# ╠═1fdedb94-6afa-4b1a-836c-c7fbdbf1b762
# ╠═85c7b430-7d36-4246-829b-4548b131892c
# ╠═3864f3ed-43a6-4df5-a258-1e5ab2d4f6c2
# ╠═cac06c6b-dd5f-44cb-b3c7-c6e0c6d2002a
# ╠═be6b2a96-7202-4d24-bd2b-95b3342305ff
# ╠═62aa8fd4-ec01-4c7b-94e0-931099fb2438
# ╠═00bda272-74c8-4f14-a2e3-df87d6e3b873
# ╠═c6ce6ae3-5f3a-4608-95f0-e0256ffc2584
# ╠═7593f656-c72b-47d2-bda0-8203aab852d6
# ╠═30fcd878-7e21-458f-af39-40249dea9dd7
# ╠═3ef9ab8d-b3dd-4de7-a9c5-f1f3affeca13
# ╠═ba812640-89f4-4c4c-b978-77f6d082e608
# ╠═b444f628-95ee-46dd-938a-aaeddbae964c
# ╠═655f9232-c85e-4094-b757-9be7901b4c61
# ╠═246223cb-df6e-4dcc-b590-f6ec5e08636b
# ╠═d4d4e674-4591-41f3-b544-4c6c4bf87fbd
# ╠═b4c15ce5-4f87-4f0a-af83-b45cd7ac9a89
# ╠═d419deee-45b4-4910-a6a7-32b9670d41e4
# ╠═1546bb8c-e34b-4a8d-89fa-979cb93b5bf9
# ╠═c10e0eac-a057-4a85-8c19-137ee8e26f63
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
Base.MainInclude.include("../common/AlarmProc.jl")
###################### # MainEntry ################################
include("../filefind/filesfind.jl")
##xxx-xxx
rootlog ="../log/xxx离线大风识别结果/"
filename1 =joinpath( rootlog, "xxx-xxx-10.14.165.41-201903.csv")
WindLogfile1 = "xxx-xxx-wind-03-31-12-S.log   xxx-xxx-wind-03-31-13-F.log
xxx-xxx-wind-03-31-14-S.log   xxx-xxx-wind-03-31-15-F.log
xxx-xxx-wind-03-31-15-S.log   xxx-xxx-wind-03-31-16-F.log
xxx-xxx-wind-03-31-16-S.log   xxx-xxx-wind-03-31-17-F.log
xxx-xxx-wind-03-31-17-S.log   xxx-xxx-wind-03-31-18-F.log
xxx-xxx-wind-03-31-19-F.log   xxx-xxx-wind-03-31-19-S.log
xxx-xxx-wind-03-31-20-F.log   xxx-xxx-wind-03-31-20-S.log
xxx-xxx-wind-03-31-21-F.log   xxx-xxx-wind-03-31-14-S.log
xxx-xxx-wind-03-31-15-F.log   xxx-xxx-wind-03-31-15-S.log
xxx-xxx-wind-03-31-16-F.log   xxx-xxx-wind-03-31-16-S.log
xxx-xxx-wind-03-31-17-F.log   xxx-xxx-wind-03-31-19-F.log
xxx-xxx-wind-03-31-19-S.log   xxx-xxx-wind-03-31-21-F.log
xxx-xxx-wind-03-31-21-S.log   xxx-xxx-wind-03-31-22-F.log
xxx-xxx-wind-03-31-23-S.log
"

function getRawFilesDT0_T1( rawDataRoot ::String )
    a=[]
    files = findAllFiles(rootdir,".conf",startTIme,endtime)
    for t in files
        Base。Filesystem.size( t)
        a += [ t,Base。Filesystem.size( t) ]
    end
    return a
end

a= getRawFilesDT0_T1( "F:/99-RAWData/xx-xxx机历史数据/" )
for a1 in a
    print(a1)
end
exit(0)
Dt0="20190310045716"
t0 =DT2UInt( DT::String,format1="yyyymmddHHMMSS"  )
deltaT =floor( UInt,2*1024*1024*1024/(526+8192+4) *0.426-1000)
t1 = t0+ deltaT
meterPerSeg1 = 5.0404
chn2SegBegIdx1 = 4097

All_alarms = Main.readInAlarm2Array(filename1, chn2SegBegIdx1,meterPerSeg1)
logfiles = split(WindLogfile1," " )
for logfile in logfiles
    WindLogfile2 =joinpath(rootlog,  logfile)
    Main.ReadInWindLog!(All_alarms,WindLogfile2 )
end

aLarmInWindoffCnt = Main.countAlarms(All_alarms, t0,t1,Main.WindOff)
aLarmInWindONCnt  = Main.countAlarms(All_alarms, t0,t1,Main.WindOn)
println("$aLarmInWindoffCnt ,$aLarmInWindONCnt " )
# exit()
Main.findWindStatusinAlarmPrintout(All_alarms,t0,t1, Main.WindOff )

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
using Plots, LaTeXStrings, JLD, Statistics,Gadfly
function RandomWalk(p,first_step)
    arr=[first_step]
    PosPerTime=[first_step]
    #condition=true
    while true
        if p>rand()
            push!(arr, 1)
        else
            push!(arr, -1)
        end
        PosPerTime=cumsum(arr,dims=1)
        if (10 ∈ PosPerTime) || (-10 ∈ PosPerTime)
            break
        end
    end
    return PosPerTime
end
# PosPerTime=RandomWalk(0.5,-1)
run_num=100000
StartingPosList=[i for i in -10:10]
AvgSteps=[]
for first_step in StartingPosList
    steps=[]
    for n in 1:run_num
        PosPerTime=RandomWalk(0.5,first_step)
        push!(steps, length(PosPerTime))
    end
    push!(AvgSteps, mean(steps))
end
save("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter5\\4.3\\chapter5_4.3_AvgSteps_p=0.5.jld",
 "data", AvgSteps)
scatter(StartingPosList, AvgSteps,legends=false, color=:black, alpha=0.5, dpi=400)
plot!(StartingPosList, AvgSteps,guidefontsize=9, color=:mediumpurple3,
    legends=false, ylabel="Average Steps Before Reaching the Boundaries", xlabel="the Starting Position",
    title="Dependance of the Random Walker's Average Lifetime to It's Starting Position",titlefontsize=9)

savefig("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter5\\4.3\\Figs\\avgstep_startpos_p=0.5.png")
# #plot!(x=StartingPosList[1:end], y=AvgSteps[1:end],guidefontsize=9, color=:mediumpurple3,
# Geom.point, Geom.smooth(method=:loess,smoothing=0.9),
# legends=false, ylabel="Average Steps Before Reaching the Boundaries", xlabel="the Starting Position",
# title="Dependance of the Random Walker's Average Lifetime to It's Starting Position",titlefontsize=9)


text:
using ReferenceTests

@testset "staticweb helpers" for file in ["test.bib"] #, "xampl.bib"] #, "ignace_ref.bib"]
   bib = Bibliography.import_bibtex("../examples/$file")
   result = ""
   for (key, val) in bib
     result *= """
         key: $(key)
           title: '$(Bibliography.xtitle(val))'
           names: '$(Bibliography.xnames(val))'
              in: '$(Bibliography.xin(val))'
            year: '$(Bibliography.xyear(val))'
            link: '$(Bibliography.xlink(val))'
            file: '$(Bibliography.xfile(val))'

         """
   end
   @test_reference "$file.txt" result
end


text:
"""
    _copy!(dst::AbstractVector, dst_sp::Integer, src::AbstractVector, src_sp::Integer, n::Integer)

Copies `n` contiguous elements from  `src` (starting at `src_sp`) into `dst` (starting at `dst_sp`).
"""
@inline function _copy!(dst::AbstractVector, dst_sp::Integer, src::AbstractVector, src_sp::Integer, n::Integer)
    @inbounds for i in 1:n
        dst[dst_sp] = src[src_sp]
        dst_sp += 1
        src_sp += 1
    end

    dst
end

"""
    _left_shift!(A::AbstractVector, sp)

Shifts `A[sp:end]` to the beggining of `A`
"""
@inline function _left_shift!(A::AbstractVector, sp)
    i = 1
    @inbounds for j in sp:length(A)
        A[i] = A[j]
        i += 1
    end

    A
end

"""
    _right_shift!(A::AbstractVector, ep)

Shifts `A[1:ep]` to the end of `A`
"""
@inline function _right_shift!(A::AbstractVector, ep)
    i = length(A)
    @inbounds for j in ep:-1:1
        A[i] = A[j]
        i -= 1
    end

    A
end
#=
@inline function _interchange!(A::AbstractVector, head, tail, n)
    @inbounds for i in 1:n
        tmp = A[head]
        A[head] = A[tail]
        A[tail] = tmp
        head += 1
        tail += 1
    end

    A
end
=#


text:
# download the NMRData repository at https://github.com/AI4DBiological-Systems/NMRData

import NMRDataSetup

using FFTW
import PyPlot
import BSON

PyPlot.close("all")
fig_num = 1

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

# creates the save path if it doesn't exist.
function loadbatchexperimentsinsubfolders(root_path::String, save_dir::String, solvent_ppm_guess::T, solvent_window_ppm::T) where T <: Real

    tmp = readdir(root_path, join = true)
    inds = findall(xx->isdir(xx), tmp)

    experiment_full_paths = tmp[inds]
    experiment_names = readdir(root_path)[inds]
    project_names = collect( "$(experiment_names[i])" for i = 1:length(experiment_names))

    isdir(save_dir) || mkdir(save_dir); # make save folder if it doesn't exist.

    for i = 1:length(experiment_names)
        project_name = project_names[i]
        experiment_full_path = experiment_full_paths[i]

        ## load.
        s_t, S, hz2ppmfunc, ppm2hzfunc, ν_0ppm, fs, SW, α_0ppm, β_0ppm, λ_0ppm, Ω_0ppm,
            results_0ppm, dic, α_solvent, β_solvent, λ_solvent, Ω_solvent,
            results_solvent = NMRDataSetup.loadspectrum(experiment_full_path;
            solvent_ppm = solvent_ppm_guess,
            solvent_window_ppm = solvent_window_ppm)

        ## store.
        save_folder_path = joinpath(save_dir, project_name)
        isdir(save_folder_path) || mkpath(save_folder_path); # make save folder if it doesn't exist.

        save_path = joinpath(save_folder_path, "experiment.bson")
        BSON.bson(save_path,
        s_t = s_t,
        fs = fs,
        SW = SW,
        ν_0ppm = ν_0ppm,
        α_0ppm = α_0ppm,
        β_0ppm = β_0ppm,
        λ_0ppm = λ_0ppm,
        dic = dic)

    end
end

### user inputs.
root_path = "/home/roy/Documents/repo/NMRData/experiments_1D1H/BMRB/similar_settings/BMRB-500-0.5mM/" # the folder that contain the experiment folders that you want to load.
save_dir = "/home/roy/MEGAsync/outputs/NMR/experiments/BMRB-500-0.5mM" # root path where the experiment BSON files will be stored.

solvent_ppm_guess = 4.7 # in units ppm.
solvent_window_ppm = 0.1 # in units ppm.
### end inputs.
loadbatchexperimentsinsubfolders(root_path, save_dir, solvent_ppm_guess, solvent_window_ppm)

### user inputs.
root_path = "/home/roy/Documents/repo/NMRData/experiments_1D1H/BMRB/similar_settings/BMRB-500-2mM/" # the folder that contain the experiment folders that you want to load.
save_dir = "/home/roy/MEGAsync/outputs/NMR/experiments/BMRB-500-2mM" # root path where the experiment BSON files will be stored.

solvent_ppm_guess = 4.7 # in units ppm.
solvent_window_ppm = 0.1 # in units ppm.
### end inputs.
loadbatchexperimentsinsubfolders(root_path, save_dir, solvent_ppm_guess, solvent_window_ppm)

### user inputs.
root_path = "/home/roy/Documents/repo/NMRData/experiments_1D1H/BMRB/similar_settings/BMRB-700-20mM/" # the folder that contain the experiment folders that you want to load.
save_dir = "/home/roy/MEGAsync/outputs/NMR/experiments/BMRB-700-20mM" # root path where the experiment BSON files will be stored.

solvent_ppm_guess = 4.7 # in units ppm.
solvent_window_ppm = 0.1 # in units ppm.
### end inputs.
loadbatchexperimentsinsubfolders(root_path, save_dir, solvent_ppm_guess, solvent_window_ppm)

### user inputs.
root_path = "/home/roy/Documents/repo/NMRData/experiments_1D1H/NRC/bioreactor/HEK293/2011-06-06" # the folder that contain the experiment folders that you want to load.
save_dir = "/home/roy/MEGAsync/outputs/NMR/experiments/NRC/bioreactor-HEK293-2011" # root path where the experiment BSON files will be stored.

solvent_ppm_guess = 4.7 # in units ppm.
solvent_window_ppm = 0.1 # in units ppm.
### end inputs.
loadbatchexperimentsinsubfolders(root_path, save_dir, solvent_ppm_guess, solvent_window_ppm)

### user inputs.
root_path = "/home/roy/Documents/repo/NMRData/experiments_1D1H/NRC/bioreactor/Feng-2015-11-30" # the folder that contain the experiment folders that you want to load.
save_dir = "/home/roy/MEGAsync/outputs/NMR/experiments/NRC/bioreactor-2015" # root path where the experiment BSON files will be stored.

solvent_ppm_guess = 4.7 # in units ppm.
solvent_window_ppm = 0.1 # in units ppm.
### end inputs.
loadbatchexperimentsinsubfolders(root_path, save_dir, solvent_ppm_guess, solvent_window_ppm)

### user inputs.
root_path = "/home/roy/Documents/repo/NMRData/experiments_1D1H/NRC/NRC_8_amino_acid_mixutre_Apr_2021" # the folder that contain the experiment folders that you want to load.
save_dir = "/home/roy/MEGAsync/outputs/NMR/experiments/NRC/NRC-8_amino_acid-Apr2021" # root path where the experiment BSON files will be stored.

solvent_ppm_guess = 4.7 # in units ppm.
solvent_window_ppm = 0.1 # in units ppm.
### end inputs.
loadbatchexperimentsinsubfolders(root_path, save_dir, solvent_ppm_guess, solvent_window_ppm)


text:
using Genie.Router
using DiscordAPI

route("/") do
    "Welcome to ZoomCord"
end

route("/start") do
    DiscordAPI.start()
end

max_stars_repo_name: MLH-Fellowship/ZoomCord

max_issues_repo_name: MLH-Fellowship/ZoomCord

max_forks_repo_name: MLH-Fellowship/ZoomCord


text:
x = rand(1:1_000_000, 100_000_000)

using SortingLab

using BenchmarkTools, SortingAlgorithms

base_sort = @belapsed sort(x)
sortingalg_sort = @belapsed sort(x, alg=RadixSort)
sortinglab_fsort = @belapsed fsort(x)

using Missings: allowmissing
x_w_missing = allowmissing(x)
x_w_missing[rand(1:length(x), 10_000_000)] .= missing

base_sort_wm = @belapsed sort(x_w_missing)
sortingalg_sort_wm = 0
sortingalg_sort_wm = @belapsed sort(x_w_missing, alg=RadixSort) # fails
sortinglab_fsort_wm = @belapsed fsort(x_w_missing)

using StatsPlots

g = groupedbar(
    repeat([" Base\n sort(x)", "SortingAlgorithms.jl\n sort(x, alg=RadixSort)", "  SortingLab.jl\n fsort(x)"], outer=2),
    [base_sort, sortingalg_sort, sortinglab_fsort,
    base_sort_wm, sortingalg_sort_wm, sortinglab_fsort_wm];
    group = repeat(["w/o missing", "with missing"], inner=3),
    ylab = "Seconds"
)

title!(g, "Sort 100m integers")

savefig(g, "benchmarks/fsort_missing_100m_int.png")

#######################
# sortperm
#######################

base_sortperm = @belapsed sortperm(x)
sortingalg_sortperm = @belapsed sortperm(x, alg=RadixSort)
sortinglab_fsortperm = @belapsed fsortperm(x)

base_sortperm_wm = @belapsed sortperm(x_w_missing)
sortingalg_sortperm_wm = 0
sortingalg_sortperm_wm = @belapsed sortperm(x_w_missing, alg=RadixSort) # fails
sortinglab_fsortperm_wm = @belapsed fsortperm(x_w_missing)

using StatsPlots

g2 = groupedbar(
    repeat([" Base\n sortperm(x)", "SortingAlgorithms.jl\n sortperm(x, alg=RadixSort)", "  SortingLab.jl\n fsortperm(x)"], outer=2),
    [base_sortperm, sortingalg_sortperm, sortinglab_fsortperm,
    base_sortperm_wm, sortingalg_sortperm_wm, sortinglab_fsortperm_wm];
    group = repeat(["w/o missing", "with missing"], inner=3),
    ylab = "Seconds"
)

title!(g2, "SortPerm (order) 100m integers")

savefig(g2, "benchmarks/fsortperm_missing_100m_int.png")


text:
include("vulkan/conventions.jl")
include("vulkan/patterns.jl")
include("vulkan/gen.jl")


text:
using MetidaNCA

include("tests.jl")


text:
export
    MultiagentNGSIMEnv,
    reset,
    step,
    observation_space_spec,
    action_space_spec,
    obs_names,
    render

#=
Description:
    Multiagent NGSIM env that plays NGSIM trajectories, allowing a variable
    number of agents to simultaneously control vehicles in the scene
=#
mutable struct MultiagentNGSIMEnv <: Env
    trajdatas::Vector{ListRecord}
    trajinfos::Vector{Dict}
    roadways::Vector{Roadway}
    roadway::Union{Nothing, Roadway} # current roadway
    scene::Scene
    rec::SceneRecord
    ext::MultiFeatureExtractor
    egoids::Vector{Int} # current ids of relevant ego vehicles
    ego_vehs::Vector{Union{Nothing, Vehicle}} # the ego vehicles
    traj_idx::Int # current index into trajdatas
    t::Int # current timestep in the trajdata
    h::Int # current maximum horizon for egoid
    H::Int # maximum horizon
    primesteps::Int # timesteps to prime the scene
    Δt::Float64

    # multiagent type members
    n_veh::Int # number of simultaneous agents
    remove_ngsim_veh::Bool # whether to remove ngsim veh from all scenes
    features::Array{Float64}
    reward::Int # reward / penalty to use in reward augmentation

    # metadata
    epid::Int # episode id
    render_params::Dict # rendering options
    infos_cache::Dict # cache for infos intermediate results
    function MultiagentNGSIMEnv(
            params::Dict;
            trajdatas::Union{Nothing, Vector{ListRecord}} = nothing,
            trajinfos::Union{Nothing, Vector{Dict}} = nothing,
            roadways::Union{Nothing, Vector{Roadway}} = nothing,
            reclength::Int = 5,
            Δt::Float64 = .1,
            primesteps::Int = 50,
            H::Int = 50,
            n_veh::Int = 20,
            remove_ngsim_veh::Bool = false,
            reward::Int = 0,
            render_params::Dict = Dict("zoom"=>5., "viz_dir"=>"/tmp"))
        param_keys = keys(params)
        @assert in("trajectory_filepaths", param_keys)

        # optionally overwrite defaults
        reclength = get(params, "reclength", reclength)
        primesteps = get(params, "primesteps", primesteps)
        H = get(params, "H", H)
        n_veh = get(params, "n_veh", n_veh)
        remove_ngsim_veh = get(params, "remove_ngsim_veh", remove_ngsim_veh)
        reward = get(params, "reward", reward)
        for (k,v) in get(params, "render_params", render_params)
            render_params[k] = v
        end

        # load trajdatas if not provided
        if trajdatas == nothing || trajinfos == nothing || roadways == nothing
            trajdatas, trajinfos, roadways = load_ngsim_trajdatas(
                params["trajectory_filepaths"],
                minlength=primesteps + H
            )
        end

        # build components
        scene_length = max_n_objects(trajdatas)
        scene = Scene(scene_length)
        rec = SceneRecord(reclength, Δt, scene_length)
        ext = build_feature_extractor(params)
        infos_cache = fill_infos_cache(ext)
        # features are stored in row-major order because they will be transferred
        # to python; this is inefficient in julia, but don't change or it will
        # break the python side of the interaction
        features = zeros(n_veh, length(ext))
        egoids = zeros(n_veh)
        ego_vehs = [nothing for _ in 1:n_veh]
        return new(
            trajdatas,
            trajinfos,
            roadways,
            nothing,
            scene,
            rec,
            ext,
            egoids, ego_vehs, 0, 0, 0, H, primesteps, Δt,
            n_veh, remove_ngsim_veh, features, reward,
            0, render_params, infos_cache
        )
    end
end

#=
Description:
    Reset the environment. Note that this environment maintains the following
    invariant attribute: at any given point, all vehicles currently being controlled
    will end their episode at the same time. This simplifies the rest of the code
    by enforcing synchronized restarts, but it does somewhat limit the sets of
    possible vehicles that can simultaneously interact. With a small enough minimum
    horizon (H <= 250 = 25 seconds) and number of vehicle (n_veh <= 100)
    this should not be a problem. If you need to run with larger numbers then those
    implement an environment with asynchronous resets.

Args:
    - env: env to reset
    - dones: bool vector indicating which indices have reached a terminal state
        these must be either all true or all false
=#
function reset(
        env::MultiagentNGSIMEnv,
        dones::Vector{Bool} = fill!(Vector{Bool}(undef, env.n_veh), true);
        offset::Int=env.H + env.primesteps,
        random_seed::Union{Nothing, Int} = nothing)
    # enforce environment invariant reset property
    # (i.e., always either all true or all false)
    @assert (all(dones) || all(.!dones))
    # first == all at this point, so if first is false, skip the reset
    if !dones[1]
        return
    end

    # sample multiple ego vehicles
    # as stated above, these will all end at the same timestep
    env.traj_idx, env.egoids, env.t, env.h = sample_multiple_trajdata_vehicle(
        env.n_veh,
        env.trajinfos,
        offset,
        rseed=random_seed
    )

    # update / reset containers
    env.epid += 1
    empty!(env.rec)
    empty!(env.scene)

    # prime
    for t in env.t:(env.t + env.primesteps)
        get!(env.scene, env.trajdatas[env.traj_idx], t)
        if env.remove_ngsim_veh
            keep_vehicle_subset!(env.scene, env.egoids)
        end
        update!(env.rec, env.scene)
    end

    # set the ego vehicle
    for (i, egoid) in enumerate(env.egoids)
        vehidx = findfirst(egoid, env.scene)
        env.ego_vehs[i] = env.scene[vehidx]
    end
    # set the roadway
    env.roadway = env.roadways[env.traj_idx]
    # env.t is the next timestep to load
    env.t += env.primesteps + 1
    # enforce a maximum horizon
    env.h = min(env.h, env.t + env.H)
    return get_features(env)
end

#=
Description:
    Propagate a single vehicle through an otherwise predeterined trajdata

Args:
    - env: environment to be stepped forward
    - action: array of floats that can be converted into an AccelTurnrate
=#
function _step!(env::MultiagentNGSIMEnv, action::Array{Float64})
    # make sure number of actions passed in equals number of vehicles
    @assert size(action, 1) == env.n_veh
    ego_states = Vector{VehicleState}(undef, env.n_veh)
    # propagate all the vehicles and get their new states
    for (i, ego_veh) in enumerate(env.ego_vehs)
        # convert action into form
	ego_action = AccelTurnrate(action[i,:]...)

        # propagate the ego vehicle
        ego_states[i] = propagate(
            ego_veh,
            ego_action,
            env.roadway,
            env.Δt
        )
        # update the ego_veh
        env.ego_vehs[i] = Entity(ego_veh, ego_states[i])
    end

    # load the actual scene, and insert the vehicles into it
    get!(env.scene, env.trajdatas[env.traj_idx], env.t)
    if env.remove_ngsim_veh
        keep_vehicle_subset!(env.scene, env.egoids)
    end
    orig_vehs = Vector{Vehicle}(undef, env.n_veh)

    for (i, egoid) in enumerate(env.egoids)
	    vehidx = findfirst(egoid, env.scene)

        # track the original vehicle for validation / infos purposes
        orig_vehs[i] = env.scene[vehidx]

    	# replace the original with the controlled vehicle
        env.scene[vehidx] = env.ego_vehs[i]

        # Raunak testing how to access lane id
        # Commented out for now as not relevant to reward augmentation
        # lane4print = env.ego_vehs[i].state.posF.roadind.tag.lane
        # println("Lane number = $lane4print")
    end

    # update rec with current scene
    update!(env.rec, env.scene)

    # Raunak adds in original vehicle properties
    step_infos = Dict{String, Vector{Float64}}(
        "rmse_pos"=>Float64[],
        "rmse_vel"=>Float64[],
        "rmse_t"=>Float64[],
        "x"=>Float64[],
        "y"=>Float64[],
        "s"=>Float64[],
        "phi"=>Float64[],
        "orig_x"=>Float64[],
        "orig_y"=> Float64[],
        "orig_theta"=>Float64[],
        "orig_length"=>Float64[],
        "orig_width"=>Float64[]
    )
    for i in 1:env.n_veh
        push!(step_infos["rmse_pos"], norm((orig_vehs[i].state.posG - env.ego_vehs[i].state.posG)))
        push!(step_infos["rmse_vel"], norm((orig_vehs[i].state.v - env.ego_vehs[i].state.v)))
        push!(step_infos["rmse_t"], norm((orig_vehs[i].state.posF.t - env.ego_vehs[i].state.posF.t)))
        push!(step_infos["x"], env.ego_vehs[i].state.posG.x)
        push!(step_infos["y"], env.ego_vehs[i].state.posG.y)
        push!(step_infos["s"], env.ego_vehs[i].state.posF.s)
        push!(step_infos["phi"], env.ego_vehs[i].state.posF.ϕ)
        push!(step_infos["orig_x"], orig_vehs[i].state.posG.x)
        push!(step_infos["orig_y"], orig_vehs[i].state.posG.y)
        push!(step_infos["orig_theta"], orig_vehs[i].state.posG.θ)
        push!(step_infos["orig_length"], orig_vehs[i].def.length)
        push!(step_infos["orig_width"], orig_vehs[i].def.width)
    end

    return step_infos
end

function _extract_rewards(env::MultiagentNGSIMEnv, infos::Dict{String, Array{Float64}})
    rewards = zeros(env.n_veh)
    for i in 1:env.n_veh
        reward_col = infos["is_colliding"][i] * env.reward
        reward_off = infos["is_offroad"][i] * env.reward
        reward_brake = infos["hard_brake"][i] * env.reward * 0.5  # braking hard is not as bad as a collision
        rewards[i] = -max(reward_col, reward_off, reward_brake)
    end
    return rewards
end

function Base.step(env::MultiagentNGSIMEnv, action::Array{Float64})
    step_infos = _step!(env, action)

    # compute features and feature_infos
    features = get_features(env)
    feature_infos = _compute_feature_infos(env, features)

    # combine infos
    infos = merge(step_infos, feature_infos)

    # update env timestep to be the next scene to load
    env.t += 1

    # compute terminal
    terminal = env.t > env.h ? true : false
    terminal = [terminal for _ in 1:env.n_veh]
    # vectorized sampler does not call reset on the environment
    # but expects the environment to handle resetting, so do that here
    # note: this mutates env.features in order to return the correct obs when resetting
    reset(env, terminal)
	rewards = _extract_rewards(env, infos)
    return deepcopy(env.features), rewards, terminal, infos
end

function _compute_feature_infos(env::MultiagentNGSIMEnv, features::Array{Float64};
                                accel_thresh_min::Float64=-2.0, accel_thresh::Float64=-3.0,
                                min_d_edge_thresh::Float64=0.5, offroad_thresh::Float64=-0.1)
    feature_infos = Dict{String, Array{Float64}}(
                "is_colliding"=>Float64[],
                "is_offroad"=>Float64[],
                "hard_brake"=>Float64[])
    for i in 1:env.n_veh
        is_colliding = features[i, env.infos_cache["is_colliding_idx"]]
        push!(feature_infos["is_colliding"], is_colliding)

        accel = features[i, env.infos_cache["accel_idx"]]
        if accel <= accel_thresh_min
            normalized_accel = abs((accel - accel_thresh_min) / (accel_thresh - accel_thresh_min))
            # linearly increase penalty up to accel_thresh
            push!(feature_infos["hard_brake"], min(1, normalized_accel))
        else
            push!(feature_infos["hard_brake"], 0)
        end

        is_offroad = features[i, env.infos_cache["out_of_lane_idx"]]
        if is_offroad < 1
            d_left = features[i, env.infos_cache["distance_road_edge_left_idx"]]
            d_right = features[i, env.infos_cache["distance_road_edge_right_idx"]]
            closest_d = min(d_left, d_right)
            if closest_d <= min_d_edge_thresh # meaning too close
                is_offroad = abs((closest_d - min_d_edge_thresh) / (offroad_thresh - min_d_edge_thresh))
            end
        end
        push!(feature_infos["is_offroad"], is_offroad)
    end
    return feature_infos
end

function AutoRisk.get_features(env::MultiagentNGSIMEnv)
    for (i, egoid) in enumerate(env.egoids)
        veh_idx = findfirst(egoid, env.scene)
        pull_features!(env.ext, env.rec, env.roadway, veh_idx)
        env.features[i, :] = deepcopy(env.ext.features)
    end
    return deepcopy(env.features)
end

function observation_space_spec(env::MultiagentNGSIMEnv)
    low = zeros(length(env.ext))
    high = zeros(length(env.ext))
    feature_infos = feature_info(env.ext)
    for (i, fn) in enumerate(feature_names(env.ext))
        low[i] = feature_infos[fn]["low"]
        high[i] = feature_infos[fn]["high"]
    end
    infos = Dict("high"=>high, "low"=>low)
    return (length(env.ext),), "Box", infos
end
action_space_spec(env::MultiagentNGSIMEnv) = (2,), "Box", Dict("high"=>[4.,.15], "low"=>[-4.,-.15])
obs_names(env::MultiagentNGSIMEnv) = feature_names(env.ext)
vectorized(env::MultiagentNGSIMEnv) = true
num_envs(env::MultiagentNGSIMEnv) = env.n_veh

#=
Description:
    Render the scene

Args:
    - env: environment to render

Returns:
    - img: returns a (height, width, channel) image to display
=#
function render(
        env::MultiagentNGSIMEnv;
        egocolor::Vector{Float64}=[0.,0.,1.],
        camtype::String="follow",
        static_camera_pos::Vector{Float64}=[0.,0.],
        camera_rotation::Float64=0.,
        canvas_height::Int=800,
        canvas_width::Int=800)
#    # define colors for all the vehicles
#    println("Default version of render")
    carcolors = Dict{Int,Colorant}()
    egocolor = ColorTypes.RGB(egocolor...)
    for veh in env.scene
        carcolors[veh.id] = in(veh.id, env.egoids) ? egocolor : colorant"green"
    end
#
#    # define a camera following the ego vehicle
    if camtype == "follow"
#        # follow the first vehicle in the scene
        cam = AutoViz.CarFollowCamera{Int}(env.egoids[1], env.render_params["zoom"])
    elseif camtype == "static"
        cam = AutoViz.StaticCamera(VecE2(static_camera_pos...), env.render_params["zoom"])
    else
        error("invalid camera type $(camtype)")
    end
    stats = [
        CarFollowingStatsOverlay(env.egoids[1], 2),
        NeighborsOverlay(env.egoids[1], textparams = TextParams(x = 600, y_start=300))
    ]

#    # rendermodel for optional rotation
#    # note that for this to work, you have to comment out a line in AutoViz
#    # src/overlays.jl:27 `clear_setup!(rendermodel)` in render
    rendermodel = RenderModel()
    camera_rotate!(rendermodel, deg2rad(camera_rotation))
#
#    # render the frame
    frame = render(
        env.scene,
        env.roadway,
        stats,
        rendermodel = rendermodel,
        cam = cam,
        car_colors = carcolors,
        canvas_height=canvas_height,
        canvas_width=canvas_width
    )
#
#    # save the frame

    if !isdir(env.render_params["viz_dir"])
        mkdir(env.render_params["viz_dir"])
    end
    ep_dir = joinpath(env.render_params["viz_dir"], "episode_$(env.epid)")
    if !isdir(ep_dir)
        mkdir(ep_dir)
    end
    filepath = joinpath(ep_dir, "step_$(env.t).png")
    AutoViz.write_to_png(frame, filepath)
#
#    # load and return the frame as an rgb array
    img = PyPlot.imread(filepath)
    return img
end


text:
module MaxHelpingHandCustomizableGlassBlockAreaController

using ..Ahorn, Maple

@mapdef Entity "MaxHelpingHand/CustomizableGlassBlockAreaController" CustomizableGlassBlockAreaController(x::Integer, y::Integer, width::Integer=Maple.defaultBlockWidth, height::Integer=Maple.defaultBlockHeight,
    starColors::String="ff7777,77ff77,7777ff,ff77ff,77ffff,ffff77", bgColor::String="302040", wavy::Bool=false)

const placements = Ahorn.PlacementDict(
    "Customizable Glass Block Area Controller (max480's Helping Hand)" => Ahorn.EntityPlacement(
        CustomizableGlassBlockAreaController
    )
)

Ahorn.minimumSize(entity::CustomizableGlassBlockAreaController) = 8, 8
Ahorn.resizable(entity::CustomizableGlassBlockAreaController) = true, true

Ahorn.selection(entity::CustomizableGlassBlockAreaController) = Ahorn.getEntityRectangle(entity)

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::CustomizableGlassBlockAreaController, room::Maple.Room)
    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    Ahorn.drawRectangle(ctx, 0, 0, width, height, (0.4, 0.4, 1.0, 0.4), (0.4, 0.4, 1.0, 1.0))
end

end

max_stars_repo_name: Kitty-Cats/MaxHelpingHand

max_issues_repo_name: Kitty-Cats/MaxHelpingHand

max_forks_repo_name: Kitty-Cats/MaxHelpingHand


text:
@testset "Isogenies" begin
  K = GF(7)
  E1 = EllipticCurve(K, [1, 2, 3, 4, 5])
  E2 = EllipticCurve(K, [1, 2, 3, 1, 1])
  phi = @inferred isogeny_from_kernel(E1, division_polynomial_univariate(E1,3)[1])
  @test @inferred domain(phi) == E1
  @test @inferred codomain(phi) == E2
  @test is_isomorphic(E1, codomain(phi))
  @test @inferred image(phi, E1([1,1])) == E2([6,2])
  @test @inferred degree(phi) == 9

  K = GF(3)
  E1 = EllipticCurve(K, [1, 0, 1, 1, 0])
  E2 = EllipticCurve(K, [1, 0, 1, 0, 2])
  Kx, x = PolynomialRing(K);
  f = x+2
  phi = @inferred isogeny_from_kernel(E1, f)
  @test @inferred domain(phi) == E1
  @test @inferred codomain(phi) == E2
  @test @inferred is_infinite(image(phi, E1([1,2])))
  phihat = @inferred dual_isogeny(phi)
  P = points_with_x(E1, 0)[1]
  @test (phi * phihat)(P) == 2*P || (phi * phihat)(P) == -2*P

  f = @inferred identity_isogeny(E1)
  f(P) == P

  E1 = EllipticCurve([1, 2, 3, 4, 5])
  E2 = EllipticCurve([1, 2, 3, 979//16, 19067//64])
  phi = @inferred isogeny_from_kernel(E1, division_polynomial_univariate(E1,2)[1])
  @test @inferred domain(phi) == E1
  @test @inferred codomain(phi) == E2
  @test is_isomorphic(E1, codomain(phi))
  P = points_with_x(E1, 1)[1]
  g = @inferred multiplication_by_m_map(E1, 32)
  g(P) == 32*P

  K= GF(2,4)
  a = gen(K)
  E1 = EllipticCurve(K,[a^2,1-a,1,0,a])
  E2 = EllipticCurve(K,[a^2,1-a,1,a^8,1])
  f = division_polynomial_univariate(E1,5)[1]
  phi = @inferred isogeny_from_kernel(E1, f)
  @test @inferred domain(phi) == E1
  @test @inferred codomain(phi) == E2
  @test @inferred image(phi, E1([a^2,a^2])) == E2([1,a^14])

  K = GF(13,2)
  E = elliptic_curve_from_j_invariant(K(2))
  div30 = division_polynomial_univariate(E,30)[1]
  L = isogeny_from_kernel_factored(E,div30)
  P = points_with_x(E, 1)[1]
  image(L, P)

  f = @inferred frobenius_map(E)
  Kx, x = PolynomialRing(base_field(E),"x")
  @test isogeny_map_phi(f) == x^13
  f(f(P)) == P

end


text:
"""
    struct DeclContext <: Any
Hold a pointer to a `clang::DeclContext` object.
"""
struct DeclContext
    ptr::CXDeclContext
end

Base.unsafe_convert(::Type{CXDeclContext}, x::DeclContext) = x.ptr
Base.cconvert(::Type{CXDeclContext}, x::DeclContext) = x


text:
# Load modules
using Plots, LaTeXStrings
pyplot()
using ThinFilmsTools
using Optim

##
function glass_reflectance(beam, incident, emergent)
    # Create a measured spectrum of reflection
    layers2 = [ LayerTMMO(incident),
                LayerTMMO(RIdb.glass(beam.λ./1e3); d=250.),
                LayerTMMO(emergent) ]
    sol = tmm_optics(beam, layers2)
    return vec(beam.p.*sol.Spectra.Rp .+ (1.0 - beam.p).*sol.Spectra.Rs)
end
##

# Wavelength range [nm]
λ = 250:900
# Angle of incidence [degrees]
θ = [5.]
# Polarisation (1.0 = p, 0.0 = s, between 0.0 and 1.0 = average)
pol = 0.5
beam = PlaneWave(λ, θ; p=pol)

# Refractive indices of incident (0) and substrate (2)
incident = RIdb.air(beam.λ)
emergent = RIdb.silicon(beam.λ)

# Define the RI model to use
layers = [
    LayerTMMO(incident),
    ModelFit(:sellmeier),
    LayerTMMO(emergent),
]

# Create reflectance spectrum to fit
Rexp = glass_reflectance(beam, incident, emergent)
# plot(beam.λ, Rexp_norm)
# gui()

options = Optim.Options(
    g_abstol=1e-8, g_reltol=1e-8, iterations=10^5, show_trace=true, store_trace=true,
);

seed = vcat(280.0, [1.0, 0.23, 1.0, 6e-3, 2e-2, 1e-2].* (rand(6).*0.2.+1.0))

solOptim = fit_tmm_optics(
    Reflectance(Rexp), [seed], beam, layers;
    options=options, alg=SAMIN(), lb=0.5.*seed, ub=1.5.*seed,
)

plot(FitSpectrum(), solOptim.beam.λ, solOptim.spectrumExp, solOptim.spectrumFit)
gui()


text:
using Dates
import DBInterface: prepare, execute, lastrowid, close!, execute
const LIBPQ_TYPE_MAPPINGS = Dict{Union{Type, Symbol}, Symbol}( # Julia / Postgres
  Char => :CHARACTER,
  String => :VARCHAR,
  Integer => :INTEGER,
  Int => :INTEGER,
  Float64 => :FLOAT,
  DateTime => :TIMESTAMP,
  Time => :TIME,
  Date => :DATE,
  Bool => :BOOLEAN,
  Dict => :JSON,

  :Serial => :SERIAL
)
primary_key_type(dbtype::Type{LibPQ.Connection}, x::Type{T}) where T<:Integer = :Serial

function insert_query(table::Table, column_names::Array, dbtype::Type{LibPQ.Connection})
    values_placeholder = join(repeat(['?'], length(column_names)), ",")
    column_names = join(column_names, ",")
    sql = """
INSERT INTO $(table.name) ($column_names)
VALUES ($values_placeholder)
RETURNING id
    """
    return sql
end

function database_column_type(dbtype::Type{LibPQ.Connection}, d::Union{Type, Symbol}) :: Symbol
    return LIBPQ_TYPE_MAPPINGS[d]
end

function clean_table_query(table::Table, dbtype::Type{LibPQ.Connection})
    return "DELETE FROM $(table.name)"
end
database_kind(c::Type{LibPQ.Connection}) = Relational
close!(db::LibPQ.Connection) = DBInterface.close!(db)
escape_value(dbtype::Type{LibPQ.Connection}, x::AbstractString) = "'$x'"

close!(db::LibPQ.Connection) = LibPQ.close(db)
execute(db::LibPQ.Connection, sql::AbstractString) = LibPQ.execute(db, sql)

struct LibPQStatement
    sql::AbstractString
    params
end
function prepare(db::LibPQ.Connection, sql::AbstractString)
    value_holders = [position[1] for position in findall("?", sql)]
    i = 1
    s = ""
    num_params = 0
    for pos in value_holders
        s *= sql[i:pos-1]
        i = pos+1
        num_params += 1
        s *= "\$$num_params"
    end
    s *= sql[i:end]
    return LibPQ.prepare(db, s)
end
execute(stmt::LibPQ.Statement, params::Array) = LibPQ.execute(stmt, params)
lastrowid(r::LibPQ.Result) =  r[1,1]
Base.getindex(row::LibPQ.Row, symbol::Symbol) = getproperty(row, symbol)

max_stars_repo_name: lucianolorenti/StructDatabaseMapping

max_issues_repo_name: lucianolorenti/StructDatabaseMapping

max_forks_repo_name: lucianolorenti/StructDatabaseMapping