text:
module Utils

using DataFrames
using DataStructures: OrderedDict
using Dates

using CodecZlib
using Serialization
using Printf

export dview, toposort, coldim, ncols, FlushedLogger
export serialize_gzip, deserialize_gzip

coldim(::AbstractDataFrame) = 2
coldim(::AbstractMatrix) = 1

"""
    ncols(d)

Return the number of columns, i.e., features of `d`.
"""
ncols(d) = size(d, coldim(d))

function dview(d::AbstractDataFrame, obs, cols)
    view(d, obs, cols)
end

dview(d::AbstractMatrix, obs, cols) = view(d, cols, obs)

"""
    isapprox(dfa::AbstractDataFrame, dfb::AbstractDataFrame)

Check if two `DataFrame` instances are approximatelly equal.
"""
function Base.isapprox(dfa::AbstractDataFrame, dfb::AbstractDataFrame)
    function isapproxifexists(a, b)
        try
            return a ≈ b
        catch e
            e isa MethodError && return a == b
            rethrow(e)
        end
    end
    names(dfa) == names(dfb) && eltypes(dfa) == eltypes(dfb) && size(dfa) == size(dfb) &&
        all(isapproxifexists(dfa[!,n], dfb[!,n]) for n in names(dfa))
end

"""
    toposort(nodes, predf, nodef = identity)

Sort graph nodes topologically.

Provide `node` an iterator over node objects, `predf` a function giving predecessor of a node and
a `nodef` function optionally applied to the node in the result. [See](https://en.wikipedia.org/wiki/Topological_sorting#Depth-first_search).

# Examples
```@meta
DocTestSetup = quote
    using DataStructures: OrderedDict
    using MLGoodies: toposort
end
```

```jldoctest
julia> nodes = OrderedDict(2 => [1], 5 => [1, 3, 4], 4 => [1, 2], 6 => [5], 1 => [], 3 => [1], 20 => [5, 10], 10 => [1]);

julia> toposort(nodes, v -> last(v), v -> first(v))
8-element Array{Any,1}:
  1
  2
  3
  4
  5
  6
 10
 20
```
"""
function toposort(nodes, predf, nodef = identity)
    ndict = OrderedDict(nodef(n) => predf(n) for n in nodes)

    tmarks = Set()
    pmarks = Set()
    l = []
    function visit(n)
        n in pmarks && return
        n in tmarks && error("not a DAG!")
        push!(tmarks, n)
        for p in ndict[n]
            visit(p)
        end
        delete!(tmarks, n)
        push!(l, n)
        push!(pmarks, n)
    end

    while length(pmarks) < length(nodes)
        for n in keys(ndict)
            n ∉ pmarks && visit(n)
        end
    end
    l
end

#-------------------------------------------------------------------------------
# Serialization
#-------------------------------------------------------------------------------
function serialize_gzip(fname, data)
    open(fname, "w") do f
        gf = GzipCompressorStream(f)
        Serialization.serialize(gf, data)
        close(gf)
    end
end

function deserialize_gzip(fname)
    open(fname) do f
        gf = GzipDecompressorStream(f)
        data = Serialization.deserialize(gf)
        close(gf)
        data
    end
end

end

max_stars_repo_name: dhonza/MLGoodies


text:
module params

using Pkg.TOML
using DSP

export MicOptics, SimulationParams, Constants, ParameterSet, Space, MicroLensArray, LightFieldSimulation

struct MicOptics

    ftl::Float64
    fobj::Float64
    M::Float64
    NA::Float64
    a0::Float64
    b0::Float64
    lambda::Float64
    n::Float64
    d::Float64

    function MicOptics(optics::Dict{String,Any}, fml::Float64)

        ftl     = pop!(optics, "ftl")
        M       = pop!(optics, "M")
        NA      = pop!(optics, "NA")
        a0      = pop!(optics, "a0")
        b0      = pop!(optics, "b0")
        lambda  = pop!(optics, "lambda")
        n       = pop!(optics, "n")
        fobj    = ftl / M

        if b0 > fml
            d = b0
        else
            d = fml
        end

        new(ftl, fobj, M, NA, a0, b0, lambda, n, d)
    end
end

struct SimulationParams

    vpix::Int64
    osr::Int64
    zmax::Float64
    zmin::Float64
    zstep::Float64
    subvpix::Int64
    pixelpitch::Float64
    subpixelpitch::Float64

    function SimulationParams(simparams::Dict{String,Any}, pitch::Float64)

        vpix    = pop!(simparams, "vpix")
        osr     = pop!(simparams, "osr")
        zmax    = pop!(simparams, "zmax")
        zmin    = pop!(simparams, "zmin")
        zstep   = pop!(simparams, "zstep")

        if mod(vpix, 2) == 0
            error("Invalid number of virtual pixels--must be odd number.")
        end

        subvpix = vpix * osr
        pixelpitch = pitch / vpix
        subpixelpitch = pixelpitch / osr
        new(vpix, osr, zmax, zmin, zstep, subvpix, pixelpitch, subpixelpitch)

    end
end

struct Constants

    k::Float64
    k0::Float64
    alpha::Float64

    function Constants(opt::MicOptics)

        k = 2 * pi * opt.n / opt.lambda
        k0 = 2 * pi * 1 / opt.lambda
        alpha = asin(opt.NA / opt.n)

        new(k, k0, alpha)
    end
end

struct Space
    x::Vector{Float64}
    y::Vector{Float64}
    z::Vector{Float64}
    xlen::Int64
    ylen::Int64
    zlen::Int64
    center::Int64
    function Space(x = [0.0], y = [0.0], z = [0.0])
        xlen = length(x)
        ylen = length(y)
        zlen = length(z)
        center = cld(xlen, 2)
        new(x, y, z, xlen, ylen, zlen, center)
    end
end
struct ParameterSet
    opt::MicOptics
    sim::SimulationParams
    con::Constants
end

struct MicroLensArray

    pitch::Float64
    fml::Float64
    array::Array{ComplexF64,2}
    space::Space

    function MicroLensArray(pitch::Float64, fml::Float64, mlaspace::Space, img::Space, par::ParameterSet)
        #center = findfirst(img.x .== 0)
        mllen = mlaspace.xlen

        allcenters = vcat(img.center:-mllen:1, (img.center + mllen):mllen:img.xlen)
        sort!(allcenters)

        a = zeros(ComplexF64, mllen, mllen)

        l2norm = mlaspace.x.^2 .+ mlaspace.y'.^2
        a .= exp.(-im * par.con.k / (2 * fml) .* l2norm);

        b = zeros(ComplexF64, img.xlen, img.ylen)
        b[allcenters, allcenters] .= 1

        #TODO: This is messy. Find a cleaner way to do it? works for now...
        c = conv(b, a)
        border = fld((size(c, 1) - img.xlen), 2) + 1
        crop = Int.(border:(img.xlen + (border-1)))

        new(pitch, fml, c[crop, crop], mlaspace)
    end
end

struct LightFieldSimulation
    par::ParameterSet
    mla::MicroLensArray
    img::Space
    obj::Space
end

end # module

ext: jl


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

using Markdown
using InteractiveUtils

# ╔═╡ 3e0c149a-3350-11eb-0b17-7bea5d1a28d4
begin
	DAY = 5
	md"# 2020 Day $DAY"
end

# ╔═╡ 18080fe2-339b-11eb-34f0-f393638506fc
md"## Part 1"

# ╔═╡ 4e084a4e-339b-11eb-3ffc-41118215f226
md"*<Insert Description>*"

# ╔═╡ 24e1613c-339b-11eb-2976-59884053d1ee
input = split(read("../data/day_$DAY.txt", String))

# ╔═╡ a9da2b56-36c7-11eb-38c5-51bf188463f9
function partition(code::AbstractString,lower::Int,upper::Int)
	# Final iteration
	if (upper - lower == 1) & length(code) == 1
		if code[1] ∈ "FL"
			return lower
		else
			return upper
		end
	# Split and go recursive
	elseif (upper - lower > 2) & (length(code) > 1)
		middle = (upper - lower) ÷ 2 + lower
		if code[1] ∈ "FL"
			return partition(code[2:end],lower,middle)
		else
			return partition(code[2:end],middle+1,upper)
		end
	else
		throw(DomainError((code, lower, upper), "this input is bad"))
	end
end

# ╔═╡ 75423650-36c8-11eb-1b63-b3a86f084690
function get_row(boarding_pass::AbstractString)
	code = boarding_pass[1:7]
	return partition(code,0,127)
end

# ╔═╡ 14af022e-36ca-11eb-3a6c-737120328af3
@assert get_row("BFFFBBFRRR") == 70

# ╔═╡ 3c6af638-36ca-11eb-3b4a-a1eb6dbbfdfc
@assert get_row("FFFBBBFRRR") == 14

# ╔═╡ 3c6b5790-36ca-11eb-325f-5b39977ccab6
@assert get_row("BBFFBBFRLL") == 102

# ╔═╡ 8bd3528a-36c8-11eb-083b-4f795b8537de
function get_col(boarding_pass::AbstractString)
	code = boarding_pass[8:end]
	return partition(code,0,7)
end

# ╔═╡ 8b722b22-36c8-11eb-0870-81aa0ca325e2
@assert get_col("BFFFBBFRRR") == 7

# ╔═╡ a017e524-36ca-11eb-3ca8-099185680195
@assert get_col("FFFBBBFRRR") == 7

# ╔═╡ a018306a-36ca-11eb-0b9e-dd5c58259bc4
@assert get_col("BBFFBBFRLL") == 4

# ╔═╡ 1c865778-36c9-11eb-2a6d-63710b4917a1
function get_id(boarding_pass::AbstractString)
	col = get_col(boarding_pass)
	row = get_row(boarding_pass)
	return row*8+col
end

# ╔═╡ cc9399cc-36ca-11eb-338e-15e78e041ad6
@assert get_id("BFFFBBFRRR") == 567

# ╔═╡ 9f987d82-36d6-11eb-04e5-d749168db73c
maximum(get_id.(input))

# ╔═╡ 9136a95e-36d5-11eb-3ab2-719c2e16c155
# alternative solution
begin
	# make it into 01110111
	pass2binary(pass::AbstractString) = join((c ∈ "BR")*1 for c in pass)
	# convert that into a base 10 representation
	pass2id(pass::AbstractString) = parse(Int, pass2binary(pass); base=2)
	maximum(pass2id.(input))
end

# ╔═╡ 09abf424-36d7-11eb-1a53-67c78def2720
maximum(pass2id.(input))

# ╔═╡ d503722e-36d6-11eb-1634-1178998b10a8
pass2seat(pass) = parse(Int, map(c -> c ∈ ('R', 'B') ? '1' : '0', pass), base = 2)

# ╔═╡ fd392980-36d6-11eb-12f8-b7e499a15059
maximum(pass2seat.(input))

# ╔═╡ 96c8893e-36df-11eb-14bf-4dcc8dd0b9a7
(
	split(read("../data/day_$DAY.txt", String)) |>
	strings -> map.(c -> c ∈ ('R', 'B') ? '1' : '0', strings) |>
	binary_strings -> parse.(Int, binary_strings; base=2) |>
	maximum
)

# ╔═╡ 2088463e-3709-11eb-0f4f-dd0ffe4ca9ce
begin
	boarding_passes = split(read("../data/day_$DAY.txt", String))
	binary_strings = map.(c -> c ∈ ('R', 'B') ? '1' : '0', boarding_passes)
	seat_ids = parse.(Int, binary_strings; base=2)
	maximum(seat_ids)
end

# ╔═╡ 0c0945e2-36e0-11eb-2953-eda4795bc3b7
maximum(pipe_seatid.(input))

# ╔═╡ 249718d2-36e0-11eb-0258-d159b87e4641
"110010010" |> s -> parse(Int, s ; base=2)

# ╔═╡ 4e2a74a0-36f9-11eb-08b2-5f7f6263d020
"110010010" |> parse(Int, _ ; base=2)

# ╔═╡ 5f797e22-339d-11eb-1ff3-dda90f946838
md"## Part 2"

# ╔═╡ e3177eae-339e-11eb-3f9f-1f7230976832
md"*Find your seat, it's not in the list but -1 and +1 will be*"

# ╔═╡ 3efb171a-36cb-11eb-1005-af68bffb30aa
function part2(input)
	ids = sort(get_id.(input))
	for i in 2:length(ids)
		if ids[i] - ids[i-1] == 2
			return ids[i]-1
		end
	end
end

# ╔═╡ c7a03d02-36cb-11eb-31c3-df3892a9fdb4
part2(input)

# ╔═╡ a81ade6a-36cb-11eb-0935-d759440ebcfb
typof(input)

# ╔═╡ a2141e96-36cb-11eb-1d9d-89c94dc1a7fb
typeof(input)

# ╔═╡ 5c846a16-3707-11eb-0339-919634714fdc

# ╔═╡ Cell order:
# ╠═3e0c149a-3350-11eb-0b17-7bea5d1a28d4
# ╟─18080fe2-339b-11eb-34f0-f393638506fc
# ╠═4e084a4e-339b-11eb-3ffc-41118215f226
# ╠═24e1613c-339b-11eb-2976-59884053d1ee
# ╠═a9da2b56-36c7-11eb-38c5-51bf188463f9
# ╠═75423650-36c8-11eb-1b63-b3a86f084690
# ╠═14af022e-36ca-11eb-3a6c-737120328af3
# ╠═3c6af638-36ca-11eb-3b4a-a1eb6dbbfdfc
# ╠═3c6b5790-36ca-11eb-325f-5b39977ccab6
# ╠═8bd3528a-36c8-11eb-083b-4f795b8537de
# ╠═8b722b22-36c8-11eb-0870-81aa0ca325e2
# ╠═a017e524-36ca-11eb-3ca8-099185680195
# ╠═a018306a-36ca-11eb-0b9e-dd5c58259bc4
# ╠═1c865778-36c9-11eb-2a6d-63710b4917a1
# ╠═cc9399cc-36ca-11eb-338e-15e78e041ad6
# ╠═9f987d82-36d6-11eb-04e5-d749168db73c
# ╠═9136a95e-36d5-11eb-3ab2-719c2e16c155
# ╠═09abf424-36d7-11eb-1a53-67c78def2720
# ╠═d503722e-36d6-11eb-1634-1178998b10a8
# ╠═fd392980-36d6-11eb-12f8-b7e499a15059
# ╠═96c8893e-36df-11eb-14bf-4dcc8dd0b9a7
# ╠═2088463e-3709-11eb-0f4f-dd0ffe4ca9ce
# ╠═0c0945e2-36e0-11eb-2953-eda4795bc3b7
# ╠═249718d2-36e0-11eb-0258-d159b87e4641
# ╠═4e2a74a0-36f9-11eb-08b2-5f7f6263d020
# ╠═5f797e22-339d-11eb-1ff3-dda90f946838
# ╠═e3177eae-339e-11eb-3f9f-1f7230976832
# ╠═3efb171a-36cb-11eb-1005-af68bffb30aa
# ╠═c7a03d02-36cb-11eb-31c3-df3892a9fdb4
# ╠═a81ade6a-36cb-11eb-0935-d759440ebcfb
# ╠═a2141e96-36cb-11eb-1d9d-89c94dc1a7fb
# ╠═5c846a16-3707-11eb-0339-919634714fdc

max_stars_repo_name: RensDimmendaal/advent-of-code-julia

max_issues_repo_name: RensDimmendaal/advent-of-code-julia

max_forks_repo_name: RensDimmendaal/advent-of-code-julia


text:
abstract type AbstractRegularization end
abstract type AbstractRegularizer end
# Regularization API #
# ------------------------------------------------------------
initialize_regularization!(lshaped::AbstractLShapedSolver) = initialize_regularization!(lshaped, lshaped.regularization)
log_regularization!(lshaped::AbstractLShapedSolver) = log_regularization!(lshaped, lshaped.regularization)
log_regularization!(lshaped::AbstractLShapedSolver, t::Integer) = log_regularization!(lshaped, t, lshaped.regularization)
take_step!(lshaped::AbstractLShapedSolver) = take_step!(lshaped, lshaped.regularization)
decision(lshaped::AbstractLShapedSolver) = decision(lshaped, lshaped.regularization)
objective(lshaped::AbstractLShapedSolver) = objective(lshaped, lshaped.regularization)
solve_regularized_master!(lshaped::AbstractLShapedSolver) = solve_regularized_master!(lshaped, lshaped.regularization)
gap(lshaped::AbstractLShapedSolver) = gap(lshaped, lshaped.regularization)
process_cut!(lshaped::AbstractLShapedSolver, cut::AbstractHyperPlane) = process_cut!(lshaped, cut, lshaped.regularization)
project!(lshaped::AbstractLShapedSolver) = project!(lshaped, lshaped.regularization)
# ------------------------------------------------------------
include("common.jl")
include("no_regularization.jl")
include("rd.jl")
include("tr.jl")
include("lv.jl")


text:
using jlmie
using Plots
pyplot()

# settings
# structure
nmat = 4  # refractive index of a sphere
radius = 100*1e-9

# environment
nenv = 1.00

# wavelenth
# lbd0 = 618*1e-9  # ED
# lbd0 = 831*1e-9  # MD
lbd0 = 918*1e-9  # 1nd Kerker
# lbd0 = 756*1e-9  # 2nd Kerker
lbdp = lbd0.*1e9

# angular range
theta = range(0,π,length=73) # per 2.5 degrees with length 73
phi = range(0,2π,length=145) # per 2.5 degrees with length 145

# other settings
nmax = -1  # -1: namx large enough (determined from x)

Isff = zeros(length(theta),length(phi))
for i = 1:length(theta)
    for j = 1:length(phi)
        Isff[i,j],_,_ = jlmie_Isff(nmat,radius,lbd0,nenv,theta[i],phi[j],nmax)
    end
end
Isff = Isff/maximum(Isff)
mat_x = Isff.* (sin.(theta) * cos.(phi)')
mat_y = Isff.* (sin.(theta) * sin.(phi)')
mat_z = Isff.* (cos.(theta) * cos.(phi*0)')

# show results
plt = plot(mat_x,mat_y,mat_z,
        fill_z=abs.(Isff),
        st=:surface,
        camera=(45,30), # azimuth, elevate
        xlabel="x",
        ylabel="y",
        zlabel="z",
        # colorbar_title="Scattering intensity (arb. units)",
        xlims=(-1,1),
        ylims=(-1,1),
        zlims=(-1,1),
)
display(plt)

# save
if true
    runningfilename = splitext(splitpath(@__FILE__)[end])[1]
    outfilename = ".\\result\\" * runningfilename
    savefig(plt,outfilename)
    println("Graph has been saved in result directory")
end

max_stars_repo_name: Hinamoooon/jlmie

max_issues_repo_name: Hinamoooon/jlmie

max_forks_repo_name: Hinamoooon/jlmie


text:
module Rater

export rate

rate(book) = "nice!"

end # module


text:
import Base: show

function show(io::IO, x::RandomStateType)
    print(io, "RandomStateType")
end

function show(io::IO, x::OpaqueType)
    print(io, "OpaqueType")
end

function show(io::IO, x::ArrayType)
    et = x.eltype
    sh = x.shape

    if isnull(et)
        print(io, "unknown(")
    else
        print(io, "$(get(et))(")
    end

    if isnull(sh)
        print(io, "unknown")
    else
        labels = map(dim -> isnull(dim) ? "unknown" : string(get(dim)), get(sh))
        join(io, labels, ", ")
    end

    print(io, ")")
end

function show(io::IO, x::Unknown)
    print(io, "U[$(x.name), $(x.typ)]")
end

function show(io::IO, x::LazyArray)
    typ = array_typeof(x.value)
    print(io, "KA[$typ]")
end

function show(io::IO, aexpr::CallNode)
    print(io, "(")
    show(io, aexpr.func)
    print(io, ", ")

    args = aexpr.args
    for (i, arg) in enumerate(args)
        show(io, arg)
        if i != length(args)
            print(io, ", ")
        end
    end
    print(io, ")")
end


text:
struct Octree{T<:AbstractFloat}
    X::Array{T}
    root::Cell        # Parent cell
    maxdepth::Int
end

struct LocationalCode
    x::String
    y::String
    z::String
end
LocationalCode() = LocationalCode("", "", "")

struct LeafData
    ids::Array{Int}
    code::LocationalCode
    depth::Int
end

LeafData(; ids, code, depth) = LeafData(ids, code, depth)


text:
bitsof(::Type{T}) where {T} = sizeof(T) * 8
bitsof(x::T) where {T} = sizeof(T) * 8

#=
    maxmag(q) whichever is the larger, abs(numerator) or abs(denominator)
    leading_zeros_maxmag(q) count of leading 0bits in maxmag(q)

    msbitidx(q) 1-based bit index of the most significant bit in maxmag(q)
                zero iff iszero(maxmag(q))
=#
maxmag(q::Rational{T}) where {T<:Integer} = max(abs(q.num), abs(q.den))
leading_zeros_maxmag(q::Rational{T}) where {T<:Integer} = leading_zeros(maxmag(q))

msbitidx(q::Rational{T}) where {T<:Integer} = bitsof(T) - leading_zeros_maxmag(q)
msbitidx(q1::Rational{T}, q2::Rational{T}) where {T<:Integer} = msbitidx(q1) + msbitidx(q2)

maxmag(q::FastRational{T}) where {T<:FastInt} = max(abs(q.num), abs(q.den))  # q.den != typemin(T)
leading_zeros_maxmag(q::FastRational{T}) where {T<:FastInt} = leading_zeros(maxmag(q))

msbitidx(q::FastRational{T}) where {T<:FastInt} = bitsof(T) - leading_zeros(maxmag(q))
msbitidx(q1::FastRational{T}, q2::FastRational{T}) where {T<:FastInt} = msbitidx(q1) + msbitidx(q2)

#=
msbitidx(q1::T) + msbitidx(q2::T)
    =  (bitsof(T) - leading_zeros(maxmag(q1))) + (bitsof(T) - leading_zeros(maxmag(q2)))
    =  (bitsof(T) + bitsof(T)) - (leading_zeros(maxmag(q1)) + leading_zeros(maxmag(q2)))
    =  2*bitsof(T) - (leading_zeros(maxmag(q1)) + leading_zeros(maxmag(q2)))
=#

mayoverflow(q1::T, q2::T) where {T} = bitsof(T) <= msbitidx(q1, q2)
#=
   = bitsof(T) <= 2*bitsof(T) - (leading_zeros(maxmag(q1)) + leading_zeros(maxmag(q2)))
   = (leading_zeros(maxmag(q1)) + leading_zeros(maxmag(q2))) <= 2*bitsof(T) - bitsof(T)
   = (leading_zeros(maxmag(q1)) + leading_zeros(maxmag(q2))) <= bitsof(T)
=#

mayoverflow(i1::T, i2::T) where {T<:Integer} =
    (leading_zeros(i1) + leading_zeros(i2)) <= bitsof(T)

#=
   this was used in deriving the faster version immediately below

   mayoverflow(q1::Rational{T}, q2::Rational{T}) where {T<:Integer} =
        bitsof(T) >= leading_zeros_maxmag(q1) + leading_zeros_maxmag(q2)
=#

mayoverflow(q1::Rational{T}, q2::Rational{T}) where {T<:Integer} =
    (bitsof(T)<<1) >= leading_zeros(q1.num) + leading_zeros(q1.den) +
                      leading_zeros(q2.num) + leading_zeros(q2.den)

mayoverflow(q1::FastRational{T}, q2::FastRational{T}) where {T<:FastInt} =
    (bitsof(T)<<1) >= leading_zeros(q1.num) + leading_zeros(q1.den) +
                      leading_zeros(q2.num) + leading_zeros(q2.den)


text:
"""
    count_nucleotides(strand)

The frequency of each nucleotide within `strand` as a dictionary.

Invalid strands raise a `DomainError`.

"""
function count_nucleotides(strand)
    counts = Dict('A' => 0, 'C' => 0, 'G' => 0, 'T' => 0)
    for nuc in strand
        if nuc == 'A' || nuc == 'C' || nuc == 'G' || nuc == 'T'
            counts[nuc] += 1
        else
            throw(DomainError(nuc, "Invalid nucleotide"))
        end
    end
    return counts
end

max_stars_repo_name: AbstractGeek/exercism-expeditions

max_issues_repo_name: AbstractGeek/exercism-expeditions

max_forks_repo_name: AbstractGeek/exercism-expeditions


text:
function calcCoilMaps(data::AbstractArray{Complex{T},3}, trj::AbstractVector{<:AbstractMatrix{T}}, U::AbstractMatrix{Complex{T}}, img_shape::NTuple{N,Int}; kernel_size = ntuple(_->6, N), calib_size =  ntuple(_->24, N), eigThresh_1=0.04, eigThresh_2=0.0, nmaps=1) where {N,T}
    Ncoils = size(data,3)
    Ndims = length(img_shape)
    imdims = ntuple(i->i, Ndims)

    # dataU = data .* U[:,1]'
    # dataU .*= dropdims(sum(abs2, combinedimsview(trj), dims=1), dims=1)
    dataU = similar(data) # size = Ncycles*Nr x Nt x Ncoils
    cU1 = conj(U[:,1])
    @batch for i ∈ CartesianIndices(dataU)
        dataU[i] = data[i] * cU1[i[2]] * sum(abs2, @view trj[i[2]][:,i[1]])
    end
    dataU = reshape(dataU, :, size(dataU, 3))

    p = plan_nfft(reduce(hcat,trj), img_shape)
    xbp = Array{Complex{T}}(undef, img_shape..., Ncoils)

    img_idx = CartesianIndices(img_shape)

    @info "BP for coils maps: "
    @time begin
        @batch for ic ∈ 1:Ncoils
            @views mul!(xbp[img_idx,ic], adjoint(copy(p)), dataU[:,ic])
        end
    end

    kbp = fftshift(xbp, imdims)
    fft!(kbp, imdims)
    kbp = fftshift(kbp, imdims)

    m = CartesianIndices(calib_size) .+ CartesianIndex((img_shape .- calib_size) .÷ 2)
    kbp = kbp[m,:]

    @info "espirit: "
    cmaps = @time MRIReco.espirit(kbp, img_shape, kernel_size, eigThresh_1=eigThresh_1, eigThresh_2=eigThresh_2, nmaps=nmaps)

    cmaps = [cmaps[img_idx,ic,1] for ic=1:Ncoils]
    xbp   = [  xbp[img_idx,ic  ] for ic=1:Ncoils]
    return cmaps, xbp
end


text:
using Test

using Random
using SparseArrays

# Ignore a statement
macro DISABLED(expr)
    return quote end
end
# Don't ignore a statement
# macro DISABLED(expr)
#     return expr
# end

const Dmax = 5

# Set reproducible random number seed
Random.seed!(0)

# Random rationals
function Base.rand(rng::AbstractRNG,
                   ::Random.SamplerType{Rational{T}}) where {T}
    return Rational{T}(T(rand(rng, -1000:1000)) // 1000)
end

# # Improved debug output for sparse matrices
# Base.show(io::IO, A::SparseMatrixCSC) = show_sparse(io, A)
# Base.show(io::IOContext, A::SparseMatrixCSC) = show_sparse(io, A)

include("test-defs.jl")
include("test-algorithms.jl")
include("test-zeroorone.jl")
include("test-sparseops.jl")
include("test-meshing.jl")
include("test-manifolds.jl")
include("test-funs.jl")
include("test-ops.jl")
include("test-manifoldops.jl")
include("test-continuum.jl")


text:
using StatsBase

abstract type Model end

layers(m::T) where {T<:Model} = error("layers(::$T) not implemented.")

inputs(m::T) where {T<:Model} = error("inputs(::$T) not implemented.")

outputs(m::T) where {T<:Model} = error("outputs(::$T) not implemented.")

Base.getindex(m::T, i::Int) where {T<:Model} = error("getindex(::$T, ::Int) not implemented")

function StatsBase.fit!(m::T, args...; kwargs...) where {T<:Model}
    error("fit(::$T, args...; kwargs...) not implemented.")
end

function evaluate(m::T, args...; kwargs...) where {T<:Model}
    error("evaluate(::$T, args...; kwargs...) not implemented.")
end

function StatsBase.predict(m::T, args...; kwargs...) where {T<:Model}
    error("predict(::$T, args...; kwargs...) not implemented.")
end

#= TODO:
- train_on_batch
- predict_on_batch
- fit_generator
- evaluate_generator
- predict_generator
- get_layer
=#

@doc PyDoc(Keras._models, :Sequential) struct Sequential <: Keras.Model
    obj::PyObject
    layers::Array{Keras.Layer}

    function Sequential()
        new(Keras._models[:Sequential](), Keras.Layer[])
    end

    function Sequential(layers::Array{Keras.Layer})
        new(
            Keras._models[:Sequential]([PyObject(l) for l in layers]),
            layers
        )
    end
end

Keras.layers(m::Sequential) = m.layers

Keras.inputs(m::Sequential) = [Tensor(i) for i in m.obj[:inputs]]

Keras.outputs(m::Sequential) = [Tensor(o) for o in m.obj[:outputs]]

Base.getindex(m::Sequential, i::Int) = m.layers[i]

function add!(m::Sequential, l::Keras.Layer)
    push!(m.layers, l)
    m.obj[:add](PyObject(l))
end

function compile!(m::Sequential, args...; kwargs...)
    m.obj[:compile](args...; kwargs...)
end

function StatsBase.fit!(m::Sequential, args...; verbose=0, kwargs...)
    m.obj[:fit](args...; verbose=verbose, kwargs...)
end

function evaluate(m::Sequential, X, y; batch_size=32, verbose=0)
    return m.obj[:evaluate](X, y, batch_size=batch_size, verbose=verbose)
end

function StatsBase.predict(m::Sequential, X; batch_size=32, verbose=0)
    return m.obj[:predict](X, batch_size=batch_size, verbose=verbose)
end

export Model, Sequential, layers, inputs, outputs, evaluate, add!, compile!


text:
@testset "Expectation" begin
    Mod.@polyvar x[1:3]
    p = x[3] - 2x[1]*x[2]^2 + 3x[3]*x[1] - 5x[1]^3
    v = (1,2,3)
    m = dirac(monomials(p), x=>v)
    @test MultivariateMoments.expectation(m, p) == p(x=>v) == MultivariateMoments.expectation(p, m)
    @test_throws ErrorException dot(x[1] * x[2] * x[3], m)
    @test dot(0.5 * x[1] * x[2]^2, m) == 2.0
    @test dot(m, x[1] * x[3]) == 3
end


text:
module CUDD

# cudd_path = joinpath(dirname(@__FILE__), "..", "deps", "cudd-3.0.0", "cudd", ".libs")
const depfile = joinpath(@__DIR__, "..", "deps", "deps.jl")

if isfile(depfile)
    include(depfile)
else
    error("CUDD not properly installed. Please run Pkg.build(\"CUDD\")")
end

using CEnum

# Files automatically generated by Clang.jl
include(joinpath(@__DIR__, "gen", "libcudd_common.jl"))
include(joinpath(@__DIR__, "gen", "libcudd_api.jl"))

# Missing C macros omitted by Clang.jl

Cudd_Not(node) =
    convert(Ptr{Nothing}, xor(convert(Int,node), 1))

Cudd_IsComplement(node) =
    isone(convert(Int,node) & 1)

Cudd_Regular(node) =
    convert(Ptr{Nothing}, convert(Int,node) & ~1)

# export everything
foreach(names(@__MODULE__, all=true)) do s
    if startswith(string(s), "Cudd_")
        @eval export $s
    end
end

# More user-friendly interface for ADD
include("ADD.jl")
include("ADD_apply.jl")
include("ADD_find.jl")
include("ADD_abs.jl")
include("ADD_inverse.jl")
include("ADD_ITE.jl")
include("ADD_negate.jl")
include("ADD_walsh.jl")

end # module


text:
module Discrete

export ConstSequence, Sequence, DFun;
export fdiff, bdiff, cdiff;
export euler, runge_kutta;
export trapezoid, simpsons;

import FastAnonymous;

immutable ConstSequence
  i_start::Int;
  values::Vector{Float64};

  Sequence(i_start::Int, values::Vector{Float64}) = new(i_start, values);
end

function call(s::ConstSequence, n::Int)
  if n < s.i_start || n > s.i_start + length(s.values)
    error("Out of domain on which sequence is defined");
  end

  return s.values[s.i_start + n];
end

function call(s::ConstSequence, nreal::Real)
  const n = convert(Int, nreal);
  if n < s.i_start || n > s.i_start + length(s.values)
    error("Out of domain on which sequence is defined");
  end

  return s.values[s.i_start + n];
end

typealias DFun     Union{Function, FastAnonymous.Fun};
typealias Sequence Union{DFun, ConstSequence};

#! Compute the forward difference of a sequence
#!
#! \param     seq     Sequence
#! \param     n       Index at which to take compute difference
#! \param     order   Order of difference
#! \return    Forward difference
function fdiff(seq::Sequence, n::Int; order::Int=1)
  @assert(order > 0, "unable to compute $order order forward difference");

  sum = 0.0;
  for i=0:order
    sum += (((order-i)%2==0) ?
              binomial(order, i) * seq(n + i) :
             -binomial(order, i) * seq(n + i));
  end

  return sum;
end

#! Compute the forward difference of a sequence
#!
#! \param     seq     Sequence
#! \param     n_range Range of index for which to compute differences
#! \param     order   Order of difference
#! \return    Forward differences
function fdiff(seq::Sequence, n_range::UnitRange{Int}; order::Int=1)
  return pmap(n -> fdiff(seq, n, order=order), n_range);
end

#! Compute the forward difference of a sequence
#!
#! \param     seq     Sequence
#! \param     n_range Range of index for which to compute differences
#! \param     order   Order of difference
#! \return    Forward differences
function fdiff(seq::Sequence, n_range::UnitRange{Float64}; order::Int=1)
  return pmap(n -> fdiff(seq, n, order=order),
              convert(UnitRange{Int}, n_range));
end

#! Compute the backward difference of a sequence
#!
#! \param     seq     Sequence
#! \param     n       Index at which to take compute difference
#! \param     order   Order of difference
#! \return    Backward difference
function bdiff(seq::Sequence, n::Int; order::Int=1)
  @assert(order > 0, "unable to compute $order order backward difference");

  sum = 0.0;
  for i=0:order
    sum += ((i % 2 == 0) ?
              binomial(order, i) * seq(n - i) :
             -binomial(order, i) * seq(n - i));
  end

  return sum;
end

#! Compute the backward difference of a sequence
#!
#! \param     seq     Sequence
#! \param     n_range Range of index for which to compute differences
#! \param     order   Order of difference
#! \return    Backward differences
function bdiff(seq::Sequence, n_range::UnitRange{Int}; order::Int=1)
  return pmap(n -> bdiff(seq, n, order=order), n_range);
end

#! Compute the backward difference of a sequence
#!
#! \param     seq     Sequence
#! \param     n_range Range of index for which to compute differences
#! \param     order   Order of difference
#! \return    Backward differences
function bdiff(seq::Sequence, n_range::UnitRange{Float64}; order::Int=1)
  return pmap(n -> bdiff(seq, n, order=order),
              convert(UnitRange{Int}, n_range));
end

#! Calculate forward difference on a continuous function
#!
#! \param   f     Continuous function
#! \param   x     Value to calculate difference at
#! \param   h     Step size
#! \param   order Order of difference
#! \retrun  Difference
function fdiff(f::DFun, x::Real, h::Real; order=1)
  sum = 0.0;
  for i=0:order
    sum += ((i % 2 == 0) ? binomial(order, i) * f(x + (n - i) * h) :
                          -binomial(order, i) * f(x + (n - i) * h));
  end
  return sum;
end

#! Calculate backward difference on a continuous function
#!
#! \param   b     Continuous function
#! \param   x     Value to calculate difference at
#! \param   h     Step size
#! \param   order Order of difference
#! \retrun  Difference
function bdiff(f::DFun, x::Real, h::Real; order=1)
  sum = 0.0;
  for i=0:order
    sum += ((i % 2 == 0) ? binomial(order, i) * f(x - i * h) :
                          -binomial(order, i) * f(x - i * h));
  end
  return sum;
end

#! Calculate central difference on a continuous function
#!
#! \param   b     Continuous function
#! \param   x     Value to calculate difference at
#! \param   h     Step size
#! \param   order Order of difference
#! \retrun  Difference
function cdiff(f::DFun, x::Real, h::Real; order=1)
  sum = 0.0;
  for i=0:order
    sum += ((i % 2 == 0) ? binomial(order, i) * f(x + (n/2 - i) * h) :
                          -binomial(order, i) * f(x + (n/2 - i) * h));
  end
  return sum;
end

#! Approximate the solution to an ODE using Euler's method
#!
#! \param   df      Differential equation, df = df(x,t)
#! \param   x0      Initial value of x
#! \param   t0      Initial value of t
#! \param   tf      Final value of t
#! \param   nsteps  Number of timesteps
#! \param   Approximate x(tf)
function euler(df::DFun, x0::Real, t0::Real, tf::Real, nsteps::Int)
  const h   =   (tf - t0) / nsteps;
  x         =   x0;
  t         =   t0;

  for i=1:nsteps
    x += h * df(x, t);
    t += h;
  end

  return x;
end

#! Approximate the solution to an ODE using Runge-Kutta method
#!
#! \param   df      Differential equation, df = df(x,t)
#! \param   x0      Initial value of x
#! \param   t0      Initial value of t
#! \param   tf      Final value of t
#! \param   nsteps  Number of timesteps
#! \param   Approximate x(tf)
function runge_kutta(df::DFun, x0::Real, t0::Real, tf::Real, nsteps::Int)
  const h   =   (tf - t0) / nsteps;
  x         =   x0;
  t         =   t0;

  for i=1:nsteps
    k1    =   h * df(x, t);
    k2    =   h * df(x + 0.5*k1, t + 0.5*h);
    k3    =   h * df(x + 0.5*k2, t + 0.5*h);
    k4    =   h * df(x + k3, t + h);

    x    +=   1.0/6.0 * (k1 + 2*k2 + 2*k3 + k4);
  end

  return x;
end

#! Integrate using trapezoidal rule
#!
#! \param   f     Function to integrate
#! \param   x_0   Starting x
#! \param   x_1   Ending x
#! \param   n     Number of discrete points (must be even)
#! \return        Approximate integral
function trapezoid(f::DFun, x_0::Real, x_1::Real, n::Int)
  const dx  =   (x_1 - x_0) / n;
  sum       =   dx * (f(x_1)+f(x_0)) / 2.0;

  for (i, x) = enumerate(linspace(x_0+dx, x_1-dx, n-1));
    sum += dx * f(x);
  end

  return sum;
end

#! Integrate using simpsons rule
#!
#! \param   f     Function to integrate
#! \param   x_0   Starting x
#! \param   x_1   Ending x
#! \param   n     Number of discrete points (must be even)
#! \return        Approximate integral
function simpsons(f::DFun, x_0::Real, x_1::Real, n::Int)
  @assert(n % 2 == 0, "Number of integration points must be even for "         *
                      "Simpson's rule. $n is not an even number.");

  const dx  =   (x_1 - x_0) / n;
  sum       =   dx * (f(x_1)+f(x_0)) / 3.0;

  for (i, x) = enumerate(linspace(x_0+dx, x_1-dx, n-1));
    sum += dx * f(x) * ((i % 2 == 0) ? 2.0/3.0 : 4.0/3.0);
  end

  return sum;
end

end # end module


text:
#=
A minor syntax issue with GtkWindowLeaf attributes #392 -  Open
pauljurczak opened this issue on Oct 9, 2018 · 4 comments

@pauljurczak pauljurczak commented on Oct 9, 2018 •
Evaluating GtkWindowLeaf outputs:
GtkWindowLeaf(name="", parent, width-request=600, height-request=-1, visible=TRUE, ...
Can it be modified to:
GtkWindowLeaf(name="", parent, width_request=600, height_request=-1, visible=true, ...
so it's syntactically correct in Julia?

Collaborator
@tknopp tknopp commented on Oct 9, 2018
not sure why booleans are printed uppercase. @vtjnash may know.

Author
@pauljurczak pauljurczak commented on Oct 9, 2018
@tknopp There is also width-request vs width_request.

Collaborator
@tknopp tknopp commented on Oct 9, 2018
There is a reason for that but I actually don't know exactly what since I did not design that part of Gtk.jl

Collaborator
@jonathanBieler jonathanBieler commented on Oct 9, 2018 •
I think the reason things are called like is that that's how they are called in Gtk.

Takeaway: Gtk.jl just calls a GObject introspection method and prints the results:

Gtk.jl/src/GLib/gvalues.jl

Lines 209 to 238 in 4190f44
=#
#Problem:
function gtkwindowleafcall(GObject, GValue)
    props = ccall((:g_object_class_list_properties, libgobject), Ptr{Ptr{GParamSpec}},
        (Ptr{Nothing}, Ptr{Cuint}), G_OBJECT_GET_CLASS(w), n)
    v = gvalue(String)
    first = true
    for i = 1:unsafe_load(n)
        param = unsafe_load(unsafe_load(props, i))
        if !first
            print(io, ", ")
        else
            first = false
        end
        print(io, GLib.bytestring(param.name))
        if (param.flags & READABLE) != 0 &&
           (param.flags & DEPRECATED) == 0 &&
           (ccall((:g_value_type_transformable, libgobject), Cint,
               (Int, Int), param.value_type, g_type(AbstractString)) != 0)
            ccall((:g_object_get_property, libgobject), Nothing,
                (Ptr{GObject}, Ptr{UInt8}, Ptr{GValue}), w, param.name, v)
            str = ccall((:g_value_get_string, libgobject), Ptr{UInt8}, (Ptr{GValue},), v)
            value = (str == C_NULL ? "NULL" : GLib.bytestring(str))
            if param.value_type == g_type(AbstractString) && str != C_NULL
                print(io, "=\"", value, '"')
            else
                print(io, '=', value)
            end
        end
    end
    print(io, ')')
    ccall((:g_value_unset, libgobject), Ptr{Nothing}, (Ptr{GValue},), v)
end
#ERROR: Invalid Access to an Invalid a Slot Number
#I guess one could "translate" those names into Julia friendly ones, since that's just for showing in the REPL it should be safe.

max_stars_repo_name: adamwillisMastery/MiniGtk

max_issues_repo_name: adamwillisMastery/MiniGtk

max_forks_repo_name: adamwillisMastery/MiniGtk


text:
######################################################################
# RuntimeEnums UTs
# -----
# Licensed under MIT License.
using Test
using RuntimeEnums

@runtime_enum OneLiner::UInt8 Byte1 = 1 Byte2 Byte3
@runtime_enum Blocked::Int begin
    SixtyNine = 69
    FourTwenty = 420
    FourtyTwo = 42
end

@runtime_enum Extendable Error Predef1 Predef2

@testset "RuntimeEnums" begin
    @testset "1-liner" begin
        @test OneLiner <: Enum{UInt8}
        @test Byte1 isa OneLiner && Byte2 isa OneLiner && Byte3 isa OneLiner
        @test issetequal(Base.Enums.instances(OneLiner), (Byte1, Byte2, Byte3))
        @test Byte1 === OneLiner(1) && Byte2 === OneLiner(2) && Byte3 === OneLiner(3)
        @test typemin(OneLiner) === Byte1 && typemax(OneLiner) === Byte3
    end

    @testset "Blocked" begin
        @test Blocked <: Enum{Int}
        @test FourtyTwo isa Blocked && SixtyNine isa Blocked && FourTwenty isa Blocked
        @test issetequal(Base.Enums.instances(Blocked), (FourtyTwo, SixtyNine, FourTwenty))
        @test FourtyTwo === Blocked(42) && SixtyNine === Blocked(69) && FourTwenty === Blocked(420)
        @test typemin(Blocked) === FourtyTwo && typemax(Blocked) === FourTwenty
    end

    @testset "Extendable" begin
        @assert issetequal(Base.Enums.instances(Extendable), (Error, Predef1, Predef2))
        @assert typemin(Extendable) === Error && typemax(Extendable) === Predef2

        @test_throws ArgumentError Extendable[:Extended1] = 0
        @test_throws ArgumentError Extendable[:Predef1] = 0

        Extendable[:NewMax] =  42
        Extendable[:NewMin] = -12
        @test issetequal(Base.Enums.instances(Extendable), (Error, Predef1, Predef2, Extendable(42), Extendable(-12)))
        @test typemin(Extendable) === Extendable(-12) && typemax(Extendable) === Extendable(42)
    end
end


text:
# record.jl
# =========
#
# Representation of a record of a VCF file.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/GeneticVariation.jl/blob/master/LICENSE

mutable struct Record
    # data and filled range
    data::Vector{UInt8}
    filled::UnitRange{Int}
    # indexes
    chrom::UnitRange{Int}
    pos::UnitRange{Int}
    id::Vector{UnitRange{Int}}
    ref::UnitRange{Int}
    alt::Vector{UnitRange{Int}}
    qual::UnitRange{Int}
    filter::Vector{UnitRange{Int}}
    infokey::Vector{UnitRange{Int}}
    format::Vector{UnitRange{Int}}
    genotype::Vector{Vector{UnitRange{Int}}}
end

"""
    VCF.Record()
Create an unfilled VCF record.
"""
function Record()
    return Record(
        # data and filled
        UInt8[], 1:0,
        # chrom-alt
        1:0, 1:0, UnitRange{Int}[], 1:0, UnitRange{Int}[],
        # qual-genotype
        1:0, UnitRange{Int}[], UnitRange{Int}[], UnitRange{Int}[], UnitRange{Int}[])
end

"""
    VCF.Record(data::Vector{UInt8})
Create a VCF object from `data` containing a VCF record.
This function verifies the format and indexes fields for accessors.
Note that the ownership of `data` is transferred to a new record.
"""
function Record(data::Vector{UInt8})
    return convert(Record, data)
end

function Base.convert(::Type{Record}, data::Vector{UInt8})
    record = Record(
        # data and filled
        data, 1:0,
        # chrom-alt
        1:0, 1:0, UnitRange{Int}[], 1:0, UnitRange{Int}[],
        # qual-genotype
        1:0, UnitRange{Int}[], UnitRange{Int}[], UnitRange{Int}[], UnitRange{Int}[])
    index!(record)
    return record
end

"""
    VCF.Record(str::AbstractString)
Create a VCF object from `str` containing a VCF record.
This function verifies the format and indexes fields for accessors.
"""
function Record(str::AbstractString)
    return convert(Record, str)
end

function Base.convert(::Type{Record}, str::AbstractString)
    return Record(Vector{UInt8}(str))
end

function initialize!(record::Record)
    record.filled = 1:0
    record.chrom = 1:0
    record.pos = 1:0
    empty!(record.id)
    record.ref = 1:0
    empty!(record.alt)
    record.qual = 1:0
    empty!(record.filter)
    empty!(record.infokey)
    empty!(record.format)
    empty!(record.genotype)
    return record
end

function isfilled(record::Record)
    return !isempty(record.filled)
end

function datarange(record::Record)
    return record.filled
end

function Base.:(==)(record1::Record, record2::Record)
    if isfilled(record1) == isfilled(record2) == true
        r1 = datarange(record1)
        r2 = datarange(record2)
        return length(r1) == length(r2) && memcmp(pointer(record1.data, first(r1)), pointer(record2.data, first(r2)), length(r1)) == 0
    else
        return isfilled(record1) == isfilled(record2) == false
    end
end

function checkfilled(record::Record)
    if !isfilled(record)
        throw(ArgumentError("unfilled VCF record"))
    end
end

function Record(base::Record;
                   chrom=nothing, pos=nothing, id=nothing,
                   ref=nothing, alt=nothing, qual=nothing,
                   filter=nothing, info=nothing, genotype=nothing)
    checkfilled(base)
    buf = IOBuffer()

    if chrom == nothing
        write(buf, base.data[base.chrom])
    else
        print(buf, string(chrom))
    end

    print(buf, '\t')
    if pos == nothing
        write(buf, base.data[base.pos])
    else
        print(buf, convert(Int, pos))
    end

    print(buf, '\t')
    if id == nothing
        if isempty(base.id)
            print(buf, '.')
        else
            for (i, r) in enumerate(base.id)
                if i != 1
                    print(buf, ';')
                end
                write(buf, base.data[r])
            end
        end
    else
        if !isa(id, Vector)
            id = [id]
        end
        if isempty(id)
            print(buf, '.')
        else
            for (i, x) in enumerate(id)
                if i != 1
                    print(buf, ';')
                end
                print(buf, string(x))
            end
        end
    end

    print(buf, '\t')
    if ref == nothing
        write(buf, base.data[base.ref])
    else
        print(buf, string(ref))
    end

    print(buf, '\t')
    if alt == nothing
        if isempty(base.alt)
            print(buf, '.')
        else
            for (i, r) in enumerate(base.alt)
                if i != 1
                    print(buf, ';')
                end
                write(buf, base.data[r])
            end
        end
    else
        if !isa(alt, Vector)
            alt = [alt]
        end
        if isempty(alt)
            print(buf, '.')
        else
            for (i, x) in enumerate(alt)
                if i != 1
                    print(buf, ';')
                end
                print(buf, string(x))
            end
        end
    end

    print(buf, '\t')
    if qual == nothing
        write(buf, base.data[base.qual])
    else
        print(buf, convert(Float64, qual))
    end

    print(buf, '\t')
    if filter == nothing
        if isempty(base.filter)
            print(buf, '.')
        else
            for (i, r) in enumerate(base.filter)
                if i != 1
                    print(buf, ';')
                end
                write(buf, base.data[r])
            end
        end
    else
        if !isa(filter, Vector)
            filter = [filter]
        end
        if isempty(filter)
            print(buf, '.')
        else
            for (i, x) in enumerate(filter)
                if i != 1
                    print(buf, ';')
                end
                print(buf, string(x))
            end
        end
    end

    print(buf, '\t')
    if info == nothing
        if isempty(base.infokey)
            print(buf, '.')
        else
            write(buf, base.data[first(base.infokey[1]):last(infovalrange(base, lastindex(base.infokey)))])
        end
    else
        if !isa(info, AbstractDict)
            throw(ArgumentError("info must be an AbstractDict"))
        elseif isempty(info)
            print(buf, '.')
        else
            for (i, (key, val)) in enumerate(info)
                if i != 1
                    print(buf, ';')
                end
                print(buf, string(key))
                if val != nothing
                    print(buf, '=', vcfformat(val))
                end
            end
        end
    end

    print(buf, '\t')
    if genotype == nothing
        if isempty(base.format)
            print(buf, '.')
        else
            write(buf, base.data[first(base.format[1]):last(base.format[end])])
        end
        if !isempty(base.genotype)
            for indiv in base.genotype
                print(buf, '\t')
                for (i, r) in enumerate(indiv)
                    if i != 1
                        print(buf, ':')
                    end
                    write(buf, base.data[r])
                end
            end
        end
    else
        if !isa(genotype, Vector)
            genotype = [genotype]
        end
        if isempty(genotype)
            print(buf, '.')
        else
            allkeys = String[]
            for indiv in genotype
                if !isa(indiv, AbstractDict)
                    throw(ArgumentError("individual must be anabstract dictionary"))
                end
                append!(allkeys, keys(indiv))
            end
            allkeys = sort!(unique(allkeys))
            if !isempty(allkeys)
                join(buf, allkeys, ':')
                for indiv in genotype
                    print(buf, '\t')
                    for (i, key) in enumerate(allkeys)
                        if i != 1
                            print(buf, ':')
                        end
                        print(buf, vcfformat(get(indiv, key, '.')))
                    end
                end
            end
        end
    end

    return Record(take!(buf))
end

function vcfformat(val)
    return string(val)
end

function vcfformat(val::Vector)
    return join(map(vcfformat, val), ',')
end

function Base.copy(record::Record)
    return Record(
        copy(record.data),
        record.filled,
        record.chrom,
        record.pos,
        copy(record.id),
        record.ref,
        copy(record.alt),
        record.qual,
        copy(record.filter),
        copy(record.infokey),
        copy(record.format),
        deepcopy(record.genotype))
end

function Base.write(io::IO, record::Record)
    checkfilled(record)
    return write(io, record.data)
end

# Accessor functions
# ------------------

"""
    chrom(record::Record)::String
Get the chromosome name of `record`.
"""
function chrom(record::Record)::String
    checkfilled(record)
    if ismissing(record, record.chrom)
        missingerror(:chrom)
    end
    return String(record.data[record.chrom])
end

function haschrom(record::Record)
    return isfilled(record) && !ismissing(record, record.chrom)
end

"""
    pos(record::Record)::Int
Get the reference position of `record`.
"""
function pos(record::Record)::Int
    checkfilled(record)
    if ismissing(record, record.pos)
        missingerror(:pos)
    end
    # TODO: no-copy accessor
    return parse(Int, String(record.data[record.pos]))
end

function haspos(record::Record)
    return isfilled(record) && !ismissing(record, record.pos)
end

"""
    id(record::Record)::Vector{String}
Get the identifiers of `record`.
"""
function id(record::Record)::Vector{String}
    checkfilled(record)
    if isempty(record.id)
        missingerror(:id)
    end
    return [String(record.data[r]) for r in record.id]
end

function hasid(record::Record)
    return isfilled(record) && !isempty(record.id)
end

"""
    ref(record::Record)::String
Get the reference bases of `record`.
"""
function ref(record::Record)::String
    checkfilled(record)
    if ismissing(record, record.ref)
        missingerror(:ref)
    end
    return String(record.data[record.ref])
end

function hasref(record::Record)
    return isfilled(record) && !ismissing(record, record.ref)
end

"""
    alt(record::Record)::Vector{String}
Get the alternate bases of `record`.
"""
function alt(record::Record)::Vector{String}
    checkfilled(record)
    if isempty(record.alt)
        missingerror(:alt)
    end
    return [String(record.data[r]) for r in record.alt]
end

function hasalt(record::Record)
    return isfilled(record) && !isempty(record.alt)
end

"""
    qual(record::Record)::Float64
Get the quality score of `record`.
"""
function qual(record::Record)::Float64
    checkfilled(record)
    if ismissing(record, record.qual)
        missingerror(:qual)
    end
    # TODO: no-copy parse
    return parse(Float64, String(record.data[record.qual]))
end

function hasqual(record::Record)
    return isfilled(record) && !ismissing(record, record.qual)
end

"""
    filter(record::Record)::Vector{String}
Get the filter status of `record`.
"""
function filter(record::Record)::Vector{String}
    checkfilled(record)
    if isempty(record.filter)
        missingerror(:filter)
    end
    return [String(record.data[r]) for r in record.filter]
end

function hasfilter(record::Record)
    return isfilled(record) && !isempty(record.filter)
end

"""
    info(record::Record)::Vector{Pair{String,String}}
Get the additional information of `record`.
"""
function info(record::Record)::Vector{Pair{String,String}}
    checkfilled(record)
    if isempty(record.infokey)
        missingerror(:info)
    end
    ret = Pair{String,String}[]
    for (i, key) in enumerate(record.infokey)
        val = infovalrange(record, i)
        push!(ret, String(record.data[key]) => String(record.data[val]))
    end
    return ret
end

function hasinfo(record::Record)
    return isfilled(record) && !isempty(record.infokey)
end

"""
    info(record::Record, key::String)::String
Get the additional information of `record` with `key`.
Keys without corresponding values return an empty string.
"""
function info(record::Record, key::String)::String
    checkfilled(record)
    i = findinfokey(record, key)
    if i == 0
        throw(KeyError(key))
    end
    val = infovalrange(record, i)
    if isempty(val)
        return ""
    else
        return String(record.data[val])
    end
end

function hasinfo(record::Record, key::String)
    return isfilled(record) && findinfokey(key) > 0
end

function findinfokey(record::Record, key::String)
    for i in 1:lastindex(record.infokey)
        if isequaldata(key, record.data, record.infokey[i])
            return i
        end
    end
    return 0
end

"""
    infokeys(record::Record)::Vector{String}
Get the keys of the additional information of `record`.
This function returns an empty vector when the INFO field is missing.
"""
function infokeys(record::Record)::Vector{String}
    checkfilled(record)
    return [String(record.data[key]) for key in record.infokey]
end

# Returns the data range of the `i`-th value.
function infovalrange(record::Record, i::Int)
    checkfilled(record)
    data = record.data
    key = record.infokey[i]
    if last(key) + 1 ≤ lastindex(data) && data[last(key)+1] == UInt8('=')
        endpos = something(findnext(isequal(0x3B), data, last(key)+1), 0) # 0x3B is byte equivalent of char ';'.
        if endpos == 0
            endpos = something(findnext(isequal(0x09), data, last(key)+1), 0) # 0x09 is byte equivalent of char '\t'
            @assert endpos != 0
        end
        return last(key)+2:endpos-1
    else
        return last(key)+1:last(key)
    end
end

"""
    format(record::Record)::Vector{String}
Get the genotype format of `record`.
"""
function format(record::Record)::Vector{String}
    checkfilled(record)
    if isempty(record.format)
        missingerror(:format)
    end
    return [String(record.data[r]) for r in record.format]
end

function hasformat(record::Record)
    return isfilled(record) && !isempty(record.format)
end

"""
    genotype(record::Record)::Vector{Vector{String}}
Get the genotypes of `record`.
"""
function genotype(record::Record)
    checkfilled(record)
    ret = Vector{String}[]
    for i in 1:lastindex(record.genotype)
        push!(ret, genotype_impl(record, i, 1:lastindex(record.format)))
    end
    return ret
end

"""
    genotype(record::Record, index::Integer)::Vector{String}
Get the genotypes of the `index`-th individual in `record`.
This is effectively equivalent to `genotype(record)[index]` but more efficient.
"""
function genotype(record::Record, index::Integer)
    checkfilled(record)
    return genotype_impl(record, index, 1:lastindex(record.format))
end

"""
    genotype(record::Record, indexes, keys)
Get the genotypes in `record` that match `indexes` and `keys`.
`indexes` and `keys` can be either a scalar or a vector value.
Trailing fields that are dropped are filled with `"."`.
"""
function genotype(record::Record, index::Integer, key::String)::String
    checkfilled(record)
    k = findgenokey(record, key)
    if k == nothing
        throw(KeyError(key))
    end
    return genotype_impl(record, index, k)
end

function genotype(record::Record, index::Integer, keys::AbstractVector{String})::Vector{String}
    checkfilled(record)
    return [genotype(record, index, key) for key in keys]
end

function genotype(record::Record, indexes::AbstractVector{T}, key::String)::Vector{String} where T<:Integer
    checkfilled(record)
    k = findgenokey(record, key)
    if k == nothing
        throw(KeyError(key))
    end
    return [genotype_impl(record, i, k) for i in indexes]
end

function genotype(record::Record, indexes::AbstractVector{T}, keys::AbstractVector{String})::Vector{Vector{String}} where T<:Integer
    checkfilled(record)
    ks = Vector{Int}(length(keys))
    for i in 1:lastindex(keys)
        key = keys[i]
        k = findgenokey(record, key)
        if k == nothing
            throw(KeyError(key))
        end
        ks[i] = k
    end
    return [genotype_impl(record, i, ks) for i in indexes]
end

function genotype(record::Record, ::Colon, key::String)::Vector{String}
    return genotype(record, 1:lastindex(record.genotype), key)
end

function findgenokey(record::Record, key::String)
    return findfirst(r -> isequaldata(key, record.data, r), record.format)
end

function genotype_impl(record::Record, index::Int, keys::AbstractVector{Int})
    return [genotype_impl(record, index, k) for k in keys]
end

function genotype_impl(record::Record, index::Int, key::Int)
    geno = record.genotype[index]
    if key > lastindex(geno)  # dropped field
        return "."
    else
        return String(record.data[geno[key]])
    end
end

function Base.show(io::IO, record::Record)
    print(io, summary(record), ':')
    if isfilled(record)
        println(io)
        println(io, "   chromosome: ", haschrom(record) ? chrom(record) : "<missing>")
        println(io, "     position: ", haspos(record) ? pos(record) : "<missing>")
        println(io, "   identifier: ", hasid(record) ? join(id(record), " ") : "<missing>")
        println(io, "    reference: ", hasref(record) ? ref(record) : "<missing>")
        println(io, "    alternate: ", hasalt(record) ? join(alt(record), " ") : "<missing>")
        println(io, "      quality: ", hasqual(record) ? qual(record) : "<missing>")
        println(io, "       filter: ", hasfilter(record) ? join(filter(record), " ") : "<missing>")
          print(io, "  information: ")
        if hasinfo(record)
            for (key, val) in info(record)
                print(io, key)
                if !isempty(val)
                    print(io, '=', val)
                end
                print(io, ' ')
            end
        else
            print(io, "<missing>")
        end
        println(io)
          print(io, "       format: ", hasformat(record) ? join(format(record), " ") : "<missing>")
        if hasformat(record)
            println(io)
            print(io, "     genotype:")
            for i in 1:lastindex(record.genotype)
                print(io, " [$(i)] ", let x = genotype(record, i); isempty(x) ? "." : join(x, " "); end)
            end
        end
    else
        print(io, " <not filled>")
    end
end

function ismissing(record::Record, range::UnitRange{Int})
    return length(range) == 1 && record.data[first(range)] == UInt8('.')
end

# Check if `str == data[range]`
function isequaldata(str::String, data::Vector{UInt8}, range::UnitRange{Int})
    rlen = length(range)
    return rlen == sizeof(str) && memcmp(pointer(data, first(range)), pointer(str), rlen) == 0
end

function memcmp(p1::Ptr, p2::Ptr, n::Integer)
    return ccall(:memcmp, Cint, (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t), p1, p2, n)
end


text:
##################################################################
# Deprecations
##################################################################

@deprecate bicgstab! bicgstab


text:
# Copyright (c) 2020-2021, Bank of Canada
# All rights reserved.

# __precompile__(false)

# CHLI library: load and bindings

using Libdl

###  check_status

export check_status

try
    include("../deps/FAMEMessages.jl")
catch
    @error "Try re-building FAME"
end

"""
    struct HLIError <: Exception … end

Exception type for errors returned by CHLI calls.
"""
struct HLIError <: Exception
    status::Int32
    msg::String
end
HLIError(st::Int32) = HLIError(st,
    st == 513 ? cfmferr() : get(chli_status_description, st, "Unknown CHLI error."))
Base.showerror(io::IO, e::HLIError) = print(io, "HLI status($(e.status)): $(e.msg)")

"""
    check_status(status)

Check the status code returned by cfmXYZ functions. If status indicates success
we do nothing, otherwise we trigger an HLIError with the error code and the
message.
"""
check_status(status::Ref{Int32}) = check_status(status[])
@inline function check_status(code::Int32)
    code == HSUCC && return
    throw(HLIError(code))
end

function Chli()
    if !haskey(Base.ENV, "FAME")
        @error "FAME environment variable is not set! Will not be able to use FAME."
        return Chli{AbstractFameDatabase}(nothing)
    end
    @static if Sys.iswindows()
        # we need only the basename of the .dll; Windows automatically
        # decorates it with extension and finds it on the PATH
        chli_lib_path = "chli"
        # ensure  chli.dll file is on the PATH so Windows can find it
        @static if Sys.ARCH == :x86_64
            Base.ENV["PATH"] *= ";" * Base.ENV["FAME"] * "\\64"
        else
            Base.ENV["PATH"] *= ";" * Base.ENV["FAME"]
        end
    elseif Sys.islinux()
        # in Linux libchli.so is located in a subdirectory hli
        # we need to specify the full path to the .so file
        @static if Sys.ARCH == :x86_64
            chli_lib_path = joinpath(Base.ENV["FAME"], "hli", "64", "libchli.so")
        else
            chli_lib_path = joinpath(Base.ENV["FAME"], "hli", "libchli.so")
        end
    else
        @error "Your OS is not supported. Will not be able to use FAME."
        return Chli{AbstractFameDatabase}(nothing)
    end
    try
        return Chli{FameDatabase}(Libdl.dlopen(chli_lib_path), FameDatabase[], Ref(FameDatabase(-1)))
    catch
        @error "Failed to load CHLI library! Will not be able to use FAME."
        return Chli{AbstractFameDatabase}(nothing)
    end
end

"""
    @cfm_global(:cfmXYZ, type)

Fetch the value of a constant by the given name from the CHLI library.

See also: [`@cfm_call_check`](@ref)
"""
macro cfm_global(symbol, type)
    sym = :($(@__MODULE__).getsym($(QuoteNode(symbol))))
    call_expr = Expr(:call, :cglobal, sym, type)
    return esc(:(unsafe_load($call_expr)))
end

"""
    @fame_call(:fame_XYZ, argTypes, args...)

Build an appropriate ccall to the given fame_XYZ function.

See also: [`@cfm_call_check`](@ref)
"""
macro fame_call(symbol, argTypes, args...)
    sym = :($(@__MODULE__).getsym($(QuoteNode(symbol))))
    call_expr = Expr(:call, :ccall, sym, :Cint, argTypes, args...)
    return esc(:($call_expr))
end

"""
    @fame_call_check(:fame_XYZ, argTypes, args...)

Build an appropriate ccall to the given fame_XYZ function. This macro also
checks the status and throws an error if not success.

See also: [`@cfm_call_check`](@ref)
"""
macro fame_call_check(symbol, argTypes, args...)
    sym = :($(@__MODULE__).getsym($(QuoteNode(symbol))))
    call_expr = Expr(:call, :ccall, sym, :Cint, argTypes, args...)
    return esc(:(check_status($call_expr)))
end

"""
    @cfm_call(:cfmXYZ, argTypes, args...)

Build an appropriate ccall to the given cfmXYZ function. When using this macro,
include the status variable and check it yourself.

See also: [`@cfm_call_check`](@ref)
"""
macro cfm_call(symbol, argTypes, args...)
    sym = :($(@__MODULE__).getsym($(QuoteNode(symbol))))
    call_expr = Expr(:call, :ccall, sym, :Cvoid, argTypes, args...)
    return esc(:($call_expr))
end

"""
    @cfm_call_check(:cfmXYZ, argTypes, args...)

Call the cfm function and check status. When using this macro, do not include
the status variable.

See also: [`@cfm_call`](@ref)
"""
macro cfm_call_check(symbol)
    sym = :($(@__MODULE__).getsym($(QuoteNode(symbol))))
    return quote
        local st = Ref{Cint}(-1)
        ccall($sym, Cvoid, (Ref{Cint},), st)
        check_status(st)
    end |> esc
end
macro cfm_call_check(symbol, argTypes, args...)
    sym = :($(@__MODULE__).getsym($(QuoteNode(symbol))))
    pushfirst!(argTypes.args, Ref{Cint})
    call_expr = Expr(:call, :ccall, sym, :Cvoid, argTypes, :st, args...)
    return quote
        local st = Ref{Cint}(-1)
        $call_expr
        check_status(st)
    end |> esc
end

cfmferr() = (
    errtxt = repeat(' ', 200);
    @cfm_call_check(cfmferr, (Cstring,), errtxt);
    errtxt
)


text:
function wave_equation(m,M,y)

# evaluates the discrete matrix equation for the waves
#
# a1 b'' + a2 b' + (a3 + a4 lambda + a5 I lambda^2)b = 0
#
# where a1-a2 are tridiagonal matrices resulting from derivatives
# and a3-a5 are diagonal matrices. The eigenvalue is lambda = C
# (see eq. 37) B&M 2019. We recover the value of frequency from
# the definition of C
#
# input
#  m -  angular order of wave
#  C -  Coriolis parameter  (dimensionless)
#  M -  magnetic parameter  (dimensionless)
#  y -  meridional coordinate
#
# output
# A0,A1,A2 - matrices used in nonlinear eigenvalue problem
#
# Form of eigenvalue problem
# (A0 + A1 lambda + A2 lambda^2 ) b = 0
#

# size of matrix system
n = length(y);

# matrices defining first and second derivatives
d1,d2 = derivatives(y);

# assemble matrix a1 (second derivative)
y2 = y.*y;
y2c = -(y2.-1);
a1 = Diagonal(y2c) * d2;

# assemble matrix a2   (first derivative)
a2 = -2*Diagonal(y) * d1;

# assemble matrix a3   (constant term)
y2c_inv = 1.0./y2c;
a3 = -m^2*Diagonal(y2c_inv);

# now combine into A0, A1 and A2
A0 = convert(Array{Complex{Float64},2},a1+a2+a3);
A1 = convert(Array{Complex{Float64},2},m*Diagonal(ones(n))/M);
A2 = convert(Array{Complex{Float64},2},Diagonal(y2)/M);

return A0,A1,A2

end

max_stars_repo_name: bbuffett/MagRossby

max_issues_repo_name: bbuffett/MagRossby

max_forks_repo_name: bbuffett/MagRossby


text:
module VectorsOfBlockArrayCooTests

using Test
using Gridap.Arrays
using BlockArrays
using FillArrays
using LinearAlgebra

l = 10
b11 = [  i*[1 2; 3 4]  for i in 1:l ]
b21 = [  i*[5 6; 7 8; 9 10] for i in 1:l]
b32 = [  i*[5 6 7 8; 7 8 1 3] for i in 1:l]
blocks = (b11,b21,b32)
blockids = [(1,1),(2,1),(3,2)]
ax = Fill((blockedrange([2,3,2]), blockedrange([2,4])),l)
al = VectorOfBlockArrayCoo(blocks,blockids,ax)

@test al[Block(1,1)] === b11
@test al[Block(2,1)] === b21
@test al[Block(2),Block(1)] === b21

@test is_zero_block(al,2,2) == true
@test is_zero_block(al,2,1) == false
@test is_nonzero_block(al,2,1) == true
@test is_zero_block(al,Block(1,2)) == true
@test is_zero_block(al,Block(1),Block(2)) == true

@test blocksize(al) == (3,2)
@test blocksize(al,1) == 3
@test blocksize(al,2) == 2

for (i,b) in enumerateblocks(al)
  @test al[i] === b
end

b1 = [  i*[3,4] for i in 1:l]
b2 = [  i*[3,4,5,6] for i in 1:l]
blocks = (b1,b2,)
blockids = [(1,),(2,)]
ax = Fill((blockedrange([2,4]),),l)
bl = VectorOfBlockArrayCoo(blocks,blockids,ax)

b2 = [  i*[3,4,5,6] for i in 1:l]
blocks = (b2,)
blockids = [(2,)]
ax = Fill((blockedrange([2,4]),),l)
bl = VectorOfBlockArrayCoo(blocks,blockids,ax)

c1 = [  i*[7,8] for i in 1:l]
c2 = [  i*[3,4,5] for i in 1:l]
blocks = (c1,c2)
blockids = [(1,),(2,)]
ax = Fill((blockedrange([2,3,2]),),l)
cl = VectorOfBlockArrayCoo(blocks,blockids,ax)

dl = apply(MulKernel(),al,bl)
test_array(dl,[ a*b for (a,b) in zip(al,bl) ])

dl = apply(MulAddKernel(2,3),al,bl,cl)
test_array(dl,[ 2*a*b + 3*c for (a,b,c) in zip(al,bl,cl) ])

dl = apply(transpose,al)
test_array(dl,transpose.(al))

dl = apply(MulKernel(),al,apply(transpose,al))
test_array(dl,[ a*transpose(a) for a in al ])

# in-homogeneous case

l1 = 4
l2 = 6
a11 = vcat([i*[1 2; 3 4] for i in 1:l1 ],[i*[1 2 3; 3 4 5; 6 7 8] for i in 1:l2])
a12 = vcat([10*i*[1 2; 3 4] for i in 1:l1 ],[10*i*[1 2 3; 3 4 5; 6 7 8] for i in 1:l2])
blocks = (a11,a12)
blockids = [(1,1),(1,2)]
ax1 = (blockedrange([2,2]),blockedrange([2,2]))
ax2 = (blockedrange([3,3]),blockedrange([3,3]))
ax = CompressedArray([ax1,ax2],vcat(fill(1,l1),fill(2,l2)))
al = VectorOfBlockArrayCoo(blocks,blockids,ax)

b2 = vcat([i*[1,2] for i in 1:l1 ],[i*[1,2,3] for i in 1:l2])
blocks = (b2,)
blockids = [(2,)]
ax1 = (blockedrange([2,2]),)
ax2 = (blockedrange([3,3]),)
ax = CompressedArray([ax1,ax2],vcat(fill(1,l1),fill(2,l2)))
bl = VectorOfBlockArrayCoo(blocks,blockids,ax)

dl = apply(MulKernel(),al,bl)
test_array(dl,[ a*b for (a,b) in zip(al,bl) ])

dl = apply(MulAddKernel(2,3),al,bl,bl)
test_array(dl,[ 2*a*b + 3*c for (a,b,c) in zip(al,bl,bl) ])

dl = apply(transpose,al)
test_array(dl,transpose.(al))

dl = apply(MulKernel(),al,apply(transpose,al))
test_array(dl,[ a*transpose(a) for a in al ])

# Blocks of blocks (in-homogeneous case)

blocks = (al,)
blockids = [(1,2),]
_ax1 = (blockedrange([2,2]),blockedrange([2,2]))
_ax2 = (blockedrange([3,3]),blockedrange([3,3]))
ax1 = (blockedrange([_ax1[1],_ax1[1]]),blockedrange([_ax1[2],_ax1[2]]))
ax2 = (blockedrange([_ax2[1],_ax2[1]]),blockedrange([_ax2[2],_ax2[2]]))
ax = CompressedArray([ax1,ax2],vcat(fill(1,l1),fill(2,l2)))
aBl = VectorOfBlockArrayCoo(blocks,blockids,ax)

blocks = (bl,)
blockids = [(1,)]
_ax1 = (blockedrange([2,2]),)
_ax2 = (blockedrange([3,3]),)
ax1 = (blockedrange([_ax1[1],_ax1[1]]),)
ax2 = (blockedrange([_ax2[1],_ax2[1]]),)
ax = CompressedArray([ax1,ax2],vcat(fill(1,l1),fill(2,l2)))
bBl = VectorOfBlockArrayCoo(blocks,blockids,ax)

dl = apply(MulKernel(),aBl,bBl)
test_array(dl,[ a*b for (a,b) in zip(aBl,bBl) ])

dl = apply(MulAddKernel(2,3),aBl,bBl,bBl)
test_array(dl,[ 2*a*b + 3*c for (a,b,c) in zip(aBl,bBl,bBl) ])

dl = apply(MulKernel(),aBl,apply(transpose,aBl))
test_array(dl,[ a*transpose(a) for a in aBl ])

end # module


text:
function MOI.set(model::JuMP.Model, attr::ForwardObjectiveFunction, func::JuMP.AbstractJuMPScalar)
    JuMP.check_belongs_to_model(func, model)
    return MOI.set(model, attr, JuMP.moi_function(func))
end
MOI.set(model::JuMP.Model, attr::ForwardObjectiveFunction, func::Number) = MOI.set(model, attr, JuMP.AffExpr(func))

function MOI.set(model::JuMP.Model, attr::ForwardConstraintFunction, con_ref::JuMP.ConstraintRef, func::JuMP.AbstractJuMPScalar)
    JuMP.check_belongs_to_model(func, model)
    return MOI.set(model, attr, con_ref, JuMP.moi_function(func))
end
MOI.set(model::JuMP.Model, attr::ForwardConstraintFunction, con_ref::JuMP.ConstraintRef, func::Number) = MOI.set(model, attr, con_ref, JuMP.AffExpr(func))

function MOI.get(model::JuMP.Model, attr::ReverseObjectiveFunction)
    func = MOI.get(JuMP.backend(model), attr)
    return JuMP.jump_function(model, func)
end

function MOI.get(model::JuMP.Model, attr::ReverseConstraintFunction, con_ref::JuMP.ConstraintRef)
    JuMP.check_belongs_to_model(con_ref, model)
    moi_func = MOI.get(JuMP.backend(model), attr, JuMP.index(con_ref))
    return JuMP.jump_function(model, moi_func)
end

# FIXME Workaround for https://github.com/jump-dev/JuMP.jl/issues/2797
function _moi_get_result(model::MOI.ModelLike, args...)
    if MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED
        throw(OptimizeNotCalled())
    end
    return MOI.get(model, args...)
end
function _moi_get_result(model::MOIU.CachingOptimizer, args...)
    if MOIU.state(model) == MOIU.NO_OPTIMIZER
        throw(NoOptimizer())
    elseif MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED
        throw(OptimizeNotCalled())
    end
    return MOI.get(model, args...)
end
function MOI.get(model::JuMP.Model, attr::ForwardVariablePrimal, var_ref::JuMP.VariableRef)
    JuMP.check_belongs_to_model(var_ref, model)
    return _moi_get_result(JuMP.backend(model), attr, JuMP.index(var_ref))
end

"""
    abstract type AbstractLazyScalarFunction <: MOI.AbstractScalarFunction end

Subtype of `MOI.AbstractScalarFunction` that is not a standard MOI scalar
function but can be converted to one using [`standard_form`](@ref).

The function can also be inspected lazily using `JuMP.coefficient` or
[`quad_sym_half`](@ref).
"""
abstract type AbstractLazyScalarFunction <: MOI.AbstractScalarFunction end

"""
    standard_form(func::AbstractLazyScalarFunction)

Converts `func` to a standard MOI scalar function.

    standard_form(func::MOItoJuMP)

Converts `func` to a standard JuMP scalar function.
"""
function standard_form end

"""
    quad_sym_half(func, vi1::MOI.VariableIndex, vi2::MOI.VariableIndex)

Return `Q[i,j] = Q[j,i]` where the quadratic terms of `func` is represented
by `x' Q x / 2` for a symmetric matrix `Q` where `x[i] = vi1` and `x[j] = vi2`.
Note that while this is equal to `JuMP.coefficient(func, vi1, vi2)` if `vi1 != vi2`,
in the case `vi1 == vi2`, it is rather equal to
`2JuMP.coefficient(func, vi1, vi2)`.
"""
function quad_sym_half end

standard_form(func::Union{MOI.VariableIndex,MOI.ScalarAffineFunction,MOI.ScalarQuadraticFunction}) = func
function Base.isapprox(func1::AbstractLazyScalarFunction, func2::MOI.AbstractScalarFunction; kws...)
    return isapprox(standard_form(func1), standard_form(func2); kws...)
end

# In the future, we could replace by https://github.com/jump-dev/MathOptInterface.jl/pull/1238
"""
    VectorScalarAffineFunction{T, VT} <: MOI.AbstractScalarFunction

Represents the function `x ⋅ terms + constant`
as an `MOI.AbstractScalarFunction` where `x[i] = MOI.VariableIndex(i)`.
Use [`standard_form`](@ref) to convert it to a `MOI.ScalarAffineFunction{T}`.
"""
struct VectorScalarAffineFunction{T, VT} <: MOI.AbstractScalarFunction
    terms::VT
    constant::T
end
MOI.constant(func::VectorScalarAffineFunction) = func.constant
function JuMP.coefficient(func::VectorScalarAffineFunction, vi::MOI.VariableIndex)
    return func.terms[vi.value]
end
function Base.convert(::Type{MOI.ScalarAffineFunction{T}}, func::VectorScalarAffineFunction) where {T}
    return MOI.ScalarAffineFunction{T}(
        # TODO we should do better if the vector is a `SparseVector`, I think
        #      I have some code working for both vector types in Polyhedra.jl
        MOI.ScalarAffineTerm{T}[
            MOI.ScalarAffineTerm{T}(func.terms[i], VI(i))
            for i in eachindex(func.terms) if !iszero(func.terms[i])
        ],
        func.constant,
    )
end
function standard_form(func::VectorScalarAffineFunction{T}) where {T}
    return convert(MOI.ScalarAffineFunction{T}, func)
end

function MOI.Utilities.operate(::typeof(-), ::Type{T}, func::VectorScalarAffineFunction{T}) where {T}
    return VectorScalarAffineFunction(LazyArrays.ApplyArray(-, func.terms), -func.constant)
end

"""
    struct MatrixScalarQuadraticFunction{T, VT, MT} <: MOI.AbstractScalarFunction
        affine::VectorScalarAffineFunction{T,VT}
        terms::MT
    end

Represents the function `x' * terms * x / 2 + affine` as an
`MOI.AbstractScalarFunction` where `x[i] = MOI.VariableIndex(i)`.
Use [`standard_form`](@ref) to convert it to a `MOI.ScalarQuadraticFunction{T}`.
"""
struct MatrixScalarQuadraticFunction{T, VT, MT} <: MOI.AbstractScalarFunction
    affine::VectorScalarAffineFunction{T,VT}
    terms::MT
end
MOI.constant(func::MatrixScalarQuadraticFunction) = MOI.constant(func.affine)
function JuMP.coefficient(func::MatrixScalarQuadraticFunction, vi::MOI.VariableIndex)
    return JuMP.coefficient(func.affine, vi)
end
function quad_sym_half(
    func::MatrixScalarQuadraticFunction,
    vi1::MOI.VariableIndex,
    vi2::MOI.VariableIndex,
)
    return func.terms[vi1.value, vi2.value]
end
function JuMP.coefficient(
    func::MatrixScalarQuadraticFunction,
    vi1::MOI.VariableIndex,
    vi2::MOI.VariableIndex,
)
    coef = quad_sym_half(func, vi1, vi2)
    if vi1 == vi2
        return coef / 2
    else
        return coef
    end
end
function Base.convert(::Type{MOI.ScalarQuadraticFunction{T}}, func::MatrixScalarQuadraticFunction) where {T}
    n = length(func.affine.terms)
    aff = convert(MOI.ScalarAffineFunction{T}, func.affine)
    quad = MOI.ScalarQuadraticTerm{T}[
        MOI.ScalarQuadraticTerm{T}(quad_sym_half(func, VI(i), VI(j)), VI(i), VI(j))
        for j in 1:n for i in 1:j if !iszero(quad_sym_half(func, VI(i), VI(j)))
    ]
    return MOI.ScalarQuadraticFunction{T}(quad, aff.terms, aff.constant)
end
function standard_form(func::MatrixScalarQuadraticFunction{T}) where {T}
    return convert(MOI.ScalarQuadraticFunction{T}, func)
end

"""
    MatrixVectorAffineFunction{T, VT} <: MOI.AbstractVectorFunction

Represents the function `terms * x + constant`
as an `MOI.AbstractVectorFunction` where `x[i] = MOI.VariableIndex(i)`.
Use [`standard_form`](@ref) to convert it to a `MOI.VectorAffineFunction{T}`.
"""
struct MatrixVectorAffineFunction{AT, VT} <: MOI.AbstractVectorFunction
    terms::AT
    constants::VT
end
MOI.constant(func::MatrixVectorAffineFunction) = func.constants
function Base.convert(::Type{MOI.VectorAffineFunction{T}}, func::MatrixVectorAffineFunction) where {T}
    return MOI.VectorAffineFunction{T}(
        MOI.VectorAffineTerm{T}[
            # TODO we should do better if the matrix is a `SparseMatrixCSC`
            MOI.VectorAffineTerm(i, MOI.ScalarAffineTerm{T}(func.terms[i, j], VI(j)))
            for i in 1:size(func.terms, 1) for j in 1:size(func.terms, 2) if !iszero(func.terms[i, j])
        ],
        func.constants,
    )
end
function standard_form(func::MatrixVectorAffineFunction{T}) where {T}
    return convert(MOI.VectorAffineFunction{T}, func)
end

# Only used for testing at the moment so performance is not critical so
# converting to standard form is ok
function MOIU.isapprox_zero(func::Union{VectorScalarAffineFunction,MatrixScalarQuadraticFunction}, tol)
    return MOIU.isapprox_zero(standard_form(func), tol)
end

"""
    IndexMappedFunction{F<:MOI.AbstractFunction} <: AbstractLazyScalarFunction

Lazily represents the function `MOI.Utilities.map_indices(index_map, DiffOpt.standard_form(func))`.
"""
struct IndexMappedFunction{F<:MOI.AbstractFunction} <: AbstractLazyScalarFunction
    func::F
    index_map::MOIU.IndexMap
end
MOI.constant(func::IndexMappedFunction) = MOI.constant(func.func)
function JuMP.coefficient(func::IndexMappedFunction, vi::MOI.VariableIndex)
    return JuMP.coefficient(func.func, func.index_map[vi])
end
function quad_sym_half(func::IndexMappedFunction, vi1::MOI.VariableIndex, vi2::MOI.VariableIndex)
    return quad_sym_half(func.func, func.index_map[vi1], func.index_map[vi2])
end
function JuMP.coefficient(func::IndexMappedFunction, vi1::MOI.VariableIndex, vi2::MOI.VariableIndex)
    return JuMP.coefficient(func.func, func.index_map[vi1], func.index_map[vi2])
end
function standard_form(func::IndexMappedFunction)
    return MOIU.map_indices(func.index_map, standard_form(func.func))
end
MOIU.isapprox_zero(func::IndexMappedFunction, tol) = MOIU.isapprox_zero(func.func, tol)

function MOIU.map_indices(index_map::MOIU.IndexMap, func::AbstractLazyScalarFunction)
    return IndexMappedFunction(func, index_map)
end

"""
    MOItoJuMP{F<:MOI.AbstractScalarFunction} <: JuMP.AbstractJuMPScalar

Lazily represents the function `JuMP.jump_function(model, DiffOpt.standard_form(func))`.
"""
struct MOItoJuMP{F<:MOI.AbstractScalarFunction} <: JuMP.AbstractJuMPScalar
    model::JuMP.Model
    func::F
end
Base.broadcastable(func::MOItoJuMP) = Ref(func)
JuMP.constant(func::MOItoJuMP) = MOI.constant(func.func)
function JuMP.coefficient(func::MOItoJuMP, var_ref::JuMP.VariableRef)
    check_belongs_to_model(var_ref, func.model)
    return JuMP.coefficient(func.func, JuMP.index(var_ref))
end
function quad_sym_half(func::MOItoJuMP, var1_ref::JuMP.VariableRef, var2_ref::JuMP.VariableRef)
    check_belongs_to_model.([var1_ref, var2_ref], Ref(func.model))
    return quad_sym_half(func.func, JuMP.index(var1_ref), JuMP.index(var2_ref))
end
function JuMP.coefficient(func::MOItoJuMP, var1_ref::JuMP.VariableRef, var2_ref::JuMP.VariableRef)
    check_belongs_to_model.([var1_ref, var2_ref], Ref(func.model))
    return JuMP.coefficient(func.func, JuMP.index(var1_ref), JuMP.index(var2_ref))
end
function Base.convert(::Type{JuMP.GenericAffExpr{T,JuMP.VariableRef}}, func::MOItoJuMP) where {T}
    return JuMP.GenericAffExpr{T,JuMP.VariableRef}(func.model, convert(MOI.ScalarAffineFunction{T}, func.func))
end
function Base.convert(::Type{JuMP.GenericQuadExpr{T,JuMP.VariableRef}}, func::MOItoJuMP) where {T}
    return JuMP.GenericQuadExpr{T,JuMP.VariableRef}(func.model, convert(MOI.ScalarQuadraticFunction{T}, func.func))
end
JuMP.moi_function(func::MOItoJuMP) = func.func
function JuMP.jump_function(model::JuMP.Model, func::AbstractLazyScalarFunction)
    return MOItoJuMP(model, func)
end
function standard_form(func::MOItoJuMP)
    return JuMP.jump_function(func.model, standard_form(func.func))
end
function JuMP.function_string(mode, func::MOItoJuMP)
    return JuMP.function_string(mode, standard_form(func))
end

# JuMP
reverse_differentiate!(model::JuMP.Model) = reverse_differentiate!(JuMP.backend(model))
forward_differentiate!(model::JuMP.Model) = forward_differentiate!(JuMP.backend(model))

# MOIU
reverse_differentiate!(model::MOI.Utilities.CachingOptimizer) = reverse_differentiate!(model.optimizer)
forward_differentiate!(model::MOI.Utilities.CachingOptimizer) = forward_differentiate!(model.optimizer)

# MOIB
reverse_differentiate!(model::MOI.Bridges.AbstractBridgeOptimizer) = reverse_differentiate!(model.model)
forward_differentiate!(model::MOI.Bridges.AbstractBridgeOptimizer) = forward_differentiate!(model.model)


text:
"""
    particle dynamics
    - 3D particle subject to contact forces

    - configuration: q = (x, y, z) ∈ R³
    - impact force (magnitude): n ∈ R₊
    - friction force: b ∈ R²
        - contact force: λ = (b, n) ∈ R² × R₊
        - friction coefficient: μ ∈ R₊

    Discrete Mechanics and Variational Integrators
        pg. 363
"""
struct Particle
    n::Int
    m::Int
    d::Int

    mass # mass
    gravity # gravity
    friction_coeff # friction coefficient

    nq # configuration dimension
    nu # control dimension
end

# mass matrix
function mass_matrix(model)
    m = model.mass

    Diagonal(@SVector [m, m, m])
end

# gravity
function gravity(model, q)
    m = model.mass
    g = model.gravity

    @SVector [0.0, 0.0, m * g]
end

# signed distance function
function signed_distance(model, q)
    q[3]
end

# contact force Jacobian
function jacobian(model, q)
    Diagonal(@SVector ones(model.nq))
end

# dynamics
function dynamics(model, q1, q2, q3, u, λ, h)
    nq = model.nq
    nu = model.nu
    SVector{nq}(mass_matrix(model) * (2.0 * q2 - q1 - q3) / h
        - h * gravity(model, q2)
        + h * jacobian(model, q3)' * (u + λ))
end


text:
export ensemblerun!
Vector_or_Tuple = Union{AbstractArray,Tuple}

"""
    ensemblerun!(models::Vector, agent_step!, model_step!, n; kwargs...)
Perform an ensemble simulation of [`run!`](@ref) for all `model ∈ models`.
Each `model` should be a (different) instance of an [`AgentBasedModel`](@ref) but probably
initialized with a different random seed or different initial agent distribution.
All models obey the same rules `agent_step!, model_step!` and are evolved for `n`.

Similarly to [`run!`](@ref) this function will collect data. It will furthermore
add one additional column to the dataframe called `:ensemble`, which has an integer
value counting the ensemble member. The function returns `agent_df, model_df, models`.

The keyword `parallel = false`, when `true`, will run the simulations in parallel using
Julia's `Distributed.pmap` (you need to have loaded `Agents` with `@everywhere`, see
docs online).

All other keywords are propagated to [`run!`](@ref) as-is.

Example usage in [Schelling's segregation model](@ref).
"""
function ensemblerun!(
    models::Vector_or_Tuple,
    agent_step!,
    model_step!,
    n;
    parallel = false,
    kwargs...,
)
    if parallel
        return parallel_ensemble(models, agent_step!, model_step!, n; kwargs...)
    else
        return series_ensemble(models, agent_step!, model_step!, n; kwargs...)
    end
end

"""
    ensemblerun!(generator, agent_step!, model_step!, n; kwargs...)
Generate many `ABM`s and propagate them into `ensemblerun!(models, ...)` using
the provided `generator` which is a one-argument function whose input is a seed.

This method has additional keywords `ensemble = 5, seeds = rand(UInt32, ensemble)`.
"""
function ensemblerun!(
    generator,
    args...;
    ensemble = 5,
    seeds = rand(UInt32, ensemble),
    kwargs...,
)
    models = [generator(seed) for seed in seeds]
    ensemblerun!(models, args...; kwargs...)
end

function series_ensemble(models, agent_step!, model_step!, n; kwargs...)
    @assert models[1] isa ABM
    df_agent, df_model = _run!(models[1], agent_step!, model_step!, n; kwargs...)
    add_ensemble_index!(df_agent, 1)
    add_ensemble_index!(df_model, 1)

    for m in 2:length(models)
        df_agentTemp, df_modelTemp =
            _run!(models[m], agent_step!, model_step!, n; kwargs...)
        add_ensemble_index!(df_agentTemp, m)
        add_ensemble_index!(df_modelTemp, m)
        append!(df_agent, df_agentTemp)
        append!(df_model, df_modelTemp)
    end
    return df_agent, df_model, models
end

function parallel_ensemble(models, agent_step!, model_step!, n; kwargs...)
    all_data = pmap(
        j -> _run!(models[j], agent_step!, model_step!, n; kwargs...),
        1:length(models),
    )

    df_agent = DataFrame()
    df_model = DataFrame()
    for (m, d) in enumerate(all_data)
        add_ensemble_index!(d[1], m)
        add_ensemble_index!(d[2], m)
        append!(df_agent, d[1])
        append!(df_model, d[2])
    end

    return df_agent, df_model, models
end

function add_ensemble_index!(df, m)
    df[!, :ensemble] = fill(m, size(df, 1))
end


text:
immutable Point{T<:FloatingPoint}
    x::T
    y::T
end

immutable Circle{T<:FloatingPoint}
    c::Point{T}
    r::T
end
Circle{T<:FloatingPoint}(a::Point{T}) = Circle(a, zero(T))

using AffineTransforms

function circlepoints{T<:FloatingPoint}(a::Point{T}, b::Point{T}, r::T)
    cp = Circle{T}[]
    r >= 0 || return (cp, "No Solution, Negative Radius")
    if a == b
        if abs(r) < 2eps(zero(T))
            return (push!(cp, Circle(a)), "Point Solution, Zero Radius")
        else
            return (cp, "Infinite Solutions, Indefinite Center")
        end
    end
    ca = Complex(a.x, a.y)
    cb = Complex(b.x, b.y)
    d = (ca + cb)/2
    tfd = tformtranslate([real(d), imag(d)])
    tfr = tformrotate(angle(cb-ca))
    tfm = tfd*tfr
    u = abs(cb-ca)/2
    r-u > -5eps(r) || return(cp, "No Solution, Radius Too Small")
    if r-u < 5eps(r)
        push!(cp, Circle(apply(Point, tfm*[0.0, 0.0]), r))
        return return (cp, "Single Solution, Degenerate Centers")
    end
    v = sqrt(r^2 - u^2)
    for w in [v, -v]
        push!(cp, Circle(apply(Point, tfm*[0.0, w]), r))
    end
    return (cp, "Two Solutions")
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
    clear_cache()

CLBlast stores binaries of compiled kernels into a cache in case the same kernel
is used later on for thesame device. This cache can be cleared to free up system
memory or it can be useful in case of debugging.
"""
function clear_cache()
    err = ccall((:CLBlastClearCache, libCLBlast), cl.CL_int, ())
    if err != cl.CL_SUCCESS
        println(stderr, "Calling function `clear_cache` failed!")
        throw(cl.CLError(err))
    end
    return err
end

function fill_cache(device::cl.Device)
    err = ccall((:CLBlastFillCache, libCLBlast), cl.CL_int, (cl.CL_device_id,), pointer(device))
    if err != cl.CL_SUCCESS
        println(stderr, "Calling function `fill_cache($device)` failed!")
        throw(cl.CLError(err))
    end
    return err
end


text:
#=
"""
    transitive_closure!(A::AbstractMatrix{Int})

Compute the transitive closure of the matrix `A`.
"""
function transitive_closure!(A::AbstractMatrix{Int})
    mustlinkto = findall(val -> val == 1, C)
    for instance in mustlinkto
        first = instance[1]
        second = instance[2]
        C[C[first, second] == && C[second, :]] == 1
    end
    return A
end

"""
"""
function pckmeans(X::AbstractMatrix{<:Real}, C::AbstractMatrix{Int},
                W::AbstractMatrix{<:Real}, k::Int;
                maxiter::Int=256, dist::SemiMetric=SqEuclidean(),
                ϵ::AbstractFloat=1.0e-6)
    pckmeans(X, C, W, M, maxiter=maxiter, dist=dist, ϵ=ϵ)
end

"""
    pckmeans(X::AbstractMatrix{<:Real}, C::AbstractMatrix{Int}, W::AbstractMatrix{<:Real}, M::AbstractMatrix{<:Real}; <keyword arguments>)

Cluster the data `X` with the pairwise constrained ``k``-means algorithm.

# Keyword arguments
- `maxiter::Int=256`: the number of maximum iterations.
- `dist::SemiMetric=SqEuclidean()`: the distance function.
- `ϵ::AbstractFloat=1.0e-6`: the absolute tolerance for convergence.
"""
function pckmeans(X::AbstractMatrix{<:Real}, C::AbstractMatrix{Int},
                W::AbstractMatrix{<:Real}, M::AbstractMatrix{<:Real};
                maxiter::Int=256, dist::SemiMetric=SqEuclidean(),
                ϵ::AbstractFloat=1.0e-6)
    mx, nx = size(X)
    nc = size(C, 2)
    nw = size(W, 2)
    mm, k = size(M)

    mx == mm || throw(DimensionMismatch("number of data features must match"))
    nc == nx || throw(DimensionMismatch("number of data instances and maximum number of constraints must match"))
    nc == nw || throw(DimensionMismatch("dimensions of constraints and weights must match"))
    k > 1 || throw(ArgumentError("number of clusters must at least be 2"))
    k <= nx || throw(ArgumentError("more clusters than data instances are not allowed"))

    n = nx
    Y = zeros(Float64, k, n)
    M = convert(Matrix{Float64}, M)

    pcs = Vector{PartitionalClustering}(undef, 0)
    pc = PartitionalClustering(copy(X), copy(C), copy(W), copy(Y), copy(M))
    push!(pcs, pc)

    # Derive neighborhood sets, and initialize the cluster centers
    transitive_closure!(C)

    distances = pairwise(dist, M, X, dims=2)
    pre_objcosts = 0
    i = 1
    while i <= maxiter
        costs = distances
        objcosts = 0
        fill!(Y, zero(eltype(Y)))
        fill!(M, zero(eltype(M)))
        # Update cluster assignments, objective function costs, and center
        # coordinates per cluster
        @inbounds for j = 1:n
            cj = view(C, :, j)
            mustlinkto = findall(val -> val == 1, cj)
            for instance in mustlinkto
                y = argmax(view(Y, :, instance))
                costs[, j] += W[, ] * 1
            end
            cannotlink = findall(val -> val == -1, cj)

            cost, y = findmin(view(costs, :, j))
            Y[y, j] = 1
            objcosts += cost
            M[:, y] += X[:, j]
        end

        # Update cluster centers
        @inbounds for j = 1:k
            M[:, j] /= sum(Y[j, :])
        end

        pairwise!(distances, dist, M, X, dims=2)

        pc = PartitionalClustering(X, C, W, Y, M)
        push!(pcs, pc)

        # Check for convergence
        isapprox(objcosts, pre_objcosts, atol=ϵ) && break

        pre_objcosts = objcosts
        i += 1
    end

    return pcs
end
=#


text:
using RandomStrainDistributions
using Test

names(RandomStrainDistributions)

@testset "PhysicalVectors.jl" begin

    @testset "In-Place Operations" begin
        A = Vector2D(1, 2)
        @test equal(A, A)

        add!(A,A)
        ans = Vector2D(2, 4)
        @test equal(A, ans)

        A = Vector2D(1, 2)
        B = Vector2D(3., 4.)
        subtract!(A, B)
        @test equal(A, Vector2D(-2., -2.))

        A = Vector2D(1, 2)
        multiply!(A, 2.)
        @test equal(A, Vector2D(2., 4.))

        A = Vector2D(2, 4)
        divide!(A, 2.)
        @test equal(A, Vector2D(1., 2.))

        A = Vector2D(-100., 0.)
        normalize!(A)
        @test equal(A, Vector2D(-1., 0.))
    end

    @testset "Normal Operations" begin
        A = Vector2D(1., 2.)
        B = Vector2D(3., 4.)
        @test equal(A + B, Vector2D(4., 6.))

        @test A ⋅ B == 11 && dot(A, B) == A ⋅ B
        @test equal(2 * A, Vector2D(2., 4.)) && equal(2 * A, A * 2)
        @test equal(A / 2, Vector2D(0.5, 1.))
        @test magnitude2(B) == 25. && magnitude(B) == 5.
        @test equal(normalize(Vector2D(0., 20.)), Vector2D(0., 1.))
    end
end

@testset "RandomStrainDistributions.jl" begin

end


text:
### Gradient-based optimization tests
# The following NLopt gradient -based algorithms are compared
grad_algorithms = [
    :LD_MMA, # Method of Moving Asymptotes
    :LD_SLSQP, # Sequential Least-Squares Quadratic Programming
    :LD_CCSAQ, # Conservative convex separable approximation
    # errors
    #:LD_LBFGS, # Low-storage BFGS
    #:LD_TNEWTON_PRECOND_RESTART, # Preconditioned truncated Newton
    #:LD_TNEWTON_PRECOND, # Same without restarting
    #:LD_TNEWTON_RESTART, # Same without preconditioning
    #:LD_TNEWTON, # Same without restarting or preconditioning
    #:LD_VAR2, # Shifted limited-memory variable-metric (rank 2)
    #:LD_VAR1  # Shifted limited-memory variable-metric (rank 1)
]

[test_alg(alg; bounds=(-1e10,1e10)) for alg in grad_algorithms]


text:
@testset "Operator Splitting" begin

    # Set parameters
    n_particles = 2
    eta_min = 0.0
    eta_max = 4π
    num_cells = 10
    delta_t = 0.1
    degree_smoother = 3

    pg = ParticleGroup{2,2}(n_particles, charge = 1.0, mass = 1.0, n_weights = 1)

    particle_info_ref = [
        11.780972450961723 5.4977871437821380
        0.78539816339744828 7.0685834705770345
        0.15731068461017067 -1.5341205443525459
        1.5341205443525459 -0.15731068461017067
        86.251495608834688 71.662174808595040
    ]

    for i_part = 1:n_particles
        pg.array[1, i_part] = particle_info_ref[1, i_part]
        pg.array[2, i_part] = particle_info_ref[2, i_part]
        pg.array[3, i_part] = particle_info_ref[3, i_part]
        pg.array[4, i_part] = particle_info_ref[4, i_part]
        pg.array[5, i_part] = particle_info_ref[5, i_part]
    end

    @test pg.array ≈ particle_info_ref

    grid = TwoDGrid(eta_min, eta_max, num_cells, eta_min, eta_max, num_cells)

    kernel = ParticleMeshCoupling2D(pg, grid, degree_smoother, :collocation)

    poisson = TwoDPoissonPeriodic(grid)

    pic = TwoDPoissonPIC(poisson, kernel)

    propagator = SplittingOperator(pic, pg)

    ParticleInCell.operator_t!(propagator, delta_t)

    particle_info_ref = [
        11.796703519422740 5.3443750893468831
        0.93881021783270291 7.0528524021160175
        0.15731068461017067 -1.5341205443525459
        1.5341205443525459 -0.15731068461017067
        86.251495608834688 71.662174808595040
    ]

    @test pg.array ≈ particle_info_ref

    charge_deposition!(propagator)
    solve_fields!(propagator)
    ParticleInCell.operator_v!(propagator, delta_t)

    particle_info_ref = [
        11.796703519422740 5.3443750893468831
        0.93881021783270291 7.0528524021160175
        0.15346588344558551 -1.5294930005799796
        1.5302579166423358 -0.15266168508042444
        86.251495608834688 71.662174808595040
    ]

    @test pg.array ≈ particle_info_ref

end


text:
module Covering1D

using ..NumericalContinuation: NumericalContinuation, Data, Problem
using ..NumericalContinuation: monitor_function
using ..NumericalContinuation: ALLVARS

function pseudoarclength(u, data)
    result = zero(eltype(u))
    for i in eachindex(u)
        result += (u[i] - data.u[i]) * data.t[i]
    end
    return result
end

struct Cover <: ProblemOwner end

function (::Cover)(::Signal{:pre_correct}, indices, problem)
    chart = get_current_chart(problem)
    u = copy(get_u(chart))
    t = copy(get_t(chart))
    tv = get_uview(chart, t)
    tv[indices[2]] = 0
    data = get_data(chart)
    data[indices[1]] = (u=u, t=t)
    return nothing
end

function covering1d(problem::Problem)
    data = Data("pseudoarclength", nothing)
    mfunc = monitor_function("pseudoarclength", pseudoarclength, (ALLVARS,), (data,))
    cover = Problem("covering", Cover(), (mfunc,))
    pass_indices(cover, data)
    return pass_indices(cover, mfunc.var[1])
end

end  # module


text:
using Distributed
using JSON3

addprocs()

@everywhere include("./src/JuliaDBQuery.jl")
@everywhere using .JuliaDBQuery
@everywhere using Dates
@everywhere include("./src/Web.jl")
@everywhere using .Web
@everywhere using OnlineStats
@everywhere using JuliaDB

data = JuliaDBQuery.load_table("pop")
sa2_filter = """{ "sa2_code":{ "\$gt" : 100100, "\$lt":100300}}"""
q = JSON3.read(sa2_filter, Dict)
a,b = JuliaDBQuery.predicate_for(q)
@time filter(a, data, select=b)

Web.start()