text:
# BED Reader
# ==========

mutable struct Reader <: BioCore.IO.AbstractReader
    state::BioCore.Ragel.State
    index::Union{GenomicFeatures.Indexes.Tabix,Nothing}

    function Reader(input::BufferedStreams.BufferedInputStream, index=nothing)
        return new(BioCore.Ragel.State(file_machine.start_state, input), index)
    end
end

"""
    BED.Reader(input::IO; index=nothing)
    BED.Reader(input::AbstractString; index=:auto)

Create a data reader of the BED file format.

The first argument specifies the data source. When it is a filepath that ends
with *.bgz*, it is considered to be block compression file format (BGZF) and the
function will try to find a tabix index file (<filename>.tbi) and read it if
any. See <http://www.htslib.org/doc/tabix.html> for bgzip and tabix tools.

Arguments
---------
- `input`: data source
- `index`: path to a tabix file
"""
function Reader(input::IO; index=nothing)
    if isa(index, AbstractString)
        index = GenomicFeatures.Indexes.Tabix(index)
    end
    return Reader(BufferedStreams.BufferedInputStream(input), index)
end

function Reader(filepath::AbstractString; index=:auto)
    if isa(index, Symbol) && index != :auto
        throw(ArgumentError("invalid index argument: ':$(index)'"))
    end
    if endswith(filepath, ".bgz")
        input = BGZFStreams.BGZFStream(filepath)
        if index == :auto
            index = GenomicFeatures.Indexes.findtabix(filepath)
        end
    else
        input = open(filepath)
    end
    return Reader(input, index=index)
end

function Base.eltype(::Type{Reader})
    return Record
end

function BioCore.IO.stream(reader::Reader)
    return reader.state.stream
end

function GenomicFeatures.eachoverlap(reader::Reader, interval::GenomicFeatures.Interval)
    if reader.index === nothing
        throw(ArgumentError("index is null"))
    end
    return GenomicFeatures.Indexes.TabixOverlapIterator(reader, interval)
end

const record_machine, file_machine = (function ()
    cat = Automa.RegExp.cat
    rep = Automa.RegExp.rep
    opt = Automa.RegExp.opt

    record = let
        chrom = re"[ -~]*"
        chrom.actions[:enter] = [:mark]
        chrom.actions[:exit]  = [:record_chrom]

        chromstart = re"[0-9]+"
        chromstart.actions[:enter] = [:mark]
        chromstart.actions[:exit]  = [:record_chromstart]

        chromend = re"[0-9]+"
        chromend.actions[:enter] = [:mark]
        chromend.actions[:exit]  = [:record_chromend]

        name = re"[ -~]*"
        name.actions[:enter] = [:mark]
        name.actions[:exit]  = [:record_name]

        score = re"[0-9]+"
        score.actions[:enter] = [:mark]
        score.actions[:exit]  = [:record_score]

        strand = re"[+\-.?]"
        strand.actions[:enter] = [:record_strand]

        thickstart = re"[0-9]+"
        thickstart.actions[:enter] = [:mark]
        thickstart.actions[:exit]  = [:record_thickstart]

        thickend = re"[0-9]+"
        thickend.actions[:enter] = [:mark]
        thickend.actions[:exit]  = [:record_thickend]

        itemrgb = cat(re"[0-9]+", opt(cat(',', re"[0-9]+", ',', re"[0-9]+")))
        itemrgb.actions[:enter] = [:mark]
        itemrgb.actions[:exit]  = [:record_itemrgb]

        blockcount = re"[0-9]+"
        blockcount.actions[:enter] = [:mark]
        blockcount.actions[:exit]  = [:record_blockcount]

        # comma-separated values
        csv(x) = cat(rep(cat(x, ',')), opt(x))

        blocksizes = let
            blocksize = re"[0-9]+"
            blocksize.actions[:enter] = [:mark]
            blocksize.actions[:exit]  = [:record_blocksizes_blocksize]

            csv(blocksize)
        end
        blocksizes.actions[:exit] = [:record_blocksizes]

        blockstarts = let
            blockstart = re"[0-9]+"
            blockstart.actions[:enter] = [:mark]
            blockstart.actions[:exit]  = [:record_blockstarts_blockstart]

            csv(blockstart)
        end
        blockstarts.actions[:exit] = [:record_blockstarts]

        cat(
            chrom, '\t',
            chromstart, '\t',
            chromend,
            opt(cat('\t', name,
            opt(cat('\t', score,
            opt(cat('\t', strand,
            opt(cat('\t', thickstart,
            opt(cat('\t', thickend,
            opt(cat('\t', itemrgb,
            opt(cat('\t', blockcount,
            opt(cat('\t', blocksizes,
            opt(cat('\t', blockstarts)))))))))))))))))))
    end
    record.actions[:enter] = [:anchor]
    record.actions[:exit]  = [:record]

    newline = let
        lf = re"\n"
        lf.actions[:enter] = [:countline]

        cat(opt('\r'), lf)
    end

    file = rep(cat(record, newline))

    return map(Automa.compile, (record, file))
end)()

#=
write("bed.dot", Automa.machine2dot(file_machine))
run(`dot -Tsvg -o bed.svg bed.dot`)
=#

const record_actions = Dict(
    :record_chrom                  => :(record.chrom      = (mark:p-1) .- offset; record.ncols += 1),
    :record_chromstart             => :(record.chromstart = (mark:p-1) .- offset; record.ncols += 1),
    :record_chromend               => :(record.chromend   = (mark:p-1) .- offset; record.ncols += 1),
    :record_name                   => :(record.name       = (mark:p-1) .- offset; record.ncols += 1),
    :record_score                  => :(record.score      = (mark:p-1) .- offset; record.ncols += 1),
    :record_strand                 => :(record.strand     =       p     - offset; record.ncols += 1),
    :record_thickstart             => :(record.thickstart = (mark:p-1) .- offset; record.ncols += 1),
    :record_thickend               => :(record.thickend   = (mark:p-1) .- offset; record.ncols += 1),
    :record_itemrgb                => :(record.itemrgb    = (mark:p-1) .- offset; record.ncols += 1),
    :record_blockcount             => :(record.blockcount = (mark:p-1) .- offset; record.ncols += 1),
    :record_blocksizes_blocksize   => :(push!(record.blocksizes, (mark:p-1) .- offset)),
    :record_blocksizes             => :(record.ncols += 1),
    :record_blockstarts_blockstart => :(push!(record.blockstarts, (mark:p-1) .- offset)),
    :record_blockstarts            => :(record.ncols += 1),
    :record => :(record.filled = 1:p-1),
    :countline => :(),
    :mark => :(mark = p),
    :anchor => :())
eval(
    BioCore.ReaderHelper.generate_index_function(
        Record,
        record_machine,
        :(offset = mark = 0),
        record_actions))
eval(
    BioCore.ReaderHelper.generate_read_function(
        Reader,
        file_machine,
        :(offset = mark = 0),
        merge(record_actions, Dict(
            :record => quote
                BioCore.ReaderHelper.resize_and_copy!(record.data, data, BioCore.ReaderHelper.upanchor!(stream):p-1)
                record.filled = (offset+1:p-1) .- offset
                found_record = true
                @escape
            end,
            :countline => :(linenum += 1),
            :anchor => :(BioCore.ReaderHelper.anchor!(stream, p); offset = p - 1)))))


text:
"""
    omega_ratio(returns, target_return)

This function calculates the Omega ratio.

# Formula

    E[max(returns - target_return, 0)] / E[max(target_return - returns, 0)]

# Arguments
- `returns`:        Vector of asset returns.
- `target_return`:  Vector or scalar value of benchmark returns having same same frequency (e.g. daily) as the provided returns.
"""
function omega_ratio(returns, target_return)
    excess = returns .- target_return
    sum1 = sum(map(x -> max(0.0, x), excess))
    sum2 = -sum(map(x -> min(0.0, x), excess))
    sum1 / sum2
end

ext: jl


text:
include("../common/PWGrid_v02.jl")
include("../common/ortho_gram_schmidt.jl")
include("../common/wrappers_fft.jl")

include("EnergiesT.jl")
include("PotentialsT.jl")
include("gen_dr.jl")
include("init_pot_harm_3d.jl")
include("op_K.jl")
include("op_V_loc.jl")
include("op_H.jl")
#include("calc_grad.jl")
include("calc_grad_v2.jl")
include("calc_rho.jl")
include("Poisson_solve.jl")
include("LDA_VWN.jl")

function printMatrix( A::Array{ComplexF64,2} )
    Nrows = size(A)[1]
    Ncols = size(A)[2]
    for ir = 1:Nrows
        for ic = 1:Ncols
            @printf("(%6.3f,%6.3f) ", real(A[ir,ic]), imag(A[ir,ic]))
        end
        @printf("\n")
    end
end

function printMatrix( A::Array{Float64,2} )
    Nrows = size(A)[1]
    Ncols = size(A)[2]
    for ir = 1:Nrows
        for ic = 1:Ncols
            @printf("%6.3f ", A[ir,ic])
        end
        @printf("\n")
    end
end

function ∇E( pw::PWGrid, Potentials, Focc, W::Array{ComplexF64,2} )
    #
    Ngwx    = size(W)[1]
    Nstates = size(W)[2]
    Ω = pw.Ω
    Ns = pw.Ns
    #
    grad = zeros( ComplexF64, Ngwx, Nstates )
    #
    F = diagm(Focc)
    HW = op_H( pw, Potentials, W )
    #U = W' * W
    #U_sqrt = sqrtm( inv(U) )
    # BbbH
    #ℍ = U_sqrt * W' * HW * U_sqrt
    ℍ = W' * HW
    #printMatrix(U)
    #println("\nU_sqrt = "); printMatrix(U_sqrt)
    HFH = ℍ*F - F*ℍ
    println("\nℍ ="); printMatrix(ℍ)
    println("\nℍ*F - F*ℍ ="); printMatrix(HFH)
    # Calculation of Q
    #mu, V = eig(U)
    #denom = sqrt(mu)*ones(1,length(mu))
    #denom = denom + denom'
    denom = ones(Nstates,Nstates)*2.0
    #ℚ = V * ( ( V' * HFH * V ) ./ denom ) * V'
    ℚ = HFH ./ denom
    #println("mu = "); println(mu)
    #println("V = "); printMatrix(V)
    #println("denom = "); printMatrix(denom)
    println("ℚ = "); printMatrix(ℚ)
    grad = (HW - W * W'*HW)*F + W*ℚ
    return grad
end

function test_main( Ns )
    #
    const LatVecs = 16.0*diagm( ones(3) )
    #
    pw = PWGrid( Ns, LatVecs )
    #
    const Ω  = pw.Ω
    const r  = pw.r
    const G  = pw.gvec.G
    const G2 = pw.gvec.G2
    const Npoints = prod(Ns)
    const Ngwx = pw.gvecw.Ngwx
    #
    # Generate array of distances
    #
    center = sum(LatVecs,2)/2
    dr = gen_dr( r, center )
    #
    # Setup potential
    #
    V_ionic = init_pot_harm_3d( pw, dr )
    #
    const Nstates = 4
    #Focc = 2.0*ones(Nstates)
    Focc = [2.0, 1, 1, 1]

    srand(1234)
    psi = randn(Ngwx,Nstates) + im*randn(Ngwx,Nstates)
    psi = ortho_gram_schmidt(psi)

    Potentials = PotentialsT( V_ionic, zeros(Npoints), zeros(Npoints) )
    rho = calc_rho( pw, Focc, psi )
    Potentials.Hartree = real( G_to_R( Ns, Poisson_solve(pw, rho) ) )
    Potentials.XC = excVWN( rho ) + rho .* excpVWN( rho )

    grad2 = calc_grad( pw, Potentials, Focc, psi )

    grad1 = ∇E( pw, Potentials, Focc, psi )

    println("sum(grad1) = ", sum(grad1))
    println("sum(grad2) = ", sum(grad2))

end

test_main( [64, 64, 64] )


text:
using PyCall;
@pyimport smtplib;
fromaddr = "malcolm.sherrington@gmail.com";
toaddrs  = "malcolm@amisllp.com";

messy = """From: $fromaddr
To: $toaddrs
Subject: Test SMTP using PyCall

Testing - 1,2,3
""";

# Note that the blank line is necessary to distinguish
# the SMTP header from the message text.

username = fromaddr;
password = "ABCDEF7890";		# Not my real password

server = pycall(smtplib.SMTP,PyAny,"smtp.gmail.com:587");
server[:ehlo]();
server[:starttls]();
server[:login](username,password);
server[:sendmail](fromaddr,toaddrs,messy);
server[:quit]();

max_stars_repo_path: Module 3/Chapter09/py-email.jl

max_issues_repo_path: Module 3/Chapter09/py-email.jl

max_forks_repo_path: Module 3/Chapter09/py-email.jl


text:
using NiLang, NiLang.AD
using Test

@testset "gvar" begin
    g1 = GVar(0.0)
    @test (~GVar)(g1) === 0.0
    @assign (g1 |> grad) 0.5
    @test g1 === GVar(0.0, 0.5)
    @test_throws InvertibilityError (~GVar)(g1)
    @test !almost_same(GVar(0.0), GVar(0.0, 1.0))
    @test zero(GVar(3.0, 2.0)) == GVar(0.0)
    @test one(GVar(3.0, 2.0)) == GVar(1.0)
    @test iszero(GVar(0.0, 2.0))
    @test zero(GVar(2, AutoBcast([1, 0, 0]))) == GVar(0, AutoBcast([0, 0, 0]))
    @test GVar(true) == true
end

@testset "assign" begin
    arg = (1,2,GVar(3.0))
    @assign (arg |> tget(3)).g 4.0
    @test arg[3].g == 4.0
    gv = GVar(1.0, GVar(0.0))
    @test gv.g.g === 0.0
    @assign gv.g.g 7.0
    @test gv.g.g === 7.0
    gv = GVar(1.0, GVar(0.0))
    @assign gv |> grad |> grad 0.0
    @test gv.g.g === 0.0
    args = (GVar(0.0, 1.0),)
    @assign (args |> tget(1) |> grad) 0.0
    @test args[1].g == 0.0
    arr = [1.0]
    arr0 = arr
    @assign arr[] 0.0
    @test arr[] == 0.0
    @test arr === arr0
end

@testset "assign tuple" begin
    x = 0.3
    @instr for i=1:length(x) GVar(x) end
    @test x === GVar(0.3)
end

@testset "assign bcast func" begin
    # vector bcast
    x = [GVar(0.1, 0.1), GVar(0.2, 0.2)]
    res = [1.0, 2.0]
    @assign (x .|> value) res
    @test x == [GVar(1.0, 0.1), GVar(2.0, 0.2)]

    # tuple bcast
    x = (GVar(0.1, 0.1), GVar(0.2, 0.2))
    res = (1.0, 2.0)
    @assign (x .|> value) res
    @test x == (GVar(1.0, 0.1), GVar(2.0, 0.2))
end

@testset "GVar over general type" begin
    struct ABC{T1, T2}
       a::T1
       b::T1
       c::T2
    end
    x = ABC(1, 2, 3.0)
    @test GVar(x) == ABC(1, 2, GVar(3.0))
    @test (~GVar)(ABC(1, 2, GVar(3.0))) == x
end


text:
function get_trajopt_data_indirect(mechanism::Mechanism,env::Environment,Δt::Real,N::Int;
                                   relax_comp=false)
    vs = VariableSelector()

    for n = 1:N
        add_var!(vs, Symbol("q", n), num_positions(mechanism))
        add_var!(vs, Symbol("v", n), num_velocities(mechanism))
        if n < N
            add_var!(vs, Symbol("u", n), num_velocities(mechanism))
            add_var!(vs, Symbol("h", n), 1)
            for i = 1:length(env.contacts)
                if relax_comp
                    add_var!(vs, Symbol("slack", i, "_", n), 3)
                end
                β_dim = size(env.contacts[i].obstacle.basis,2)
                add_var!(vs, Symbol("β", i, "_", n), β_dim)
                add_var!(vs, Symbol("λ", i, "_", n), 1)
                add_var!(vs, Symbol("c_n", i, "_", n), 1)
            end
        end
    end

    cs = ConstraintSelector()

    for n = 1:N-1
        add_eq!(cs, Symbol("kin", n), num_positions(mechanism))
        add_eq!(cs, Symbol("dyn", n), num_velocities(mechanism))
        add_ineq!(cs, Symbol("h_pos", n), 1)
        for i = 1:length(env.contacts)
            β_dim = size(env.contacts[i].obstacle.basis,2)
            add_ineq!(cs, Symbol("β_pos", i, "_", n), β_dim)
            add_ineq!(cs, Symbol("λ_pos", i, "_", n), 1)
            add_ineq!(cs, Symbol("c_n_pos", i, "_", n), 1)
            add_ineq!(cs, Symbol("ϕ_pos", i, "_", n), 1)
            add_ineq!(cs, Symbol("fric_pos", i, "_", n), β_dim)
            add_ineq!(cs, Symbol("cone_pos", i, "_", n), 1)
            if relax_comp
                add_ineq!(cs, Symbol("slack_pos", i, "_", n), 3)
                add_ineq!(cs, Symbol("ϕ_c_n_comp", i, "_", n), 1)
                add_ineq!(cs, Symbol("fric_β_comp", i, "_", n), β_dim)
                add_ineq!(cs, Symbol("cone_λ_comp", i, "_", n), 1)
            else
                add_eq!(cs, Symbol("ϕ_c_n_comp", i, "_", n), 1)
                add_eq!(cs, Symbol("fric_β_comp", i, "_", n), β_dim)
                add_eq!(cs, Symbol("cone_λ_comp", i, "_", n), 1)
            end
        end
    end

    state_cache = [StateCache(mechanism) for n = 1:2]
    envj_cache = [EnvironmentJacobianCache(env) for n = 1:2]

    generate_solver_fn = :generate_solver_fn_trajopt_indirect
    extract_sol = :extract_sol_trajopt_indirect

    sim_data = SimData(mechanism,env,
                       state_cache,envj_cache,
                       Δt,vs,cs,generate_solver_fn,extract_sol,
                       [],[],[],[],[],[],[],[],[],N,[],[],[],[])

    sim_data
end

function extract_sol_trajopt_indirect(sim_data::SimData, xopt::AbstractArray{T}) where T
    N = sim_data.N
    vs = sim_data.vs
    relax_comp = haskey(vs.vars, :slack1_1)
    env = sim_data.env

    qtraj = Array{Array{Float64,1},1}(undef, 0)
    vtraj = Array{Array{Float64,1},1}(undef, 0)
    utraj = Array{Array{Float64,1},1}(undef, 0)
    htraj = Array{Float64,1}(undef, 0)
    contact_traj = Array{Array{Float64,1},1}(undef, 0)
    slack_traj = Array{Array{Float64,1},1}(undef, 0)
    x = MechanismState(sim_data.mechanism)
    for n = 1:N
        set_configuration!(x, vs(xopt, Symbol("q", n)))
        normalize_configuration!(x)
        q = configuration(x)
        push!(qtraj, q)
        push!(vtraj, vs(xopt, Symbol("v", n)))
        if n < N
            push!(utraj, vs(xopt, Symbol("u", n)))
            push!(htraj, vs(xopt, Symbol("h", n))[1])
            contact_sol = []
            for i = 1:length(env.contacts)
                if relax_comp
                    push!(slack_traj, vs(xopt, Symbol("slack", i, "_", n)))
                end
                contact_sol = vcat(contact_sol,
                                   vs(xopt, Symbol("c_n", i, "_", n)),
                                   vs(xopt, Symbol("β", i, "_", n)),
                                   vs(xopt, Symbol("λ", i, "_", n)))
            end
            push!(contact_traj, contact_sol)
        end
    end

    # some other useful vectors
    ttraj = vcat(0., cumsum(max.(0.,htraj))...)
    qv_mat = vcat(hcat(qtraj...), hcat(vtraj...))

    qtraj, vtraj, utraj, htraj, contact_traj, slack_traj, ttraj, qv_mat, xopt
end

function generate_solver_fn_trajopt_indirect(sim_data::SimData)
    N = sim_data.N
    vs = sim_data.vs
    cs = sim_data.cs

    relax_comp = haskey(vs.vars, :slack1_1)
    num_contacts = length(sim_data.env.contacts)
    num_vel = num_velocities(sim_data.mechanism)
    world = root_body(sim_data.mechanism)
    world_frame = default_frame(world)

    function eval_obj(x::AbstractArray{T}) where T
        f = 0.

        if relax_comp
            for n = 1:N-1
                for i = 1:num_contacts
                    slack = vs(x, Symbol("slack", i, "_", n))
                    f += sum(slack)
                end
            end
        end

        # extra user-defined objective
        for i = 1:length(sim_data.obj_fns)
            obj_name, obj_fn = sim_data.obj_fns[i]
            f += obj_fn(x)
        end

        f
    end

    function eval_cons(x::AbstractArray{T}) where T
        g = Vector{T}(undef, cs.num_eqs + cs.num_ineqs) # TODO preallocate

#         @threads for n = 1:(N-1)
        for n = 1:(N-1)
            q0 = vs(x, Symbol("q", n))
            v0 = vs(x, Symbol("v", n))
            u0 = vs(x, Symbol("u", n))
            h = vs(x, Symbol("h", n))
            qnext = vs(x, Symbol("q", n+1))
            vnext = vs(x, Symbol("v", n+1))

            x0 = sim_data.state_cache[1][T]
            xn = sim_data.state_cache[2][T]
            envj = sim_data.envj_cache[2][T]

            set_configuration!(x0, q0)
            set_velocity!(x0, v0)
            setdirty!(x0)
            set_configuration!(xn, qnext)
            set_velocity!(xn, vnext)
            setdirty!(xn)

            H = mass_matrix(x0)
            config_derivative = configuration_derivative(xn)
            dyn_bias = dynamics_bias(xn)

            contact_bias = zeros(T, num_vel)
            if (num_contacts > 0)
                contact_jacobian!(envj, xn)
                contact_τ_indirect!(contact_bias, sim_data, envj, x, n)
            end

            g[cs(Symbol("kin", n))] .= qnext .- q0 .- h .* config_derivative
            g[cs(Symbol("dyn", n))] .= H * (vnext - v0) .- h .* (u0 .- dyn_bias .- contact_bias)
            g[cs(Symbol("h_pos", n))] .= -h

            for i = 1:num_contacts
                cj = envj.contact_jacobians[i]
                c = cj.contact

                # TODO preallocate
                contact_v = cj.contact_rot' * point_jacobian(xn, path(sim_data.mechanism, world, c.body), transform(xn, c.point, world_frame)).J * vnext

                β = vs(x, Symbol("β", i, "_", n))
                λ = vs(x, Symbol("λ", i, "_", n))
                c_n = vs(x, Symbol("c_n", i, "_", n))

                Dtv = c.obstacle.basis' * contact_v

                g[cs(Symbol("β_pos", i, "_", n))] .= -β
                g[cs(Symbol("λ_pos", i, "_", n))] .= -λ
                g[cs(Symbol("c_n_pos", i, "_", n))] .= -c_n
                g[cs(Symbol("ϕ_pos", i, "_", n))] .= -envj.contact_jacobians[i].ϕ
                # λe + D'*v >= 0
                g[cs(Symbol("fric_pos", i, "_", n))] .= -(λ .+ Dtv)
                # μ*c_n - sum(β) >= 0
                g[cs(Symbol("cone_pos", i, "_", n))] .= -(c.obstacle.μ .* c_n .- sum(β))
                if !relax_comp
                    # dist * c_n = 0
                    g[cs(Symbol("ϕ_c_n_comp", i, "_", n))] .= envj.contact_jacobians[i].ϕ .* c_n
                    # (λe + Dtv)' * β = 0
                    g[cs(Symbol("fric_β_comp", i, "_", n))] .= (λ .+ Dtv) .* β
                    # (μ * c_n - sum(β)) * λ = 0
                    g[cs(Symbol("cone_λ_comp", i, "_", n))] .= (c.obstacle.μ .* c_n .- sum(β)) .* λ
                else
                    slack = vs(x, Symbol("slack", i, "_", n))
                    g[cs(Symbol("slack_pos", i, "_", n))] .= -slack
#                     g[cs(Symbol("ϕ_c_n_comp", i, "_", n))] .= envj.contact_jacobians[i].ϕ .* c_n .- slack[1]
                    g[cs(Symbol("ϕ_c_n_comp", i, "_", n))] .= envj.contact_jacobians[i].ϕ .* c_n
                    g[cs(Symbol("fric_β_comp", i, "_", n))] .= (λ .+ Dtv) .* β .- slack[2]
                    g[cs(Symbol("cone_λ_comp", i, "_", n))] .= (c.obstacle.μ .* c_n .- sum(β)) .* λ .- slack[3]
                end
            end
        end

        # extra user-defined constraints
        for i = 1:length(sim_data.con_fns)
            con_name, con_fn = sim_data.con_fns[i]
            g[cs(con_name)] .= con_fn(x)
        end

        g
    end

    generate_autodiff_solver_fn(eval_obj,eval_cons,cs.eqs,cs.ineqs,vs.num_vars)
end


text:
using JSON

const data = JSON.parse(readline())

# Remove ternary for part 1
JSON_value(d::Dict) = "red" ∉ values(d) ? sum(JSON_value.(values(d))) : 0

JSON_value(a::Vector) = sum(JSON_value.(a))

JSON_value(n::Number) = n

JSON_value(s::String) = 0

JSON_value(data) |> println

max_stars_repo_name: pLOPeGG/JuliaOfCode

max_issues_repo_name: pLOPeGG/JuliaOfCode

max_forks_repo_name: pLOPeGG/JuliaOfCode


text:
export QuellerGoodnight, Ritland, Lynch, LynchRitland, LynchLi, LiHorvitz, Moran, Blouin, Loiselle #, Wang

"""
    Blouin(ind1::GenoArray, ind2::GenoArray, locus_names::Vector{Symbol}; alleles::NamedTuple)
Allele sharing index described by Blouin (1996)

- Single Locus Equation: The number of alleles shared between individuals over ploidy.
    - If both allele positions are shared (e.g. AA x AA or AB x AB) then 1
    - If one allele position is shared (e.g. AB x AC) then 0.5
    - If neither allele position is shared (e.g. AB x CD) then 0
- How to combine multiple loci: Single locus estimates are simply averaged together
- Assumes no inbreeding

Blouin, M. S., Parsons, M., Lacaille, V., & Lotz, S. (1996). Use of microsatellite loci to classify individuals by relatedness. Molecular ecology, 5(3), 393-401.
"""
function Blouin(ind1::T, ind2::T, locus_names::Vector{Symbol}, alleles::U; kwargs...) where T <: GenoArray where U <: NamedTuple
    isempty(locus_names) && return missing

    Mxy = Vector{Float64}(undef, length(locus_names))
    loc_id = 0

    for (loc,gen1,gen2) in zip(locus_names, ind1, ind2)
        loc_id += 1
        i,j = gen1
        k,l = gen2

        Mxy[loc_id] = (((i ∈ gen2) & (k ∈ gen1)) + ((j ∈ gen2) & (l ∈ gen1))) / 2
    end
    return mean(Mxy)
end

"""
    LiHorvitz(ind1::GenoArray, ind2::GenoArray, locus_names::Vector{Symbol}; alleles::NamedTuple)
Allele sharing index described by Li and Horvitz (1953)

-Single Locus Equation: If all alleles are the same between individuals (eg. AA x AA) then 1.
    - If two alleles are shared between individuals (eg.  AA x AB or AB x AB) then 0.5.
    - If only one allele is shared between individuals (eg. AB x AC) then 0.25.
    - If no alleles are shared (eg. AB x CD) then 0.
- How to combine multiple loci: Single locus estimates are simply averaged together
- Assumes no inbreeding

Li, C. C., & Horvitz, D. G. (1953). Some methods of estimating the inbreeding coefficient. American journal of human genetics, 5(2), 107.
"""
function LiHorvitz(ind1::T, ind2::T, locus_names::Vector{Symbol}, alleles::U; kwargs...) where T <: GenoArray where U <: NamedTuple
    isempty(locus_names) && return missing

    Bxy = Vector{Float64}(undef, length(locus_names))

    loc_id = 0
    for (loc,gen1,gen2) in zip(locus_names, ind1, ind2)
        loc_id += 1
        i,j = gen1
        k,l = gen2

        Bxy[loc_id] = sum([i, j] .∈ [k,l]') / 4
    end
    return mean(Bxy)
end

"""
    Loiselle(ind1::GenoArray, ind2::GenoArray, locus_names::Vector{Symbol}; alleles::NamedTuple)
Calculates the moments based estimator of pairwise relatedness using the estimator propsed by
Loiselle et al (1995) and modified to individual dyads by Heuertz et al. (2003).

- Multiple Locus Equation:
- Assumes no inbreeding

See equations 22 in: Wang(2017) for variant of estimator used

Loiselle, B. A., Sork, V. L., Nason, J., & Graham, C. (1995). Spatial genetic structure of a tropical understory shrub, <i>Psychotria officinalis</i> (Rubiaceae). American journal of botany, 82(11), 1420-1425.
Heuertz, M., Vekemans, X., Hausman, J. F., Palada, M., & Hardy, O. J. (2003). Estimating seed vs. pollen dispersal from spatial genetic structure in the common ash. Molecular Ecology, 12(9), 2483-2495.
Wang, J. (2017). Estimating pairwise relatedness in a small sample of individuals. Heredity, 119(5), 302-313.
"""
function Loiselle(ind1::T, ind2::T, locus_names::Vector{Symbol}, alleles::U; kwargs...) where T <: GenoArray where U <: NamedTuple
    isempty(locus_names) && return missing
    d_kw = Dict(kwargs...)
    numerator1 = 0.0
    denominator1 = 0.0

    for (loc,gen1,gen2) in zip(locus_names, ind1, ind2)
        for allele in keys(alleles[loc])
            fq = alleles[loc][allele]
            numerator1 += ((sum(gen1 .== allele) / 2.0) - fq) * ((sum(gen2 .== allele) / 2.0) - fq)
            denominator1 += fq * (1.0 - fq)
        end
    end
    return numerator1 / denominator1 + 2.0 / (2 * d_kw[:n_samples] - 1)
end

"""
    Lynch(ind1::GenoArray, ind2::GenoArray, locus_names::Vector{Symbol}; alleles::NamedTuple)
Allele sharing index described by Lynch (1988)

- Single Locus Equation: If all alleles are the same between individuals (eg. AA x AA) then 1.
    - If both individuals are heterozygous with the same alleles or one is homozygous for the shared allele (eg. AB x AB or AA x AB) then 0.75.
    - If only one allele is shared between individuals (eg. AB x AC) then 0.5.
    - If no alleles are shared (eg. AB x CD) then 0.
- How to combine multiple loci: Single locus estimates are simply averaged together
- Assumes no inbreeding

Lynch, M. (1988). Estimation of relatedness by DNA fingerprinting. Molecular biology and evolution, 5(5), 584-599.
"""
function Lynch(ind1::T, ind2::T, locus_names::Vector{Symbol}, alleles::U; kwargs...) where T <: GenoArray where U <: NamedTuple
    isempty(locus_names) && return missing

    Sxy = Vector{Float64}(undef, length(locus_names))
    loc_id = 0

    for (loc,gen1,gen2) in zip(locus_names, ind1, ind2)
        loc_id += 1
        i,j = gen1
        k,l = gen2

        Sxy[loc_id] = ((i ∈ gen2) + (j ∈ gen2) + (k ∈ gen1) + (l ∈ gen1)) / 4
    end
    return mean(Sxy)
end

"""
    LynchLi(ind1::GenoArray, ind2::GenoArray, locus_names::Vector{Symbol}; alleles::NamedTuple)
Calculates the moments based estimator of pairwise relatedness by Lynch (1988) & improved by Li et al. (1993).

- Single Locus Equation:
- How to combine multiple loci: Sum the difference between observed and expected similarity across all loci and then divide by the sum of 1 - the expected similarity
- Assumes no inbreeding

See equations 13 - 16 in Wang (2017) for variant of estimator used

Li, C. C., Weeks, D. E., & Chakravarti, A. (1993). Similarity of DNA fingerprints due to chance and relatedness. Human heredity, 43(1), 45-52.
Wang, J. (2017). Estimating pairwise relatedness in a small sample of individuals. Heredity, 119(5), 302-313.
"""
function LynchLi(ind1::T, ind2::T, locus_names::Vector{Symbol}, alleles::U; kwargs...) where T <: GenoArray where U <: NamedTuple
    isempty(locus_names) && return missing

    numerator1 = 0.0
    denominator1 = 0.0

    for (loc,gen1,gen2) in zip(locus_names, ind1, ind2)
        a,b = gen1
        c,d = gen2

        Sxy = (0.5) * (((a == c) + (a == d) + (b == c) + (b == d)) / (2.0 * (1.0 + (a == b))) + ((a == c) + (a == d) + (b == c) + (b == d)) / (2.0 * (1.0 + (c == d))))
        #TODO Change to unbiased formulation (eq 25)
        S0 = 2.0 * sum(values(alleles[loc]) .^ 2) - sum(values(alleles[loc]) .^ 3)

        numerator1 += Sxy - S0
        denominator1 += 1.0 - S0
    end
    return numerator1 / denominator1
end

"""
    LynchRitland(ind1::GenoArray, ind2::GenoArray, locus_names::Vector{Symbol}; alleles::NamedTuple)
Calculates the moments based estimator of pairwise relatedness by Lynch and Ritland (1999).

- Single Locus Equation:
- How to combine multiple loci: Weighted average of each term seperately weighted by the sample variance (assuming zero relatedness) and subsequently divided by the average sampling variance
- Assumes no inbreeding

See equation 10 in Wang (2017) for variant of estimator used

Lynch, M., & Ritland, K. (1999). Estimation of pairwise relatedness with molecular markers. Genetics, 152(4), 1753-1766.
Wang, J. (2017). Estimating pairwise relatedness in a small sample of individuals. Heredity, 119(5), 302-313.
"""
function LynchRitland(ind1::T, ind2::T, locus_names::Vector{Symbol}, alleles::U; kwargs...) where T <: GenoArray where U <: NamedTuple
    isempty(locus_names) && return missing

    numerator1 = 0.0
    denominator1 = 0.0

    for (loc,gen1,gen2) in zip(locus_names, ind1, ind2)
        a,b = gen1
        c,d = gen2
        fq_a, fq_b, fq_c, fq_d = map(i -> alleles[loc][i], (a,b,c,d))

        n1 = fq_a * ((b == c) + (b == d)) + fq_b * ((a == c) + (a == d)) - 4.0 * fq_a * fq_b
        n2 = fq_c * ((d == a) + (d == b)) + fq_d * ((c == a) + (c == b)) - 4.0 * fq_c * fq_d

        d1 = 2.0 * (1.0 + (a == b)) * (fq_a + fq_b) - 8.0 * fq_a * fq_b
        d2 = 2.0 * (1.0 + (c == d)) * (fq_c + fq_d) - 8.0 * fq_c * fq_d

        WL1 = ((1 + (a == b)) * (fq_a + fq_b) - 4 * fq_a * fq_b) / (2 * fq_a * fq_b)
        WL2 = ((1 + (c == d)) * (fq_c + fq_d) - 4 * fq_c * fq_d) / (2 * fq_c * fq_d)

        numerator1 += ((n1 / d1) * WL1 + (n2 / d2) * WL2)
        denominator1 += (WL1 + WL2) / 2.0
    end
    return numerator1 / denominator1
end

"""
    Moran(ind1::GenoArray, ind2::GenoArray, locus_names::Vector{Symbol}; alleles::NamedTuple)
Reinterpretation of Moran's I (commonly used for spatial autocorrelation) to estimate genetic relatedness
by Hardy and Vekemans (1999)

- Multiple Locus Equation:
- Assumes no inbreeding

Hardy, O. J., & Vekemans, X. (1999). Isolation by distance in a continuous population: reconciliation between spatial autocorrelation analysis and population genetics models. Heredity, 83(2), 145-154.
"""
function Moran(ind1::T, ind2::T, locus_names::Vector{Symbol}, alleles::U; kwargs...) where T <: GenoArray where U <: NamedTuple
    #TODO NEED TO CHECK TO CONFIRM EQUATIONS
    isempty(locus_names) && return missing

    numerator1 = 0.0
    denominator1 = 0.0

    for (loc,gen1,gen2) in zip(locus_names, ind1, ind2)
        for allele in keys(alleles[loc])
            fq = alleles[loc][allele]
            numerator1 += ((sum(gen1 .== allele) / 2.0) - fq) * ((sum(gen2 .== allele) / 2.0) - fq)
            #denominator1 += ((sum(gen1 .== allele) / 2.0) - fq)^2
            denominator1 += (((sum(gen1 .== allele) / 2.0) - fq)^2 + ((sum(gen2 .== allele) / 2.0) - fq)^2) / 2.0
        end
        #denominator1 += (1 / (length(alleles[loc]) - 1))
    end
    return (numerator1 / denominator1)
end

function Moran_experimental(ind1::T, ind2::T, locus_names::Vector{Symbol}, alleles::U; kwargs...) where T <: GenoArray where U <: NamedTuple
    #TODO NEED TO CHECK TO CONFIRM EQUATIONS
    isempty(locus_names) && return missing

    numerator1 = Vector{Float64}(undef, length(locus_names))
    denominator1 = similar(numerator1)

    numerator1 = numerator1 .* 0.0
    denominator1 = denominator1 .* 0.0

    idx = 0
    for (loc,gen1,gen2) in zip(locus_names, ind1, ind2)
        idx += 1
        for allele in keys(alleles[loc])
            fq = alleles[loc][allele]
            numerator1[idx] += ((sum(gen1 .== allele) / 2.0) - fq) * ((sum(gen2 .== allele) / 2.0) - fq)
            denominator1[idx] += (((sum(gen1 .== allele) / 2.0) - fq)^2) #+ ((sum(gen2 .== allele) / 2.0) - fq)^2)  / 2
        end
        denominator1[idx] += (1 / (length(alleles[loc]) - 1))
    end

    return mean(numerator1 ./ denominator1)
end

"""
    QuellerGoodnight(ind1::GenoArray, ind2::GenoArray, locus_names::Vector{Symbol}; alleles::NamedTuple)
Calculates the moments based estimator of pairwise relatedness developed by Queller & Goodnight (1989).

- Single Locus Equation:
- How to combine multiple loci:
    - Multiple loci are combined by independently summing the two numerator and two denominator terms before performing the final division and averaging the two components.
- Assumes no inbreeding
See equation 3 in Wang(2017) for variant of estimator used.

Queller, D. C., & Goodnight, K. F. (1989). Estimating relatedness using genetic markers. Evolution, 43(2), 258-275.
Wang, J. (2017). Estimating pairwise relatedness in a small sample of individuals. Heredity, 119(5), 302-313.
"""
function QuellerGoodnight(ind1::T, ind2::T, locus_names::Vector{Symbol}, alleles::U; kwargs...) where T <: GenoArray where U <: NamedTuple
    isempty(locus_names) && return missing

    numerator1 = 0.0
    numerator2 = 0.0
    denominator1 = 0.0
    denominator2 = 0.0

    for (loc,gen1,gen2) in zip(locus_names, ind1, ind2)
        a,b = gen1
        c,d = gen2
        ident = ((a == c) + (a == d) + (b == c) + (b == d))
        fq_a, fq_b, fq_c, fq_d = map(i -> alleles[loc][i], (a,b,c,d))

        numerator1 += ident - 2.0 * (fq_a + fq_b)
        numerator2 += ident - 2.0 * (fq_c + fq_d)

        denominator1 += (2.0 * (1.0 + (a==b) - fq_a - fq_b))
        denominator2 += (2.0 * (1.0 + (c==d) - fq_c - fq_d))
    end
    return (numerator1/denominator1 + numerator2/denominator2)/2.0
end

"""
    Ritland(ind1::GenoArray, ind2::GenoArray, locus_names::Vector{Symbol}; alleles::NamedTuple)
Calculates the moments based estimator of pairwise relatedness proposed by Li and Horvitz (1953) and implemented/made popular by Ritland (1996).

- Single Locus Equation:
- How to combine multiple loci: A weighted average of individual locus specific estimates weighted by sampling variance
- Assumes no inbreeding

See equation 7 in: Wang (2017) for variant of estimator used

Ritland, K. (1996). Estimators for pairwise relatedness and individual inbreeding coefficients. Genetics Research, 67(2), 175-185.
Wang, J. (2017). Estimating pairwise relatedness in a small sample of individuals. Heredity, 119(5), 302-313.
"""
function Ritland(ind1::T, ind2::T, locus_names::Vector{Symbol}, alleles::U; kwargs...) where T <: GenoArray where U <: NamedTuple
    isempty(locus_names) && return missing

    numerator1 = 0.0
    denominator1 = 0.0

    for (loc,gen1,gen2) in zip(locus_names, ind1, ind2)
        a,b = gen1
        c,d = gen2

        A = ((alleles[loc] |> length) - 1)

        R = 0.0
        for i in unique((a,b,c,d))
            # Individual locus relatedness value (eq 7 in paper)
            R += ((((a == i) + (b == i)) * ((c == i) + (d == i))) / (4.0 * alleles[loc][i]))
        end
        R = (2.0 / A) * (R - 1.0)
        # numerator for weighted combination of loci
        numerator1 += (R * A)
        # denominator for weighted combination of loci
        denominator1 += A
    end
    return numerator1 / denominator1
end

### Wang 2002 helper functions ###
function a_wang_base(m::Int, alleles::Dict)
    sum(values(alleles) .^ m)
end

function a_wang(N::Int, alleles::Dict)
    #unbiased estimator
    a = 0.0

    b = (N * a_wang_base(2, alleles) - 1) / (N - 1)

    c = (N^2 * a_wang_base(3, alleles) - 3 * (N - 1) * b - 1) / ((N - 1) * (N - 2))

    d = (N^3 * a_wang_base(4, alleles) - 6 * (N - 1) * (N - 2) * c - 7 * (N - 1) * b - 1) / (N^3 - 6 * N^2 + 11 * N - 6)

    return [a, b, c, d]
end

"""
Wang(ind1::GenoArray, ind2::GenoArray, locus_names::Vector{Symbol}; alleles::NamedTuple)
Calculates the moments based estimator of pairwise relatedness by Wang (2002).

-Single Locus Equation:
-How to combine multiple loci: Each individual locus subcomponent (b-g) and each genotypic state (P1-P3) is averaged weighted by the average similarity of unrelated dyads at each locus. Then the values of V, Φ, Δ, and r are calculated

-Assumes no inbreeding
-Corrected for sampling bias in allele frequencies to get an unbiased estimator

Wang, J. (2002). An estimator for pairwise relatedness using molecular markers. Genetics, 160(3), 1203-1215.
"""
function Wang(ind1::T, ind2::T, locus_names::Vector{Symbol}, alleles::U; kwargs...) where T <: GenoArray where U <: NamedTuple
    #TODO NEED TO CHECK TO CONFIRM EQUATIONS
    isempty(locus_names) && return missing
    kw_dict = Dict(kwargs...)
    P1 = Vector{Float64}(undef, length(locus_names))
    P2, P3, P4, u, b, c, d, e, f, g = map(i -> similar(P1), 1:10)
    loc_id = 0

    for (loc,gen1,gen2, N) in zip(locus_names, ind1, ind2, kw_dict[:loc_n])
        loc_id += 1
        i,j = gen1
        k,l = gen2

        #N = nonmissing(data.loci[data.loci.locus .== string(loc), :genotype])

        a = a_wang(2 * N, alleles[loc])
        a2_sq = a[2] ^ 2

        u[loc_id] = 2 * a[2] - a[3]

        # Which category of dyad
        Sxy = ((i ∈ gen2) + (j ∈ gen2) + (k ∈ gen1) + (l ∈ gen1)) / 4

        # Both alleles shared between individuals either the same or different
        P1[loc_id] = 1.0 * (Sxy == 1)
        # One allele shared between individuals and one is homozygous for that allele
        P2[loc_id] = 1.0 * (Sxy == (3/4))
        # One allele shared with the other two being unique
        P3[loc_id] = 1.0 * (Sxy == (1/2))
        P4[loc_id] = 1.0 * ((P1 + P2 + P3) == 0)

        b[loc_id] = (2.0 * a2_sq - a[4])
        c[loc_id] = (a[2] - 2.0 * a2_sq + a[4])
        d[loc_id] = (4.0 * (a[3] - a[4]))
        e[loc_id] = (2.0 * (a[2] - 3.0 * a[3] + 2.0 * a[4]))
        f[loc_id] = (4.0 * (a[2] - a2_sq - 2.0 * a[3] + 2.0 * a[4]))
        g[loc_id] = (1.0 - 7.0 * a[2] + 4.0 * a2_sq + 10.0 * a[3] - 8.0 * a[4])

    end
    #return (1 / (sum(1/u) * u)) * r
    w = (1 / (sum(1/u) * u))

    P1 = w * P1
    P2 = w * P2
    P3 = w * P3

    b = w * b
    c = w * c
    d = w * d
    e = w * e
    f = w * f
    g = w * g

    #Eq 11
    V = (1.0 - b)^2 * (e^2 * f + d * g^2) -
        (1.0 - b) * (e * f - d * g)^2 +
        2.0 * c * d * f * (1.0 - b) * (g + e) +
        c^2 * d * f * (d + f)

    #Eq 9
    Φ = (d * f * ((e + g) * (1.0 - b) + c * (d + f)) * (P1 - 1.0) +
        d * (1.0 - b) * (g * (1.0 - b - d) + f * (c + e)) * P3 +
        f * (1.0 - b) * (e * (1.0 - b - f) + d * (c + g)) * P2) / V

    #Eq 10
    Δ = (c * d * f * (e + g) * (P1 + 1.0 - 2 * b) +
        ((1.0 - b) * (f * e^2 + d * g^2) - (e * f - d * g)^2) * (P1 - b) +
        c * (d * g - e * f) * (d * P3 - f * P2) - c^2 * d * f * (P3 + P2 - d - f) -
        c * (1.0 - b) * (d * g * P3 + e * f * P2)) / V

    r = (Φ/2.0 + Δ)
    return r
end


text:
import StatsPlots.scatter
import StatsPlots.mm

using LocalizationMicroscopy

function showpairs(molecules1, molecules2, pairedmolecules1, pairedmolecules2, distances)
    if length(pairedmolecules1) != length(pairedmolecules2) != length(distances)
        throw(ErrorException("Lengths of paired inputs must be identical."))
    end

    p = plot()
    pairedcoordinates1 = extractcoordinates(pairedmolecules1)
    pairedcoordinates2 = extractcoordinates(pairedmolecules2)
    coordinates1 = extractcoordinates(molecules1)
    coordinates2 = extractcoordinates(molecules2)
    scatter!(coordinates1[1,:], coordinates1[2,:], markersize = 2, markercolor = :white)
    scatter!(coordinates2[1,:], coordinates2[2,:], markersize = 2, markercolor = :white)
    scatter!(pairedcoordinates1[1,:], pairedcoordinates1[2,:], markersize = 1, markerstrokecolor=nothing, markercolor = :magenta)
    scatter!(pairedcoordinates2[1,:], pairedcoordinates2[2,:], markersize = 1, markerstrokecolor=nothing, markercolor = :green)
end

using Printf

distanceprobabilityplot(result::Result) = distanceprobabilityplot(result.distances, result.percentileranks)
function distanceprobabilityplot(distances::Vector{AbstractFloat}, percentileranks::Vector{AbstractFloat})
    lessthan10 = count((percentileranks .< 0.1) .& (distances .< 200)) / length(percentileranks)
    p1 = scatter(distances, percentileranks,
                 xaxis=("distance (nm)", (0,1000), 45),
                 yaxis=(:hide),
                 marker=(2, stroke(0)),
                 legend=:none,
                 left_margin=-30mm,
                 top_margin=-30mm,
                 annotations=(190,0,text(Printf.@sprintf("%.3f", lessthan10), 8, :right, :bottom)),
                 ywiden=false,
                 tick_direction=:out)
    line = plot!(p1, [0; 1000], [0.1, 0.1], line=(:gray))
    line = plot!(p1, [200; 200], [0, 1.0], line=(:gray))
    p2 = scatter(distances, percentileranks,
                 xaxis=((1000,15000), (1000:7000:15000), 45),
                 xformatter=:plain,
                 yaxis=(:hide),
                 marker=(2, stroke(0)),
                 legend=:none,
                 top_margin=-30mm,
                 left_margin=-20mm,
                 ywiden=false,
                 tick_direction=:out)

    p3 = histogram(distances,
                   bins=0:25:1000,
                   left_margin=-10mm,
                   xaxis=((0,1000), :hide),
                   linealpha=0,
                   legend=:none,
                   tick_direction=:out)
    p4 = histogram(distances,
                   bins=1000:25:15000,
                   xaxis=((1000,15000), (1000:7000:15000), :hide),
                   legend=:none,
                   linealpha=0,
                   left_margin=-30mm,
                   yaxis=(:hide))
    p5 = histogram(percentileranks,
                   bins=0:0.05:1.01,
                   orientation = :horizontal,
                   linealpha=0,
                   xaxis=(:flip, 45),
                   yaxis=("p(chance association)"),
                   legend=:none,
                   top_margin=-30mm,
                   ywiden=false,
                   tick_direction=:out)
    plot(Plots.Plot(), p3, p4, p5, p1, p2, layout = grid(2, 3, widths=[0.15,0.65,0.15], heights=[0.1,0.9]), link=:y)
end

moleculesplot(result::Result) = moleculesplot(result.channels[1].molecules, result.channels[2].molecules)
function moleculesplot(molecules1::Vector{Molecule}, molecules2::Vector{Molecule})
    mol1coords = molecules1 |> extractcoordinates
    mol2coords = molecules2 |> extractcoordinates
    scatter(mol1coords[1,:], mol1coords[2,:], marker=(2, stroke(0), :magenta))
    plot!([1000; 6000], [1000; 1000], line=(3, :black), annotations=(1000,1250,text("5 \\mum", 10, :left, :top)))

    scatter!(mol2coords[1,:], mol2coords[2,:], marker=(2, stroke(0), :green))
    plot!(aspect_ratio=:equal, xlims=(0, 40960), yaxis=((0, 40960), :flip), legend=:none, grid=:hide, ticks=(0),
          framestyle=:box) # check limit
end

moleculesplot_sim(result::Result) = moleculesplot_sim(result.channels[1].molecules, result.channels[2].molecules)
function moleculesplot_sim(molecules1::Vector{Molecule}, molecules2::Vector{Molecule})
    mol1coords = molecules1 |> extractcoordinates
    mol2coords = molecules2 |> extractcoordinates
    scatter(mol1coords[1,:], mol1coords[2,:], marker=(4, stroke(0), :magenta))
    scatter!(mol2coords[1,:], mol2coords[2,:], marker=(4, stroke(0), :green))
    plot!(aspect_ratio=:equal, xlims=(-2821, 2821), yaxis=((-2821, 2821), :flip), legend=:none, grid=:hide, ticks=(0),
          framestyle=:box) # check limit
end

neighborsplot(result::Result) = neighborsplot(result.channels[1].molecules, result.channels[2].molecules, result.channels[1].neighbormolecules, result.channels[2].neighbormolecules, result.percentileranks, result.distances)
function neighborsplot(molecules1::Vector{Molecule}, molecules2::Vector{Molecule}, neighbormolecules1::Vector{Molecule}, neighbormolecules2::Vector{Molecule}, percentileranks::Vector{AbstractFloat}, distances::Vector{AbstractFloat})
    if length(neighbormolecules1) != length(neighbormolecules2) != length(percentileranks) != length(distances)
        throw(ErrorException("Lengths of paired inputs must be identical."))
    end
    mol1coords = molecules1 |> extractcoordinates
    mol2coords = molecules2 |> extractcoordinates
    scatter(mol1coords[1,:], mol1coords[2,:], marker=(2, stroke(0), :magenta, 0.25))
    scatter!(mol2coords[1,:], mol2coords[2,:], marker=(2, stroke(0), :green, 0.25))

    mol1coords = neighbormolecules1 |> extractcoordinates
    mol2coords = neighbormolecules2 |> extractcoordinates
    scatter!(mol1coords[1,:], mol1coords[2,:], marker=(2, stroke(0), :magenta, 0.25))
    plot!([1000; 6000], [1000; 1000], line=(3, :black), annotations=(1000,1250,text("5 \\mum", 10, :left, :top)))
    scatter!(mol2coords[1,:], mol2coords[2,:], marker=(2, stroke(0), :green, 0.25))

    isbound = (percentileranks .< 0.1) .& (distances .< 200)
    boundmolecules1 = neighbormolecules1[isbound]
    boundmolecules2 = neighbormolecules2[isbound]
    mol1coords = boundmolecules1 |> extractcoordinates
    mol2coords = boundmolecules2 |> extractcoordinates
    scatter!(mol1coords[1,:], mol1coords[2,:], marker=(2, stroke(0), :black, 0.5))
    scatter!(mol2coords[1,:], mol2coords[2,:], marker=(2, stroke(0), :black, 0.5))

    plot!(aspect_ratio=:equal, xlims=(0, 40960), yaxis=((0, 40960), :flip), legend=:none, grid=:hide, ticks=(0),
          framestyle=:box) # check limit
end

neighborsplot_forprint(result::Result, include_scalebar=true) = neighborsplot_forprint(result.channels[1].molecules, result.channels[2].molecules, result.channels[1].neighbormolecules, result.channels[2].neighbormolecules, result.percentileranks, result.distances, true)
function neighborsplot_forprint(molecules1::Vector{Molecule}, molecules2::Vector{Molecule}, neighbormolecules1::Vector{Molecule}, neighbormolecules2::Vector{Molecule}, percentileranks::Vector{AbstractFloat}, distances::Vector{AbstractFloat}, include_scalebar=true)
    if length(neighbormolecules1) != length(neighbormolecules2) != length(percentileranks) != length(distances)
        throw(ErrorException("Lengths of paired inputs must be identical."))
    end
    mol1coords = molecules1 |> extractcoordinates
    mol2coords = molecules2 |> extractcoordinates
    scatter(mol1coords[1,:], mol1coords[2,:], marker=(8, stroke(0), :magenta, 0.25))
    scatter!(mol2coords[1,:], mol2coords[2,:], marker=(8, stroke(0), :green, 0.25))
    mol1coords = neighbormolecules1 |> extractcoordinates
    mol2coords = neighbormolecules2 |> extractcoordinates

    scatter!(mol1coords[1,:], mol1coords[2,:], marker=(8, stroke(0), :magenta, 0.25))
    plot!([1000; 6000], [1000; 1000], line=(3, :black), annotations=(1000,1250,text("5 \\mum", 10, :left, :top)))
    scatter!(mol2coords[1,:], mol2coords[2,:], marker=(8, stroke(0), :green, 0.25))

    isbound = (percentileranks .< 0.1) .& (distances .< 200)
    boundmolecules1 = neighbormolecules1[isbound]
    boundmolecules2 = neighbormolecules2[isbound]
    mol1coords = boundmolecules1 |> extractcoordinates
    mol2coords = boundmolecules2 |> extractcoordinates
    scatter!(mol1coords[1,:], mol1coords[2,:], marker=(8, stroke(0), :black, 0.5))
    scatter!(mol2coords[1,:], mol2coords[2,:], marker=(8, stroke(0), :black, 0.5))

    plot!(aspect_ratio=:equal, xlims=(0, 40960), yaxis=((0, 40960), :flip), legend=:none, grid=:hide, ticks=(0),
          framestyle=:box, size=(2048,2048)) # check limit
end

function localizationsplot(result::Result; color1 = :magenta, color2 = :green, insetbox = [[0,0], [0,0]], include_scalebar = false, forprint = false)
    if forprint
        locmarkersize = 8
        boxlinewidth = 12
        scalewidth = 12
        textsize = 48
        size = (2048, 2048)
    else
        locmarkersize = 2
        boxlinewidth = 3
        scalewidth = 6
        textsize = 20
        size = (512, 512)
    end

    loc1coords = mapreduce(x -> x.group.localizations, vcat, result.channels[1].molecules) |> extractcoordinates
    loc2coords = mapreduce(x -> x.group.localizations, vcat, result.channels[2].molecules) |> extractcoordinates

    scatter(loc1coords[1,:], loc1coords[2,:], marker=(locmarkersize, stroke(0), 0.75, color1))
    scatter!(loc2coords[1,:], loc2coords[2,:], marker=(locmarkersize, stroke(0), 0.75, color2))

    if insetbox != [[0,0], [0,0]]
        xlims, ylims = first(insetbox), last(insetbox)
        plot!([xlims[1], xlims[2], xlims[2], xlims[1], xlims[1]], [ylims[1], ylims[1], ylims[2], ylims[2], ylims[1]], line = (boxlinewidth, :black))
    end

    include_scalebar && plot!([1000; 6000], [1800; 1800], line=(scalewidth, :black),
          annotations=(1000, 2000, text("5 \\mum", textsize, :left, :top)))

    plot!(aspect_ratio=:equal, xlims=(0, 40960), yaxis=((0, 40960), :flip), legend=:none, grid=:hide, ticks=(0), framestyle=:none, size=size)
end

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

function insetplot(result::Result, xlims, ylims; color1 = :magenta, color2 = :green, include_scalebar = false, forprint = false)
    if forprint
        locmarkersize = 16
        molmarkersize = 64
		molstrokesize = 2
		scalewidth = 12
		textsize = 48
        size = (2048, 2048)
    else
        locmarkersize = 4
        molmarkersize = 16
		molstrokesize = 4
		scalewidth = 6
		textsize = 20
        size = (512, 512)
    end

    mol1coords = result.channels[1].molecules |> extractcoordinates
    mol2coords = result.channels[2].molecules |> extractcoordinates

    scatter(mol1coords[1,:], mol1coords[2,:], marker=(molmarkersize, stroke(molstrokesize, color1), 0.75, :white))
    scatter!(mol2coords[1,:], mol2coords[2,:], marker=(molmarkersize, stroke(molstrokesize, color2), 0.75, :white))
	scatter!(mol1coords[1,:], mol1coords[2,:], marker=(molmarkersize, stroke(0), 0.25, color1))
    scatter!(mol2coords[1,:], mol2coords[2,:], marker=(molmarkersize, stroke(0), 0.25, color2))

    loc1coords = mapreduce(x -> x.group.localizations, vcat, result.channels[1].molecules) |> extractcoordinates
    loc2coords = mapreduce(x -> x.group.localizations, vcat, result.channels[2].molecules) |> extractcoordinates

    scatter!(loc1coords[1,:], loc1coords[2,:], marker=(locmarkersize, stroke(0), 0.75, color1))
    scatter!(loc2coords[1,:], loc2coords[2,:], marker=(locmarkersize, stroke(0), 0.75, color2))

    include_scalebar && plot!(first(xlims) .+ [10; 210], first(ylims) .+ [10; 10], line=(scalewidth, :black),
          annotations=(first(xlims) + 10, first(ylims) + 20, text("200 nm", textsize, :left, :top)))

    plot!(aspect_ratio=:equal, xlims=xlims, yaxis=(ylims, :flip), legend=:none, grid=:hide, ticks=(0), framestyle=:none, size=size)
end

function yzframelocalizationplot(result::Result)
    loc1coords = mapreduce(x -> x.group.localizations, vcat, result.channels[1].molecules) |> extractcoordinates
    loc1frames = [y.frame for y ∈ mapreduce(x -> x.group.localizations, vcat, result.channels[1].molecules)]
    loc2coords = mapreduce(x -> x.group.localizations, vcat, result.channels[2].molecules) |> extractcoordinates
    loc2frames = [y.frame for y ∈ mapreduce(x -> x.group.localizations, vcat, result.channels[2].molecules)]

    scatter(loc1frames, loc1coords[2,:], marker=(2, :magenta, stroke(0)), xaxis=("frames", (0,22000)),
            yaxis=("y", (0,40960), :flip, 90), legend=:none)
    scatter!(loc2frames, loc2coords[2,:], marker=(2, :green, stroke(0)))
    plot!(tick_direction=:out, size=(1024,512))
end


text:
module Quadrant_Data
using MAT
using Plots
using DataFrames
using StatsBase

export calculatepvalue
export openmatdata
export midpoint
export totaltime
export testdata
export createquadcarlo

include("utils.jl")
# Write your package code here.

### Currently Needs F and G to be defined ###
function midpoint(pos1, pos2)
    x = (pos1[1] + pos2[1])/2
    y = (pos1[2]+pos2[2])/2
    return [x, y]
end

 function slope(pos1, pos2)
   m = (pos1[2]-pos2[2])/(pos1[1]-pos1[1])
   return m
end

function openmatdata(file_name, center)
    file = matread(file_name)
    file = file["trajectories"]
    quadrantdf = DataFrame(x = file[:,1,1], y = file[:,1,2])
    filter!(isfinite ∘ sum, quadrantdf)
    array = quadrantarray(quadrantdf, center)
    quadrantdf[!, :quadrant] = array
    return quadrantdf
end

function quadrantarray(dataframe, center)
    array = []
    for z in 1:length(dataframe[!,1])
        if(dataframe[z,1] > center[1] && dataframe[z,2] > center[2])
            push!(array, "quad 1")
        elseif(dataframe[z,1] <= center[1] && dataframe[z,2] > center[2])
            push!(array, "quad 2")
        elseif(dataframe[z,1] <= center[1] && dataframe[z,2] <= center[2])
            push!(array, "quad 3")
        elseif(dataframe[z,1] > center[1] && dataframe[z,2] <= center[2])
            push!(array, "quad 4")
        else
            throw(DomainError("Not in a Quadrant"))
        end
    end
    return array
end

function createquadcarlo(dataframe) #creates a data frame with the time spent in each quadrant before moving to the next quadrant
    quadcarlo = DataFrame(quadrant = String[], time = Float64[])
    current = 0
    previous = dataframe[1,3]
    time = 0
    for z in 1:length(dataframe[!,1])
        current = dataframe[z,3]
        if current == previous
            time = time + 1
        else
            push!(quadcarlo, (previous, time))
            time = 1
        end
        previous = dataframe[z,3]
    end
    push!(quadcarlo, (previous, time))
    return quadcarlo
end

function randomizequadcarlo(dataframe)  #randomizes the amount of time spent in each quadrant
    randorder = sample(1:nrow(dataframe), nrow(dataframe), replace = false)
    temp = DataFrame(time = dataframe[!,2], randorder = randorder)
    sort!(temp, 2)
    randquadcarlo = DataFrame(quadrant = dataframe[!,1], time = temp[!,1])
    return randquadcarlo
end

function testdata(quadcarlo, quad) #looks at a dataframe and returns a test statistic based on difference between the overall average and the chosen quadrant's average
    quad1, quad2, quad3, quad4 = groupby(quadcarlo, [:quadrant], sort = true)
    quad1sum, quad2sum, quad3sum, quad4sum = sum(quad1[!,2]), sum(quad2[!,2]), sum(quad3[!,2]), sum(quad4[!,2])
    if(quad=="quad 1")
        teststat = quad1sum - (quad1sum + quad2sum + quad3sum + quad4sum)/4
    elseif(quad=="quad 2")
        teststat = quad2sum - (quad1sum + quad2sum + quad3sum + quad4sum)/4
    elseif(quad=="quad 3")
        teststat = quad3sum - (quad1sum + quad2sum + quad3sum + quad4sum)/4
    elseif(quad=="quad 4")
        teststat = quad4sum - (quad1sum + quad2sum + quad3sum + quad4sum)/4
    else
        throw(DomainError("Invalid Quadrant"))
    end
    return teststat
end

function createdist(quadcarlo, quad)  #creates a 10000 point distribution random test statistics for a certain quadrant
    distribution = []
    for z in 1:10000
        randquadcarlo = randomizequadcarlo(quadcarlo)
        push!(distribution, testdata(randquadcarlo, quad))
    end
    return distribution
end

function totaltime(quadcarlo, quad) #looks at a dataframe and returns the total time spent in a specified quadrant
    quad1, quad2, quad3, quad4 = groupby(quadcarlo, [:quadrant], sort = true)
    quad1sum, quad2sum, quad3sum, quad4sum = sum(quad1[!,2]), sum(quad2[!,2]), sum(quad3[!,2]), sum(quad4[!,2])
    if(quad=="quad 1")
        return quad1sum
    elseif(quad=="quad 2")
        return quad2sum
    elseif(quad=="quad 3")
        return quad3sum
    elseif(quad=="quad 4")
        return quad4sum
    else
        throw(DomainError("Invalid Quadrant"))
    end
end

function calculatepvalue(quadrantdf, quad)
    quadcarlo = createquadcarlo(quadrantdf)
    expvalue = testdata(quadcarlo, quad)
    dist = createdist(quadcarlo, quad)
    function twosided(x)
        if(x > 0)
            return x >= abs(expvalue)
        elseif(x <= 0)
            return abs(x) >= abs(expvalue)
        else
            throw(DomainError("Not a valid input"))
        end
    end
    expormore = filter(twosided, dist)
    pvalue = length(expormore)/length(dist)
    return pvalue
end

end


text:
using DeferredFutures
using Base.Test

@testset "Comparison" begin
    f = Future()

    @test DeferredFuture(f) == DeferredFuture(f)
    @test hash(DeferredFuture(f)) == hash(DeferredFuture(f))
end

@testset "Distributed" begin
    top = myid()
    bottom = addprocs(1)[1]
    @everywhere using DeferredFutures

    try
        val = "hello"
        df = DeferredFuture(top)
        @test !isready(df)

        fut = remotecall_wait(bottom, df) do dfr
            put!(dfr, val)
        end
        @test fetch(fut) == df
        @test isready(df)
        @test fetch(df) == val
        @test wait(df) == df
        @test_throws ErrorException put!(df, val)

        @test df.outer.where == top
        @test fetch(df.outer).where == bottom
    finally
        rmprocs(bottom)
    end
end

@testset "Allocation" begin
    rand_size = 800000000  # sizeof(rand(10000, 10000))
    gc()
    main_size = Base.summarysize(Main)

    top = myid()
    bottom = addprocs(1)[1]
    @everywhere using DeferredFutures

    try
        df = DeferredFuture(top)

        remote_size = remotecall_fetch(bottom, df) do dfr
            gc()
            main_size = Base.summarysize(Main)
            put!(dfr, rand(10000, 10000))
            main_size
        end

        gc()
        @test Base.summarysize(Main) < main_size + rand_size

        remote_size_new = remotecall_fetch(bottom) do
            gc()
            Base.summarysize(Main)
        end

        @test remote_size_new >= remote_size + rand_size
    finally
        rmprocs(bottom)
    end
end

@testset "Transfer" begin
    rand_size = 800000000  # sizeof(rand(10000, 10000))
    gc()
    main_size = Base.summarysize(Main)

    top = myid()
    left, right = addprocs(2)
    @everywhere using DeferredFutures

    try
        df = DeferredFuture(top)

        left_remote_size = remotecall_fetch(left, df) do dfr
            gc()
            main_size = Base.summarysize(Main)
            put!(dfr, rand(10000, 10000))
            main_size
        end

        right_remote_size = remotecall_fetch(right, df) do dfr
            gc()
            main_size = Base.summarysize(Main)
            global data = fetch(dfr)
            main_size
        end

        gc()
        @test Base.summarysize(Main) < main_size + rand_size

        right_remote_size_new = remotecall_fetch(right) do
            gc()
            Base.summarysize(Main)
        end

        @test right_remote_size_new >= right_remote_size + rand_size
    finally
        rmprocs([left, right])
    end
end


text:
#############################################################
# This file was automatically generated on 2021-12-15.      #
#                                                           #
# JuliaPy Bindings Version 2.1.29                           #
#                                                           #
# If you have a bugfix for this file and want to commit it, #
# please fix the bug in the generator. You can find a link  #
# to the generators git repository on tinkerforge.com       #
#############################################################

include("brick_dc.jl")
include("brick_hat.jl")
include("brick_hat_zero.jl")
include("brick_imu.jl")
include("brick_imu_v2.jl")
include("brick_master.jl")
include("brick_red.jl")
include("brick_servo.jl")
include("brick_silent_stepper.jl")
include("brick_stepper.jl")
include("bricklet_accelerometer.jl")
include("bricklet_accelerometer_v2.jl")
include("bricklet_air_quality.jl")
include("bricklet_ambient_light.jl")
include("bricklet_ambient_light_v2.jl")
include("bricklet_ambient_light_v3.jl")
include("bricklet_analog_in.jl")
include("bricklet_analog_in_v2.jl")
include("bricklet_analog_in_v3.jl")
include("bricklet_analog_out.jl")
include("bricklet_analog_out_v2.jl")
include("bricklet_analog_out_v3.jl")
include("bricklet_barometer.jl")
include("bricklet_barometer_v2.jl")
include("bricklet_can.jl")
include("bricklet_can_v2.jl")
include("bricklet_co2.jl")
include("bricklet_co2_v2.jl")
include("bricklet_color.jl")
include("bricklet_color_v2.jl")
include("bricklet_compass.jl")
include("bricklet_current12.jl")
include("bricklet_current25.jl")
include("bricklet_dc_v2.jl")
include("bricklet_distance_ir.jl")
include("bricklet_distance_ir_v2.jl")
include("bricklet_distance_us.jl")
include("bricklet_distance_us_v2.jl")
include("bricklet_dmx.jl")
include("bricklet_dual_button.jl")
include("bricklet_dual_button_v2.jl")
include("bricklet_dual_relay.jl")
include("bricklet_dust_detector.jl")
include("bricklet_e_paper_296x128.jl")
include("bricklet_energy_monitor.jl")
include("bricklet_gps.jl")
include("bricklet_gps_v2.jl")
include("bricklet_hall_effect.jl")
include("bricklet_hall_effect_v2.jl")
include("bricklet_humidity.jl")
include("bricklet_humidity_v2.jl")
include("bricklet_imu_v3.jl")
include("bricklet_industrial_analog_out.jl")
include("bricklet_industrial_analog_out_v2.jl")
include("bricklet_industrial_counter.jl")
include("bricklet_industrial_digital_in_4.jl")
include("bricklet_industrial_digital_in_4_v2.jl")
include("bricklet_industrial_digital_out_4.jl")
include("bricklet_industrial_digital_out_4_v2.jl")
include("bricklet_industrial_dual_0_20ma.jl")
include("bricklet_industrial_dual_0_20ma_v2.jl")
include("bricklet_industrial_dual_ac_relay.jl")
include("bricklet_industrial_dual_analog_in.jl")
include("bricklet_industrial_dual_analog_in_v2.jl")
include("bricklet_industrial_dual_relay.jl")
include("bricklet_industrial_ptc.jl")
include("bricklet_industrial_quad_relay.jl")
include("bricklet_industrial_quad_relay_v2.jl")
include("bricklet_io16.jl")
include("bricklet_io16_v2.jl")
include("bricklet_io4.jl")
include("bricklet_io4_v2.jl")
include("bricklet_isolator.jl")
include("bricklet_joystick.jl")
include("bricklet_joystick_v2.jl")
include("bricklet_laser_range_finder.jl")
include("bricklet_laser_range_finder_v2.jl")
include("bricklet_lcd_128x64.jl")
include("bricklet_lcd_16x2.jl")
include("bricklet_lcd_20x4.jl")
include("bricklet_led_strip.jl")
include("bricklet_led_strip_v2.jl")
include("bricklet_line.jl")
include("bricklet_linear_poti.jl")
include("bricklet_linear_poti_v2.jl")
include("bricklet_load_cell.jl")
include("bricklet_load_cell_v2.jl")
include("bricklet_moisture.jl")
include("bricklet_motion_detector.jl")
include("bricklet_motion_detector_v2.jl")
include("bricklet_motorized_linear_poti.jl")
include("bricklet_multi_touch.jl")
include("bricklet_multi_touch_v2.jl")
include("bricklet_nfc.jl")
include("bricklet_nfc_rfid.jl")
include("bricklet_oled_128x64.jl")
include("bricklet_oled_128x64_v2.jl")
include("bricklet_oled_64x48.jl")
include("bricklet_one_wire.jl")
include("bricklet_outdoor_weather.jl")
include("bricklet_particulate_matter.jl")
include("bricklet_performance_dc.jl")
include("bricklet_piezo_buzzer.jl")
include("bricklet_piezo_speaker.jl")
include("bricklet_piezo_speaker_v2.jl")
include("bricklet_ptc.jl")
include("bricklet_ptc_v2.jl")
include("bricklet_real_time_clock.jl")
include("bricklet_real_time_clock_v2.jl")
include("bricklet_remote_switch.jl")
include("bricklet_remote_switch_v2.jl")
include("bricklet_rgb_led.jl")
include("bricklet_rgb_led_button.jl")
include("bricklet_rgb_led_matrix.jl")
include("bricklet_rgb_led_v2.jl")
include("bricklet_rotary_encoder.jl")
include("bricklet_rotary_encoder_v2.jl")
include("bricklet_rotary_poti.jl")
include("bricklet_rotary_poti_v2.jl")
include("bricklet_rs232.jl")
include("bricklet_rs232_v2.jl")
include("bricklet_rs485.jl")
include("bricklet_segment_display_4x7.jl")
include("bricklet_segment_display_4x7_v2.jl")
include("bricklet_servo_v2.jl")
include("bricklet_silent_stepper_v2.jl")
include("bricklet_solid_state_relay.jl")
include("bricklet_solid_state_relay_v2.jl")
include("bricklet_sound_intensity.jl")
include("bricklet_sound_pressure_level.jl")
include("bricklet_temperature.jl")
include("bricklet_temperature_ir.jl")
include("bricklet_temperature_ir_v2.jl")
include("bricklet_temperature_v2.jl")
include("bricklet_thermal_imaging.jl")
include("bricklet_thermocouple.jl")
include("bricklet_thermocouple_v2.jl")
include("bricklet_tilt.jl")
include("bricklet_uv_light.jl")
include("bricklet_uv_light_v2.jl")
include("bricklet_voltage.jl")
include("bricklet_voltage_current.jl")
include("bricklet_voltage_current_v2.jl")
include("bricklet_xmc1400_breakout.jl")

export get_device_type
function get_device_type(device_identifier::Integer)
    device_types = Dict{Integer, String}(
        11 => "BrickDC",
        111 => "BrickHAT",
        112 => "BrickHATZero",
        16 => "BrickIMU",
        18 => "BrickIMUV2",
        13 => "BrickMaster",
        17 => "BrickRED",
        14 => "BrickServo",
        19 => "BrickSilentStepper",
        15 => "BrickStepper",
        250 => "BrickletAccelerometer",
        2130 => "BrickletAccelerometerV2",
        297 => "BrickletAirQuality",
        21 => "BrickletAmbientLight",
        259 => "BrickletAmbientLightV2",
        2131 => "BrickletAmbientLightV3",
        219 => "BrickletAnalogIn",
        251 => "BrickletAnalogInV2",
        295 => "BrickletAnalogInV3",
        220 => "BrickletAnalogOut",
        256 => "BrickletAnalogOutV2",
        2115 => "BrickletAnalogOutV3",
        221 => "BrickletBarometer",
        2117 => "BrickletBarometerV2",
        270 => "BrickletCAN",
        2107 => "BrickletCANV2",
        262 => "BrickletCO2",
        2147 => "BrickletCO2V2",
        243 => "BrickletColor",
        2128 => "BrickletColorV2",
        2153 => "BrickletCompass",
        23 => "BrickletCurrent12",
        24 => "BrickletCurrent25",
        2165 => "BrickletDCV2",
        25 => "BrickletDistanceIR",
        2125 => "BrickletDistanceIRV2",
        229 => "BrickletDistanceUS",
        299 => "BrickletDistanceUSV2",
        285 => "BrickletDMX",
        230 => "BrickletDualButton",
        2119 => "BrickletDualButtonV2",
        26 => "BrickletDualRelay",
        260 => "BrickletDustDetector",
        2146 => "BrickletEPaper296x128",
        2152 => "BrickletEnergyMonitor",
        222 => "BrickletGPS",
        276 => "BrickletGPSV2",
        240 => "BrickletHallEffect",
        2132 => "BrickletHallEffectV2",
        27 => "BrickletHumidity",
        283 => "BrickletHumidityV2",
        2161 => "BrickletIMUV3",
        258 => "BrickletIndustrialAnalogOut",
        2116 => "BrickletIndustrialAnalogOutV2",
        293 => "BrickletIndustrialCounter",
        223 => "BrickletIndustrialDigitalIn4",
        2100 => "BrickletIndustrialDigitalIn4V2",
        224 => "BrickletIndustrialDigitalOut4",
        2124 => "BrickletIndustrialDigitalOut4V2",
        228 => "BrickletIndustrialDual020mA",
        2120 => "BrickletIndustrialDual020mAV2",
        2162 => "BrickletIndustrialDualACRelay",
        249 => "BrickletIndustrialDualAnalogIn",
        2121 => "BrickletIndustrialDualAnalogInV2",
        284 => "BrickletIndustrialDualRelay",
        2164 => "BrickletIndustrialPTC",
        225 => "BrickletIndustrialQuadRelay",
        2102 => "BrickletIndustrialQuadRelayV2",
        28 => "BrickletIO16",
        2114 => "BrickletIO16V2",
        29 => "BrickletIO4",
        2111 => "BrickletIO4V2",
        2122 => "BrickletIsolator",
        210 => "BrickletJoystick",
        2138 => "BrickletJoystickV2",
        255 => "BrickletLaserRangeFinder",
        2144 => "BrickletLaserRangeFinderV2",
        298 => "BrickletLCD128x64",
        211 => "BrickletLCD16x2",
        212 => "BrickletLCD20x4",
        231 => "BrickletLEDStrip",
        2103 => "BrickletLEDStripV2",
        241 => "BrickletLine",
        213 => "BrickletLinearPoti",
        2139 => "BrickletLinearPotiV2",
        253 => "BrickletLoadCell",
        2104 => "BrickletLoadCellV2",
        232 => "BrickletMoisture",
        233 => "BrickletMotionDetector",
        292 => "BrickletMotionDetectorV2",
        267 => "BrickletMotorizedLinearPoti",
        234 => "BrickletMultiTouch",
        2129 => "BrickletMultiTouchV2",
        286 => "BrickletNFC",
        246 => "BrickletNFCRFID",
        263 => "BrickletOLED128x64",
        2112 => "BrickletOLED128x64V2",
        264 => "BrickletOLED64x48",
        2123 => "BrickletOneWire",
        288 => "BrickletOutdoorWeather",
        2110 => "BrickletParticulateMatter",
        2156 => "BrickletPerformanceDC",
        214 => "BrickletPiezoBuzzer",
        242 => "BrickletPiezoSpeaker",
        2145 => "BrickletPiezoSpeakerV2",
        226 => "BrickletPTC",
        2101 => "BrickletPTCV2",
        268 => "BrickletRealTimeClock",
        2106 => "BrickletRealTimeClockV2",
        235 => "BrickletRemoteSwitch",
        289 => "BrickletRemoteSwitchV2",
        271 => "BrickletRGBLED",
        282 => "BrickletRGBLEDButton",
        272 => "BrickletRGBLEDMatrix",
        2127 => "BrickletRGBLEDV2",
        236 => "BrickletRotaryEncoder",
        294 => "BrickletRotaryEncoderV2",
        215 => "BrickletRotaryPoti",
        2140 => "BrickletRotaryPotiV2",
        254 => "BrickletRS232",
        2108 => "BrickletRS232V2",
        277 => "BrickletRS485",
        237 => "BrickletSegmentDisplay4x7",
        2137 => "BrickletSegmentDisplay4x7V2",
        2157 => "BrickletServoV2",
        2166 => "BrickletSilentStepperV2",
        244 => "BrickletSolidStateRelay",
        296 => "BrickletSolidStateRelayV2",
        238 => "BrickletSoundIntensity",
        290 => "BrickletSoundPressureLevel",
        216 => "BrickletTemperature",
        217 => "BrickletTemperatureIR",
        291 => "BrickletTemperatureIRV2",
        2113 => "BrickletTemperatureV2",
        278 => "BrickletThermalImaging",
        266 => "BrickletThermocouple",
        2109 => "BrickletThermocoupleV2",
        239 => "BrickletTilt",
        265 => "BrickletUVLight",
        2118 => "BrickletUVLightV2",
        218 => "BrickletVoltage",
        227 => "BrickletVoltageCurrent",
        2105 => "BrickletVoltageCurrentV2",
        279 => "BrickletXMC1400Breakout",
    )

    return device_types[device_identifier]
end

export create_device
function create_device(device_identifier::Integer, uid::String, ipcon::IPConnection)
    return get_device_type(device_identifier)(uid, ipcon)
end


text:
using BiodiversityObservationNetworks
using NeutralLandscapes
using Plots

dims = (200, 200)
uncert = rand(MidpointDisplacement(0.95), dims)
uncertplt = heatmap(uncert', frame=:box, axis=:none, colorbartitle="SDM Uncertainty")

npts = 30

weight = 10
noweighting = rand(BalancedAcceptance(numpoints=npts, α=-weight), uncert)
strongweighting = rand(BalancedAcceptance(numpoints=npts, α=weight), uncert)

lowplt = heatmap(uncert',titlefontsize=8, colorbar=:none, title="Weighted toward low uncertainty", size=(350,300), frame=:box, axis=:none, colorbartitle="SDM Uncertainty")
scatter!(lowplt, [x[1] for x in noweighting], [x[2] for x in noweighting], ms=2.5, mc=:white, label="")

hiplt = heatmap(uncert', colorbar=:none, titlefontsize=8, title="Weighted toward high uncertainty", size=(350,300), frame=:box, axis=:none, colorbartitle="SDM Uncertainty")
scatter!(hiplt, [x[1] for x in strongweighting], [x[2] for x in strongweighting], ms=2.5, mc=:white, label="")

plot(uncertplt, lowplt, hiplt, size=(1000, 250), dpi=300, layout=(1,3) )

savefig("balanced_acceptance_with_weights.png")


text:
@testset "ASCII Conversion" begin
    n = readnimrod("test_hexed_composite_2BI.dat")
    test_ascii_file = "test_ascii_tmp.nim"
    toasciigrid(n, test_ascii_file)

    expectedFileContents = """nrows 8
ncols 8
xllcenter 5.0
yllcenter -75.0
cellsize 10.0
nodata_value -1
5 5 -1 -1 -1 -1 -1 -1
-1 5 -1 -1 -1 -1 -1 -1
-1 -1 5 -1 -1 -1 -1 -1
-1 -1 -1 5 -1 -1 -1 -1
-1 -1 -1 -1 5 -1 -1 -1
-1 -1 -1 -1 -1 5 -1 -1
-1 -1 -1 -1 -1 -1 5 -1
-1 -1 -1 -1 -1 -1 -1 5
"""
    @test readstring(test_ascii_file) == expectedFileContents
    rm(test_ascii_file)
end


text:
#Symmetric and Hermitian matrices
immutable Symmetric{T} <: AbstractMatrix{T}
    S::Matrix{T}
    uplo::Char
end
Symmetric(A::Matrix, uplo::Symbol=:U) = (chksquare(A);Symmetric(A, string(uplo)[1]))
immutable Hermitian{T} <: AbstractMatrix{T}
    S::Matrix{T}
    uplo::Char
end
Hermitian(A::Matrix, uplo::Symbol=:U) = (chksquare(A);Hermitian(A, string(uplo)[1]))
typealias HermOrSym{T} Union(Hermitian{T}, Symmetric{T})

size(A::HermOrSym, args...) = size(A.S, args...)
getindex(A::HermOrSym, i::Integer, j::Integer) = (A.uplo == 'U') == (i < j) ? getindex(A.S, i, j) : conj(getindex(A.S, j, i))
full(A::Symmetric) = copytri!(A.S, A.uplo)
full(A::Hermitian) = copytri!(A.S, A.uplo, true)
convert{T}(::Type{Symmetric{T}},A::Symmetric) = Symmetric(convert(Matrix{T},A.S),A.uplo)
convert{T}(::Type{Hermitian{T}},A::Hermitian) = Hermitian(convert(Matrix{T},A.S),A.uplo)
copy(A::Symmetric) = Symmetric(copy(A.S),A.uplo)
copy(A::Hermitian) = Hermitian(copy(A.S),A.uplo)
ishermitian(A::Hermitian) = true
ishermitian{T<:Real}(A::Symmetric{T}) = true
ishermitian{T<:Complex}(A::Symmetric{T}) = all(imag(A.S) .== 0)
issym{T<:Real}(A::Hermitian{T}) = true
issym{T<:Complex}(A::Hermitian{T}) = all(imag(A.S) .== 0)
issym(A::Symmetric) = true
transpose(A::Symmetric) = A
ctranspose(A::Hermitian) = A

*(A::HermOrSym, B::HermOrSym) = full(A)*full(B)
*(A::HermOrSym, B::StridedMatrix) = full(A)*B
*(A::StridedMatrix, B::HermOrSym) = A*full(B)

factorize(A::HermOrSym) = bkfact(A.S, symbol(A.uplo), issym(A))
\(A::HermOrSym, B::StridedVecOrMat) = \(bkfact(A.S, symbol(A.uplo), issym(A)), B)

eigfact!{T<:BlasReal}(A::Symmetric{T}) = Eigen(LAPACK.syevr!('V', 'A', A.uplo, A.S, 0.0, 0.0, 0, 0, -1.0)...)
eigfact!{T<:BlasComplex}(A::Hermitian{T}) = Eigen(LAPACK.syevr!('V', 'A', A.uplo, A.S, 0.0, 0.0, 0, 0, -1.0)...)
eigfact{T<:BlasFloat}(A::HermOrSym{T}) = eigfact!(copy(A))
eigfact{T}(A::HermOrSym{T}) = (S = promote_type(Float32,typeof(one(T)/norm(one(T)))); S != T ? eigfact!(convert(typeof(A).name.primary{S}, A)) : eigfact!(copy(A)))
eigvals!{T<:BlasReal}(A::Symmetric{T}, il::Int=1, ih::Int=size(A,1)) = LAPACK.syevr!('N', 'I', A.uplo, A.S, 0.0, 0.0, il, ih, -1.0)[1]
eigvals!{T<:BlasReal}(A::Symmetric{T}, vl::Real, vh::Real) = LAPACK.syevr!('N', 'V', A.uplo, A.S, vl, vh, 0, 0, -1.0)[1]
eigvals!{T<:BlasComplex}(A::Hermitian{T}, il::Int=1, ih::Int=size(A,1)) = LAPACK.syevr!('N', 'I', A.uplo, A.S, 0.0, 0.0, il, ih, -1.0)[1]
eigvals!{T<:BlasComplex}(A::Hermitian{T}, vl::Real, vh::Real) = LAPACK.syevr!('N', 'V', A.uplo, A.S, vl, vh, 0, 0, -1.0)[1]
eigvals{T<:BlasFloat}(A::HermOrSym{T},l::Real=1,h::Real=size(A,1)) = eigvals!(copy(A),l,h)
eigvals{T}(A::HermOrSym{T},l::Real=1,h::Real=size(A,1)) = (S = promote_type(Float32,typeof(one(T)/norm(one(T)))); S != T ? eigvals!(convert(typeof(A).name.primary{S}, A, l, h)) : eigvals!(copy(A), l, h))
eigmax(A::HermOrSym) = eigvals(A, size(A, 1), size(A, 1))[1]
eigmin(A::HermOrSym) = eigvals(A, 1, 1)[1]

function eigfact!{T<:BlasReal}(A::Symmetric{T}, B::Symmetric{T})
    vals, vecs, _ = LAPACK.sygvd!(1, 'V', A.uplo, A.S, B.uplo == A.uplo ? B.S : B.S')
    GeneralizedEigen(vals, vecs)
end
function eigfact!{T<:BlasComplex}(A::Hermitian{T}, B::Hermitian{T})
    vals, vecs, _ = LAPACK.sygvd!(1, 'V', A.uplo, A.S, B.uplo == A.uplo ? B.S : B.S')
    GeneralizedEigen(vals, vecs)
end
eigvals!{T<:BlasReal}(A::Symmetric{T}, B::Symmetric{T}) = LAPACK.sygvd!(1, 'N', A.uplo, A.S, B.uplo == A.uplo ? B.S : B.S')[1]
eigvals!{T<:BlasComplex}(A::Hermitian{T}, B::Hermitian{T}) = LAPACK.sygvd!(1, 'N', A.uplo, A.S, B.uplo == A.uplo ? B.S : B.S')[1]

#Matrix-valued functions
expm(A::HermOrSym) = (F = eigfact(A); F.vectors*Diagonal(exp(F.values))*F.vectors')
function sqrtm(A::HermOrSym)
    F = eigfact(A)
    isposdef(F) && return F.vectors*Diagonal(sqrt(F.values))*F.vectors'
    return F.vectors*Diagonal(sqrt(complex(F.values)))*F.vectors'
end

max_stars_repo_name: alyst/julia

max_issues_repo_name: alyst/julia

max_forks_repo_name: alyst/julia


text:
#
# Function that computes the fit (least square sense) of a SIRD simulation to the
# trajectory data
#
# This function modifies the sird_traj data of input, which is to be taken as an
# auxiliary vector only
#

# If the sum of all data is used (do not specify ipop):

function fiteval!(x0 :: Vector{Float64}, md_traj :: MDTraj,
                  sird_input :: SIRDInput, sird_traj :: SIRDTraj)

  # perform SIRD simulation
  sird!(x0,sird_input,sird_traj)

  mU = md_traj.U
  mS = md_traj.S
  mD = md_traj.D
  mI = md_traj.I

  sU = sird_traj.U
  sS = sird_traj.S
  sD = sird_traj.D
  sI = sird_traj.I

  fit = 0.
  for i in 1:md_traj.nsteps
     fit = fit + (mU[i]-sU[i])^2 +
                 (mS[i]-sS[i])^2 +
                 (mD[i]-sD[i])^2 +
                 (mI[i]-sI[i])^2
  end
  fit = sqrt(fit/(4*md_traj.nsteps))

  return fit

end

# If a single set is used (defined by ipop):

function fiteval!(x0 :: Vector{Float64}, md_traj :: MDTraj,
                  sird_input :: SIRDInput, sird_traj :: SIRDTraj,
                  ipop :: Int64)

  # perform SIRD simulation
  sird!(x0,sird_input,sird_traj)

  # Which data will be used
  if ipop == 0  # Susceptible
    md = md_traj.U
    sird = sird_traj.U
  elseif ipop == 1 # Sick
    md = md_traj.S
    sird = sird_traj.S
  elseif ipop == 2 # Dead
    md = md_traj.D
    sird = sird_traj.D
  elseif ipop == 3 # Immune
    md = md_traj.I
    sird = sird_traj.I
  else
    error(" Invalid value for ipop = ", ipop)
  end

  fit = 0.
  for i in 1:md_traj.nsteps
    fit = fit + (md[i]-sird[i])^2
  end
  fit = sqrt(fit/md_traj.nsteps)

  return fit

end


text:
ado_id_to_replace="PEDESTRIAN/236"
ego_pos_init_vec=nothing;
ego_pos_goal_vec=nothing;
target_speed=nothing;


text:
using Optim: AbstractOptimizer, Options, MultivariateOptimizationResults

abstract type ObjectiveFunction end
abstract type PenaltyFunction end

struct GrapeResult{D,T<:AbstractMatrix,P<:MultivariateOptimizationResults}
    Ut::Operator{T,D}
    Ui::Operator{Matrix{ComplexF64},D}
    ui::Matrix{Float64}
    fi::Float64
    Uf::Operator{Matrix{ComplexF64},D}
    uf::Matrix{Float64}
    ff::Float64
    t::Float64
    optim_res::P
end

grape(O::ObjectiveFunction,ui::Array,opt::Options) = grape(O,ui,ConjugateGradient(),opt)

grape(O::ObjectiveFunction,P::PenaltyFunction,ui::Array,opt::Options) = grape(O,P,ui,ConjugateGradient(),opt)

function grape(O::ObjectiveFunction,ui::Array,method::AbstractOptimizer=ConjugateGradient(),opt::Options=Options())
    f(u) = objective(O,u)
    g!(fp,u) = gradient!(O,fp,u)
    return grape(f,g!,ui,method,opt)
end

function grape(O::ObjectiveFunction,P::PenaltyFunction,ui::Array,method::AbstractOptimizer=ConjugateGradient(),opt::Options=Options())
    f(u) = objective(O,u) + objective(P,u)
    g!(fp,u) = (gradient!(O,fp,u); gradient!(P,fp,u))
    return grape(f,g!,ui,method,opt)
end

function grape(f::Function,g!::Function,ui::Array,method::AbstractOptimizer=ConjugateGradient(),opt::Options=Options())
    # Build OnceDifferentiable object from objective and gradient functions
    od = OnceDifferentiable(f,g!,vec(ui))
    # Calculate initial fidelity and initial propagator
    fi = 1 - objective(f.O,ui)
    Ui = current_propagator(f.O)
    # Run optimization
    optim_res = optimize(od,vec(ui),method,opt)
    # Grab optimized control amplitudes
    uf = reshape(Optim.minimizer(optim_res),size(ui,1),size(ui,2))
    # Calculate final fidelity (cheap since propagators are already calculated)
    ff = 1 - objective(f.O,uf)
    # Grab final evolution propagator
    Uf = current_propagator(f.O)
    # Target evolution operator
    Ut = target_propagator(f.O)
    # Return result
    return GrapeResult(Ut,Ui,ui[:,:],fi,Uf,uf,ff,f.O.δt*size(ui,1),optim_res)
end

calc_fprops!(O,u) = calc_fprops!(O.U,O.X,O.D,O.V,u,O.δt,O.Hd,O.Hc,O.H,O.u_last)

function calc_fprops!(U,X,D,V,u,δt,Hd,Hc,H,u_last)
    # If the control amplitudes u did not change, return
    u == u_last && return nothing
    # Otherwise calculate each individual propagator
    n = length(U); m = length(Hc)
    for j = 1:n # loop over timeslices
        H.data .= -δt.*Hd # copy drift (constant) Hamiltonian
        for k = 1:m # loop over controls
            H.data .+= -δt*u[(k-1)*n+j].*Hc[k] # add in time dependent parts
        end
        # Calculate U and store eigenvectors and eigenvalues for this timeslice
        expim!(U[j],H,D[j],V[j],X[1])
    end
    # Calculate forward propagators (cumulative product of U)
    copyto!(X[1],U[1])
    for i = 2:n
        mul!(X[i],U[i],X[i-1])
    end
    # Save control amplitudes
    copyto!(u_last,u)
    return nothing
end

calc_bprops!(O) = calc_bprops!(O.P,O.U,O.Ut)

function calc_bprops!(P::Vector{<:AbstractMatrix},U,Ut)
    # Calculate backward propagators
    for i = length(P)-1:-1:1
        mul!(P[i],adjoint(U[i+1]),P[i+1])
    end
    return nothing
end

function calc_bprops!(P::Vector{<:NTuple{M,AbstractMatrix}},U,Ut) where M
    # Calculate backward propagators
    for i = length(P)-1:-1:1, m = 1:M
        mul!(P[i][m],adjoint(U[i+1]),P[i+1][m])
    end
    return nothing
end

Jmat!(O,k,j) = Jmat!(O.Jkj,O.Hc[k],O.cisDj,O.D[j],O.V[j],O.δt,O.A)

function Jmat!(Jkj,Hck,cisDj,Dj,Vj,δt,A)
    # Jₖⱼ = Vⱼ*((Vⱼ'*Hₖ*Vⱼ).*Λⱼ)*Vⱼ'
    # Λⱼ[l,m] = λl≈λm ? -1im*δt*cis(λl) : -δt*(cis(λl)-cis(λm))/(λl-λm)
    mul!(A,adjoint(Vj),Hck)
    mul!(Jkj,A,Vj)
    _Jmathermprod!(Jkj,cisDj,Dj,δt)
    mul!(A,Vj,Jkj)
    mul!(Jkj,A,adjoint(Vj))
    return nothing
end

function _Jmathermprod!(J,cisDj,Dj,δt)
    for m = 1:length(Dj), l = 1:length(Dj)
        λl, λm = Dj[l], Dj[m]
        cisλl, cisλm = cisDj[l], cisDj[m]
        J[l,m] *= abs(λl-λm)<1E-10 ? -1im*δt*cisλl : -δt*(cisλl-cisλm)/(λl-λm)
    end
    return nothing
end

function current_propagator(O)
    # Return the last propagator calculated
    return Operator(copy(O.X[end]),O.dims)
end

function target_propagator(O)
    # Return the target propagator
    return Operator(O.Ut isa NTuple ? sum(O.Ut) : copy(O.Ut), O.dims)
end

function step_sample(fun,params,tspan,steps)
    t₁, t₂ = tspan
    Δt = (t₂-t₁)/steps
    u = Vector{typeof(fun(t₁,params))}(steps)
    for i = 1:steps
        tᵢ = t₁ + (i-1)*Δt
        u[i] = (fun(tᵢ,params) + fun(tᵢ+0.5Δt,params) + fun(tᵢ+Δt,params))/3
    end
    return u
end


text:
using MultivariateStats
using Base.Test

srand(15678)

## icagfun

f = icagfun(:tanh)
u, v = evaluate(f, 1.5)
@test_approx_eq u 0.905148253644866438242
@test_approx_eq v 0.180706638923648530597

f = icagfun(:tanh, 1.5)
u, v = evaluate(f, 1.2)
@test_approx_eq u 0.946806012846268289646
@test_approx_eq v 0.155337561057228069719

f = icagfun(:gaus)
u, v = evaluate(f, 1.5)
@test_approx_eq u 0.486978701037524594696
@test_approx_eq v -0.405815584197937162246

## data

# sources
n = 1000
k = 3
m = 8

t = linspace(0.0, 10.0, n)
s1 = sin(t * 2)
s2 = s2 = 1.0 - 2.0 * Bool[isodd(ifloor(_ / 3)) for _ in t]
s3 = Float64[mod(_, 5.0) for _ in t]

s1 += 0.1 * randn(n)
s2 += 0.1 * randn(n)
s3 += 0.1 * randn(n)

S = hcat(s1, s2, s3)'
@assert size(S) == (k, n)

A = randn(m, k)

X = A * S
mv = vec(mean(X,2))
@assert size(X) == (m, n)
C = cov(X; vardim=2)

# FastICA

M = fit(ICA, X, k; do_whiten=false)
@test isa(M, ICA)
@test indim(M) == m
@test outdim(M) == k
@test mean(M) == mv
W = M.W
@test_approx_eq transform(M, X) W' * (X .- mv)
@test_approx_eq W'W eye(k)

M = fit(ICA, X, k; do_whiten=true)
@test isa(M, ICA)
@test indim(M) == m
@test outdim(M) == k
@test mean(M) == mv
W = M.W
@test_approx_eq W'C * W eye(k)


text:
#=
author: mingchuno
smart permutation inspired by python-high-perf
TODO: parallel & mutil core, compile?, README and output
=#

BASE = 22
PPP = sum(map(x -> BASE^x,0:2)) # generic PPP
println("PPP is:$(PPP)")
array = sort(collect(setdiff(Set(0:BASE-1), Set(1))))
sarray = Set(array)
count = 0

@time for last4_arr in combinations(array, 4), efgh_arr in permutations(last4_arr)
  @inbounds e = efgh_arr[1]
  if e == 0 continue end # fail fast

  @inbounds f = efgh_arr[2]

  @inbounds g = efgh_arr[3]
  if g == 0 continue end # fail fast

  @inbounds h = efgh_arr[4]

  ef = e*BASE + f
  gh = g*BASE + h

  if ef + gh == PPP
    remain_arr = collect(setdiff(sarray, Set(last4_arr)))
    for first4_arr in combinations(remain_arr, 4), abcd_arr in permutations(first4_arr)
      @inbounds a = abcd_arr[1]
      if a == 0 continue end # fail fast

      @inbounds b = abcd_arr[2]

      @inbounds c = abcd_arr[3]
      if c == 0 continue end # fail fast

      @inbounds d = abcd_arr[4]
      ab = a*BASE + b
      cd = c*BASE + d
      if ab - cd == ef
        println("($(base(BASE,ab)) - $(base(BASE,cd)) = $(base(BASE, ef))) + $(base(BASE,gh)) = 111")
        count += 1
      end
    end
  end
end

println("Total number of soln is:$(count)")

max_stars_repo_name: jackhftang/ABCDEFGHPPP

max_issues_repo_name: jackhftang/ABCDEFGHPPP

max_forks_repo_name: jackhftang/ABCDEFGHPPP


text:
module CountClick

using Gtk

function julia_main()::Cint
    win = GtkWindow("Count Click")
    v = GtkBox(:v)
    l = GtkLabel("You clicked 0 times.")
    b = GtkButton("Click!")
    push!(win, v)
    push!(v, l)
    push!(v, b)
    set_gtk_property!(v, :expand, l, true)

    # callback function
    global nclick = 0
    function click(x)
        nclick += 1
        set_gtk_property!(l, :label, "You clicked $(nclick) times.")
    end

    # connect button and function
    signal_connect(click, b, "clicked")
    showall(win)

    if !isinteractive()
        c = Condition()
        signal_connect(win, :destroy) do widget
            notify(c)
        end
        wait(c)
    end

    return 0
end

end

max_stars_repo_name: goropikari/JuliaGUIAppExample

max_issues_repo_name: goropikari/JuliaGUIAppExample

max_forks_repo_name: goropikari/JuliaGUIAppExample


text:
using oneAPI.oneL0

# constructors
voidptr_a = ZePtr{Cvoid}(Int(0xDEADBEEF))
@test reinterpret(Ptr{Cvoid}, voidptr_a) == Ptr{Cvoid}(Int(0xDEADBEEF))

# getters
@test eltype(voidptr_a) == Cvoid

# comparisons
voidptr_b = ZePtr{Cvoid}(Int(0xCAFEBABE))
@test voidptr_a != voidptr_b

@testset "conversions" begin

# between regular and CUDA pointers
@test_throws ArgumentError convert(Ptr{Cvoid}, voidptr_a)

# between CUDA pointers
intptr_a = ZePtr{Int}(Int(0xDEADBEEF))
@test convert(typeof(intptr_a), voidptr_a) == intptr_a

# convert back and forth from UInt
intptr_b = ZePtr{Int}(Int(0xDEADBEEF))
@test convert(UInt, intptr_b) == 0xDEADBEEF
@test convert(ZePtr{Int}, Int(0xDEADBEEF)) == intptr_b
@test Int(intptr_b) == Int(0xDEADBEEF)

# pointer arithmetic
intptr_c = ZePtr{Int}(Int(0xDEADBEEF))
intptr_d = 2 + intptr_c
@test isless(intptr_c, intptr_d)
@test intptr_d - intptr_c == 2
@test intptr_d - 2 == intptr_c
end

@testset "GPU or CPU integration" begin

a = [1]
ccall(:clock, Nothing, (Ptr{Int},), a)
@test_throws Exception ccall(:clock, Nothing, (ZePtr{Int},), a)
ccall(:clock, Nothing, (PtrOrZePtr{Int},), a)

b = oneArray{eltype(a), ndims(a)}(undef, size(a))
ccall(:clock, Nothing, (ZePtr{Int},), b)
@test_throws Exception ccall(:clock, Nothing, (Ptr{Int},), b)
ccall(:clock, Nothing, (PtrOrZePtr{Int},), b)

end

@testset "reference values" begin
    # Ref

    @test typeof(Base.cconvert(Ref{Int}, 1)) == typeof(Ref(1))
    @test Base.unsafe_convert(Ref{Int}, Base.cconvert(Ref{Int}, 1)) isa Ptr{Int}

    ptr = reinterpret(Ptr{Int}, C_NULL)
    @test Base.cconvert(Ref{Int}, ptr) == ptr
    @test Base.unsafe_convert(Ref{Int}, Base.cconvert(Ref{Int}, ptr)) == ptr

    arr = [1]
    @test Base.cconvert(Ref{Int}, arr) isa Base.RefArray{Int, typeof(arr)}
    @test Base.unsafe_convert(Ref{Int}, Base.cconvert(Ref{Int}, arr)) == pointer(arr)

    # ZeRef

    @test typeof(Base.cconvert(ZeRef{Int}, 1)) == typeof(ZeRef(1))
    @test Base.unsafe_convert(ZeRef{Int}, Base.cconvert(ZeRef{Int}, 1)) isa ZeRef{Int}

    zeptr = reinterpret(ZePtr{Int}, C_NULL)
    @test Base.cconvert(ZeRef{Int}, zeptr) == zeptr
    @test Base.unsafe_convert(ZeRef{Int}, Base.cconvert(ZeRef{Int}, zeptr)) == Base.bitcast(ZeRef{Int}, zeptr)

    zearr = oneAPI.oneArray([1])
    @test Base.cconvert(ZeRef{Int}, zearr) isa oneL0.ZeRefArray{Int, typeof(zearr)}
    @test Base.unsafe_convert(ZeRef{Int}, Base.cconvert(ZeRef{Int}, zearr)) == Base.bitcast(ZeRef{Int}, pointer(zearr))

    # RefOrZeRef

    @test typeof(Base.cconvert(RefOrZeRef{Int}, 1)) == typeof(Ref(1))
    @test Base.unsafe_convert(RefOrZeRef{Int}, Base.cconvert(RefOrZeRef{Int}, 1)) isa RefOrZeRef{Int}

    @test Base.cconvert(RefOrZeRef{Int}, ptr) == ptr
    @test Base.unsafe_convert(RefOrZeRef{Int}, Base.cconvert(RefOrZeRef{Int}, ptr)) == Base.bitcast(RefOrZeRef{Int}, ptr)

    @test Base.cconvert(RefOrZeRef{Int}, zeptr) == zeptr
    @test Base.unsafe_convert(RefOrZeRef{Int}, Base.cconvert(RefOrZeRef{Int}, zeptr)) == Base.bitcast(RefOrZeRef{Int}, zeptr)

    @test Base.cconvert(RefOrZeRef{Int}, arr) isa Base.RefArray{Int, typeof(arr)}
    @test Base.unsafe_convert(RefOrZeRef{Int}, Base.cconvert(RefOrZeRef{Int}, arr)) == Base.bitcast(RefOrZeRef{Int}, pointer(arr))

    @test Base.cconvert(RefOrZeRef{Int}, zearr) isa oneL0.ZeRefArray{Int, typeof(zearr)}
    @test Base.unsafe_convert(RefOrZeRef{Int}, Base.cconvert(RefOrZeRef{Int}, zearr)) == Base.bitcast(RefOrZeRef{Int}, pointer(zearr))
end


text:
using S3Dicts

as3  = S3Dict("s3://seunglab/jpwu/test/")

a = rand(UInt8, 50)

as3["test"] = a

b = as3["test"]

b = reinterpret(UInt8, b)

@assert all(a .== b)

info("delete the file in s3")
delete!(as3, "test")


text:
using ArrayLayouts, BlockBandedMatrices, BlockArrays, LinearAlgebra, Random, Test
import BlockBandedMatrices: colblockbandwidths, BroadcastStyle, BlockSkylineStyle, blockisequal

Random.seed!(0)

@testset "BlockSkylineMatrix" begin
    @testset "@jagot lmul! bug" begin
        rows = rand(1:10, 5)
        l = rand(-2:2, 5)
        u = rand(-2:2, 5)

        m = sum(rows)

        A = BlockSkylineMatrix(Zeros(m,m), rows,rows, (l,u))
        @test @inferred(size(A)) == (sum(rows),sum(rows))

        @test Array(A) == Matrix(A)

        A.data .= rand(size(A.data)...)

        V = zeros(m,2)
        V[:,1] .= rand(m)
        reference = A*V[:,1]

        @view(V[:,2]) .= MulAdd(A, @view(V[:,1]))
        @test V[:,2] ≈ reference

        V[:,2] .= NaN
        @view(V[:,2]) .= MulAdd(A, @view(V[:,1]))
        @test V[:,2] ≈ reference
    end

    @testset "BlockSkylineMatrix multiplication" begin
        rows = [3, 1, 2, 1, 2, 1, 2, 1, 2, 1, 3]
        l,u = [1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1], [1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1]

        M = BlockSkylineMatrix{Float64}(undef, rows,rows, (l,u))
        M.data .= 1

        d = Diagonal(1.0:size(M,2))
        D = BandedBlockBandedMatrix(d, rows,rows, (0,0), (0,0))

        MD = M*D
        @test MD isa BlockSkylineMatrix
        @test MD == Matrix(M)*d
        @test colblockbandwidths(MD) == (l,u)

        MM = M*M
        @test MM isa BlockSkylineMatrix
        @test MM == Matrix(M)^2
        # Ensure correct (minimal) bandedness of product
        MMl,MMu = colblockbandwidths(MM)
        @test MMl[1:7] == [3,4,3,4,3,4,3]
        @test all(MMl[8:10] .≥ [3,2,1])
        @test all(MMu[2:4] .≥ [1,2,3])
        @test MMu[5:11] == [3,4,3,4,3,4,3]

        N = BlockBandedMatrix{Float64}(undef, rows,rows, (1,1))
        N.data .= 1
        NN = N*N
        # We don't want a BlockBandedMatrix^2 to become a general
        # BlockSkylineMatrix
        @test NN isa BlockBandedMatrix
        @test NN == Matrix(N)^2

        rows = [9, 4, 1, 10, 6]
        O = BlockSkylineMatrix{Int64}(undef, rows,rows, ([-2, 2, 0, 2, -1],[-1, 2, 1, 0, -1]))
        O.data .= 1
        OO = O*O
        @test OO isa BlockSkylineMatrix
        @test OO == Matrix(O)^2
    end

    @testset "Broadcasting" begin
        rows = [3, 1, 2, 1, 2, 1, 2, 1, 2, 1, 3]
        l,u = [1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1], [1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1]

        M = BlockSkylineMatrix{Float64}(undef, rows,rows, (l,u))
        M.data .= 1

        @test BroadcastStyle(BlockSkylineMatrix) == BlockSkylineStyle()
        @test blocksize(M + M) == blocksize(M)

        # A is a BlockSkylineMatrix with varying bandwidths
        rows = [3,4,5]
        A = BlockSkylineMatrix{Float64}(undef, rows,rows, ([1,0,1],[0,1,0]))
        A.data .= 1

        # B is diagonal
        B = BandedBlockBandedMatrix{Float64}(undef, rows, rows, (0,0), (0,0))
        B.data .= -1

        # AB should have the same structure as A
        AB = A+B
        @test Matrix(AB) == Matrix(A) + Matrix(B)
        @test blockisequal(axes(AB), axes(A))

        # C has larger bandwidths than A
        C = BandedBlockBandedMatrix{Float64}(undef, rows, rows, (1,1), (1,1))
        C.data .= -1

        # Thus, AC should be a BlockSkylineMatrix with C's bandwidths.
        AC = A+C
        @test Matrix(AC) == Matrix(A) + Matrix(C)
        @test blockisequal(axes(AC),axes(C))
    end

    @testset "BlockSkylineMatrix arithmetic" begin
        t = Int
        rows = cols = 1:4
        N = sum(rows)
        D = Diagonal(ones(t, N))
        Bu = Bidiagonal(ones(t, N), 2ones(t, N-1), 'U')
        Bl = Bidiagonal(ones(t, N), 2ones(t, N-1), 'L')
        T = Tridiagonal(ones(t, N-1), 2ones(t, N), 3ones(t, N-1))
        ST = SymTridiagonal(ones(t, N), 2ones(t, N-1))

        function check_blockskylinematrix_arithmetic(A, op, B, (l,u))
            C = op(A, B)
            Cdense = op((A isa BlockSkylineMatrix ? Matrix(A) : A),
                        (B isa BlockSkylineMatrix ? Matrix(B) : B))
            @test C == Cdense

            bs = (A isa BlockSkylineMatrix ? A : B).block_sizes
            expl = max.(l, bs.l)
            expu = max.(u, bs.u)

            Cbs = C.block_sizes
            @test all(Cbs.l .== expl)
            @test all(Cbs.u .== expu)
        end

        for (l,u) = [([0,0,0,0],[0,0,0,0]),
                     ([-1,-1,-1,-1],[1,2,2,2]),
                     ([-1,-2,-2,-2],[1,1,2,2]),
                     ([2,2,1,1],[-1,-1,-1,-1]),
                     ([2,2,1,1],[-2,-2,-2,-1])]
            M = BlockSkylineMatrix{t}(undef, rows, cols, (l,u))
            M.data .= 1
            for (A,expected_bws) = [(3I, ([0,0,0,0],[0,0,0,0])),
                                    (D, ([0,0,0,0],[0,0,0,0])),
                                    (Bu, ([0,0,0,0],[0,1,1,1])),
                                    (Bl, ([1,1,1,0],[0,0,0,0])),
                                    (T, ([1,1,1,0],[0,1,1,1])),
                                    (ST, ([1,1,1,0],[0,1,1,1]))]
                for op = (+, -)
                    check_blockskylinematrix_arithmetic(M, op, A, expected_bws)
                    check_blockskylinematrix_arithmetic(A, op, M, expected_bws)
                end
            end
        end
    end
end


text:
# Functions used to initiate ycreator type, for adjusting phenotypes with gensel outputs
function build_ycreator()
    return YCreator(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
end
#
function add_markers_gensel(ycreator, markersfileprefix)
    info = read_bim_file(markersfileprefix .* ".bim")
    iid = read_fam_file(markersfileprefix .* ".fam")
    bed = read_bed_file(markersfileprefix .* ".bed")
    genselvalues = (bedvals[bed] .- 1).*-1  # convert to -1, 0, 1 # *-1 as plinks ref/alt coding is differnt for bed files vs raw files. gensel alphas are calculated wrt to raw files.
    ycreator.genselmarkers = GenselMarkers(iid, info, genselvalues)
end
#
function add_alphas(ycreator, alphasfilename)
    read_filtered_alphas(ycreator, alphasfilename)
end
# rsplit acts in reverse, and limit means it just splits on the 1st occurence. (This prevents removal of missingIDs (CHRN_POS) format)
function read_filtered_alphas(ycreator, alphasfilename)
    alphas = readdlm(alphasfilename, '\t', Any, '\n'; header = false)
    ycreator.alphas = Alphas(alphas[:,1], alphas[:,2:end], 0 ,0)
    ycreator.alphas.m = size(ycreator.alphas.markerNames, 1)
    ycreator.alphas.k = size(ycreator.alphas.alphas, 2)
    ycreator.alphas.markerNames = map(x -> String(x[1]), rsplit.(ycreator.alphas.markerNames, '_'; limit = 2))
end

function add_phenotypes(ycreator, phenotypesfilename)
    ycreator.phenotypes = read_y(phenotypesfilename)

end

function read_y(phenotypesfilename)
    y = readdlm(phenotypesfilename, ' ', Float64, '\n'; header = true)
    return Phenotypes(map(x->string(Int64(x)), y[1][:,1]), y[1][:,2:2])
end
# marker and alpha orders should be the same
function sort_ycreator(ycreator)
    markerOrderM = collect(1:ycreator.alphas.m)
    markerOrderN, phenoOrderN = get_orders(ycreator.genselmarkers.IID, ycreator.phenotypes.IID)
    get_genselmarker_intersect!(ycreator.genselmarkers, markerOrderN, markerOrderM)
    get_phenotype_intersect!(ycreator.phenotypes, phenoOrderN)
end

function get_genselmarker_intersect!(genselmarkers, new_orderN, new_orderM)
    genselmarkers.IID = genselmarkers.IID[new_orderN]
    genselmarkers.SNPID = genselmarkers.SNPID[new_orderM, :]
    genselmarkers.markers = genselmarkers.markers[new_orderN, new_orderM]
end

function get_alpha_intersect!(alphas, new_order)
    alphas.markerNames = alphas.markerNames[new_order]
    alphas.alphas = alphas.alphas[new_order, :]
    alphas.m = size(alphas.markerNames, 1)
end

function initiate_ycreator(ycreator, interval, chr, outputprefix, region, code)
    ycreator.n = size(ycreator.phenotypes.IID, 1)
    ycreator.m = size(ycreator.genselmarkers.SNPID, 1)
    ycreator.k = size(ycreator.alphas.k)
    ycreator.outputprefix = outputprefix
    ycreator.interval = interval
    ycreator.region = region
    ycreator.chr = chr
    ycreator.code = code
    if ycreator.code == "LOCO"
        ycreator.chrlength = 0
    else
        ycreator.chrlength = get_chr_length(ycreator.genselmarkers.SNPID, ycreator.chr)
    end
    ycreator.nlayer1 = ceil(ycreator.chrlength / ycreator.interval)
    ycreator.nlayer2 = ycreator.nlayer1 - 1
end

function get_chr_length(markerinfo, chr)
    chrindex = findall(x -> x == chr, markerinfo[:,1])
    chrview = view(markerinfo, chrindex, 2)
    chrviewfloat = map(x -> parse(Float64, x), chrview)
    chrlength = maximum(chrviewfloat)
    return chrlength
end

function fill_missing_genotypes!(ycreator)
    if any(ismissing.(ycreator.genselmarkers.markers))
        ycreator.genselmarkers.markers = Array{Union{Missing, Float64},2}(ycreator.genselmarkers.markers)
        for i in 1:ycreator.m
            l = length(ycreator.genselmarkers.markers[ismissing.(ycreator.genselmarkers.markers[:,i]),i])
            if l != 0
                ycreator.genselmarkers.markers[ismissing.(ycreator.genselmarkers.markers[:,i]),i] = fill(mean(skipmissing(ycreator.genselmarkers.markers[:,i])),l)
            end
        end
    end
end

max_stars_repo_name: egmreynolds/MANA

max_issues_repo_name: egmreynolds/MANA

max_forks_repo_name: egmreynolds/MANA