text:
# ## model
include("trajopt_model.jl")

# ## horizon

freq = 100
h = 1.0 / freq
T = Int(floor(0.5 * 0.65 / h)) + 1
Tm = 17 # T / 2

# ## centroidal_quadruped
s = get_simulation("centroidal_quadruped", "flat_3D_lc", "flat")
model = s.model
env = s.env

nx = 2 * model.nq
nc = 4 #model.nc
nu = model.nu + nc + 4 * nc + nc + 4 * nc + 1
nθ = 5

# ## model
d1 = DTO.Dynamics((y, x, u, w) -> centroidal_quadruped_dyn1(model, env, [h], y, x, u, w), nx + nθ + nx + model.nu, nx, nu)
dt = DTO.Dynamics((y, x, u, w) -> centroidal_quadruped_dynt(model, env, [h], y, x, u, w), nx + nθ + nx + model.nu, nx + nθ + nx + model.nu, nu)

dyn = [d1, [dt for t = 2:T-1]...]

# ## initial conditions
mode = :left
body_height = 0.3
foot_x = 0.17
foot_y = 0.15
foot_height = 0.08

q1 = zeros(model.nq)

# body position
q1[1:3] = [0.0; 0.0; body_height]
q1[4:6] = [0.0; 0.0; 0.0]

# foot1
q1[7:9] = [foot_x; foot_y; 0]

# foot2
q1[10:12] = [foot_x; -foot_y; 0]

# foot3
q1[13:15] = [-foot_x; foot_y; 0]

# foot4
q1[16:18] = [-foot_x; -foot_y; 0]

qM = deepcopy(q1)

if mode == :left
    qM[6 + 3] += foot_height # front left
    qM[15 + 3] += foot_height # back right
elseif mode == :right
    qM[9 + 3] += foot_height # front right
    qM[12 + 3] += foot_height # back left
elseif mode == :none
end

qT = copy(q1)

# visualize!(vis, model, [qM], Δt=h);

q_ref = [linear_interpolation(q1, qM, Tm)..., linear_interpolation(qM, qT, Tm)...]
x_ref = [[q_ref[t]; q_ref[t+1]] for t = 1:T]
x1 = x_ref[1]
xM = x_ref[Tm]
xT = x_ref[T]

visualize!(vis, model, q_ref, Δt=h);

# ## objective
obj = DTO.Cost{Float64}[]
for t = 1:T
    if t == T
        function objT(x, u, w)
            J = 0.0
            v = (x[model.nq .+ (1:model.nq)] - x[1:model.nq]) ./ h
            J += 0.5 * 1.0e-1 * dot(v, v)
            J += 0.5 * transpose(x[1:nx] - x_ref[t]) * Diagonal(1000.0 * ones(nx)) * (x[1:nx] - x_ref[t])
            return J
        end
        push!(obj, DTO.Cost(objT, nx + nθ + nx + model.nu, 0))
    elseif t == 1
        function obj1(x, u, w)
            J = 0.0
            v = (x[model.nq .+ (1:model.nq)] - x[1:model.nq]) ./ h
            J += 0.5 * 1.0e-3 * dot(v, v)
            J += 0.5 * transpose(x[1:nx] - x_ref[t]) * Diagonal(1000.0 * ones(nx)) * (x[1:nx] - x_ref[t])
            J += 0.5 * transpose(u[1:model.nu]) * Diagonal(1.0e-3 * ones(model.nu)) * u[1:model.nu]
            J += 1000.0 * u[end] # slack
            return J
        end
        push!(obj, DTO.Cost(obj1, nx, nu))
    else
        function objt(x, u, w)
            J = 0.0
            v = (x[model.nq .+ (1:model.nq)] - x[1:model.nq]) ./ h
            J += 0.5 * 1.0e-3 * dot(v, v)
            u_previous = x[nx + 5 + nx .+ (1:model.nu)]
            u_control = u[1:model.nu]
            w = (u_control - u_previous) ./ h
            J += 0.5 * 1.0e-2 * dot(w, w)
            J += 0.5 * transpose(x[1:nx] - x_ref[t]) * Diagonal(1000.0 * ones(nx)) * (x[1:nx] - x_ref[t])
            J += 0.5 * transpose(u[1:model.nu]) * Diagonal(1.0e-3 * ones(model.nu)) * u[1:model.nu]
            J += 1000.0 * u[end] # slack
            return J
        end
        push!(obj, DTO.Cost(objt, nx + nθ + nx + model.nu, nu))
    end
end

# ## constraints
# initial condition
xl1 = [q1; q1]
xu1 = [q1; q1]
xlt = [-Inf * ones(nx); -Inf * ones(nθ); -Inf * ones(nx); -Inf * ones(model.nu)]
xut = [Inf * ones(nx); Inf * ones(nθ); Inf * ones(nx); Inf * ones(model.nu)]

# final condition
xlT = [qT; qT; -Inf * ones(nθ); -Inf * ones(nx); -Inf * ones(model.nu)]
xuT = [qT; qT; Inf * ones(nθ); Inf * ones(nx); Inf * ones(model.nu)]

ul = [-Inf * ones(model.nu); zeros(nu - model.nu)]
uu = [Inf * ones(model.nu); Inf * ones(nu - model.nu)]

bnd1 = DTO.Bound(nx, nu, xl=xl1, xu=xu1, ul=ul, uu=uu)
bndt = DTO.Bound(nx + nθ + nx + model.nu, nu, xl=xlt, xu=xut, ul=ul, uu=uu)
bndT = DTO.Bound(nx + nθ + nx + model.nu, 0, xl=xlT, xu=xuT)
bnds = [bnd1, [bndt for t = 2:T-1]..., bndT];

cons = DTO.Constraint{Float64}[]
for t = 1:T
    if t == 1
        function constraints_1(x, u, w)
            [
            # equality (16)
            contact_constraints_equality(model, env, h, x, u, w);
            # inequality (28)
            contact_constraints_inequality_1(model, env, h, x, u, w);

            # body/feet constraints
            # x[3] - x_ref[t][3]; # body height
            # x[model.nq + 3] - x_ref[t][model.nq + 3]; # body height
            ]
        end
        push!(cons, DTO.Constraint(constraints_1, nx, nu, idx_ineq=collect(16 .+ (1:28))))
    elseif t == T
        function constraints_T(x, u, w)
            [
            # inequality (8)
            contact_constraints_inequality_T(model, env, h, x, u, w);

            # body/feet constraints
            # x[3] - x_ref[t][3]; # body height
            # x[model.nq + 3] - x_ref[t][model.nq + 3]; # body height
            ]
        end
        push!(cons, DTO.Constraint(constraints_T, nx + nθ + nx + model.nu, nu, idx_ineq=collect(0 .+ (1:8))))
    else
        function constraints_t(x, u, w)
            [
            # equality (16)
            contact_constraints_equality(model, env, h, x, u, w);
            # inequality (32)
            contact_constraints_inequality_t(model, env, h, x, u, w);

            # body/feet constraints
            # x[3] - x_ref[t][3]; # body height
            # x[model.nq + 3] - x_ref[t][model.nq + 3]; # body height
            ]
        end
        push!(cons, DTO.Constraint(constraints_t, nx + nθ + nx + model.nu, nu, idx_ineq=collect(16 .+ (1:32))) )
    end
end

# ## problem
tolerance = 1.0e-3
p = DTO.solver(dyn, obj, cons, bnds,
    options=DTO.Options(
        max_iter=2000,
        tol=tolerance,
        constr_viol_tol=tolerance,
        ))

# ## initialize
x_interpolation = [x_ref[1], [[x_ref[t]; zeros(nθ); zeros(nx); zeros(model.nu)] for t = 2:T]...]
u_guess = [1.0e-4 * rand(nu) for t = 1:T-1] # may need to run more than once to get good trajectory
DTO.initialize_states!(p, x_interpolation)
DTO.initialize_controls!(p, u_guess)

# ## solve
@time DTO.solve!(p)

# ## solution
x_sol, u_sol = DTO.get_trajectory(p)
@show x_sol[1]
@show x_sol[T]
sum([u[end] for u in u_sol[1:end-1]])

# ## visualize
vis = Visualizer()
render(vis)
visualize!(vis, model, x_sol, Δt=h);

q_opt = [x_sol[1][1:model.nq], [x[model.nq .+ (1:model.nq)] for x in x_sol]...]
v_opt = [(x[model.nq .+ (1:model.nq)] - x[0 .+ (1:model.nq)]) ./ h for x in x_sol]
u_opt = [u[1:model.nu] for u in u_sol]
λ_opt = [u[model.nu .+ (1:4)] for u in u_sol]
b_opt = [u[model.nu + 4 .+ (1:16)] for u in u_sol]

# # mirror
# q_mirror = Vector{Float64}[]
# u_mirror = Vector{Float64}[]
# reflection = Diagonal([1.0; -1.0; 1.0])

# for (t, q) in enumerate(q_opt)
#     qfl = q[6 .+ (1:3)]  # front left
#     qfr = q[9 .+ (1:3)]  # front right
#     qbl = q[12 .+ (1:3)] # back left
#     qbr = q[15 .+ (1:3)] # back right

#     push!(q_mirror, [q[1:6]; reflection * qfr; reflection * qfl; reflection * qbr; reflection * qbl])

#     if t < T
#         # set left using right
#         ufl = u_opt[t][0 .+ (1:3)] # front left
#         ufr = u_opt[t][3 .+ (1:3)] # front right
#         ubl = u_opt[t][6 .+ (1:3)] # back left
#         ubr = u_opt[t][9 .+ (1:3)] # back right
#         if t < T
#             push!(u_mirror, [reflection * ufr; reflection * ufl; reflection * ubr; reflection * ubl])
#         end
#     end
# end

# using Plots
# plot(hcat(q_opt..., q_mirror...)', xlabel="time", ylabel="configuration", label="")
# plot(hcat(u_opt..., u_mirror...)', xlabel="time", ylabel="control", label="")

# RoboDojo.visualize!(vis, RoboDojo.centroidal_quadruped, [q_opt..., q_mirror...], Δt=h);

# function mirror_gait(q, u, γ, b, ψ, η, T)
# 	qm = [deepcopy(q)...]
# 	um = [deepcopy(u)...]
# 	γm = [deepcopy(γ)...]
# 	bm = [deepcopy(b)...]
# 	ψm = [deepcopy(ψ)...]
# 	ηm = [deepcopy(η)...]

#     # stride = zero(qm[1])
# 	@show strd = q[T+1][1] - q[2][1]

# 	p10 = q[2][6 .+ (1:3)]
# 	p20 = q[2][9 .+ (1:3)]
# 	p30 = q[2][12 .+ (1:3)]
# 	p40 = q[2][15 .+ (1:3)]

# 	for t = 1:T-1
# 		# configuration
# 		pt = q[t+2][1:3]
# 		a = q[t+2][3 .+ (1:3)]
# 		p1 = q[t+2][6 .+ (1:3)]
# 		# p2 = q[t+2][9 .+ (1:3)]
# 		# p3 = q[t+2][12 .+ (1:3)]
# 		p4 = q[t+2][15 .+ (1:3)]

# 		p2_diff = q[t+2][9 .+ (1:3)] - p20
# 		p3_diff = q[t+2][12 .+ (1:3)] - p30

# 		push!(qm, [pt + [0.5 * strd; 0.0; 0.0]; a;
# 			p10[1] + p2_diff[1]; q[t+2][8]; p10[3] + p2_diff[3]
# 			q[end][9 .+ (1:3)];
# 			q[end][12 .+ (1:3)];
# 			p40[1] + p3_diff[1]; q[t+2][17]; p40[3] + p3_diff[3];
# 			])
# 		# push!(qm, q[t+2])

# 		# control
# 		u1 = u[t][1:3]
# 		u2 = u[t][3 .+ (1:3)]
# 		u3 = u[t][6 .+ (1:3)]
# 		u4 = u[t][9 .+ (1:3)]
# 		push!(um, [u2; u1; u4; u3])

# 		# impact
# 		γ1 = γ[t][1]
# 		γ2 = γ[t][2]
# 		γ3 = γ[t][3]
# 		γ4 = γ[t][4]
# 		push!(γm, [γ2; γ1; γ4; γ3])
# 		# push!(γm, γ[t])

# 		# friction
# 		b1 = b[t][1:4]
# 		b2 = b[t][4 .+ (1:4)]
# 		b3 = b[t][8 .+ (1:4)]
# 		b4 = b[t][12 .+ (1:4)]
# 		push!(bm, [b2; b1; b4; b3])
# 		# push!(bm, b[t])

# 		# dual
# 		ψ1 = ψ[t][1]
# 		ψ2 = ψ[t][2]
# 		ψ3 = ψ[t][3]
# 		ψ4 = ψ[t][4]
# 		push!(ψm, [ψ2; ψ1; ψ4; ψ3])
# 		# push!(ψm, ψ[t])

# 		# dual
# 		η1 = η[t][1:4]
# 		η2 = η[t][4 .+ (1:4)]
# 		η3 = η[t][8 .+ (1:4)]
# 		η4 = η[t][12 .+ (1:4)]
# 		push!(ηm, [η2; η1; η4; η3])
# 		# push!(ηm, η[t])
# 	end

# 	return qm, um, γm, bm, ψm, ηm
# end

q_opt = [x_sol[1][1:model.nq], [x[model.nq .+ (1:model.nq)] for x in x_sol]...]
v_opt = [(x[model.nq .+ (1:model.nq)] - x[0 .+ (1:model.nq)]) ./ h for x in x_sol]
u_opt = [u[1:model.nu] for u in u_sol]
γ_opt = [u[model.nu .+ (1:4)] for u in u_sol]
b_opt = [u[model.nu + 4 .+ (1:16)] for u in u_sol]
ψ_opt = [u[model.nu + 4 + 16 .+ (1:4)] for u in u_sol]
η_opt = [u[model.nu + 4 + 16 + 4 .+ (1:16)] for u in u_sol]

if mode == :right
    qr = q_opt
    vr = v_opt
    ur = u_opt
    γr = γ_opt
    br = b_opt
    ψr = ψ_opt
    ηr = η_opt

    using JLD2
    @save joinpath(@__DIR__, "inplace_trot_right.jld2") qr ur γr br ψr ηr
    @load joinpath(@__DIR__, "inplace_trot_right.jld2") qr ur γr br ψr ηr
end

if mode == :left
    ql = q_opt
    vl = v_opt
    ul = u_opt
    γl = γ_opt
    bl = b_opt
    ψl = ψ_opt
    ηl = η_opt

    using JLD2
    @save joinpath(@__DIR__, "inplace_trot_left.jld2") ql ul γl bl ψl ηl
    @load joinpath(@__DIR__, "inplace_trot_left.jld2") ql ul γl bl ψl ηl
end

if mode == :none
    ql = q_opt
    vl = v_opt
    ul = u_opt
    γl = γ_opt
    bl = b_opt
    ψl = ψ_opt
    ηl = η_opt

    using JLD2
    @save joinpath(@__DIR__, "inplace_trot_none.jld2") ql ul γl bl ψl ηl
    @load joinpath(@__DIR__, "inplace_trot_none.jld2") ql ul γl bl ψl ηl
end

@load joinpath(@__DIR__, "inplace_trot_left.jld2") ql ul γl bl ψl ηl

qm = [qr[2:end]..., ql[2:end]...]
um = [ur..., ul...]
γm = [γr..., γl...]
bm = [br..., bl...]
ψm = [ψr..., ψl...]
ηm = [ηr..., ηl...]
μm = model.μ_world
hm = h
timesteps = range(0.0, stop=(h * (length(qm) - 2)), length=(length(qm) - 2))
plot(timesteps, hcat(qm[2:end-1]...)', labels="")
plot(timesteps, hcat(um...)', labels="")
plot(timesteps, hcat(γm...)', labels="")
plot(timesteps, hcat(bm...)', labels="")
plot(timesteps, hcat(ψm...)', labels="")
plot(timesteps, hcat(ηm...)', labels="")

using JLD2
@save joinpath(@__DIR__, "inplace_trot.jld2") qm um γm bm ψm ηm μm hm
@load joinpath(@__DIR__, "inplace_trot.jld2") qm um γm bm ψm ηm μm hm


text:
#
# this is from https://docs.julialang.org/en/v1/manual/performance-tips/#Avoid-changing-the-type-of-a-variable
#
# conclusion:
#  1. add concrete types to local variables
#  2. or initialise to something explicit
#  3. abstract types can work as well in some cases (abstract float)
#  4. adding return type doesn't seem to matter for performance purposes
#
module Stability

    using BenchmarkTools
    using InteractiveUtils

    const NR = 10_000

    function unstable()
        x = 1
        for i in 1:NR
            r = rand()
            if( r < 0.5 )
                x += r
            end
        end
        return x
    end

    function stable1()::Float64
        x = 1
        for i in 1:NR
            r = rand()
            if( r < 0.5 )
                x += r
            end
        end
        return x
    end

    function stable2()
        x = 1.0
        for i in 1:NR
            r = rand()
            if( r < 0.5 )
                x += r
            end
        end
        return x
    end

    function stable3()::Float64
        x = 1.0
        for i::Int in 1:10_000
            r = rand()
            if( r < 0.5 )
                x += r
            end
        end
        return x
    end

    function stable4()::Float64
        x :: AbstractFloat = 1
        for i in 1:NR
            r = rand()
            if( r < 0.5 )
                x += r
            end
        end
        return x
   end

    function do_all()::BenchmarkGroup
        stables = BenchmarkGroup()
        stables[:stable1] = @benchmark stable1()
        stables[:stable2] = @benchmark stable2()
        stables[:stable3] = @benchmark stable3()
        stables[:stable4] = @benchmark stable4()
        stables[:unstable] = @benchmark unstable()
        stables
   end

end # module

@code_warntype( Stability.unstable() )
@code_warntype( Stability.stable3() )
@code_warntype( Stability.stable4() )

Stability.do_all()

ext: jl

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
using Test
using Suppressor
using LambdaTAME

#src: https://github.com/charlescolley/DistributedTensorConstruction.jl.git
using DistributedTensorConstruction

import LinearAlgebra: norm
import Random: seed!

#= Potential Kwargs:
   tol - (Float)
   iterations - (Int)
   profile- (bool)
   alphas - (Array{Float64,1})
   betas - (Array{Float64,1})
   rank - Int
=#

#seed
seed!(54321)

NORM_CHECK_TOL = 1e-15
    # norm checks may be scaled by n*m to account for round off error
    #    please see 2.7.6 of 4th edition
    #    of Matrix Computations for more.
function round_off_bound(X::Matrix{T}) where T
    return size(X,1)*size(X,2)*NORM_CHECK_TOL
end
function round_off_bound(x::Vector{T}) where T
    return length(x)*NORM_CHECK_TOL
end

tensor_A_file="test_tensors/test_tensorB.ssten"
tensor_B_file="test_tensors/test_tensorA.ssten"
matrix_A_file="test_tensors/testA.smat"
matrix_B_file="test_tensors/testB.smat"
    #NOTE Flipped to address m >= n assertion error in bipartite_matching_primal_dual

# Third Order Symmetric Tensor
A_TOST = load_ThirdOrderSymTensor(tensor_A_file)
B_TOST = load_ThirdOrderSymTensor(tensor_B_file)

#TODO: rename to A_TOSTU
# Unweighted Third Order Symmetric Tensors
A_UTOST  = load_UnweightedThirdOrderSymTensor(tensor_A_file)
B_UTOST = load_UnweightedThirdOrderSymTensor(tensor_B_file)

# Unweighted Symmetric Tensors
#     src: DistributedTensorConstruction
A_UST = load_SymTensorUnweighted(tensor_A_file,Clique(),'\t')
B_UST = load_SymTensorUnweighted(tensor_B_file,Clique(),'\t')

A_UST_Cycle = load_SymTensorUnweighted(tensor_A_file,Cycle(),'\t')
B_UST_Cycle = load_SymTensorUnweighted(tensor_B_file,Cycle(),'\t')

#for multi-motif routines
A_UST_MM = [A_UST]#, A_UST]
B_UST_MM = [B_UST]#,] B_UST]

A_UST_Cycle_MM = [A_UST_Cycle]#, A_UST]
B_UST_Cycle_MM = [B_UST_Cycle]#,] B_UST]

# -- test_noerror solution -- #
# src: https://github.com/JuliaLang/julia/issues/18780#issuecomment-863505347
struct NoException <: Exception end
macro test_nothrow(ex)
    esc(:(@test_throws NoException ($(ex); throw(NoException()))))
end

include("Experiments_test.jl")
include("fileio_test.jl")
include("Contraction_test.jl")
include("Matchings_test.jl")
include("AlignmentDrivers_test.jl")
include("PostProcessing_test.jl")
include("TAME_Implementation_test.jl")
include("LowRankTAME_Implementation_test.jl")
include("LambdaTAME_Implementation_test.jl")

max_stars_repo_name: chuckcol/LambdaTAME

max_issues_repo_name: chuckcol/LambdaTAME

max_forks_repo_name: chuckcol/LambdaTAME


text:
# from #605

# using Revise

using IncrementalInference
using Test
using TensorCast

##

@testset "test mixture prior" begin

##

# init graph
fg = initfg()
N = 100
getSolverParams(fg).N = N;

# add first variable
addVariable!(fg, :x0, ContinuousScalar)

# also test AliasingScalingSampler
v = rand(50)
v[20:29] .+= 5*rand(10)
v ./= sum(v)
bss = AliasingScalarSampler(collect(1:50), v)

# add bi-modal mixture prior
Prior0 = Mixture(Prior,(a=Normal(-5.0,1.0), b=Uniform(0.0,1.0)), (0.5,0.5))
Prior0 = Mixture(Prior,(Normal(-5.0,1.0), Normal(0.0,1.0)), (0.5,0.5))
Prior0 = Mixture(Prior,[Normal(-5.0,1.0), Normal(0.0,1.0)], (0.5,0.5))
Prior0 = Mixture(Prior,(Normal(-5.0,1.0), Normal(0.0,1.0)), [0.5;0.5])
Prior0 = Mixture(Prior,[Normal(-5.0,1.0), Normal(0.0,1.0)], [0.5;0.5])
Prior0 = Mixture(Prior,(Normal(-5.0,1.0), bss), Categorical([0.5;0.5]))
f1 = addFactor!(fg, [:x0], Prior0)

# also test serialization of AliasingScalarSampler
saveDFG("/tmp/test_fg_bss", fg)

# check numerics -- replaced by CalcFactor{<:Mixture},
smpls_ = approxConv(fg, :x0f1, :x0)
# smpls, = getSample(Prior0, N) # ,lb=

# lazy
@cast smpls[i,j] := smpls_[j][i]

# should be a balance of particles
# @test sum(lb .== 1) - sum(lb .== 2) |> abs < 0.3*N
@test sum(smpls .< -2.5) - sum(-2.5 .< smpls) |> abs < 0.3*N

# solve
solveTree!(fg);

marginalPts_ = getBelief(fg, :x0) |> getPoints

# lazy
@cast marginalPts[i,j] := marginalPts_[j][i]

# check solver solution consistent too
@test sum(marginalPts .< -2.5) - sum(-2.5 .< marginalPts) |> abs < 0.3*N

##

end

@testset "Serialization of Mixture(Prior,..) including a AliasingScalarSampler" begin

##

fg_ = loadDFG("/tmp/test_fg_bss")

N = getSolverParams(fg_).N

solveTree!(fg_);

marginalPts_ = getBelief(fg_, :x0) |> getPoints

# lazy
@cast marginalPts[i,j] := marginalPts_[j][i]

# check solver solution consistent too
@test sum(marginalPts .< -2.5) - sum(-2.5 .< marginalPts) |> abs < 0.3*N

# cleanup
Base.rm("/tmp/test_fg_bss.tar.gz")

##

end

# using RoMEPlotting
# Gadfly.set_default_plot_size(35cm,20cm)

# # plot the results
# plotKDE(fg, :x0)

#


text:
"""Lateral flow in a HF135 nitrocellulose membrane."""
module ExampleHF135

using Fronts
using Fronts.PorousModels
using Plots

ϵ = 1e-7

# Wetting of an HF135 membrane, Van Genuchten model
# Data from Buser (PhD thesis, 2016)
# http://hdl.handle.net/1773/38064
k = 5.50e-13  # m^2
α = 0.2555  # 1/m
n = 2.3521
θi = 0.102755  # Computed from P0

θr = 0.0473
θs = 0.945

model = VanGenuchten(n=n, α=α, k=k, θr=θr, θs=θs)

prob = DirichletProblem(model, i=θi, b=θs-ϵ)

θ = solve(prob)

r = range(0, 0.05, length=500)

plt = plot(r, r -> θ(r,60), label="t=60 s", xguide="r [m]", yguide="θ")
display(plt)

plt = plot(r, r -> flux(θ, r, 60), label="t=60 s", xguide="r [m]", yguide="U [m/s]")
display(plt)

end


text:
using MuZero
using GridGames

include("game.jl")
include("muzero_config.jl")
config = getconfig(Tac)

learner = Learner(
    config;
    uuid="ad24fc40-a07b-4418-9fbc-adb0278dd4e1",
    dir="~/models/MuZeroModels.jl/tictactoe/",
)
l = learner

games = newgames(learner, 2)

is_terminal(games[1])
is_terminal(games)

movesmaxedout(learner, games[1])
movesmaxedout(learner, games)

get_state(games[1]) |> typeof
get_state.(games) |> typeof

root = Node.(0, get_state.(games))
net_output = initial_inference(net, games) |> NetworkOutput

add_exploration_noise(learner.config, root)

store_search_statistics.(games, root)

action = [1,2]
num_moves = length.(history.(game))

run_mcts(net, learner.config, root[1], action_history(games[1]))

action = select_action(
    root[1];
    num_moves=num_moves[1],
    num_steps=net.step,
    temp_fn=learner.config.visit_softmax_temperature_fn,
    preplanned=nothing
)

action = select_action(
    root;
    num_moves=num_moves,
    num_steps=net.step,
    temp_fn=learner.config.visit_softmax_temperature_fn,
    preplanned=nothing
)

apply!.(games, action)

#
# legacy
#
net = getplaynet(learner)
game = new_game(learner.config, net.step)
config = learner.config
while !is_terminal(game) && !movesmaxedout(learner, game)
    root = Node(0, get_state(game))
    net_output = initial_inference(net, game) |> SingleNetworkOutput
    expand_node!(root, legal_actions(game), net_output, Float32(0))
    add_exploration_noise(config, root)
    run_mcts(net, config, root, action_history(game))
    store_search_statistics(game, root)
    action = select_action(
        root;
        num_moves=length(history(game)),
        num_steps=net.step,
        temp_fn=config.visit_softmax_temperature_fn,
        preplanned=nothing
    )
    apply!(game, action)
end

#
# novel
#
isactive(game) = !is_terminal(game) && !movesmaxedout(l, game)
net = getplaynet(l)
games = newgames(l, 100)
config = l.config
active = view(games, :)

@time begin
    while !isempty(active)
        root = Node.(0, get_state.(active))
        net_output = initial_inference(net, active) |> NetworkOutput
        expand_node!(root, legal_actions.(active), net_output, zeros(Float32, length(root)))
        add_exploration_noise(config, root)
        run_mcts(net, config, root, action_history.(active))
        store_search_statistics.(active, root)
        action = select_action.(
            root;
            num_moves=length.(history.(active)),
            num_steps=net.step,
            temp_fn=config.visit_softmax_temperature_fn,
            preplanned=nothing
        )
        apply!.(active, action)
        active = view(active, isactive.(active))
    end
end

@time begin
    for _ in 1:100
        play_game(l)
    end
end


text:
# A.*B

FAhadamardprodB(A,B)=(A.*B,nothing)

function FAhadamardprodB_inplace(value,aux,A,B)
    copy!(value,A.*B)
end

function DAhadamardprodB(derivativeIDX,f_c,faux_c,grad_c,grad_n,A,B)
    if derivativeIDX==1
        axpy!(1.0,grad_c.*B,grad_n)
    elseif derivativeIDX==2
        axpy!(1.0,grad_c.*A,grad_n)
    end
end

if (PROC=="GPU") || (PROC=="GPU32")

    function FAhadamardprodB(A::CudaArray,B::CudaArray)
        tmp=CudaArray(Float64,size(A))
        vmult!(1.0,A,B,tmp)
        return (tmp,nothing)
    end

    function FAhadamardprodB_inplace(value::CudaArray,aux,A::CudaArray,B::CudaArray)
        vmult!(1.0,A,B,value)
    end

    function DAhadamardprodB(derivativeIDX,f_c,faux_c,grad_c,grad_n,A::CudaArray,B::CudaArray)
    if derivativeIDX==1
        vmultupdate!(1.0,grad_c,B,grad_n)
    elseif derivativeIDX==2
        vmultupdate!(1.0,grad_c,A,grad_n)
    end
end

end

Derivative[FAhadamardprodB]=DAhadamardprodB
Inplace[FAhadamardprodB]=FAhadamardprodB_inplace

import Base. .*
.*(A::ADnode,B::ADnode)=ADnode(FAhadamardprodB,[A B])

export .*

export FAhadamardprodB

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
if VERSION < v"1.3"
    include("build_p7zip.v16.2.0.jl")
end


text:
"""
Δdielectric, update_locations = get_update_maps(background_dielectric, desired_dielectric)\n
Computes the difference between two dielectric maps.\n
Useful when comparing with other simulation techniques like FDTD.\n
From the resulting outcomes the mapping matrix S and the dielectric update vector C can be computed.\n
The inputs can be obtained from:\n
\n
background_dielectric = cell_to_yee_dielectric(vσ, vϵ, x, y, z)\n
desired_dielectric = cell_to_yee_dielectric(vσ, vϵ, x, y, z)\n
\n
Where the vσ and vϵ are the conductivity and permittivity for the different simulations.\n
"""
function get_update_maps(background_dielectric::NamedTuple, desired_dielectric::NamedTuple)
    dimensionsˣ = size(desired_dielectric[:σˣ])
    dimensionsʸ = size(desired_dielectric[:σʸ])
    dimensionsᶻ = size(desired_dielectric[:σᶻ])

    # compute the update dielectric
    σˣ = round.(desired_dielectric[:σˣ] - background_dielectric[:σˣ],digits = 3)
    σʸ = round.(desired_dielectric[:σʸ] - background_dielectric[:σʸ],digits = 3)
    σᶻ = round.(desired_dielectric[:σᶻ] - background_dielectric[:σᶻ],digits = 3)
    ϵˣ = round.(desired_dielectric[:ϵˣ] - background_dielectric[:ϵˣ],digits = 3)
    ϵʸ = round.(desired_dielectric[:ϵʸ] - background_dielectric[:ϵʸ],digits = 3)
    ϵᶻ = round.(desired_dielectric[:ϵᶻ] - background_dielectric[:ϵᶻ],digits = 3)

    #alloc update map
    updateLocˣ = zeros(Bool, dimensionsˣ...)
    updateLocʸ = zeros(Bool, dimensionsʸ...)
    updateLocᶻ = zeros(Bool, dimensionsᶻ...)

    #compute the update maps
    tmp = findall(x -> abs(x) > 0.0, σˣ)
    updateLocˣ[tmp] .= true
    tmp = findall(x -> abs(x) > 0.0, ϵˣ)
    updateLocˣ[tmp] .= true

    tmp = findall(x -> abs(x) > 0.0, σʸ)
    updateLocʸ[tmp] .= true
    tmp = findall(x -> abs(x) > 0.0, ϵʸ)
    updateLocʸ[tmp] .= true

    tmp = findall(x -> abs(x) > 0.0, σᶻ)
    updateLocᶻ[tmp] .= true
    tmp = findall(x -> abs(x) > 0.0, ϵᶻ)
    updateLocᶻ[tmp] .= true

    return (;σˣ, ϵˣ, σʸ, ϵʸ, σᶻ, ϵᶻ),(;updateLocˣ, updateLocʸ, updateLocᶻ)
end


text:
import GeoInterface
using Base.Test, LibGEOS

@testset "LibGEOS" begin
    include("test_geos_functions.jl")
    include("test_geos_operations.jl")
    include("test_geo_interface.jl")
    include("test_regressions.jl")
    include("test_invalid_geometry.jl")
end


text:
"""
    add_obj_fun!(ata_model::Union{MaximinModel,MinimaxModel}; kwargs...)

# Description

It adds the objective function as specified in the `settings_file`. It requires the [`start_ata`](#ATA.start_ata) build step.

Computes the IIFs at predefined ability points using the item parameter estimates in the item bank.
For each pair of item, and ability point \$(i, θ_k)\$, an IIF\$(θ_k)_{i}\$ is computed.

# Arguments

- **`ata_model::Union{MaximinModel, MinimaxModel}`** : Required.
"""
function add_obj_fun!(ata_model::Union{MaximinModel,MinimaxModel}; kwargs...)

    try
        compute_estimated_iif!(ata_model)
        success!(ata_model, "Estimated IIFs computed.")

    catch e
        error!(ata_model, string(sprint(showerror, e)))

    end
    return nothing
end


text:
module EMLcstr
	using DanaTypes
	using DotPlusInheritance
	using Reexport
	@reexport using ...streams.EMLstreams
	import EMLtypes.length
	include("cstr/cstr_basic.jl")
	include("cstr/cstr_yield.jl")
	include("cstr/cstr.jl")
end


text:
module GEMDataReader

import JSON
using FileOps
export convert_file, convert_files, GEMDataFile, read_run

const MAX_TAPPERS = 0x04
# ---------------------------------------------------------------------------- #
struct GEMDataFile
    path::String
    hdr::Dict
    idx_map::Vector{UInt64}
    nrun::Int
    nwin::Int
end
# ---------------------------------------------------------------------------- #
struct GEMPacket
    dtp::UInt8
    window::UInt16
    click_time::UInt32
    async::Vector{Int16}
    adjust::Int16
end
# ---------------------------------------------------------------------------- #
function GEMDataFile(filepath::String)
    open(filepath, "r") do io
        hdr_length = read(io, UInt64)
        println("Header length: ", hdr_length)

        # hdr = JSON.parse(String(read(io, UInt8, hdr_length)))
        hdr = JSON.parse(String(read(io, hdr_length)))
        println("Header: ", hdr)

        nrun = length(hdr["metronome_alpha"]) * Int64(hdr["repeats"])
        println("#runs: ", nrun)

        # idx_map = read(io, UInt64, nrun)
        idx_map = reinterpret(UInt64, read(io, nrun*sizeof(UInt64)))
        println("Run index map: ", idx_map)

        return GEMDataFile(filepath, hdr, idx_map, nrun, Int(hdr["windows"]))
    end
end
# ---------------------------------------------------------------------------- #
function swapext(file::String, ext::String)
    ext = ext[1] == "." ? ext : "." * ext
    return splitext(file)[1] * ext
end
# ---------------------------------------------------------------------------- #
"""
    ofiles = convert_files(dir::String, force::Bool=false, depth::Integer=1)

# In:
- `dir::String`: the path to a base directory to search within
- `force::Bool`: true to force reconversion of gdf files where the resulting csv file already exists
- `depth::Integer`: the depth of sub-directories to search

# Out:
- `ofiles::Vector{String}`: resulting csv files paths (successful conversions only)

# Examples:

```julia
ofiles = convert_files("../Data/")
```
"""
function convert_files(dir::String, force::Bool=false, depth::Integer=1)
    dirs = find_directories(dir, r"\d{8}")
    for k in eachindex(dirs)
        for j in 1:depth-1
            # only folders named as yyyymmdd
            append!(dirs, find_directories(dirs[j], r"\d{8}"))
        end
    end

    files = Vector{String}()
    for dir in dirs
        append!(files, find_files(dir, r"\.gdf$"))
    end

    ofiles = Vector{String}()

    for k in eachindex(files)
        ofile = swapext(files[k], "csv")
        if !isfile(ofile) || force
            try
                convert_file(files[k], ofile)
                push!(ofiles, ofile)
            catch err
                warn("Failed to convert file: \"$(files[k])\" - $(err)")
            end
        end
    end
    return ofiles
end
# ---------------------------------------------------------------------------- #
function convert_file(ifile::String, ofile::String)

    df = GEMDataFile(ifile)
    hdrs = Vector{Dict}()

    open(ofile, "w") do io
        labels = join(["\"async_$(x)\"" for x in 1:MAX_TAPPERS], ", ")
        write(io,
            """
            \"run\", \"alpha\", \"window\", \"click_time\", $(labels), \"next_adjust\"
            """
        )
        for k = 1:df.nrun
            hdr, data = read_run(df, k)
            if isempty(data)
                break
            end
            push!(hdrs, hdr)
            write_run(io, hdr, data, k)
        end

        # add run headers to the file header dict for output
        df.hdr["run_headers"] = hdrs
    end

    open(swapext(ofile, "json"), "w") do io
        write(io, JSON.json(df.hdr))
    end

    return nothing
end
# ---------------------------------------------------------------------------- #
# Columns: run #, alpha, window #, scheduled click, async 1-4, next_adjust
function write_run(io::IOStream, hdr::Dict, data::Vector{GEMPacket},
    run::Integer)
    for k in eachindex(data)
        write(io, "$(run), $(hdr["alpha"]), $(k), " * tostring(data[k]) * "\n")
    end
end
# ---------------------------------------------------------------------------- #
function tostring(x::GEMPacket)
    return "$(x.click_time), " * join(x.async, ", ") * ", $(x.adjust)"
end
# ---------------------------------------------------------------------------- #
function read_run(file::GEMDataFile, run::Integer=1)
    hdr, offset = read_run_header(file, run)
    if offset > 0
        data = Vector{GEMPacket}(undef,file.nwin)
        open(file.path, "r") do io
            seek(io, offset)

            for k = 1:file.nwin
                data[k] = read_packet(io)
            end
            return hdr, data
        end
    else
        return hdr, Vector{GEMPacket}()
    end
end
# ---------------------------------------------------------------------------- #
"""
returns the run header and the position of the start of the data block for that
run
"""
function read_run_header(file::GEMDataFile, k::Integer=1)
    @assert(k <= file.nrun, "Invalid run number!")

    offset = file.idx_map[k]
    println("Reading header for run ", k, " at offset ", offset)

    if offset > 0
        open(file.path, "r") do io
            seek(io, offset)
            hdr_length = read(io, UInt64)
            println("Header length: ", hdr_length)

            hdr = JSON.parse(String(read(io, hdr_length)))
            return hdr, position(io)
        end
    else
        return Dict(), -1
    end

end
# ---------------------------------------------------------------------------- #
function read_packet(io::IOStream)
    return GEMPacket(
        read(io, UInt8),            #data-transfer-protocol id
        read(io, UInt16),           #window number
        read(io, UInt32),           #click_time
        read!(io, Vector{Int16}(undef, MAX_TAPPERS)),#asynchronies
        read(io, Int16)             #next adjustment factor
    )
end
# ---------------------------------------------------------------------------- #
end # end module

max_stars_repo_name: janatalab/GEM

max_issues_repo_name: janatalab/GEM

max_forks_repo_name: janatalab/GEM


text:
# make sure docs is working directory before calling
push!(LOAD_PATH,"..")
using Pkg; Pkg.activate("."); Pkg.instantiate()
using Documenter, LTspice, Printf

# Build docs.
# ===========

const PAGES = [
  "Home" => "index.md",
  "Installation" => "install.md",
  "Quickstart" => "quickstart.md",
  "Examples" => "examples.md",
  "Public API" => "public_api.md",
  "Contents" => "contents.md",
]

makedocs(
  modules = [LTspice],
  clean = true,
  doctest   = "doctest" in ARGS,
  linkcheck = "linkcheck" in ARGS,
  format    = "pdf" in ARGS ? :latex : :html,
  build = "site",
  sitename = "LTspice",
  authors = "Chris Stook",
  pages = PAGES,
  html_prettyurls = "deploy" in ARGS,
)

if "deploy" in ARGS
  fake_travis = "C:/Users/Chris/fake_travis_LTspice.jl"
  if isfile(fake_travis)
    include(fake_travis)
  end
  deploydocs(
    repo = "github.com/cstook/LTspice.jl.git",
    target = "site",
    branch = "gh-pages",
    latest = "master",
    osname = "linux",
    julia  = "1.0",
    deps = nothing,
    make = nothing,
  )
end


text:
#function _deepcopy(vtud::VTUDataField{T}) where {T<:Number}
#	len = length(vtud.dat)
#	dat = Vector{T}(undef,len)
#	for i = 1:len
#		dat[i] = vtud.dat[i]
#	end
#	return VTUDataField(dat)
#end

function Base.deepcopy(vtud::VTUData)
	ipdat = vtud.interp_data
	#ip = GC.@preserve ipdat map(deepcopy,ipdat)
	ip = map(deepcopy,ipdat)
	return VTUData(vtud.names,vtud.header,deepcopy(vtud.data),ip,vtud.idat)
end

function fill_zeros(arr::Base.ReinterpretArray{T,N,S,A,IsReshaped}) where {T,N,S,A,IsReshaped}
	fill!(arr.parent,zero(S))
	return nothing
end

function fill_zeros(arr::VTUDataField{T}) where {T}
	fill!(arr,zero(T))
	return nothing
end

#function fill_zeros(arr::Vector{T}) where {T<:Number}
#	fill!(arr,zero(T))
#	return nothing
#end

function Base.similar(vtud::VTUData)
	ret = deepcopy(vtud)
	for dat in ret.interp_data
		fill_zeros(dat)
	end
	return ret
end

function Base.fill!(ret::VTUDataField{T}, c::T) where {T}
	@inbounds for i in 1:length(ret.dat)
		ret.dat[i] = c
	end
	return nothing
end

function Base.fill!(ret::VTUData, c::Float64)
	for dat in ret.interp_data
		fill!(dat,c)
	end
	return nothing
end

function Base.zero(vtud::VTUData)
	ret = similar(vtud)
	fill!(ret,0.0)
	return nothing
end

function update_data!(vtud::VTUData)
	for (i,j) in enumerate(vtud.idat)
		vtud.data[j] = vtud.interp_data[i]
	end
end

max_stars_repo_name: nagelt/VTUFileHandler

max_issues_repo_name: nagelt/VTUFileHandler

max_forks_repo_name: nagelt/VTUFileHandler


text:
## functions to indentify lowest and highest values by boolean
##
## useful together with TimeData.getVars

function islowest(x::Array{Float64, 1}, n::Integer = 1)
    nVals = length(x)
    inds = sortperm(x)
    logics = falses(nVals)
    logics[inds[1:n]] = true
    return logics
end

function ishighest(x::Array{Float64, 1}, n::Integer = 1)
    nVals = length(x)
    inds = sortperm(x, rev=true)
    logics = falses(nVals)
    logics[inds[1:n]] = true
    return logics
end

function createTestData()
    # create test data

    # get test data from EconDatasets
    intRates = []
    try
        intRates = dataset("Treasuries")
    catch
        getDataset("Treasuries")
        intRates = dataset("Treasuries")
    end

    # create subset
    dats = Date(1991, 1, 2):Date(2013, 12, 31)
    intData = intRates[dats, 4]
    names!(intData.vals, [:Y1])

    # save subset as test data to make it available to R
    fname = joinpath(Pkg.dir("Econometrics"), "test/data/no_git_intRates.csv")
    writeTimedata(fname, intData)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiCoreV1ContainerPort <: SwaggerModel
    containerPort::Any # spec type: Union{ Nothing, Int32 } # spec name: containerPort
    hostIP::Any # spec type: Union{ Nothing, String } # spec name: hostIP
    hostPort::Any # spec type: Union{ Nothing, Int32 } # spec name: hostPort
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    protocol::Any # spec type: Union{ Nothing, String } # spec name: protocol

    function IoK8sApiCoreV1ContainerPort(;containerPort=nothing, hostIP=nothing, hostPort=nothing, name=nothing, protocol=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1ContainerPort, Symbol("containerPort"), containerPort)
        setfield!(o, Symbol("containerPort"), containerPort)
        validate_property(IoK8sApiCoreV1ContainerPort, Symbol("hostIP"), hostIP)
        setfield!(o, Symbol("hostIP"), hostIP)
        validate_property(IoK8sApiCoreV1ContainerPort, Symbol("hostPort"), hostPort)
        setfield!(o, Symbol("hostPort"), hostPort)
        validate_property(IoK8sApiCoreV1ContainerPort, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(IoK8sApiCoreV1ContainerPort, Symbol("protocol"), protocol)
        setfield!(o, Symbol("protocol"), protocol)
        o
    end
end # type IoK8sApiCoreV1ContainerPort

const _property_map_IoK8sApiCoreV1ContainerPort = Dict{Symbol,Symbol}(Symbol("containerPort")=>Symbol("containerPort"), Symbol("hostIP")=>Symbol("hostIP"), Symbol("hostPort")=>Symbol("hostPort"), Symbol("name")=>Symbol("name"), Symbol("protocol")=>Symbol("protocol"))
const _property_types_IoK8sApiCoreV1ContainerPort = Dict{Symbol,String}(Symbol("containerPort")=>"Int32", Symbol("hostIP")=>"String", Symbol("hostPort")=>"Int32", Symbol("name")=>"String", Symbol("protocol")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1ContainerPort }) = collect(keys(_property_map_IoK8sApiCoreV1ContainerPort))
Swagger.property_type(::Type{ IoK8sApiCoreV1ContainerPort }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1ContainerPort[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1ContainerPort }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1ContainerPort[property_name]

function check_required(o::IoK8sApiCoreV1ContainerPort)
    (getproperty(o, Symbol("containerPort")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1ContainerPort }, name::Symbol, val)
end


text:
## Optimization Model for dtu42380a1
## Mar 11, 2020

f = open("prod_plan.csv","w");

println(f, ",,Demand,,,Production,,,Inventory,,,Capacity,,,,Unmet demand,,,");
println(f, "DC,Period,Product 1,Product 2,Product 3,Product 1,Product 2,Product 3,Product 1,Product 2,Product 3,Blend (RT),Blend (OT), Pack(RT), Pack(OT),Product 1,Product 2,Product 3");

for i = 1:L
    for t = 1:T
        print(f,i,",",t,",");
        for p = 1:P
            print(f, D[i,p,t],",");
        end
        for p = 1:P
            print(f,value(x[i,p,t]),",");
        end
        for p = 1:P
            print(f,value(I[i,p,t]),",");
        end
        print(f, sum(value(mBR[i,p,t]) for p = 1:P) ,",");
        print(f, sum(value(mBO[i,p,t]) for p = 1:P) ,",");
        print(f, sum(value(mPR[i,p,t]) for p = 1:P) ,",");
        print(f, sum(value(mPO[i,p,t]) for p = 1:P) ,",");
        for p = 1:P
            print(f,value(delta[i,p,t]),",");
        end
        println(f,"");
    end
end

close(f);


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    TaitBryanIntr

Intrinsic right-handed rotation by the ZXY axes.
"""
struct TaitBryanIntr <: RotationConvention end
axesseq(::Type{TaitBryanIntr}) = :ZXY
orientation(::Type{TaitBryanIntr}) = (:CCW,:CCW,:CCW)
angleunits(::Type{TaitBryanIntr}) = :RAD
mainaxis(::Type{TaitBryanIntr}) = :X
isextrinsic(::Type{TaitBryanIntr}) = false


text:
using SPRING2
using Test


text:
@testset "496.next-greater-element-i.jl" begin
    @test next_greater_element([4,1,2], [1,3,4,2]) == [-1,3,-1]
    @test next_greater_element([2, 4], [1, 2, 3, 4]) == [3, -1]
end


text:
#=
Treeview Functions
=#

#=
Callback
=#
function unit_select_cb(w::Ptr,p1,p2,user_data::Tuple{Single_Channel})

    sc, = user_data
    select_unit(sc)
end

function select_unit(sc::Single_Channel)
    clus=get_cluster_id(sc)

    old_clus=sc.buf.selected_clus

    sc.buf.selected_clus=clus
    if clus>0
        mytol=sc.temp.tol[clus]
        setproperty!(sc.adj_sort, :value, mytol)
    end

    ctx=Gtk.getgc(sc.c3)

    if old_clus>0
        (x1_i,x2_i,y1_i,y2_i)=get_template_dims(sc,old_clus)
        draw_box(x1_i,y1_i,x2_i,y2_i,(0.0,0.0,0.0),2.0,ctx)
        draw_box(x1_i,y1_i,x2_i,y2_i,(1.0,1.0,1.0),1.0,ctx)
    end

    if sc.buf.selected_clus>0
        (x1_f,x2_f,y1_f,y2_f)=get_template_dims(sc,sc.buf.selected_clus)
        draw_box(x1_f,y1_f,x2_f,y2_f,(1.0,0.0,1.0),1.0,ctx)
    end

    nothing
end

function get_cluster_id(sc::Single_Channel)
    selmodel=Gtk.GAccessor.selection(sc.sort_tv)
    iter=Gtk.selected(selmodel)

    myind=parse(Int64,Gtk.get_string_from_iter(TreeModel(sc.sort_list), iter))
end

function update_treeview(sc::Single_Channel)

    for i=length(sc.sort_list):-1:2
        deleteat!(sc.sort_list,i)
    end

    for i=1:sc.total_clus
        push!(sc.sort_list,(i,))
    end

    selmodel=Gtk.GAccessor.selection(sc.sort_tv)
    Gtk.select!(selmodel, Gtk.iter_from_index(sc.sort_list,1))

    nothing
end

function is_selected(store,tv,ind)
    iter=Gtk.iter_from_string_index(store,string(ind))
    selection=Gtk.GAccessor.selection(tv)
    ccall((:gtk_tree_selection_iter_is_selected, Gtk.libgtk),Bool,
    (Ptr{Gtk.GObject}, Ptr{Gtk.GtkTreeIter}),selection, Gtk.mutable(iter))
end


text:
module ParticleSwarmAlignment

export PSO_MSA, score_sequences, generate_sequences, global_align, progressive_alignment_inorder

using Random
using Printf
using Dates

# Random module must be seeded with a constant in order to have
# reproducible test results
Random.seed!(1)

# Debug printf function which can be turned on/off with debug flag
debug = false
if debug
    function dprintf(s::String, args...)
        @eval @printf($s, $(args...))
    end
else
    function dprintf(s, args...)
        # nothing
    end
end

# Needleman-Wunsch Global Alignment using fixed penalties.
# Extending with keyword arguments for optional preallocated s, b matrices
# will drastically improve performance during its use in progressive alignment
function global_align(v, w, match_penalty=1, mismatch_penalty=-1, deletion_penalty=-1)
    n1 = length(v)
    n2 = length(w)
    #if !use_preallocated_matrices
    s = zeros(Float64, n1+1, n2+1)
    b = zeros(Float64, n1+1, n2+1)
    #end

    for i in 1:(n1+1)
        s[i,1] = (i-1) * deletion_penalty
        b[i,1] = 2
    end
    for j in 1:(n2+1)
        s[1,j] = (j-1) * deletion_penalty
        b[1,j] = 3
    end

    for i in 2:(n1+1)
        for j in 2:(n2+1)
            if v[i-1] == w[j-1]
                ms = s[i-1,j-1] + match_penalty
            else
                # ignore cases where a letter is paired with a gap
                # do not consider this a mismatch
                # if v[i-1] != '-' && w[j-1] != '-'
                #     ms = s[i-1,j-1] + mismatch_penalty
                # else
                #     # if a letter is paired with a gap, add no penalty
                #     ms = s[i-1,j-1] #+ match_penalty # + 0.5 * mismatch_penalty
                # end
                ms = s[i-1,j-1] + mismatch_penalty
            end
            test = [ms, s[i-1,j] + deletion_penalty, s[i,j-1] + deletion_penalty]
            p = argmax(test)
            s[i,j] = test[p]
            b[i,j] = p
        end
    end

    i = n1+1
    j = n2+1
    sv = []
    sw = []
    while(i > 1 || j > 1)
        p = b[i,j]
        if (p == 1)
            i = i-1
            j = j-1
            push!(sv, v[i])
            push!(sw, w[j])
        elseif p == 2
            i=i-1
            push!(sv, v[i])
            push!(sw, "-")
        elseif p == 3
            j = j-1
            push!(sv, "-")
            push!(sw, w[j])
        else
            break
        end
    end

    return (s[n1+1,n2+1], join(reverse(sv)), join(reverse(sw)))
end

# Generates t sequences all of length l, and returns an array of strings
function generate_sequences(t::Int64, l::Int64)
    # t is the number of sequences to create
    # l is the length of the sequences
    DNA = Array{String,1}(undef,0)
    base_arr = ["A", "T", "G", "C"]

    for t_index in 1:t
        push!(DNA, "")
        for l_value in 1:l
            r = convert(Int64, floor(Random.rand() * 4) + 1)
            DNA[t_index] = string(DNA[t_index], base_arr[r])
        end
    end
    return DNA
end

# Performs an *ordered* equals comparison of two ordered arrays
function array_equals_ordered(A::Array, B::Array)
    if length(A) != length(B)
        return false
    end

    for i in 1:length(A)
        if A[i] != B[i]
            return false
        end
    end
    return true
end

# Determines whether Any element e is contained within Any array A.
# Returns true if so.
function array_contains(A::Array, e)
    for i in 1:length(A)
        if array_equals_ordered(A[i], e)
            return true
        end
    end
    return false
end

# Creates num random permutations of the elements of A. If must_be_unique is left
# as true, then it ensures that the permutations returned are unique. This can fail
# if it takes too long to find a unique permutation, as the method randomly searches
# for permutations instead of iterating through them in order.
function random_permutations(A::Array{String,1}, num::Int64, must_be_unique::Bool=true)::Array{Array{String,1},1}
    # Return num random permutations of the elements of A
    results = Array{Array{String,1},1}(undef,0)
    if num >= factorial(big(length(A)))
        throw(ErrorException("num is too high"))
    end
    max_check = 100
    try_count = 0
    for i in 1:num
        localA = deepcopy(A)
        randomA = shuffle!(localA)
        if must_be_unique
            while array_contains(results, randomA)
                randomA = shuffle!(A)
                try_count += 1
                if try_count > max_check
                    msg = @sprintf("Exceeded max tries (%d) while trying to find unique permutations\n", max_check)
                    throw(ErrorException(msg))
                end
            end
        end
        push!(results, randomA)
    end
    return results
end

# Returns a swap sequence [(idx1, idx2), ...] to turn A into B
# Not necessarily the optimal swap sequence, of the least number of swaps
function get_swap_sequence(A_in::Array{String,1}, B_in::Array{String,1})
    if length(A_in) != length(B_in)
        throw(ErrorException("Length of A and B must be the same"))
    end
    dprintf("Getting swaps between:\n%s,\n%s\n", string(A_in), string(B_in))
    A = Array{String,1}(undef,length(A_in))
    B = Array{String,1}(undef,length(B_in))
    for i in 1:length(A_in)
        A[i] = A_in[i]
        B[i] = B_in[i]
    end

    swap_sequence = Tuple{Int,Int}[]

    for i in 1:length(A)
        if A[i] != B[i]
            # We want to turn A into B
            # so find where B[i] is in A[i:end], and swap that element in A with the element at A[i]
            a_index = findfirst(x -> x == B[i], A[i:end])
            if a_index === nothing
                msg = string("Could not find ", string(B[i], string(" in ", string(A[i:end]))))
                throw(ErrorException(msg))
            end
            a_index += i - 1

            temp = A[i]
            A[i] = A[a_index]
            A[a_index] = temp

            push!(swap_sequence, (i, a_index))
        end
    end
    equals_check = array_equals_ordered(A, B)
    if !equals_check
        throw(ErrorException("A and B were not the same:\n%s\n%s\n", string(A), string(B)))
    end
    return swap_sequence
end

# Calculates a new "velocity" for the particle position given global best position,
# and the local best position. Alpha and Beta are tuning parameters which determine
# how quickly the particle trends towards the local/global best.
function pso_particle_velocity(particle::Array{String,1}, global_best::Array{String,1}, local_best::Array{String,1}, alpha::Float64, beta::Float64)
    dprintf("current_position: %s\n", string(particle))
    dprintf("local_best: %s\n", string(local_best))
    dprintf("global_best: %s\n", string(global_best))
    swaps_to_local_best = get_swap_sequence(particle, local_best)
    swaps_to_global_best = get_swap_sequence(particle, global_best)
    dprintf("swaps_to_global_best: %s\n", swaps_to_global_best)
    dprintf("swaps_to_local_best: %s\n", swaps_to_local_best)
    velocity = Tuple{Int,Int}[]

    for lswap in swaps_to_local_best
        r = Random.rand()
        if r < alpha
            push!(velocity, lswap)
        end
    end
    for gswap in swaps_to_global_best
        r = Random.rand()
        if r < beta
            push!(velocity, gswap)
        end
    end
    return velocity
end

# This takes a particle position, and performs the swaps specified by
# swap_sequence_velocity, in order. Returns the resulting position array.
function apply_velocity_to_particle(position_in::Array{String,1}, swap_sequence_velocity::Array)::Array{String,1}
    position = deepcopy(position_in)
    for (idxA, idxB) in swap_sequence_velocity
        temp = position[idxA]
        position[idxA] = position[idxB]
        position[idxB] = temp
    end
    return position
end

# Determines where gaps need to be inserted into old to obtain the new string value.
# Returns the indexes.
function get_new_gap_indexes(old::String, new::String)::Array{Int64,1}
    old_local = old
    indexes = Int[]
    for i in 1:length(new)
        if (new[i] == '-' && i == (length(old_local) + 1)) || new[i] == '-' && old_local[i] != '-'
            #@printf("get_new_gap_indexes, new gap at %d, old: %s, new: %s\n", i, old_local, new)
            push!(indexes, i)
            old_local = string(string(old_local[1:i-1], '-'), old_local[i:end])
        elseif new[i] == old_local[i]
            continue
        else
            throw(ErrorException("Sequences ignoring gaps were not the same"))
        end
    end
    return indexes
end

# Determines where gaps need to be inserted into old to obtain the new string value.
# Returns the indexes.
# This is an optimization of get_new_gap_indexes which does incur large memory copies
# due to string copying.
function get_new_gap_indexes2(old::String, new::String)::Array{Int64,1}
    indexes = Int[]
    old_i = 1
    new_i = 1
    while new_i <= length(new)
        # if new string had a gap either matched with an old non-gap, or past the end of the old string
        if (new[new_i] == '-' && (old_i > length(old) || old[old_i] != '-'))
            # new gap
            push!(indexes, new_i)
            new_i += 1
        elseif (new[new_i] == old[old_i])
            new_i += 1
            old_i += 1
        else
            @printf("Unknown error old[%d] = %s, new[%d] = %s\n", old_i, string(old[old_i]), new_i, string(new[new_i]))
            throw(ErrorException("unknown error!"))
        end
    end
    return indexes
end

# Inserts gap characters into the sequences specified by the indexes in seq_indexes.
# The gaps are inserted at the indexes within each specified sequence at positions gap_indexes.
function insert_gaps_at(sequences::Array, seq_indexes::Array{Int64,1}, gap_indexes::Array{Int64,1})
    if length(seq_indexes) > 0 && length(gap_indexes) > 0
        gap_indexes = sort(gap_indexes)

        for sidx in 1:length(seq_indexes)
            seq_index = seq_indexes[sidx]
            s = sequences[seq_index]

            for i in gap_indexes
                # insert a '-' character at index i
                s = string(string(s[1:i-1], '-'), s[i:end])
            end
            # put the modified string back in the array
            sequences[seq_index] = s
        end
    end
end

# k is the length of the sequences, t is the number of sequences, and A is the sequences.
# The sequences of A must be the same length, or k must be the length of the longest sequence.
function Make_Profile(k::Int, t::Int, A::Array{String,1})::Dict{String,Array{Int64,1}}
    idxA = zeros(Int64, k)
    idxC = zeros(Int64, k)
    idxG = zeros(Int64, k)
    idxT = zeros(Int64, k)
    idxGap = zeros(Int64, k)

    for index in range(1, length=k)
        for seq in A
            if seq[index] == 'A'
                idxA[index] += 1
            elseif seq[index] == 'C'
                idxC[index] += 1
            elseif seq[index] == 'G'
                idxG[index] += 1
            elseif seq[index] == 'T'
                idxT[index] += 1
            elseif seq[index] == '-'
                idxGap[index] += 1
            else
                msg = @sprintf("character [%s] not allowed in DNA sequence", string(seq[index]))
                throw(ErrorException(msg))
            end
        end
    end

    return Dict(
        "A" => idxA,
        "C" => idxC,
        "G" => idxG,
        "T" => idxT,
        "-" => idxGap
    )
end

# Score the sequences in A by obtaining the consensus score using Make_Profile
function score_sequences(A::Array{String,1})::Int
    k = length(A[1])    # how long are the sequences
    t = length(A)       # how many sequences
    profile = Make_Profile(k, t, A)
    # score is the sum of the maximum occurring letter (not including gaps) in each position
    score = 0
    for i in 1:k
        a = profile["A"][i]
        t = profile["T"][i]
        g = profile["G"][i]
        c = profile["C"][i]
        score += max(a, t, g, c)
    end
    return score
end

# Performs a comparison of A and B, but ignoring the gap characters in each.
function sequence_equals_ignore_gaps(A::String, B::String)
    A = filter(x -> x != '-', A)
    B = filter(x -> x != '-', B)
    return A == B
end

# This performs a progressive alignment of sequences, using the elements of the edges array
# to determine the order with which to perform neighbor alignments.
# The edges array must contain the pairwise edges and weights for all consecutive
# elements of the sequences array.
function progressive_alignment_inorder(sequences::Array, edges::Array{Tuple{String,String,Int},1})::Array{String,1}
    if length(edges) != length(sequences) - 1
        msg = @sprintf("Expected %d edge weights, got %d", length(sequences) - 1, length(edges))
        throw(ErrorException(msg))
    end

    aligned_sequences = Array{String,1}(undef,length(sequences))

    sets = Array{Array{Int64,1},1}(undef,0)

    for i in 1:length(sequences)
        #seq in sequences
        seq = sequences[i]
        aligned_sequences[i] = seq
        push!(sets, [i])
    end
    dprintf("Edges: %s\n", string(edges))

    local_edges = sort(edges, by=x -> x[3], rev=true) # sort desc by weight

    for i in 1:length(edges)
        min_edge = pop!(local_edges)

        dprintf("min_edge: %s\n", string(min_edge))

        # use aligned sequences here instead of original sequences
        #idxA = findfirst(s -> sequence_equals_ignore_gaps(s, min_edge[1]), sequences)
        idxA = findfirst(s -> s == min_edge[1], sequences)
        #idxB = findfirst(s -> sequence_equals_ignore_gaps(s, min_edge[2]), sequences)
        idxB = findfirst(s -> s == min_edge[2], sequences)
        A = aligned_sequences[idxA]
        B = aligned_sequences[idxB]
        dprintf("\n\nPerforming global alignment of %d %d\n%s\nand\n%s\nedge_weight: %d\n", idxA, idxB, A, B, min_edge[3])
        score, alignedA, alignedB = global_align(A, B)
        dprintf("Got alignment\n%s\nand\n%s\n", alignedA, alignedB)

        # take new_gap_indexes_A, and insert same gaps into all current
        # predecessors of A.  And same for B
        new_gap_indexes_A = get_new_gap_indexes2(A, alignedA)
        seqs_to_update = [x for x in sets[idxA] if x != idxA]
        insert_gaps_at(aligned_sequences, seqs_to_update, new_gap_indexes_A)

        new_gap_indexes_B = get_new_gap_indexes2(B, alignedB)
        seqs_to_update = [x for x in sets[idxB] if x != idxB]
        insert_gaps_at(aligned_sequences, seqs_to_update, new_gap_indexes_B)

        aligned_sequences[idxA] = alignedA
        aligned_sequences[idxB] = alignedB

        # Update predecessors for A and B
        # TODO avoid the loop somehow by reusing the sets array from one of the existing idxA or idxB elements
        newset = Int[]#[idxA, idxB]
        push!(newset, sets[idxA]...)
        push!(newset, sets[idxB]...)
        #@printf("newset: %s\n", string(newset))
        sets[idxA] = newset
        sets[idxB] = newset
        for p in newset
            sets[p] = newset
        end
    end
    return aligned_sequences
end

# Searches the edges array to find the edge between nodeA and nodeB.
# Returns the edge Tuple{String,String,Int}, as well as the index it was found at.
function get_edge_between(nodeA::String, nodeB::String, edges::Array)
    #nodeA = filter(x -> x != '-', nodeA)
    #nodeB = filter(x -> x != '-', nodeB)
    for idx in 1:length(edges)
        edge = edges[idx]
        if edge[1] == nodeA && edge[2] == nodeB
            return (edge, idx)
        elseif edge[1] == nodeB && edge[2] == nodeA
            return (edge, idx)
        end
    end
    return nothing
end

# Holds all relevant information for a particle in the swarm
mutable struct Particle
    sequences::Array{String,1}              # initial order/position
    aligned_sequences::Array{String,1}      # hold the aligned sequences to avoid unnecessary recomputes
    position::Array{String,1} # ordering    # current order/position
    best_sequences::Array{String,1}         # best order/position
    best_score::Int                         # best score
    pidx::Int                               # particle index (for reverse mapping into the particle array)
    edges::Array{Tuple{String,String,Int},1}# edges for the current position
end

# Utility function for printing the more relevant information about a particle for debugging.
function print_particle(p::Particle)
    print("Particle: ")
    println(p.pidx)
    print("  position:\n")
    for seq in p.position
        print("    ")
        println(seq)
    end
    print("  best_sequences:\n")
    for seq in p.best_sequences
        print("    ")
        println(seq)
    end
    print("  best_score: ")
    println(p.best_score)
    print("  edges:")
    for edge in p.edges
        print(" ")
        print(edge[1])
    end
    println()
end

# The primary function of the module, which takes the sequences, the number of swarm iterations to run, and
# optionally the number of particles to use. If no num_particles is provided, then a reasonable value is
# calculated based on the number of sequences.
# The length of the sequences do not all need to be the same, however the performance may reduce if they
# are not because additional work must be done to initially align neighbors before computing edge weights
function PSO_MSA(sequences::Array{String,1}, iterations::Int, num_particles::Int=0)
    # Number of sequences
    t = length(sequences)
    # Number of permutations ignoring complete reversals
    solution_space = factorial(big(t)) / 2

    if num_particles == 0
        if t >= 8
            num_particles = 100
        else
            #convert(Int64, t*(t-1)/2) # num edges over 2
            num_particles = t
        end
    end
    @printf("Number of Particles: %d\n", num_particles)
    search_space = num_particles * iterations
    search_solution_ratio = search_space / solution_space

    @printf("Solution Space: %.02f\n", solution_space)
    @printf("Search Space: %.02f\n", search_space)
    @printf("Search/Solution space: %.010f\n", search_solution_ratio)

    nodes = Array{String,1}(undef,0)
    for i in 1:length(sequences)
        push!(nodes, sequences[i])
    end

    permutations = random_permutations(nodes, num_particles)

    # Initialize particles list from the generated permutations
    println("Initializing particles")
    particles = Array{Particle,1}(undef,length(permutations))
    for pidx in 1:length(permutations)
        sequences = copy(permutations[pidx])
        aligned_sequences = copy(sequences)
        @printf("Initializing particle %d\n", pidx)
        # Initial edges
        p_edges = Array{Tuple{String,String,Int},1}(undef,0)
        for sidx in 2:length(sequences)
            seqA_orig = sequences[sidx-1]
            seqB_orig = sequences[sidx]
            #scoreAB = score_sequences([seqA, seqB])
            # if length(seqA_orig) != length(seqB_orig)
            #     (align_score, seqA, seqB), t, bytes, gct, memalloc = @timed global_align(seqA_orig, seqB_orig)
            #     println(t)
            # else
            #     seqA, seqB = seqA_orig, seqB_orig
            # end
            # scoreAB = score_sequences([seqA, seqB])

            # Here the initial score for each neighbor distance is set to zero
            # to avoid having to align them to obtain a valid distance score
            scoreAB = 0
            push!(p_edges, (seqA_orig, seqB_orig, scoreAB))
        end

        @time aligned_sequences = progressive_alignment_inorder(sequences, p_edges)

        initial_best_score = score_sequences(aligned_sequences)
        p = Particle(
            permutations[pidx],         # original sequences/position
            aligned_sequences,          # aligned_sequences
            copy(permutations[pidx]),   # current_position
            copy(permutations[pidx]),   # best_sequences
            initial_best_score,         # initial best score
            pidx,                       # index
            p_edges)                    # edges for the current position
        particles[pidx] = p
        if debug
            print_particle(particles[pidx])
        end
    end

    #alpha = Random.rand()
    #beta = Random.rand()

    # Here are two important tuning values. These determine how rapidly the particles trend
    # towards their own past best and the swarm's best solution. With some non-systematic testing
    # it was determined that 0.2 is a reasonable value for both. Too high of a value results
    # in particles not searching the area around the local optimum enough. Too low of a value results
    # in not enough motion happening.
    alpha = 0.2
    beta = 0.2

    # keep track of the local best score+ordering for each particle
    local_best_scores = [p.best_score for p in particles]
    # keep track of the global best particle
    global_best_particle = particles[argmax(local_best_scores)]
    println("Initial Best")
    @printf("Global best (score = %d, pidx = %d):\n", global_best_particle.best_score, global_best_particle.pidx)
    for seq in global_best_particle.sequences
        println(seq)
    end
    @printf("Initial global best score: %d\n", global_best_particle.best_score)

    iteration_metrics = Millisecond[]

    # These take a potential solution to swarm stagnation. Every so often, we can reassign
    # a certain number of particles in the swarm to a random position.
    # Setting random_subgroup_prob to -1 disables this. We have disabled it due to
    # it introducing too much randomness into the system.
    random_subgroup_move_max = 10
    random_subgroup_move_count = 0
    random_subgroup_prob = -1

    # Percentage of particles which if they don't move, move them randomly
    # This is also a significant tuning value. Setting it to one means that any time a particle
    # doesn't move for the past iteration, immediately pick it up and give it a new position.
    #reshuffle_threshold = length(particles) * 0.5
    reshuffle_threshold = 1

    for iteration in 1:iterations

        unmoved_particles = []

        @printf("Beginning iteration %d\n", iteration)
        start = Dates.now()

        if random_subgroup_move_count >= random_subgroup_move_max
            do_random_subgroup_move = true
            random_subgroup_move_count = 0
        else
            random_subgroup_move_count += 1
            do_random_subgroup_move = false
        end

        for pidx in 1:length(particles)
            particle = particles[pidx]
            if debug
                println("Looking at particle:")
                print_particle(particles[pidx])
            end

            new_velocity = nothing

            if do_random_subgroup_move
                if Random.rand() < random_subgroup_prob
                    random_position = random_permutations(particle.position, 1)[1]
                    dprintf("Randomly moving particle from %s to %s\n", string(particle.position), string(random_position))
                    new_velocity = pso_particle_velocity(
                        random_position,
                        global_best_particle.best_sequences,
                        particle.best_sequences,
                        alpha, beta)
                end
            end
            if new_velocity === nothing
                new_velocity = pso_particle_velocity(
                    particle.position,
                    global_best_particle.best_sequences,
                    particle.best_sequences,
                    alpha, beta)
            end

            # This particle did not move in this iteration
            if length(new_velocity) == 0
                push!(unmoved_particles, pidx)
            end

            if length(new_velocity) > 0
                # Apply the velocity to the position of the particle
                new_position = apply_velocity_to_particle(particle.position, new_velocity)
                dprintf("Applied velocity %s to particle %d, moved from\n%s, to %s\n",
                    string(new_velocity),
                    particle.pidx,
                    particle.position,
                    new_position)
                particle.position = new_position

                # Perform the same swaps to the aligned sequences
                particle.aligned_sequences = apply_velocity_to_particle(particle.aligned_sequences, new_velocity)

                # Calculate the new edges for the current position (sequence ordering)
                p_edges = Array{Tuple{String,String,Int},1}(undef,length(particle.position)-1)
                for sidx in 2:length(particle.position)
                    seqA_orig = particle.position[sidx-1]
                    seqB_orig = particle.position[sidx]
                    seqA_aligned = particle.aligned_sequences[sidx-1]
                    seqB_aligned = particle.aligned_sequences[sidx]

                    scoreAB = score_sequences([seqA_aligned, seqB_aligned])
                    p_edges[sidx-1] = (seqA_orig, seqB_orig, scoreAB)
                end
                particle.edges = p_edges

                # Perform a progressive alignment of the new position with the new edges
                aligned_sequences = progressive_alignment_inorder(particle.position, particle.edges)
                score = score_sequences(aligned_sequences)

                particle.aligned_sequences = aligned_sequences

                # if the new position score is better, update the local best
                if score > particle.best_score
                    dprintf("Updating particle best score to %d\n", score)
                    particle.best_score = score
                    particle.best_sequences = copy(particle.position)
                end

                # if particle.best_score > global_best_particle.best_score
                #     @printf("Updating global best score to %d\n", particle.best_score)
                #     global_best_particle = deepcopy(particle)
                # end
            end

            # Update the global best at the end of the round
            particle_best_scores = [p.best_score for p in particles]
            best_particle = particles[argmax(particle_best_scores)]
            if best_particle.best_score > global_best_particle.best_score
                @printf("Updating global best score to %d\n", best_particle.best_score)
                global_best_particle = deepcopy(best_particle)
            end
        end

        if length(unmoved_particles) >= reshuffle_threshold
            dprintf("Reshuffling particles: %s\n", string(unmoved_particles))
            for reshuffle_pidx in unmoved_particles
                p = particles[reshuffle_pidx]
                p.position = random_permutations(p.position, 1)[1]
            end
        end

        end_time = Dates.now()
        duration = end_time - start
        push!(iteration_metrics, duration)
        @printf("Iteration %d took %d milliseconds\n", iteration, duration.value)
        if debug
            @printf("Global best (score = %d, pidx = %d):\n",
                global_best_particle.best_score,
                global_best_particle.pidx)
            for seq in global_best_particle.aligned_sequences
                println(seq)
            end
            println()
        end
    end

    println("FINISHED PARTICLE SWARM")

    # print the final results
    @printf("Global best (score = %d, pidx = %d):\n",
        global_best_particle.best_score,
        global_best_particle.pidx)
    for seq in global_best_particle.aligned_sequences
        println(seq)
    end
    println()

    return global_best_particle.best_score, global_best_particle.aligned_sequences
end

end # end module


text:
using Game2048
using BenchmarkTools
using StatsBase

board = initboard()
move!(board, left)

@time res = simulate_games(10_000);
res
@time mean(res)
@time extrema(res)

@benchmark simulate_game()


text:
using Plots
using KernelAbstractions: CPU
using ClimateMachine.MPIStateArrays: array_device
using ClimateMachine.BalanceLaws: Prognostic, Auxiliary, GradientFlux

"""
    plot_friendly_name(ϕ)

Get plot-friendly string, since many Unicode
characters do not render in plot labels.
"""
function plot_friendly_name(ϕ)
    s = ϕ
    s = replace(s, "ρ" => "rho")
    s = replace(s, "α" => "alpha")
    s = replace(s, "∂" => "partial")
    s = replace(s, "∇" => "nabla")
    return s
end

"""
    export_plot(
        z,
        time_data,
        dons_arr::Array,
        ϕ_all,
        filename;
        xlabel,
        ylabel,
        time_units = "[s]",
        round_digits = 2,
        horiz_layout = false,
        kwargs...
    )

Export plot of all variables, or all
available time-steps in `dons_arr`.
"""
function export_plot(
    z,
    time_data,
    dons_arr::Array,
    ϕ_all,
    filename;
    xlabel,
    ylabel,
    time_units = "[s]",
    round_digits = 2,
    sample_rate = 1,
    horiz_layout = false,
    kwargs...,
)
    ϕ_all isa Tuple || (ϕ_all = (ϕ_all,))
    single_var = ϕ_all[1] == xlabel || length(ϕ_all) == 1
    p = plot()
    sample = 1:sample_rate:length(time_data)
    for (t, data) in zip(time_data[sample], dons_arr[sample])
        for ϕ in ϕ_all
            ϕ_string = String(ϕ)
            ϕ_data = data[ϕ_string][:]
            t_label = "t=$(round(t, digits=round_digits)) $time_units"
            label = single_var ? t_label : "$ϕ_string, $t_label"
            args = horiz_layout ? (z, ϕ_data) : (ϕ_data, z)
            plot!(
                args...;
                xlabel = xlabel,
                ylabel = ylabel,
                label = label,
                kwargs...,
            )
        end
    end
    savefig(filename)
end

"""
    export_contour(
        z,
        time_data,
        dons_arr::Array,
        ϕ,
        filename;
        xlabel = "time [s]",
        ylabel = "z [m]",
        label = String(ϕ),
        kwargs...
    )

Export contour plots given
 - `z` Array of altitude. Note: this must not include duplicate nodal points.
 - `time_data` array of time data
 - `dons_arr` an array whose elements are populated by `dict_of_nodal_states`
 - `ϕ` variable to contour
 - `filename` file name to export to.
 - `xlabel` x-label
 - `ylabel` y-label
 - `label` contour labels
"""
function export_contour(
    z,
    time_data,
    dons_arr::Array,
    ϕ,
    filename;
    xlabel = "time [s]",
    ylabel = "z [m]",
    label = String(ϕ),
    kwargs...,
)
    ϕ_string = String(ϕ)
    ϕ_data = hcat([data[ϕ_string][:] for data in dons_arr]...)
    args = (time_data, z, ϕ_data)
    try
        contourf(
            args...;
            xlabel = xlabel,
            ylabel = ylabel,
            label = label,
            c = :viridis,
            kwargs...,
        )
        savefig(filename)
    catch
        @warn "Contour plot $label failed. Perhaps the field is all zeros"
    end
end

function save_binned_surface_plots(
    x,
    y,
    z,
    title,
    filename,
    n_plots = (3, 3),
    z_label_prefix = "z",
    n_digits = 5,
)
    n_z_partitions = prod(n_plots)
    z_min_global = min(z...)
    z_max_global = max(z...)
    Δz = (z_max_global - z_min_global) / n_z_partitions
    z_min = ntuple(i -> z_min_global + (i - 1) * Δz, n_z_partitions)
    z_max = ntuple(i -> z_min_global + (i) * Δz, n_z_partitions)
    p = []
    for i in 1:n_z_partitions
        mask = z_min[i] .<= z .<= z_max[i]
        x_i = x[mask]
        y_i = y[mask]
        sz_min = string(z_min[i])[1:min(n_digits, length(string(z_min[i])))]
        sz_max = string(z_max[i])[1:min(n_digits, length(string(z_max[i])))]
        p_i = plot(
            x_i,
            y_i,
            title = "$(title), in ($sz_min, $sz_max)",
            seriestype = :scatter,
            markersize = 5,
        )
        push!(p, p_i)
    end
    plot(p..., layout = n_plots, legend = false)
    savefig(filename)
end;

state_prefix(::Prognostic) = "prog_"
state_prefix(::Auxiliary) = "aux_"
state_prefix(::GradientFlux) = "grad_flux_"

"""
    export_state_plots(
        solver_config,
        dons_arr,
        time_data,
        output_dir;
        state_types = (Prognostic(), Auxiliary()),
        z = Array(get_z(solver_config.dg.grid)),
        sample_rate = 1,
        time_units = "[s]",
        ylabel = "z [m]",
        kwargs...
    )

Export line plots of states given
 - `solver_config` a `SolverConfiguration`
 - `dons_arr` an array of dictionaries, returned from `dict_of_nodal_states`
 - `time_data` an array of time values
 - `output_dir` output directory
"""
function export_state_plots(
    solver_config,
    dons_arr,
    time_data,
    output_dir;
    state_types = (Prognostic(), Auxiliary()),
    z = Array(get_z(solver_config.dg.grid)),
    sample_rate = 1,
    time_units = "[s]",
    ylabel = "z [m]",
    kwargs...,
)
    FT = eltype(solver_config.Q)
    mkpath(output_dir)
    for st in state_types
        vs = vars_state(solver_config.dg.balance_law, st, FT)
        for fn in flattenednames(vs)
            base_name = state_prefix(st) * replace(fn, "." => "_")
            file_name = joinpath(output_dir, "$(base_name).png")
            export_plot(
                z,
                time_data,
                dons_arr,
                (fn,),
                file_name;
                xlabel = fn,
                sample_rate = sample_rate,
                ylabel = ylabel,
                time_units = time_units,
                round_digits = 5,
                kwargs...,
            )
        end
    end
end

"""
    export_state_contours(
        solver_config,
        dons_arr,
        time_data,
        output_dir;
        state_types = (Prognostic(),),
        xlabel = "time [s]",
        ylabel = "z [m]",
        z = Array(get_z(solver_config.dg.grid; rm_dupes=true)),
        kwargs...
    )

Call `export_contour` for every
state variable given `state_types`.
"""
function export_state_contours(
    solver_config,
    dons_arr,
    time_data,
    output_dir;
    state_types = (Prognostic(),),
    xlabel = "time [s]",
    ylabel = "z [m]",
    z = Array(get_z(solver_config.dg.grid; rm_dupes = true)),
    kwargs...,
)
    FT = eltype(solver_config.Q)
    mkpath(output_dir)
    for st in state_types
        vs = vars_state(solver_config.dg.balance_law, st, FT)
        for fn in flattenednames(vs)
            base_name = state_prefix(st) * replace(fn, "." => "_")
            filename = joinpath(output_dir, "cnt_$(base_name).png")
            label = string(replace(fn, "." => "_"))
            args = (z, time_data, dons_arr, fn, filename)
            export_contour(
                args...;
                xlabel = xlabel,
                ylabel = ylabel,
                label = label,
                kwargs...,
            )
        end
    end
end


text:
videos = Dict(
              "Julia的安装" => ("468292801","BV1K5411m7tn","578865923"),
              "VSCode的安装" => ("213279400","BV17a411v7vt","579486498"),
              "WindowsTerminal的安装" => ("595763705","BV1rB4y1U7Qz","579887100"),
              "msys2的安装" => ("510800726","BV1xu411y7je","580583233"),
              "CLion安装及Fortran插件设置" => ("425850842","BV123411T798","581173295"),
              "IJulia-Jupyter安装" => ("683428323","BV1VS4y1c7yh","582013046"),
)


text:
@service S3

BUCKET_NAME = "aws-jl-test-issues---" * _now_formatted()

try
    S3.create_bucket(BUCKET_NAME)

    @testset "issue 223" begin
        # https://github.com/JuliaCloud/AWS.jl/issues/223
        body = "Hello World!"
        file_name = "contains spaces"

        try
            S3.put_object(BUCKET_NAME, file_name, Dict("body" => body))
            resp = S3.get_object(BUCKET_NAME, file_name)

            @test String(resp) == body
        finally
            S3.delete_object(BUCKET_NAME, file_name)
        end
    end

    @testset "issue 227" begin
        @testset "s3 public bucket" begin
            # https://github.com/JuliaCloud/AWS.jl/issues/227
            config = AWSConfig(; creds=nothing)
            resp = S3.get_object("www.invenia.ca", "index.html"; aws_config=config)

            @test !isempty(resp)
        end

        @testset "s3 private bucket" begin
            bucket_name = "aws-jl-test-issues---" * _now_formatted()
            file_name = "hello_world"

            try
                S3.create_bucket(bucket_name)
                S3.put_object(bucket_name, file_name)

                @test_throws AWSException S3.get_object(
                    bucket_name, file_name; aws_config=AWSConfig(; creds=nothing)
                )
            finally
                S3.delete_object(bucket_name, file_name)
                S3.delete_bucket(bucket_name)
            end
        end

        @testset "lambda" begin
            @service Lambda

            @test_throws NoCredentials Lambda.list_functions(;
                aws_config=AWSConfig(; creds=nothing)
            )
        end
    end

    @testset "issue 324" begin
        body = "Hello World!"
        file_name = "streaming.bin"

        try
            S3.put_object(BUCKET_NAME, file_name, Dict("body" => body))
            resp = S3.get_object(BUCKET_NAME, file_name)
            @test String(resp) == body

            # ERROR: MethodError: no method matching iterate(::Base.BufferStream)
            #   => BUG: header `response_stream` is pushed into the query...
            io = Base.BufferStream()
            S3.get_object(
                BUCKET_NAME,
                file_name,
                Dict("response_stream" => io, "return_stream" => true),
            )
            if bytesavailable(io) > 0
                @test String(readavailable(io)) == body
            else
                @test "no body data was available" == body
            end

        finally
            S3.delete_object(BUCKET_NAME, file_name)
        end
    end

    @testset "issue 466" begin
        file_name = "hang.txt"

        try
            S3.put_object(BUCKET_NAME, file_name)

            # The tests below validate the current behavior of how streams are handled.
            # Note: Avoid using `eof` for these tests can hang when using an unclosed `Base.BufferStream`

            stream = S3.get_object(BUCKET_NAME, file_name, Dict("return_stream" => true))
            if AWS.DEFAULT_BACKEND[] isa AWS.HTTPBackend
                @test !isopen(stream)
            else
                @test isopen(stream)
            end

            stream = Base.BufferStream()
            S3.get_object(BUCKET_NAME, file_name, Dict("response_stream" => stream))
            if AWS.DEFAULT_BACKEND[] isa AWS.HTTPBackend
                @test !isopen(stream)
            else
                # See: https://github.com/JuliaCloud/AWS.jl/issues/471
                @test_broken isopen(stream)
            end

            stream = Base.BufferStream()
            S3.get_object(
                BUCKET_NAME,
                file_name,
                Dict("response_stream" => stream, "return_stream" => true),
            )
            if AWS.DEFAULT_BACKEND[] isa AWS.HTTPBackend
                @test !isopen(stream)
            else
                @test isopen(stream)
            end
        finally
            S3.delete_object(BUCKET_NAME, file_name)
        end
    end

    @testset "issue 474" begin
        body = "foo\0bar"
        expected = Vector{UInt8}(body)
        file_name = "null.txt"

        try
            S3.put_object(BUCKET_NAME, file_name, Dict("body" => body))

            raw = S3.get_object(BUCKET_NAME, file_name, Dict("return_raw" => true))
            @test raw == expected

            stream = S3.get_object(BUCKET_NAME, file_name, Dict("return_stream" => true))
            if AWS.DEFAULT_BACKEND[] isa AWS.HTTPBackend
                @test stream isa Base.BufferStream
                @test !isopen(stream)

                if !isopen(stream)
                    @test read(stream) == expected
                end
            else
                @test stream isa IOBuffer
                @test isopen(stream)

                seekstart(stream)
                @test read(stream) == expected
            end
        finally
            S3.delete_object(BUCKET_NAME, file_name)
        end
    end

    # https://github.com/JuliaCloud/AWS.jl/issues/515
    @testset "issue 515" begin
        function _incomplete_patch(; data, num_attempts_to_fail=4)
            attempt_num = 0
            n = length(data)

            function _downloads_response(content_length)
                headers = ["content-length" => string(content_length)]
                return Downloads.Response("http", "", 200, "HTTP/1.1 200 OK", headers)
            end

            patch = if AWS.DEFAULT_BACKEND[] isa AWS.HTTPBackend
                @patch function HTTP.request(args...; response_stream, kwargs...)
                    attempt_num += 1
                    if attempt_num <= num_attempts_to_fail
                        write(response_stream, data[1:(n - 1)]) # an incomplete stream that shouldn't be retained
                        throw(HTTP.IOError(EOFError(), "msg"))
                    else
                        write(response_stream, data)
                        return HTTP.Response(200, "{\"Location\": \"us-east-1\"}")
                    end
                end
            elseif AWS.DEFAULT_BACKEND[] isa AWS.DownloadsBackend
                @patch function Downloads.request(args...; output, kwargs...)
                    attempt_num += 1
                    if attempt_num <= num_attempts_to_fail
                        write(output, data[1:(n - 1)])  # an incomplete stream that shouldn't be retained
                        message = "transfer closed with 1 bytes remaining to read"
                        e = Downloads.RequestError("", 18, message, _downloads_response(n))
                        throw(e)
                    else
                        write(output, data)
                        return _downloads_response(n)
                    end
                end
            end

            return patch
        end

        n = 100
        data = rand(UInt8, n)
        bucket = "www.invenia.ca"  # use public bucket as dummy
        key = "index.html"
        config = AWSConfig(; creds=nothing)

        @testset "Fail 2 attempts then succeed" begin
            apply(_incomplete_patch(; data=data, num_attempts_to_fail=2)) do
                retrieved = S3.get_object(bucket, key; aws_config=config)

                @test length(retrieved) == n
                @test retrieved == data
            end
        end

        @testset "Fail all 4 attempts then throw" begin
            err_t = if AWS.DEFAULT_BACKEND[] isa AWS.HTTPBackend
                HTTP.IOError
            else
                Downloads.RequestError
            end
            io = IOBuffer()

            apply(_incomplete_patch(; data=data, num_attempts_to_fail=4)) do
                params = Dict("response_stream" => io)
                @test_throws err_t S3.get_object(bucket, key, params; aws_config=config)

                seekstart(io)
                retrieved = read(io)
                @test length(retrieved) == n - 1
                @test retrieved == data[1:(n - 1)]
            end
        end
    end

finally
    S3.delete_bucket(BUCKET_NAME)
end


text:
const ptrdiff_t = Cptrdiff_t
const FILE = Libc.FILE
const SHRT_MAX = 32767
const USHRT_MAX = 65535
const CHAR_BIT = 8
const INT_MAX = typemax(Cint)
const INT_MIN = typemin(Cint)
const UINT_MAX = typemax(Cuint)
const UINT_MIN = typemin(Cuint)
const LONG_MAX = typemax(Clong)
const LONG_MIN = typemin(Clong)
const ULONG_MAX = typemax(Culong)


text:
import Acorn.EditorConfig

ed = Acorn.Editor()

Acorn.commandBind(ed, "h help")
@test EditorConfig.isKeyBound('h') == true
@test EditorConfig.getKeyBinding('h') == "help"

Acorn.commandBind(ed, "h ~")
@test EditorConfig.isKeyBound('h') == false

EditorConfig.rmKeyBinding('c')
Acorn.commandBind(ed, "c")
@test EditorConfig.isKeyBound('c') == false

Acorn.commandBind(ed, "a echo a b c d")
@test EditorConfig.isKeyBound('a') == true


text:
#
# accessors.jl --
#
# Accessor functions for OI-FITS objects.  Most of these are kept for backward
# compatibility but superseded by the `obj.key` syntax.
#
# This file must be included after OI-fITS formats definitions.
#
#------------------------------------------------------------------------------

# Automatically define getters from all fields of a data-block.
for extname in EXTNAMES
    let T = get_datablock_type(extname)
        for symb in Builder.FIELDS[extname]
            eval(Meta.parse("get_$symb(db::$T) = db.$symb"))
        end
    end
end

# Define getters which rely on indirections.
get_eff_wave(db::Union{OIVis,OIVis2,OIT3}) = db.instr.eff_wave
get_eff_band(db::Union{OIVis,OIVis2,OIT3}) = db.instr.eff_band

"""
    OIFITS.get_target_id(obj)

yields the "TARGET_ID" field of OI-FITS data-block object `obj` (an instance of
`OITarget`, `OIVis`, `OIVis2`, `OIT3`, `OIFlux` or `OIPolarization`).

!!! note
    This method is to be deprecated, use `obj.target_id` instead.

""" get_target_id

"""
    OIFITS.get_target(obj)

yields the `OITarget` data-block associated with OI-FITS data-block object
`obj`, if any, `nothing` otherwise.

!!! note
    This method is to be deprecated, use `obj.target` or `obj.target[arrname]`
    instead.

"""
get_target(obj::OIMaster{T}) where {T} =
    obj.target :: Union{Nothing,OITarget{T}}
#get_target(obj::OITarget) = obj
get_target(obj::OIDataBlock) = nothing
get_target(obj::Union{OIVis,OIVis2,OIT3,OIFlux,OIPolarization}) =
    get_target(obj.owner)

"""
    OIFITS.get_array(obj)

yields the `OIArray` data-block associated with OI-FITS data-block object
`obj`, if any, `nothing` otherwise.

    OIFITS.get_array(obj, arrname)

yields the `OIArray` data-block whose name is `arrname` in OI-FITS matser
object `obj` or `nothing` if not found.

!!! note
    This method is to be deprecated, use `obj.array` or `obj.array[arrname]`
    instead.

"""
get_array(obj::OIMaster{T}, arrname::AbstractString) where {T} =
    get(get_array(obj), fix_name(arrname),
        nothing) :: Union{Nothing,OIArray{T}}
get_array(obj::OIDataBlock) = nothing
get_array(obj::OIArray) = obj
get_array(obj::Union{OIMaster,OIVis,OIVis2,OIT3,OIFlux,OIPolarization}) =
    obj.array

"""
    OIFITS.get_instrument(obj)

yields the `OIWavelength` data-block associated with OI-FITS data-block object
`obj`, if any, `nothing` otherwise.

    OIFITS.get_instrument(obj, arrname)

yields the `OIWavelength` data-block whose name is `arrname` in OI-FITS matser
object `obj` or `nothing` if not found.

!!! note
    This method is to be deprecated, use `obj.instr` or `obj.instr[arrname]`
    instead.

"""
get_instrument(obj::OIMaster{T}, insname::AbstractString) where {T} =
    get(get_instrument(obj), fix_name(insname),
        nothing) :: Union{Nothing,OIWavelength{T}}
get_instrument(obj::OIDataBlock) = nothing
get_instrument(obj::OIWavelength) = obj
get_instrument(obj::Union{OIMaster,OIVis,OIVis2,OIT3,OIFlux,OIPolarization}) =
    obj.instr

"""
    OIFITS.get_targets(obj)

yields the array of names of the targets defined in OI-FITS master object `obj`.

"""
get_targets(obj::OIMaster) =
    (tgt = obj.target) === nothing ? String[] : tgt.target

"""
    OIFITS.get_arrays(obj)

yields the names of the interferometric arrays defined in OI-FITS master object
`obj`.

"""
get_arrays(obj::OIMaster) = collect(keys(obj.array))

"""
    OIFITS.get_instruments(obj)

yields the names of the instruments defined in OI-FITS master object `obj`.

"""
get_instruments(obj::OIMaster) = collect(keys(obj.instr))


text:
using Flux

struct Embedding{T}
    E::Flux.TrackedArray{T, 2}
end

Embedding(vocab::Integer, dims::Integer) = Embedding(Flux.param(Flux.gpu(randn(Float32, vocab, dims))))
SequenceTransform(x::AbstractArray{T, 1} where T, vsize::Integer) = permutedims(Float32.(Flux.onehotbatch(x, 1:vsize)), [2, 1]) |> Flux.gpu
(m::Embedding)(x::AbstractArray{T, 1} where T) = map(seq -> SequenceTransform(seq, size(m.E, 1)) * m.E, x)
Flux.@treelike Embedding


text:
function pretty_print(t::Trial)
    for i in 1:5
        if t.result.codes[i] == PRESENT_NOWHERE
            print('⬜')
        elseif t.result.codes[i] == PRESENT_ELSEWHERE
            print('🟨')
        elseif t.result.codes[i] == PRESENT_HERE
            print('🟩')
        else
            error("Invalid code encountered")
        end
    end
    println()
end

pretty_print(trials::Vector{Trial}) = foreach(pretty_print, trials)


text:
using DynamicalSystemsBase
using Test

println("\nTesting continuous system evolution...")

@testset "Lorenz63 oop" begin

    lo11 = Systems.lorenz() #with Jac
    ts1 = trajectory(lo11, 1.0, Δt = 0.01)
    @test size(ts1) == (101, 3)
    @test ts1[1, :] == SVector{3}(lo11.u0)
    @test !any(x -> abs(x) > 1e3, ts1[end])

    ts1 = trajectory(lo11, 1.0; Δt=0.1, diffeq = (abstol=1e-8, reltol=1e-8))
    @test size(ts1) == (11, 3)
    @test ts1[1, :] == SVector{3}(lo11.u0)

    data2 = trajectory(lo11, 100; save_idxs = 1:2)
    @test size(data2)[2] == 2
end

@testset "Lorenz96 iin" begin
    u = ones(5)
    lo11 = Systems.lorenz96(5, u)
    ts1 = trajectory(lo11, 2.0, Δt = 0.01)
    @test size(ts1) == (201, 5)
    @test ts1[1, :] == SVector{5}(u)

    ts1 = trajectory(lo11, 2.0; Δt=0.1, diffeq = (abstol=1e-8, reltol=1e-8))
    @test size(ts1) == (21, 5)
    @test ts1[1, :] == SVector{5}(u)
    @test !any(x -> abs(x) > 1e3, ts1[end])
end